# Thermodynamics#

Burnman has a number of functions and classes which deal with the thermodynamics of single phases and aggregates.

## Lattice Vibrations#

### Debye model#

burnman.eos.debye.debye_fn(x)[source]#

Evaluate the Debye function. Takes the parameter xi = Debye_T/T

burnman.eos.debye.debye_fn_cheb(x)[source]#

Evaluate the Debye function using a Chebyshev series expansion coupled with asymptotic solutions of the function. Shamelessly adapted from the GSL implementation of the same function (Itself adapted from Collected Algorithms from ACM). Should give the same result as debye_fn(x) to near machine-precision.

burnman.eos.debye.thermal_energy(T, debye_T, n)[source]#

calculate the thermal energy of a substance. Takes the temperature, the Debye temperature, and n, the number of atoms per molecule. Returns thermal energy in J/mol

burnman.eos.debye.molar_heat_capacity_v(T, debye_T, n)[source]#

Heat capacity at constant volume. In J/K/mol

burnman.eos.debye.helmholtz_free_energy(T, debye_T, n)[source]#

Helmholtz free energy of lattice vibrations in the Debye model [J]. It is important to note that this does NOT include the zero point energy for the lattice. As long as you are calculating relative differences in F, this should cancel anyways.

burnman.eos.debye.entropy(T, debye_T, n)[source]#

Entropy due to lattice vibrations in the Debye model [J/K].

burnman.eos.debye.dmolar_heat_capacity_v_dT(T, debye_T, n)[source]#

First temperature derivative of the heat capacity at constant volume [J/K^2/mol].

### Einstein model#

burnman.eos.einstein.thermal_energy(T, einstein_T, n)[source]#

calculate the thermal energy of a substance. Takes the temperature, the Einstein temperature, and n, the number of atoms per molecule. Returns thermal energy in J/mol

burnman.eos.einstein.molar_heat_capacity_v(T, einstein_T, n)[source]#

Heat capacity at constant volume. In J/K/mol

burnman.eos.einstein.helmholtz_free_energy(T, einstein_T, n)[source]#

Helmholtz free energy of lattice vibrations in the Einstein model [J]. It is important to note that this does NOT include the zero point energy for the lattice. As long as you are calculating relative differences in F, this should cancel anyway.

burnman.eos.einstein.entropy(T, einstein_T, n)[source]#

Entropy due to lattice vibrations in the Einstein model [J/K]

burnman.eos.einstein.dmolar_heat_capacity_v_dT(T, einstein_T, n)[source]#

First temperature derivative of the heat capacity at constant volume according to the Einstein model [J/K^2/mol].

## Chemistry parsing and thermodynamics#

A simple function to read a file with a two column list of elements and their masses into a dictionary

burnman.utils.chemistry.atomic_masses = {'Ag': 0.107868, 'Al': 0.0269815, 'Ar': 0.039948, 'As': 0.0749216, 'Au': 0.196967, 'B': 0.010811, 'Ba': 0.137327, 'Be': 0.00901218, 'Bi': 0.20898, 'Br': 0.079904, 'C': 0.0120107, 'Ca': 0.040078, 'Cd': 0.112411, 'Ce': 0.140116, 'Cl': 0.035453, 'Co': 0.0589332, 'Cr': 0.0519961, 'Cs': 0.132905, 'Cu': 0.063546, 'Dy': 0.1625, 'Er': 0.167259, 'Eu': 0.151964, 'F': 0.0189984, 'Fe': 0.055845, 'Ga': 0.069723, 'Gd': 0.15725, 'Ge': 0.07264, 'H': 0.00100794, 'He': 0.0040026, 'Hf': 0.17849, 'Hg': 0.20059, 'Ho': 0.16493, 'I': 0.126904, 'In': 0.114818, 'Ir': 0.192217, 'K': 0.0390983, 'Kr': 0.083798, 'La': 0.138905, 'Li': 0.006941, 'Lu': 0.174967, 'Mg': 0.024305, 'Mn': 0.054938, 'Mo': 0.09596, 'N': 0.0140067, 'Na': 0.0229898, 'Nb': 0.0929064, 'Nd': 0.144242, 'Ne': 0.0201797, 'Ni': 0.0586934, 'O': 0.0159994, 'Os': 0.19023, 'P': 0.0309738, 'Pa': 0.231036, 'Pb': 0.2072, 'Pd': 0.10642, 'Pr': 0.140908, 'Pt': 0.195084, 'Rb': 0.0854678, 'Re': 0.186207, 'Rh': 0.102905, 'Ru': 0.10107, 'S': 0.032065, 'Sb': 0.12176, 'Sc': 0.0449559, 'Se': 0.07896, 'Si': 0.0280855, 'Sm': 0.15036, 'Sn': 0.11871, 'Sr': 0.08762, 'Ta': 0.180948, 'Tb': 0.158925, 'Te': 0.1276, 'Th': 0.232038, 'Ti': 0.047867, 'Tl': 0.204383, 'Tm': 0.168934, 'U': 0.238029, 'V': 0.0509415, 'Vc': 0.0, 'W': 0.18384, 'Xe': 0.131293, 'Y': 0.0889058, 'Yb': 0.173054, 'Zn': 0.06538, 'Zr': 0.091224}#

