# 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 .material import material_property, cached_property
from .mineral import Mineral
from .solutionmodel import SolutionModel
from .solutionmodel import MechanicalSolution, IdealSolution
from .solutionmodel import SymmetricRegularSolution, AsymmetricRegularSolution
from .solutionmodel import SubregularSolution, FunctionSolution
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.
Parameters
----------
name : string
Name of the solution
solution_type : string
String determining which SolutionModel to use. One of 'mechanical',
'ideal', 'symmetric', 'asymmetric' or 'subregular'.
endmembers : list of lists
List of endmembers in this solution. The first item of each
list should be a :class:`burnman.Mineral` object. The second item
should be a string with the site formula of the endmember.
molar_fractions : numpy array (optional)
The molar fractions of each endmember in the solution.
Can be reset using the set_composition() method.
"""
def __init__(self,
name=None,
solution_type=None,
endmembers=None,
energy_interaction=None,
volume_interaction=None,
entropy_interaction=None,
energy_ternary_terms=None,
volume_ternary_terms=None,
entropy_ternary_terms=None,
alphas=None,
excess_gibbs_function=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_type is not None:
self.solution_type = solution_type
if endmembers is not None:
self.endmembers = endmembers
if energy_interaction is not None:
self.energy_interaction = energy_interaction
if volume_interaction is not None:
self.volume_interaction = volume_interaction
if entropy_interaction is not None:
self.entropy_interaction = entropy_interaction
if energy_ternary_terms is not None:
self.energy_ternary_terms = energy_ternary_terms
if volume_ternary_terms is not None:
self.volume_ternary_terms = volume_ternary_terms
if entropy_ternary_terms is not None:
self.entropy_ternary_terms = entropy_ternary_terms
if alphas is not None:
self.alphas = alphas
if endmembers is not None:
self.endmembers = endmembers
if excess_gibbs_function is not None:
self.excess_gibbs_function = excess_gibbs_function
if hasattr(self, 'endmembers') is False:
raise Exception("'endmembers' attribute missing "
"from solution")
# Set default solution model type
if hasattr(self, 'solution_type'):
if self.solution_type == 'mechanical':
self.solution_model = MechanicalSolution(self.endmembers)
elif self.solution_type == 'ideal':
self.solution_model = IdealSolution(self.endmembers)
else:
if hasattr(self, 'energy_interaction') is False:
self.energy_interaction = None
if hasattr(self, 'volume_interaction') is False:
self.volume_interaction = None
if hasattr(self, 'entropy_interaction') is False:
self.entropy_interaction = None
if self.solution_type == 'symmetric':
self.solution_model = SymmetricRegularSolution(
self.endmembers, self.energy_interaction,
self.volume_interaction, self.entropy_interaction)
elif self.solution_type == 'asymmetric':
if hasattr(self, 'alphas') is False:
raise Exception(
"'alphas' attribute missing from solution")
self.solution_model = AsymmetricRegularSolution(
self.endmembers, self.alphas, self.energy_interaction,
self.volume_interaction, self.entropy_interaction)
elif self.solution_type == 'subregular':
if hasattr(self, 'energy_ternary_terms') is False:
self.energy_ternary_terms = None
if hasattr(self, 'volume_ternary_terms') is False:
self.volume_ternary_terms = None
if hasattr(self, 'entropy_ternary_terms') is False:
self.entropy_ternary_terms = None
self.solution_model = SubregularSolution(
self.endmembers,
self.energy_interaction, self.volume_interaction,
self.entropy_interaction,
self.energy_ternary_terms, self.volume_ternary_terms,
self.entropy_ternary_terms)
elif self.solution_type == 'function':
self.solution_model = FunctionSolution(self.endmembers,
self.excess_gibbs_function)
else:
raise Exception("Solution model type "
+ self.solution_type + "not recognised.")
else:
self.solution_model = SolutionModel()
# Equation of state
for i in range(self.n_endmembers):
self.endmembers[i][0].set_method(
self.endmembers[i][0].params['equation_of_state'])
# Molar fractions
if molar_fractions is not None:
self.set_composition(molar_fractions)
[docs] def get_endmembers(self):
return self.endmembers
[docs] def set_composition(self, molar_fractions):
"""
Set the composition for this solution.
Resets cached properties.
Parameters
----------
molar_fractions: list of float
molar abundance for each endmember, needs to sum to one.
"""
assert(len(self.endmembers) == len(molar_fractions))
if self.solution_type != 'mechanical':
assert(sum(molar_fractions) > 0.9999)
assert(sum(molar_fractions) < 1.0001)
self.reset()
self.molar_fractions = np.array(molar_fractions)
[docs] def set_method(self, method):
for i in range(self.n_endmembers):
self.endmembers[i][0].set_method(method)
# note: do not set self.method here!
self.reset()
[docs] def set_state(self, pressure, temperature):
Mineral.set_state(self, pressure, temperature)
for i in range(self.n_endmembers):
self.endmembers[i][0].set_state(pressure, temperature)
@material_property
def formula(self):
"""
Returns molar chemical formula of the solution.
"""
return sum_formulae(self.endmember_formulae, self.molar_fractions)
@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.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.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.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.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.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.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.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.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. / (sum([self.endmembers[i][0].V
/ (self.endmembers[i][0].K_T)
* self.molar_fractions[i]
for i in range(self.n_endmembers)]))
@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. / 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. / 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.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. / 3. * 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. / self.V) * sum([self.endmembers[i][0].alpha
* self.endmembers[i][0].V
* self.molar_fractions[i]
for i in range(self.n_endmembers)])
@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.endmembers[i][0].molar_heat_capacity_p
* self.molar_fractions[i]
for i in range(self.n_endmembers)])
[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.endmembers]
[docs] @cached_property
def n_endmembers(self):
"""
The number of endmembers in the solution.
"""
return len(self.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