# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit for
# the Earth and Planetary Sciences
# Copyright (C) 2012 - 2022 by the BurnMan team, released under the GNU
# GPL v2 or later.
from __future__ import absolute_import
import numpy as np
from sympy import Matrix, nsimplify
from collections import OrderedDict
from .material import material_property, cached_property
from .mineral import Mineral
from .solutionmodel import MechanicalSolution
from .solutionmodel import PolynomialSolution
from .averaging_schemes import reuss_average_function
from ..utils.reductions import independent_row_indices
from ..utils.chemistry import sum_formulae, sort_element_list_to_IUPAC_order
[docs]class Solution(Mineral):
"""
This is the base class for all solutions.
Site occupancies, endmember activities and the constant
and pressure and temperature dependencies of the excess
properties can be queried after using set_composition().
States of the solution can only be queried after setting
the pressure, temperature and composition using set_state().
This class is available as :class:`burnman.Solution`.
It uses an instance of :class:`burnman.SolutionModel` to
calculate interaction terms between endmembers.
All the solution parameters are expected to be in SI units. This
means that the interaction parameters should be in J/mol, with the T
and P derivatives in J/K/mol and m^3/mol.
The parameters are relevant to all solution models. Please
see the documentation for individual models for details about
other parameters.
:param name: Name of the solution.
:type name: string
:param solution_model: The SolutionModel object defining the properties
of the solution.
:type solution_model: :class:`burnman.SolutionModel`
:param molar_fractions: The molar fractions of each endmember in the solution.
Can be reset using the set_composition() method.
:type molar_fractions: numpy.array
"""
def __init__(self, name=None, solution_model=None, molar_fractions=None):
"""
Set up matrices to speed up calculations for when P, T, X is defined.
"""
Mineral.__init__(self)
# Solution needs a method attribute to call Mineral.set_state().
# Note that set_method() below will not change self.method
self.method = "SolutionMethod"
if name is not None:
self.name = name
if solution_model is not None:
self.solution_model = solution_model
# Equation of state
for i in range(self.n_endmembers):
self.solution_model.endmembers[i][0].set_method(
self.solution_model.endmembers[i][0].params["equation_of_state"]
)
# Molar fractions
if molar_fractions is not None:
self.set_composition(molar_fractions)
[docs] @cached_property
def endmembers(self):
return self.solution_model.endmembers
[docs] def set_composition(self, molar_fractions):
"""
Set the composition for this solution.
Resets cached properties.
:param molar_fractions: Molar abundance for each endmember, needs to sum to one.
:type molar_fractions: list of float
"""
assert len(self.solution_model.endmembers) == len(molar_fractions)
if type(self.solution_model) != MechanicalSolution:
assert sum(molar_fractions) > 0.9999
assert sum(molar_fractions) < 1.0001
if type(self.solution_model) == PolynomialSolution:
self.solution_model.set_composition(molar_fractions)
self.reset()
self.molar_fractions = np.array(molar_fractions)
[docs] def set_method(self, method):
for i in range(self.n_endmembers):
self.solution_model.endmembers[i][0].set_method(method)
# note: do not set self.method here!
self.reset()
[docs] def set_state(self, pressure, temperature):
if type(self.solution_model) == PolynomialSolution:
self.solution_model.set_state(pressure, temperature)
Mineral.set_state(self, pressure, temperature)
for i in range(self.n_endmembers):
self.solution_model.endmembers[i][0].set_state(pressure, temperature)
@material_property
def formula(self):
"""
Returns molar chemical formula of the solution.
:rtype: Counter
"""
return sum_formulae(self.endmember_formulae, self.molar_fractions)
@material_property
def site_occupancies(self):
"""
:returns: The fractional occupancies of species on each site.
:rtype: list of OrderedDicts
"""
occs = np.einsum(
"ij, i", self.solution_model.endmember_occupancies, self.molar_fractions
)
site_occs = []
k = 0
for i in range(self.solution_model.n_sites):
site_occs.append(OrderedDict())
for j in range(len(self.solution_model.sites[i])):
site_occs[-1][self.solution_model.sites[i][j]] = occs[k]
k += 1
return site_occs
@material_property
def activities(self):
"""
Returns a list of endmember activities [unitless].
"""
return self.solution_model.activities(
self.pressure, self.temperature, self.molar_fractions
)
@material_property
def activity_coefficients(self):
"""
Returns a list of endmember activity coefficients
(gamma = activity / ideal activity) [unitless].
"""
return self.solution_model.activity_coefficients(
self.pressure, self.temperature, self.molar_fractions
)
@material_property
def molar_internal_energy(self):
"""
Returns molar internal energy of the mineral [J/mol].
Aliased with self.energy
"""
return self.molar_helmholtz + self.temperature * self.molar_entropy
@material_property
def excess_partial_gibbs(self):
"""
Returns excess partial molar gibbs free energy [J/mol].
Property specific to solutions.
"""
return self.solution_model.excess_partial_gibbs_free_energies(
self.pressure, self.temperature, self.molar_fractions
)
@material_property
def excess_partial_volumes(self):
"""
Returns excess partial volumes [m^3].
Property specific to solutions.
"""
return self.solution_model.excess_partial_volumes(
self.pressure, self.temperature, self.molar_fractions
)
@material_property
def excess_partial_entropies(self):
"""
Returns excess partial entropies [J/K].
Property specific to solutions.
"""
return self.solution_model.excess_partial_entropies(
self.pressure, self.temperature, self.molar_fractions
)
@material_property
def partial_gibbs(self):
"""
Returns endmember partial molar gibbs free energy [J/mol].
Property specific to solutions.
"""
return (
np.array(
[
self.solution_model.endmembers[i][0].gibbs
for i in range(self.n_endmembers)
]
)
+ self.excess_partial_gibbs
)
@material_property
def partial_volumes(self):
"""
Returns endmember partial volumes [m^3].
Property specific to solutions.
"""
return (
np.array(
[
self.solution_model.endmembers[i][0].molar_volume
for i in range(self.n_endmembers)
]
)
+ self.excess_partial_volumes
)
@material_property
def partial_entropies(self):
"""
Returns endmember partial entropies [J/K].
Property specific to solutions.
"""
return (
np.array(
[
self.solution_model.endmembers[i][0].molar_entropy
for i in range(self.n_endmembers)
]
)
+ self.excess_partial_entropies
)
@material_property
def excess_gibbs(self):
"""
Returns molar excess gibbs free energy [J/mol].
Property specific to solutions.
"""
return self.solution_model.excess_gibbs_free_energy(
self.pressure, self.temperature, self.molar_fractions
)
@material_property
def gibbs_hessian(self):
"""
Returns an array containing the second compositional derivative
of the Gibbs free energy [J]. Property specific to solutions.
"""
return self.solution_model.gibbs_hessian(
self.pressure, self.temperature, self.molar_fractions
)
@material_property
def entropy_hessian(self):
"""
Returns an array containing the second compositional derivative
of the entropy [J/K]. Property specific to solutions.
"""
return self.solution_model.entropy_hessian(
self.pressure, self.temperature, self.molar_fractions
)
@material_property
def volume_hessian(self):
"""
Returns an array containing the second compositional derivative
of the volume [m^3]. Property specific to solutions.
"""
return self.solution_model.volume_hessian(
self.pressure, self.temperature, self.molar_fractions
)
@material_property
def molar_gibbs(self):
"""
Returns molar Gibbs free energy of the solution [J/mol].
Aliased with self.gibbs.
"""
return (
sum(
[
self.solution_model.endmembers[i][0].gibbs * self.molar_fractions[i]
for i in range(self.n_endmembers)
]
)
+ self.excess_gibbs
)
@material_property
def molar_helmholtz(self):
"""
Returns molar Helmholtz free energy of the solution [J/mol].
Aliased with self.helmholtz.
"""
return self.molar_gibbs - self.pressure * self.molar_volume
@material_property
def molar_mass(self):
"""
Returns molar mass of the solution [kg/mol].
"""
return sum(
[
self.solution_model.endmembers[i][0].molar_mass
* self.molar_fractions[i]
for i in range(self.n_endmembers)
]
)
@material_property
def excess_volume(self):
"""
Returns excess molar volume of the solution [m^3/mol].
Specific property for solutions.
"""
return self.solution_model.excess_volume(
self.pressure, self.temperature, self.molar_fractions
)
@material_property
def molar_volume(self):
"""
Returns molar volume of the solution [m^3/mol].
Aliased with self.V.
"""
return (
sum(
[
self.solution_model.endmembers[i][0].molar_volume
* self.molar_fractions[i]
for i in range(self.n_endmembers)
]
)
+ self.excess_volume
)
@material_property
def density(self):
"""
Returns density of the solution [kg/m^3].
Aliased with self.rho.
"""
return self.molar_mass / self.molar_volume
@material_property
def excess_entropy(self):
"""
Returns excess molar entropy [J/K/mol].
Property specific to solutions.
"""
return self.solution_model.excess_entropy(
self.pressure, self.temperature, self.molar_fractions
)
@material_property
def molar_entropy(self):
"""
Returns molar entropy of the solution [J/K/mol].
Aliased with self.S.
"""
return (
sum(
[
self.solution_model.endmembers[i][0].S * self.molar_fractions[i]
for i in range(self.n_endmembers)
]
)
+ self.excess_entropy
)
@material_property
def excess_enthalpy(self):
"""
Returns excess molar enthalpy [J/mol].
Property specific to solutions.
"""
return self.solution_model.excess_enthalpy(
self.pressure, self.temperature, self.molar_fractions
)
@material_property
def molar_enthalpy(self):
"""
Returns molar enthalpy of the solution [J/mol].
Aliased with self.H.
"""
return (
sum(
[
self.solution_model.endmembers[i][0].H * self.molar_fractions[i]
for i in range(self.n_endmembers)
]
)
+ self.excess_enthalpy
)
@material_property
def isothermal_bulk_modulus(self):
"""
Returns isothermal bulk modulus of the solution [Pa].
Aliased with self.K_T.
"""
return (
self.V
* 1.0
/ (
sum(
[
self.solution_model.endmembers[i][0].V
/ (self.solution_model.endmembers[i][0].K_T)
* self.molar_fractions[i]
for i in range(self.n_endmembers)
]
)
+ self.solution_model.VoverKT_excess()
)
)
@material_property
def adiabatic_bulk_modulus(self):
"""
Returns adiabatic bulk modulus of the solution [Pa].
Aliased with self.K_S.
"""
if self.temperature < 1e-10:
return self.isothermal_bulk_modulus
else:
return (
self.isothermal_bulk_modulus
* self.molar_heat_capacity_p
/ self.molar_heat_capacity_v
)
@material_property
def isothermal_compressibility(self):
"""
Returns isothermal compressibility of the solution.
(or inverse isothermal bulk modulus) [1/Pa].
Aliased with self.K_T.
"""
return 1.0 / self.isothermal_bulk_modulus
@material_property
def adiabatic_compressibility(self):
"""
Returns adiabatic compressibility of the solution.
(or inverse adiabatic bulk modulus) [1/Pa].
Aliased with self.K_S.
"""
return 1.0 / self.adiabatic_bulk_modulus
@material_property
def shear_modulus(self):
"""
Returns shear modulus of the solution [Pa].
Aliased with self.G.
"""
G_list = np.fromiter(
(e[0].G for e in self.solution_model.endmembers),
dtype=float,
count=self.n_endmembers,
)
return reuss_average_function(self.molar_fractions, G_list)
@material_property
def p_wave_velocity(self):
"""
Returns P wave speed of the solution [m/s].
Aliased with self.v_p.
"""
return np.sqrt(
(self.adiabatic_bulk_modulus + 4.0 / 3.0 * self.shear_modulus)
/ self.density
)
@material_property
def bulk_sound_velocity(self):
"""
Returns bulk sound speed of the solution [m/s].
Aliased with self.v_phi.
"""
return np.sqrt(self.adiabatic_bulk_modulus / self.density)
@material_property
def shear_wave_velocity(self):
"""
Returns shear wave speed of the solution [m/s].
Aliased with self.v_s.
"""
return np.sqrt(self.shear_modulus / self.density)
@material_property
def grueneisen_parameter(self):
"""
Returns grueneisen parameter of the solution [unitless].
Aliased with self.gr.
"""
if self.temperature < 1e-10:
return float("nan")
else:
return (
self.thermal_expansivity
* self.isothermal_bulk_modulus
* self.molar_volume
/ self.molar_heat_capacity_v
)
@material_property
def thermal_expansivity(self):
"""
Returns thermal expansion coefficient (alpha)
of the solution [1/K].
Aliased with self.alpha.
"""
return (1.0 / self.V) * (
sum(
[
self.solution_model.endmembers[i][0].alpha
* self.solution_model.endmembers[i][0].V
* self.molar_fractions[i]
for i in range(self.n_endmembers)
]
)
+ self.solution_model.alphaV_excess()
)
@material_property
def molar_heat_capacity_v(self):
"""
Returns molar heat capacity at constant volume of the
solution [J/K/mol].
Aliased with self.C_v.
"""
return (
self.molar_heat_capacity_p
- self.molar_volume
* self.temperature
* self.thermal_expansivity
* self.thermal_expansivity
* self.isothermal_bulk_modulus
)
@material_property
def molar_heat_capacity_p(self):
"""
Returns molar heat capacity at constant pressure
of the solution [J/K/mol].
Aliased with self.C_p.
"""
return (
sum(
[
self.solution_model.endmembers[i][0].molar_heat_capacity_p
* self.molar_fractions[i]
for i in range(self.n_endmembers)
]
)
+ self.solution_model.Cp_excess()
)
[docs] @cached_property
def stoichiometric_matrix(self):
"""
A sympy Matrix where each element M[i,j] corresponds
to the number of atoms of element[j] in endmember[i].
"""
def f(i, j):
e = self.elements[j]
if e in self.endmember_formulae[i]:
return nsimplify(self.endmember_formulae[i][e])
else:
return 0
return Matrix(len(self.endmember_formulae), len(self.elements), f)
[docs] @cached_property
def stoichiometric_array(self):
"""
An array where each element arr[i,j] corresponds
to the number of atoms of element[j] in endmember[i].
"""
return np.array(self.stoichiometric_matrix)
[docs] @cached_property
def reaction_basis(self):
"""
An array where each element arr[i,j] corresponds
to the number of moles of endmember[j] involved in reaction[i].
"""
reaction_basis = np.array(
[v[:] for v in self.stoichiometric_matrix.T.nullspace()]
)
if len(reaction_basis) == 0:
reaction_basis = np.empty((0, len(self.endmember_names)))
return reaction_basis
[docs] @cached_property
def n_reactions(self):
"""
The number of reactions in reaction_basis.
"""
return len(self.reaction_basis[:, 0])
[docs] @cached_property
def independent_element_indices(self):
"""
A list of an independent set of element indices. If the amounts of
these elements are known (element_amounts),
the amounts of the other elements can be inferred by
-compositional_null_basis[independent_element_indices].dot(element_amounts).
"""
return sorted(independent_row_indices(self.stoichiometric_matrix.T))
[docs] @cached_property
def dependent_element_indices(self):
"""
The element indices not included in the independent list.
"""
return [
i
for i in range(len(self.elements))
if i not in self.independent_element_indices
]
[docs] @cached_property
def compositional_null_basis(self):
"""
An array N such that N.b = 0 for all bulk compositions that can
be produced with a linear sum of the endmembers in the solution.
"""
null_basis = np.array([v[:] for v in self.stoichiometric_matrix.nullspace()])
M = null_basis[:, self.dependent_element_indices]
assert (M.shape[0] == M.shape[1]) and (M == np.eye(M.shape[0])).all()
return null_basis
[docs] @cached_property
def endmember_names(self):
"""
A list of names for all the endmember in the solution.
"""
return [mbr[0].name for mbr in self.solution_model.endmembers]
[docs] @cached_property
def n_endmembers(self):
"""
The number of endmembers in the solution.
"""
return len(self.solution_model.endmembers)
[docs] @cached_property
def elements(self):
"""
A list of the elements which could be contained in the solution,
returned in the IUPAC element order.
"""
keys = []
for f in self.endmember_formulae:
keys.extend(f.keys())
return sort_element_list_to_IUPAC_order(set(keys))
SolidSolution = Solution