IUPAC_element_order provides a list of all the elements. Element order is based loosely on electronegativity, following the scheme suggested by IUPAC, except that H comes after the Group 16 elements, not before them.

burnman.utils.chemistry.dictionarize_formula(formula)[source]#

A function to read a chemical formula string and convert it into a dictionary

Parameters:

formula (str) – Chemical formula, written in the XnYm format, where the formula has n atoms of element X and m atoms of element Y

Returns:

The same chemical formula, but expressed as a dictionary.

Return type:

dict

burnman.utils.chemistry.sum_formulae(formulae, amounts=None)[source]#

Adds together a set of formulae.

Parameters:
• formulae (list of dictionary or counter objects) – List of chemical formulae.

• amounts (list of floats) – List of amounts of each formula.

Returns:

The sum of the user-provided formulae

Return type:

Counter object

burnman.utils.chemistry.formula_mass(formula)[source]#

A function to take a chemical formula and compute the formula mass.

Parameters:

formula (dict or Counter object) – A chemical formula

Returns:

The mass per mole of formula [kg]

Return type:

float

burnman.utils.chemistry.convert_formula(formula, to_type='mass', normalize=False)[source]#

Converts a chemical formula from one type (mass or molar) into the other. Renormalises amounts if normalize=True.

Parameters:
• formula (dict or Counter object) – A chemical formula.

• to_type (str) – Conversion type, one of ‘mass’ or ‘molar’.

• normalize (bool) – Whether or not to normalize the converted formula to 1.

Returns:

The converted formula.

Return type:

dict

burnman.utils.chemistry.process_solution_chemistry(solution_model)[source]#

This function parses a class instance with a “formulas” attribute containing site information, e.g.

[ ‘[Mg]3[Al]2Si3O12’, ‘[Mg]3[Mg1/2Si1/2]2Si3O12’ ]

It outputs the bulk composition of each endmember (removing the site information), and also a set of variables and arrays which contain the site information. These are output in a format that can easily be used to calculate activities and gibbs free energies, given molar fractions of the phases and pressure and temperature where necessary.

Parameters:

solution_model – Class must have a “formulas” attribute, containing a list of chemical formulae with site information

Return type:

None

Note

Nothing is returned from this function, but the solution_model object gains the following attributes:

• solution_formulae [list of dictionaries]

List of endmember formulae in dictionary form.

• empty_formula [string]

Abbreviated chemical formula with sites denoted by empty square brackets.

• general_formula [string]

General chemical formula with sites denoted by square brackets filled with a comma-separated list of species

• n_sites [integer]

Number of sites in the solution. Should be the same for all endmembers.

• sites [list of lists of strings]

A list of species for each site in the solution.

• site_names [list of strings]

A list of species_site pairs in the solution, where each distinct site is given by a unique uppercase letter e.g. [‘Mg_A’, ‘Fe_A’, ‘Al_A’, ‘Al_B’, ‘Si_B’].

• n_occupancies [integer]

Sum of the number of possible species on each of the sites in the solution. Example: A binary solution [[A][B],[B][C1/2D1/2]] would have n_occupancies = 5, with two possible species on Site 1 and three on Site 2.

• site_multiplicities [2D array of floats]

A 1D array for each endmember in the solution, containing the multiplicities of each site per formula unit. To simplify computations later, the multiplicities are repeated for each species on each site, so the shape of this attribute is (n_endmembers, n_site_species).

• endmember_occupancies [2d array of floats]

A 1D array for each endmember in the solution, containing the fraction of atoms of each species on each site.

• endmember_noccupancies [2d array of floats]

A 1D array for each endmember in the solution, containing the number of atoms of each species on each site per mole of endmember.

burnman.utils.chemistry.site_occupancies_to_strings(site_species_names, site_multiplicities, endmember_occupancies)[source]#

Converts a list of endmember site occupancies into a list of string representations of those occupancies.

Parameters:
• site_species_names (2D list of strings) – A list of list of strings, giving the names of the species which reside on each site. List of sites, each of which contains a list of the species occupying each site.

• site_multiplicities (1D or 2D numpy array of floats) – List of floats giving the multiplicity of each site. If 2D, must have the same shape as endmember_occupancies. If 1D, must be either the same length as the number of sites, or the same length as site_species_names (with an implied repetition of the same number for each species on a given site).

• endmember_occupancies (2D numpy array of floats) – A list of site-species occupancies for each endmember. The first dimension loops over the endmembers, and the second dimension loops over the site-species occupancies for that endmember. The total number and order of occupancies must be the same as the strings in site_species_names.

Returns:

A list of strings in standard burnman format. For example, [Mg]3[Al]2 would correspond to the classic two-site pyrope garnet.

Return type:

list of strings

burnman.utils.chemistry.compositional_array(formulae)[source]#
Parameters:

formulae (list of dicts) – List of chemical formulae

Returns:

Array of endmember formulae and a list of elements.

Return type:

2D numpy.array of floats and a list of strs

burnman.utils.chemistry.ordered_compositional_array(formulae, elements)[source]#
Parameters:

formulae (list of dicts) – List of chemical formulae

:param elements : List of elements :type elements: list of strings

Returns:

Array of endmember formulae

Return type:

2D array of floats

burnman.utils.chemistry.formula_to_string(formula)[source]#
Parameters:

formula (dict or Counter) – Chemical formula

Returns:

A formula string, with element order as given in the list IUPAC_element_order. If one or more keys in the dictionary are not one of the elements in the periodic table, then they are added at the end of the string.

Return type:

str

burnman.utils.chemistry.sort_element_list_to_IUPAC_order(element_list)[source]#

:param element_list : List of elements. :type element_list: list

Returns:

List of elements sorted into IUPAC order

Return type:

list

burnman.utils.chemistry.convert_fractions(composite, phase_fractions, input_type, output_type)[source]#

Takes a composite with a set of user defined molar, volume or mass fractions (which do not have to be the fractions currently associated with the composite) and converts the fractions to molar, mass or volume.

Conversions to and from mass require a molar mass to be defined for all phases. Conversions to and from volume require set_state to have been called for the composite.

Parameters:
• composite (Composite) – Composite for which fractions are to be defined.

• phase_fractions (list of floats) – List of input phase fractions (of type input_type).

• input_type (str) – Input fraction type. One of ‘molar’, ‘mass’ or ‘volume’.

• output_type (str) – Output fraction type. One of ‘molar’, ‘mass’ or ‘volume’.

Returns:

List of output phase fractions (of type output_type)

Return type:

list of floats

burnman.utils.chemistry.reaction_matrix_as_strings(reaction_matrix, compound_names)[source]#

Returns a list of string representations of all the reactions in reaction_matrix.

Parameters:
• reaction_matrix (2D numpy array) – Matrix of stoichiometric amounts of each compound j in reaction i.

• compound_names (list of strings) – List of compound names.

Returns:

List of strings corresponding to each reaction.

Return type:

list of strings

burnman.tools.chemistry.fugacity(standard_material, assemblage)[source]#

Calculates the fugacity of a standard material in another assemblage.

Note

set_method and set_state should already have been used on both assemblages.

Parameters:
• standard_material – Standard material for which to calculate the fugacity. The material must have a formula as a dictionary parameter.

• assemblage (burnman.Composite) – Assemblage for which to calculate the fugacity.

Returns:

Value of the fugacity of the component with respect to the standard material.

Return type:

float

burnman.tools.chemistry.relative_fugacity(component_formula, assemblage, reference_assemblage)[source]#

Calculates the fugacity of a chemical component in one assemblage relative to another one.

Note

set_method and set_state should already have been used on both assemblages.

Parameters:
• component_formula (dictionary) – Chemical formula for which to compute the relative fugacity.

• assemblage (burnman.Composite) – Assemblage for which to calculate the fugacity.

• reference_assemblage (burnman.Composite) – Reference assemblage against which to measure the fugacity.

Returns:

Value of the fugacity of the component in the assemblage with respect to the reference_assemblage.

Return type:

float

burnman.tools.chemistry.equilibrium_pressure(minerals, stoichiometry, temperature, pressure_initial_guess=100000.0)[source]#

Given a list of minerals, their reaction stoichiometries and a temperature of interest, compute the equilibrium pressure of the reaction.

Parameters:
• minerals (list of burnman.Mineral) – List of minerals involved in the reaction.

• stoichiometry (list of floats) – Reaction stoichiometry for the minerals provided. Reactants and products should have the opposite signs [mol].

• temperature (float) – Temperature of interest [K].

• pressure_initial_guess (float) – Initial pressure guess [Pa].

Returns:

The equilibrium pressure of the reaction [Pa].

Return type:

float

burnman.tools.chemistry.equilibrium_temperature(minerals, stoichiometry, pressure, temperature_initial_guess=1000.0)[source]#

Given a list of minerals, their reaction stoichiometries and a pressure of interest, compute the equilibrium temperature of the reaction.

Parameters:
• minerals (list of burnman.Mineral) – List of minerals involved in the reaction.

• stoichiometry (list of floats) – Reaction stoichiometry for the minerals provided. Reactants and products should have the opposite signs [mol].

• pressure (float) – Pressure of interest [Pa].

• temperature_initial_guess (float) – Initial temperature guess [K].

Returns:

The equilibrium temperature of the reaction [K].

Return type:

float

burnman.tools.chemistry.invariant_point(minerals_r1, stoichiometry_r1, minerals_r2, stoichiometry_r2, pressure_temperature_initial_guess=[1000000000.0, 1000.0])[source]#

Given a list of minerals, their reaction stoichiometries and a pressure of interest, compute the equilibrium temperature of the reaction.

Parameters:
• minerals (list of burnman.Mineral) – List of minerals involved in the reaction.

• stoichiometry (list of floats) – Reaction stoichiometry for the minerals provided. Reactants and products should have the opposite signs [mol].

• pressure (float) – Pressure of interest [Pa].

• temperature_initial_guess (float) – Initial temperature guess [K].

Returns:

The equilibrium temperature of the reaction [K].

Return type:

float

burnman.tools.chemistry.hugoniot(mineral, P_ref, T_ref, pressures, reference_mineral=None)[source]#

Calculates the temperatures (and volumes) along a Hugoniot as a function of pressure according to the Hugoniot equation U2-U1 = 0.5*(p2 - p1)(V1 - V2) where U and V are the internal energies and volumes (mass or molar) and U = F + TS

Parameters:
• mineral (burnman.Mineral) – Mineral for which the Hugoniot is to be calculated.

• P_ref (float) – Reference pressure [Pa]

• T_ref (float) – Reference temperature [K]

• pressures (numpy.array of floats) – Set of pressures [Pa] for which the Hugoniot temperature and volume should be calculated.

• reference_mineral (burnman.Mineral) – Mineral which is stable at the reference conditions Provides an alternative U_0 and V_0 when the reference mineral transforms to the mineral of interest at some (unspecified) pressure.

Returns:

The Hugoniot temperatures and volumes at the given pressures.

Return type:

tuple of numpy.arrays

burnman.tools.chemistry.reactions_from_stoichiometric_matrix(stoichiometric_matrix)[source]#

Returns a list of all the balanced reactions between compounds of fixed chemical composition. Includes both the forward and reverse reactions (so there will always be an even number of reactions).

Parameters:

stoichiometric_matrix (2D numpy array) – An array of the stoichiometric (molar) amounts of component j in compound i.

Returns:

An array of the stoichiometric (molar) amounts of compound j in reaction i.

Return type:

2D numpy array

burnman.tools.chemistry.reactions_from_formulae(formulae, compound_names, return_strings=True)[source]#

Returns a list of all the balanced reactions between compounds of fixed chemical composition. Includes both the forward and reverse reactions (so there will always be an even number of reactions).

Parameters:
• formulae (list of dictionaries or list of strings) – List of the chemical formulae, either as strings or as a list of dictionaries of elements.

• compound_names (list of strings) – List of the compound names in the formula list.

• return_strings (bool) – Whether to return the reactions as strings or array.

Returns:

Either a 2D array of the stoichiometric (molar) amounts of compound j in reaction i, or a list of strings. The parameter compound_names is only used if strings are requested.

Return type:

2D numpy array or list of strings

# Equilibrium Thermodynamics#

burnman.tools.equilibration.equilibrate(composition, assemblage, equality_constraints, free_compositional_vectors=[], tol=0.001, store_iterates=False, store_assemblage=True, max_iterations=100.0, verbose=False)[source]#

A function that finds the thermodynamic equilibrium state of an assemblage subject to given equality constraints by solving a set of nonlinear equations related to the chemical potentials and other state variables of the system.

The user chooses an assemblage (e.g. olivine, garnet and orthopyroxene) and $$2+n_c$$ equality constraints, where $$n_c$$ is the number of bulk compositional degrees of freedom. The equilibrate function attempts to find the remaining unknowns that satisfy those constraints.

There are a number of equality constraints implemented in burnman. These are given as a list of lists. Each constraint should have the form: [<constraint type>, <constraint>], where <constraint type> is one of ‘P’, ‘T’, ‘S’, ‘V’, ‘X’, ‘PT_ellipse’, ‘phase_fraction’, or ‘phase_composition’. The format of the <constraint> object depends on the constraint type:

• P: float or numpy.array of

pressures [Pa]

• T: float or numpy.array of

temperatures [K]

• S: float or numpy.array of

entropies [J/K]

• V: float or numpy.array of

volumes [m:math:^3]

• PT_ellipse: list of two floats or numpy.arrays, where the equality

satifies the equation norm(([P, T] - arr[0])/arr[1]) = 1

• phase_fraction: tuple with the form (<phase> <fraction(s)>),

where <phase> is one of the phase objects in the assemblage and <fraction(s)> is a float or numpy.array corresponding to the desired phase fractions.

• phase_composition: tuple with the form (<phase> <constraint>),

where <phase> is one of the phase objects in the assemblage and <constraint> has the form (site_names, n, d, v), where $$(nx)/(dx) = v$$, n and d are constant vectors of site coefficients, and v is a float or numpy.array. For example, a constraint of the form ([Mg_A, Fe_A], [1., 0.], [1., 1.], [0.5]) would correspond to equal amounts Mg and Fe on the A site (Mg_A / (Mg_A + Fe_A) = 0.5).

• X: list of a numpy.array and a float or numpy.array,

where the equality satisfies the linear equation arr[0] x = arr[1]. The first numpy.array is fixed, and the second is to be looped over by the equilibrate function. This is a generic compositional equality constraint.

Parameters:
• composition (dict) – The bulk composition that the assemblage must satisfy.

• assemblage (burnman.Composite) – The assemblage to be equilibrated.

• equality_constraints (list of list) – The list of equality constraints. See above for valid formats.

• free_compositional_vectors (list of dict) – A list of dictionaries containing the compositional freedom of the solution. For example, if the list contains the vector {‘Mg’: 1., ‘Fe’: -1}, that implies that the bulk composition is equal to composition + $$a$$ (n_Mg - n_Fe), where the value of $$a$$ is to be determined by the solve. Vector given in atomic (molar) units of elements.

• tol (float) – The tolerance for the nonlinear solver.

• store_iterates (bool) – Whether to store the parameter values for each iteration in each solution object.

• store_assemblage (bool) – Whether to store a copy of the assemblage object in each solution object.

• max_iterations (int) – The maximum number of iterations for the nonlinear solver.

• verbose (bool) – Whether to print output updating the user on the status of equilibration.

Returns:

Solver solution object (or a list, or a 2D list of solution objects) created by burnman.optimize.nonlinear_solvers.damped_newton_solve(), and a namedtuple object created by burnman.tools.equilibration.get_equilibration_parameters(). See documentation of these functions for the return types. If store_assemblage is True, each solution object also has an attribute called assemblage, which contains a copy of the input assemblage with the equilibrated properties. So, for a 2D grid of solution objects, one could call either sols[0][1].x[0] or sols[0][1].assemblage.pressure to get the pressure.

Return type:

tuple

burnman.tools.equilibration.get_equilibration_parameters(assemblage, composition, free_compositional_vectors)[source]#

Builds a named tuple containing the parameter names and various other parameters needed by the equilibrium solve.

Parameters:
• assemblage (burnman.Composite) – The assemblage to be equilibrated.

• composition (dict) – The bulk composition for the equilibrium problem.

• free_compositional_vectors (list of dictionaries) – The bulk compositional degrees of freedom for the equilibrium problem.

Returns:

A tuple with attributes n_parameters (the number of parameters for the current equilibrium problem) and phase_amount_indices (the indices of the parameters that correspond to phase amounts).

Return type:

namedtuple