BurnMan 1.1.0: a thermodynamic and geophysics toolkit for the Earth and planetary sciences¶
BurnMan is a Python library for computing the thermodynamic and thermoelastic properties of geological materials from simple mineral endmembers to complex multilayered planetary interiors.
BurnMan is released under the GNU GPL v2 or newer. It relies heavily on numpy, scipy, and matplotlib.
Homepage: https://geodynamics.github.io/burnman/
Documentation: http://burnman.readthedocs.io
Source code: https://github.com/geodynamics/burnman
If you haven’t yet installed BurnMan, you can go straight to Requirements for detailed instructions. After that, you might want to try out The BurnMan Tutorial or the other Examples. Finally, and most importantly, have fun!
Citing BurnMan¶
If you use BurnMan in your work, we ask that you cite the following publications:
Myhill, R., Cottaar, S., Heister, T., Rose, I., and Unterborn, C. (2021): BurnMan v1.0.1 [Software]. Computational Infrastructure for Geodynamics. Zenodo. (https://doi.org/10.5281/zenodo.7080174)
Cottaar S., Heister, T., Rose, I., and Unterborn, C., 2014, BurnMan: A lower mantle mineral physics toolkit, Geochemistry, Geophysics, and Geosystems, 15(4), 11641179 (https://doi.org/10.1002/2013GC005122)
Contributing to BurnMan¶
If you would like to contribute bug fixes, new functions or new modules to the existing codebase, please contact us at info@burnman.org or make a pull request at https://github.com/geodynamics/burnman.
BurnMan also includes a contrib directory that contains python and ipython scripts used to reproduce published results. We welcome the submission of new contributions to this directory. As with the contribution of code, please contact us at info@burnman.org or make a pull request at https://github.com/geodynamics/burnman.
Acknowledgement and Support¶
This project was initiated at, and followup research support was received through, Cooperative Institute of Deep Earth Research, CIDER (NSF FESD grant 1135452) – see www.deepearth.org
We thank all the members of the CIDER Mg/Si team for their input: Valentina Magni, Yu Huang, JiaChao Liu, Marc Hirschmann, and Barbara Romanowicz. We also thank Lars Stixrude for providing benchmarking calculations and Zack Geballe, Motohiko Murakami, Bill McDonough, Quentin Williams, Wendy Panero, and Wolfgang Bangerth for helpful discussions.
We thank CIG (www.geodynamics.org) for support and accepting our donation of BurnMan as an official project.
Introducing BurnMan 1.1.0¶
Overview¶
BurnMan is an open source mineral physics and seismological toolkit written in Python which can enable a user to calculate (or fit) the physical and chemical properties of endmember minerals, fluids/melts, solutions, and composite assemblages.
Properties which BurnMan can calculate include:
the thermodynamic free energies, allowing phase equilibrium calculations, endmember activities, chemical potentials and oxygen (and other) fugacities.
entropy, enabling the user to calculate isentropes for a given assemblage.
volume, to allow the user to create density profiles.
seismic velocities, including VoigtReussHill and HashinStrikman bounds and averages.
The toolkit itself comes with a large set of classes and functions which are designed to allow the user to easily combine mineral physics with geophysics, and geodynamics. The features of BurnMan include:
the full codebase, which includes implementations of many static and thermal equations of state (including Vinet, Birch Murnaghan, MieDebyeGrueneisen, Modified Tait), and solution models (ideal, symmetric, asymmetric, subregular).
popular endmember and solution datasets already coded into burnmanusable format (including [HollandPowell11], [SLB05] and [SLB11])
Optimal least squares fitting routines for multivariate data with (potentially correlated) errors in pressure and temperature. As an example, such functions can be used to simultaneously fit volumes, seismic velocities and enthalpies.
a “Planet” class, which selfconsistently calculates gravity profiles, mass, moment of inertia of planets given the chemical and temperature structure of a planet
published geotherms
a tutorial on the basic use of BurnMan
a large collection of annotated examples
a set of highlevel functions which create files readable by seismological and geodynamic software, including: Mineos [MWF11], AxiSEM [NissenMeyervanDrielStahler+14] and ASPECT
an extensive suite of unit tests to ensure code functions as intended
a series of benchmarks comparing BurnMan output with published data
a directory containing usercontributed code from published papers
BurnMan makes extensive use of SciPy, NumPy and SymPy which are widely used Python libraries for scientific computation. Matplotlib is used to display results and produce publication quality figures. The computations are consistently formulated in terms of SI units.
The code documentation including class and function descriptions can be found online at http://burnman.readthedocs.io.
This software has been designed to allow the enduser a great deal of freedom to do whatever calculations they may wish and to add their own modules. The underlying Python classes have been designed to make new endmember, solution and composite models easy to read and create. We have endeavoured to provide examples and benchmarks which cover the most popular uses of the software, some of which are included in the figure below. This list is certainly not exhaustive, and we will definitely have missed interesting applications. We will be very happy to accept contributions in form of corrections, examples, or new features.
Structure¶
Requirements¶
Python 3.7+
Python modules: NumPy, SciPy, SymPy, Matplotlib
Optional modules: cvxpy, pycddlib
Installation¶
Installation of BurnMan is mostly platform independent. As long as you know how to use a terminal, the process should be straightforward. The following instructions should help, but let us know if you have any problems.
Dependencies¶
First, make sure you have a sufficiently recent version of python installed on your machine (see above for the latest requirements). To check your version of python, type the following in a terminal:
python –version
If your version is not recent enough, visit https://www.python.org/ to find out how to install a newer version.
Once you have checked your version of python, you should make sure you have installed the python module pip. We will use this module to install BurnMan. If you don’t have it already, you can install it by opening a terminal window and typing:
python m ensurepip –upgrade
Mac users will also need to install Xcode, which can be found in the MacStore.
Stable version¶
If you are only interested in using BurnMan (rather than developing the software), and you aren’t interested in any of the latest changes, you can install the stable version by typing the following into a terminal window:
python m pip install burnman
This method of installation does not give you easy access to all the examples, or the test suite. These can be found in the latest release package which can be downloaded from https://github.com/geodynamics/burnman/releases.
Development version¶
If you want to install the development version of BurnMan (with all the latest features), you will first need to download the source code. The best way to do this is by using git (a version control system). To install git, follow the instructions at https://gitscm.com/downloads.
Then, using a terminal, navigate to the directory into which you want to clone the BurnMan repository, and type
(If you don’t want to use git, you can download the current master branch from https://github.com/geodynamics/burnman/archive/master.zip.)
Once the repository is cloned, navigate to the toplevel directory by typing cd burnman in the terminal, and then install BurnMan, either in static mode: python m pip install . or in development mode (if you want to develop or modify the code): python m pip install e ..
Checking that the installation worked¶
To check that the installation has worked, you can run the test suite (./test.sh). This takes a few minutes to run.
A more basic check that BurnMan is installed is to navigate to the Burnman examples directory and type:
python example_beginner.py
If figures show up, BurnMan has been installed.
Citing BurnMan¶
If you use BurnMan in your work, we ask that you cite the following publications:
Myhill, R., Cottaar, S., Heister, T., Rose, I., and Unterborn, C. (2022): BurnMan v1.1.0 [Software]. Computational Infrastructure for Geodynamics. Zenodo. (https://doi.org/10.5281/zenodo.7080174)
Cottaar S., Heister, T., Rose, I., and Unterborn, C., 2014, BurnMan: A lower mantle mineral physics toolkit, Geochemistry, Geophysics, and Geosystems, 15(4), 11641179 (https://doi.org/10.1002/2013GC005122)
Contributing to BurnMan¶
If you would like to contribute bug fixes, new functions or new modules to the existing codebase, please contact us at info@burnman.org or make a pull request at https://github.com/geodynamics/burnman.
BurnMan also includes a contrib directory that contains python and ipython scripts used to reproduce published results. We welcome the submission of new contributions to this directory. As with the contribution of code, please contact us at info@burnman.org or make a pull request at https://github.com/geodynamics/burnman.
Acknowledgement and Support¶
This project was initiated at, and followup research support was received through, Cooperative Institute of Deep Earth Research, CIDER (NSF FESD grant 1135452) – see www.deepearth.org
We thank all the members of the CIDER Mg/Si team for their input: Valentina Magni, Yu Huang, JiaChao Liu, Marc Hirschmann, and Barbara Romanowicz. We also thank Lars Stixrude for providing benchmarking calculations and Zack Geballe, Motohiko Murakami, Bill McDonough, Quentin Williams, Wendy Panero, and Wolfgang Bangerth for helpful discussions.
We thank CIG (www.geodynamics.org) for support and accepting our donation of BurnMan as an official project.
Mathematical Background¶
Here is a bit of background on the methods used to calculate thermoelastic and thermodynamic properties in BurnMan. More detail can be found in the cited papers.
Endmember Properties¶
Calculating Thermoelastic Properties¶
To calculate the bulk (\(K\)) modulus, shear modulus (\(G\)) and density (\(\rho\)) of a material at a given pressure (\(P\)) and temperature (\(T\)), optionally defined by a geotherm) and determine the seismic velocities (\(V_S, V_P, V_\Phi\)), one uses an Equation of State (EoS). Currently the following EoSs are supported in BurnMan:
BirchMurnaghan finitestrain EoS (excludes temperature effects, [Poi91]),
BirchMurnaghan finitestrain EoS with a MieGrüneisenDebye thermal correction, as formulated by [SLB05].
BirchMurnaghan finitestrain EoS with a MieGrüneisenDebye thermal correction, as formulated by [MBR+07].
Modified Tait EoS (excludes temperature effects, [HuangChow74]),
Modified Tait EoS with a pseudoEinstein model for thermal corrections, as formulated by [HollandPowell11].
CompensatedRedlichKwong for fluids, as formulated by [HP91].
To calculate these thermoelastic parameters, the EoS requires the user to input the pressure, temperature, and the phases and their molar fractions. These inputs and outputs are further discussed in User input.
BirchMurnaghan (isothermal)¶
The BirchMurnaghan equation is an isothermal Eulerian finitestrain EoS relating pressure and volume. The negative finitestrain (or compression) is defined as
where \(V\) is the volume at a given pressure and \(V_0\) is the volume at a reference state (\(P = 10^5\) Pa , \(T\) = 300 K). The pressure and elastic moduli are derived from a thirdorder Taylor expansion of Helmholtz free energy in \(f\) and evaluating the appropriate volume and strain derivatives (e.g., [Poi91]). For an isotropic material one obtains for the pressure, isothermal bulk modulus, and shear modulus:
Here \(K_0\) and \(G_0\) are the reference bulk modulus and shear modulus and \(K_0^\prime\) and \({G}^\prime_{0}\) are the derivative of the respective moduli with respect to pressure.
BurnMan has the option to use the secondorder expansion for shear modulus by dropping the \(f^2\) terms in these equations (as is sometimes done for experimental fits or EoS modeling).
Modified Tait (isothermal)¶
The Modified Tait equation of state was developed by [HuangChow74]. It has the considerable benefit of allowing volume to be expressed as a function of pressure. It performs very well to pressures and temperatures relevant to the deep Earth [HollandPowell11].
MieGrüneisenDebye (thermal correction to BirchMurnaghan)¶
The Debye model for the Helmholtz free energy can be written as follows [MBR+07]
where \(\theta\) is the Debye temperature and \(\gamma\) is the Grüneisen parameter.
Using thermodynamic relations we can derive equations for the thermal pressure and bulk modulus
The thermal shear correction used in BurnMan was developed by [HamaSuito98]
The total pressure, bulk and shear moduli can be calculated from the following sums
This equation of state is substantially the same as that in SLB2005 (see below). The primary differences are in the thermal correction to the shear modulus and in the volume dependences of the Debye temperature and the Gruneisen parameter.
HP2011 (thermal correction to Modified Tait)¶
The thermal pressure can be incorporated into the Modified Tait equation of state, replacing \(P\) with \(P\left(P_{\textrm{th}}  P_{\textrm{th0}}\right)\) in Equation (5) [HollandPowell11]. Thermal pressure is calculated using a MieGrüneisen equation of state and an Einstein model for heat capacity, even though the Einstein model is not actually used for the heat capacity when calculating the enthalpy and entropy (see following section).
\(\Theta\) is the Einstein temperature of the crystal in Kelvin, approximated for a substance \(i\) with \(n_i\) atoms in the unit formula and a molar entropy \(S_i\) using the empirical formula
SLB2005 (for solids, thermal)¶
Thermal corrections for pressure, and isothermal bulk modulus and shear modulus are derived from the MieGrüneisenDebye EoS with the quasiharmonic approximation. Here we adopt the formalism of [SLB05] where these corrections are added to equations (2)–(4):
The \(\Delta\) refers to the difference in the relevant quantity from the reference temperature (300 K). \(\gamma\) is the Grüneisen parameter, \(q\) is the logarithmic volume derivative of the Grüneisen parameter, \(\eta_{S}\) is the shear strain derivative of the Grüneisen parameter, \(C_V\) is the heat capacity at constant volume, and \(\mathcal{U}\) is the internal energy at temperature \(T\). \(C_V\) and \(\mathcal{U}\) are calculated using the Debye model for vibrational energy of a lattice. These quantities are calculated as follows:
where \(\theta\) is the Debye temperature of the mineral, \(\nu\) is the frequency of vibrational modes for the mineral, \(n\) is the number of atoms per formula unit (e.g. 2 for periclase, 5 for perovskite), and \(R\) is the gas constant. Under the approximation that the vibrational frequencies behave the same under strain, we may identify \(\nu/\nu_0 = \theta/\theta_0\). The quantities \(\gamma_0\), \(\eta_{S0}\) \(q_0\), and \(\theta_0\) are the experimentally determined values for those parameters at the reference state.
Due to the fact that a planetary mantle is rarely isothermal along a geotherm, it is more appropriate to use the adiabatic bulk modulus \(K_S\) instead of \(K_T\), which is calculated using
where \(\alpha\) is the coefficient of thermal expansion:
There is no difference between the isothermal and adiabatic shear moduli for an isotropic solid. All together this makes an eleven parameter EoS model, which is summarized in the Table below. For more details on the EoS, we refer readers to [SLB05].
User Input 
Symbol 
Definition 
Units 

V_0 
\(V_{0}\) 

m \(^{3}\) mol \(^{1}\) 
K_0 
\(K_{0}\) 
Isothermal bulk modulus at P=10^5 Pa, T = 300 K 
Pa 
Kprime_0 
\(K^\prime_0\) 
Pressure derivative of \(K_{0}\) 

G_0 
\(G_{0}\) 
Shear modulus at P = \(10^5\) Pa, T = 300 K 
Pa 
Gprime_0 
\(G^\prime_0\) 
Pressure derivative of \(G_{0}\) 

molar_mass 
\(\mu\) 
mass per mole formula unit 
kg \(\mathrm{mol}^{1}\) 
n 
n 
number of atoms per formula unit 

Debye_0 
\(\theta_{0}\) 
Debye Temperature 
K 
grueneisen_0 
\(\gamma_{0}\) 
Grüneisen parameter at P = \(10^5\) Pa, T = 300 K 

q0 
\(q_{0}\) 
Logarithmic volume derivative of the Grüneisen parameter 

eta_s_0 
\(\eta_{S0}\) 
Shear strain derivative of the Grüneisen parameter 
This equation of state is substantially the same as that of the MieGruneisenDebye (see above). The primary differences are in the thermal correction to the shear modulus and in the volume dependences of the Debye temperature and the Gruneisen parameter.
CompensatedRedlichKwong (for fluids, thermal)¶
The CORK equation of state [HP91] is a simple virialtype extension to the modified RedlichKwong (MRK) equation of state. It was designed to compensate for the tendency of the MRK equation of state to overestimate volumes at high pressures and accommodate the volume behaviour of coexisting gas and liquid phases along the saturation curve.
Calculating Thermodynamic Properties¶
So far, we have concentrated on the thermoelastic properties of minerals. There are, however, additional thermodynamic properties which are required to describe the thermal properties such as the energy, entropy and heat capacity. These properties are related by the following expressions:
where \(P\) is the pressure, \(T\) is the temperature and \(\mathcal{E}\), \(\mathcal{F}\), \(\mathcal{H}\), \(\mathcal{S}\) and \(V\) are the molar internal energy, Helmholtz free energy, enthalpy, entropy and volume respectively.
HP2011¶
The heat capacity at one bar is given by an empirical polynomial fit to experimental data
The entropy at high pressure and temperature can be calculated by differentiating the expression for \(\mathcal{G}\) with respect to temperature
Finally, the enthalpy at high pressure and temperature can be calculated
SLB2005¶
The Debye model yields the Helmholtz free energy and entropy due to lattice vibrations
Property modifiers¶
The thermodynamic models above consider the effects of strain and quasiharmonic lattice vibrations on the free energies of minerals at given temperatures and pressures. There are a number of additional processes, such as isochemical orderdisorder and magnetic effects which also contribute to the total free energy of a phase. Corrections for these additional processes can be applied in a number of different ways. Burnman currently includes implementations of the following:
Linear excesses (useful for DQF modifications for [HollandPowell11])
Tricritical Landau model (two formulations)
BraggWilliams model
Magnetic excesses
In all cases, the excess Gibbs free energy \(\mathcal{G}\) and first and second partial derivatives with respect to pressure and temperature are calculated. The thermodynamic properties of each phase are then modified in a consistent manner; specifically:
Subscripts \(_o\) and \(_m\) indicate original properties and modifiers respectively. Importantly, this allows us to stack modifications such as multiple Landau transitions in a simple and straightforward manner. In the burnman code, we add property modifiers as an attribute to each mineral as a list. For example:
from burnman.minerals import SLB_2011
stv = SLB_2011.stishovite()
stv.property_modifiers = [
['landau',
{'Tc_0': 4250.0, 'S_D': 0.012, 'V_D': 1e09}]]
['linear',
{'delta_E': 1.e3, 'delta_S': 0., 'delta_V': 0.}]]
Each modifier is a list with two elements, first the name of the modifier type, and second a dictionary with the required parameters for that model. A list of parameters for each model is given in the following sections.
Linear excesses (linear)¶
A simple linear correction in pressure and temperature. Parameters are ‘delta_E’, ‘delta_S’ and ‘delta_V’.
Tricritical Landau model (landau)¶
Applies a tricritical Landau correction to the properties of an endmember which undergoes a displacive phase transition. These transitions are not associated with an activation energy, and therefore they occur rapidly compared with seismic wave propagation. Parameters are ‘Tc_0’, ‘S_D’ and ‘V_D’.
This correction follows [Putnis92], and is done relative to the completely ordered state (at 0 K). It therefore differs in implementation from both [SLB11] and [HollandPowell11], who compute properties relative to the completely disordered state and standard states respectively. The current implementation is preferred, as the excess entropy (and heat capacity) terms are equal to zero at 0 K.
If the temperature is above the critical temperature, Q (the order parameter) is equal to zero, and the Gibbs free energy is simply that of the disordered phase:
If temperature is below the critical temperature, Q is between 0 and 1. The gibbs free energy can be described thus:
Tricritical Landau model (landau_hp)¶
Applies a tricritical Landau correction similar to that described above. However, this implementation follows [HollandPowell11], who compute properties relative to the standard state. Parameters are ‘P_0’, ‘T_0’, ‘Tc_0’, ‘S_D’ and ‘V_D’.
It is worth noting that the correction described by [HollandPowell11] has been incorrectly used throughout the geological literature, particularly in studies involving magnetite (which includes studies comparing oxygen fugacities to the FMQ buffer (due to an incorrect calculation of the properties of magnetite). Note that even if the implementation is correct, it still allows the order parameter Q to be greater than one, which is physically impossible.
We include this implementation in order to reproduce the dataset of [HollandPowell11]. If you are creating your own minerals, we recommend using the standard implementation.
If the temperature is above the critical temperature, Q (the order parameter) is equal to zero. Otherwise
The second derivatives of the Gibbs free energy are only nonzero if the order parameter exceeds zero. Then
BraggWilliams model (bragg_williams)¶
The BraggWilliams model is a symmetric solution model between endmembers with an excess configurational entropy term determined by the specifics of orderdisorder in the mineral, multiplied by some empirical factor. Expressions for the excess Gibbs free energy can be found in [HP96]. Parameters are ‘deltaH’, ‘deltaV’, ‘Wh’, ‘Wv’, ‘n’ and ‘factor’.
Magnetic model (magnetic_chs)¶
This model approximates the excess energy due to magnetic ordering. It was originally described in [CHS87]. The expressions used by BurnMan can be found in [Sun91]. Parameters are ‘structural_parameter’, ‘curie_temperature’[2] (zero pressure value and pressure dependence) and ‘magnetic_moment’[2] (zero pressure value and pressure dependence).
Calculating Solution Properties¶
Many phases (whether minerals, melts, fluids or gases) can exist over a finite region of composition space. These spaces are bounded by endmembers (which may themselves not be stable), and each phase can then be described as a solution of those endmembers. In a solid solution, different elements substitute for one another on distinct crystallographic sites in the structure. For example, low pressure silicate garnets have two distinct sites on which mixing takes place; a dodecahedral site (of which there are three per unit cell on an eightcation basis) and octahedral site (of which there are two per unit cell). A third tetrahedral cation site (three per unit cell) is usually assumed to be occupied solely by silicon, and therefore can be ignored in solution calculations. The chemical formula of many low pressure garnets exist within the solution:
We typically calculate solution properties by appropriate differentiation of the Gibbs energy, where
Implemented models¶
Ideal solutions¶
A solution is not simply a mechanical mixture of its constituent endmembers. The mixing of different species results in an excess configurational entropy \(S\). In BraggWilliamstype solutions, the entropy only depends on the amounts of species on sites, and the site multiplicities.
where \(s\) is a site in the lattice, \(c\) are the species mixing on site \(s\). \(x_c^s\) is the absolute number of species \(c\) on site \(s\) in the lattice; it is calculated by multiplying the proportion of the species on the site by the multiplicity of the site per formula unit and the number of moles of formula units.
Solutions where this configurational entropy is the only deviation from a mechanical mixture are termed ideal, because the enthalpy of mixing is zero. In BurnMan, the multiplicities of each site are allowed to vary linearly between endmembers. This is known as a Temkin model [Tem45].
Symmetric solutions¶
Many real phases are not well approximated as ideal solutions. Deviations are the result of elastic and chemical interactions between ions with different physical and chemical characteristics. Regular (symmetric) solution models account for the simplest form of deviations from ideality by incorporating terms describing excess enthalpies, nonconfigurational entropies and volumes relative to the ideal solution model. These excess terms have the matrix form [DPWH07]
where \(p\) is a vector of molar fractions of each of the \(n\) endmembers and \(W\) is a strictly uppertriangular matrix of interaction terms between endmembers. Excesses within binary systems (\(i\)\(j\)) have a quadratic form and a maximum of \(W_{ij}/4\) halfway between the two endmembers.
Asymmetric solutions¶
Some solutions exhibit asymmetric excess terms. These can be accounted for with an asymmetric solution [DPWH07]
\(\alpha\) is a vector of “van Laar parameters” governing asymmetry in the excess properties.
The \(w_{ij}\) terms are a set of interaction terms between endmembers \(i\) and \(j\). If all the \(\alpha\) terms are equal to unity, a nonzero \(w\) yields an excess with a quadratic form and a maximum of \(w/4\) halfway between the two endmembers.
Subregular solutions¶
An alternative way to create asymmetric solution models is to expand each binary term as a cubic expression [HW89]. In this case,
Note the similarity with the symmetric solution model, the primary difference being that there are two interaction terms for each binary and also additional ternary terms.
Arbitrary solutions¶
BurnMan also allows the user to define their own excess Gibbs energy function for a solution model.
Thermodynamic and thermoelastic properties¶
From the preceeding equations, we can define the thermodynamic potentials of solutions:
We can also define the derivatives of volume with respect to pressure and temperature
Making the approximation that the excess entropy has no temperature dependence
Including orderdisorder¶
Orderdisorder can be treated trivially with solutions. The only difference between mixing between ordered and disordered endmembers is that disordered endmembers have a nonzero configurational entropy, which must be accounted for when calculating the excess entropy within a solution.
Including spin transitions¶
The regular solution formalism should provide an elegant way to model spin transitions in phases such as periclase and bridgmanite. High and low spin iron can be treated as different elements, providing distinct endmembers and an excess configurational entropy. Further excess terms can be added as necessary.
Calculating Multiphase Composite Properties¶
Averaging schemes¶
After the thermoelastic parameters (\(K_S\), \(G\), \(\rho\)) of each phase are determined at each pressure and/or temperature step, these values must be combined to determine the seismic velocity of a multiphase assemblage. We define the volume fraction of the individual minerals in an assemblage:
where \(V_i\) and \(n_i\) are the molar volume and the molar fractions of the \(i\) th individual phase, and \(V\) is the total molar volume of the assemblage:
The density of the multiphase assemblage is then
where \(\rho_i\) is the density and \(\mu_i\) is the molar mass of the \(i\) th phase.
Unlike density and volume, there is no straightforward way to average the bulk and shear moduli of a multiphase rock, as it depends on the specific distribution and orientation of the constituent minerals. BurnMan allows several schemes for averaging the elastic moduli: the Voigt and Reuss bounds, the HashinShtrikman bounds, the VoigtReussHill average, and the HashinShtrikman average [WDOConnell76].
The Voigt average, assuming constant strain across all phases, is defined as
where \(X_i\) is the bulk or shear modulus for the \(i\) th phase. The Reuss average, assuming constant stress across all phases, is defined as
The VoigtReussHill average is the arithmetic mean of Voigt and Reuss bounds:
The HashinShtrikman bounds make an additional assumption that the distribution of the phases is statistically isotropic and are usually much narrower than the Voigt and Reuss bounds [WDOConnell76]. This may be a poor assumption in regions of Earth with high anisotropy, such as the lowermost mantle, however these bounds are more physically motivated than the commonlyused VoigtReussHill average. In most instances, the VoigtReussHill average and the arithmetic mean of the HashinShtrikman bounds are quite similar with the pure arithmetic mean (linear averaging) being well outside of both.
It is worth noting that each of the above bounding methods are derived from mechanical models of a linear elastic composite. It is thus only appropriate to apply them to elastic moduli, and not to other thermoelastic properties, such as wave speeds or density.
Computing seismic velocities¶
Once the moduli for the multiphase assemblage are computed, the compressional (\(P\)), shear (\(S\)) and bulk sound (\(\Phi\)) velocities are then result from the equations:
To correctly compare to observed seismic velocities one needs to correct for the frequency sensitivity of attenuation. Moduli parameters are obtained from experiments that are done at high frequencies (MHzGHz) compared to seismic frequencies (mHzHz). The frequency sensitivity of attenuation causes slightly lower velocities for seismic waves than they would be for high frequency waves. In BurnMan one can correct the calculated acoustic velocity values to those for long period seismic tomography [MA81]:
Similar to [MBR+07], we use a \(\beta\) value of 0.3, which falls in the range of values of \(0.2\) to \(0.4\) proposed for the lower mantle (e.g. [KS90]). The correction is implemented for \(Q\) values of PREM for the lower mantle. As \(Q_S\) is smaller than \(Q_P\), the correction is more significant for S waves. In both cases, though, the correction is minor compared to, for example, uncertainties in the temperature (corrections) and mineral physical parameters. More involved models of relaxation mechanisms can be implemented, but lead to the inclusion of more poorly constrained parameters, [MB07]. While attenuation can be ignored in many applications [TVV01], it might play a significant role in explaining strong variations in seismic velocities in the lowermost mantle [DGD+12].
Thermodynamic Equilibration¶
For a composite with fixed phases at a given pressure, temperature and composition, equilibrium is reached when the following relationships are satisfied:
where \(\mu_j\) are the chemical potentials of all of the endmembers in all of the phases, and \(R_{ij}\) is an independent set of balanced reactions between endmembers.
It is generally true that at a fixed composition, one can choose two equilibrium constraints (such as fixed temperature, pressure, entropy, volume, phase proportion or some composition constraint) and solve for the remaining unknowns. In BurnMan, this can be achieved using the equilibrate function (see Equilibrium Thermodynamics).
User input¶
Mineralogical composition¶
A number of predefined minerals are included in the mineral library and users can create their own. The library includes wrapper functions to include a transition from the highspin mineral to the lowspin mineral [LSMM13] or to combine minerals for a given iron number.
Standard minerals – The ‘standard’ mineral format includes a list of parameters given in the above table. Each mineral includes a suggested EoS with which the mineral parameters are derived. For some minerals the parameters for the thermal corrections are not yet measured or calculated, and therefore the corrections can not be applied. An occasional mineral will not have a measured or calculated shear moduli, and therefore can only be used to compute densities and bulk sound velocities. The mineral library is subdivided by citation. BurnMan includes the option to produce a LaTeX; table of the mineral parameters used. BurnMan can be easily setup to incorporate uncertainties for these parameters.
Minerals with a spin transition – A standard mineral for the high spin and low spin must be defined separately. These minerals are “wrapped,” so as to switch from the high spin to high spin mineral at a give pressure. While not realistic, for the sake of simplicity, the spin transitions are considered to be sharp at a given pressure.
Minerals depending on Fe partitioning – The wrapper function can partition iron, for example between ferropericlase, fp, and perovskite, pv. It requires the input of the iron mol fraction with regards to Mg, \(X_\mathrm{fp}\) and \(X_\mathrm{pv}\), which then defines the chemistry of an MgFe solid solution according to (\(\mathrm{Mg}_{1X_{\mathrm{Fe}}^{\mathrm{fp}}},\mathrm{Fe}_{X_{\mathrm{Fe}}^{\mathrm{fp}}})\mathrm{O}\) or \((\mathrm{Mg}_{1X_{\mathrm{Fe}}^{\mathrm{pv}}},\mathrm{Fe}_{X_{\mathrm{Fe}}^{\mathrm{pv}}})\mathrm{SiO_3}\). The iron mol fractions can be set to be constant or varying with P and T as needed. Alternatively one can calculate the iron mol fraction from the distribution coefficient \(K_D\) defined as
We adopt the formalism of [NFR12] choosing a reference distribution coefficient \(K_{D0}\) and standard state volume change (\(\Delta \upsilon^{0}\)) for the FeMg exchange between perovskite and ferropericlase
where \(R\) is the gas constant and \(P_0\) the reference pressure. As a default, we adopt the average \(\Delta \upsilon^{0}\) of [NFR12] of \(2\cdot10^{7}\) \(m^3 mol^{1}\) and suggest using their \({K_D}_0\) value of \(0.5\).
The multiphase mixture of these minerals can be built by the user in three ways:
1. Molar fractions of an arbitrary number of predefined minerals, for example mixing standard minerals mg_perovskite (\(\mathrm{MgSiO_3}\)), fe_perovskite (\(\mathrm{FeSiO_3}\)), periclase (\(\mathrm{MgO}\)) and wüstite (\(\mathrm{FeO}\)).
2. A twophase mixture with constant or (\(P,T\)) varying Fe partitioning using the minerals that include Fedependency, for example mixing \(\mathrm{(Mg,Fe)SiO_3}\) and \(\mathrm{(Mg,Fe)O}\) with a predefined distribution coefficient.
3. Weight percents (wt%) of (Mg, Si, Fe) and distribution coefficient (includes (P,T)dependent Fe partitioning). This calculation assumes that each element is completely oxidized into its corresponding oxide mineral (\(\mathrm{MgO}\), \(\mathrm{FeO}\), \(\mathrm{SiO_2}\)) and then combined to form ironbearing perovskite and ferropericlase taking into account some Fe partition coefficient.
Geotherm¶
Unlike the pressure, the temperature of the lower mantle is relatively unconstrained. As elsewhere, BurnMan provides a number of builtin geotherms, as well as the ability to use userdefined temperaturedepth relationships. A geotherm in BurnMan is an object that returns temperature as a function of pressure. Alternatively, the user could ignore the geothermal and compute elastic velocities for a range of temperatures at any give pressure.
Currently, we include geotherms published by [BS81] and [And82]. Alternatively one can use an adiabatic gradient defined by the thermoelastic properties of a given mineralogical model. For a homogeneous material, the adiabatic temperature profile is given by integrating the ordinary differential equation (ODE)
This equation can be extended to multiphase composite using the first law of thermodynamics to arrive at
where the subscripts correspond to the \(i\) th phase, \(C_P\) is the heat capacity at constant pressure of a phase, and the other symbols are as defined above. Integrating this ODE requires a choice in anchor temperature (\(T_0\)) at the top of the lower mantle (or including this as a parameter in an inversion). As the adiabatic geotherm is dependent on the thermoelastic parameters at high pressures and temperatures, it is dependent on the equation of state used.
Seismic Models¶
BurnMan allows for direct visual and quantitative comparison with seismic velocity models. Various ways of plotting can be found in the examples. Quantitative misfits between two profiles include an L2norm and a chisquared misfit, but user defined norms can be implemented. A seismic model in BurnMan is an object that provides pressure, density, and seismic velocities (\(V_P, V_\Phi, V_S\)) as a function of depth.
To compare to seismically constrained profiles, BurnMan provides the 1D seismic velocity model PREM [DA81]. One can choose to evaluate \(V_P, V_\Phi, V_S, \rho, K_S\) and/or \(G\). The user can input their own seismic profile, an example of which is included using AK135 [KEB95].
Besides standardized 1D radial profiles, one can also compare to regionalized average profiles for the lower mantle. This option accommodates the observation that the lowermost mantle can be clustered into two regions, a ‘slow’ region, which represents the socalled Large Low Shear Velocity Provinces, and ‘fast’ region, the continuous surrounding region where slabs might subduct [LCDR12]. This clustering as well as the averaging of the 1D model occurs over five tomographic S wave velocity models (SAW24B16: [MegninR00]; HMSLS: [HMSL08]; S362ANI: [KED08]; GyPSuM: [SFBG10]; S40RTS: [RDvHW11]). The strongest deviations from PREM occur in the lowermost 1000 km. Using the ‘fast’ and ‘slow’ S wave velocity profiles is therefore most important when interpreting the lowermost mantle. Suggestion of compositional variation between these regions comes from seismology [HW12, TRCT05] as well as geochemistry [DCT12, JCK+10]. Based on thermochemical convection models, [SDG11] also show that averaging profiles in thermal boundary layers may cause problems for seismic interpretation.
We additionally apply cluster analysis to and provide models for P wave velocity based on two tomographic models (MITP08: [LvdH08]; GyPSuM: [SMJM12]). The clustering results correlate well with the fast and slow regions for S wave velocities; this could well be due to the fact that the initial model for the P wave velocity models is scaled from S wave tomographic velocity models. Additionally, the variations in P wave velocities are a lot smaller than for S waves. For this reason using these adapted models is most important when comparing the S wave velocities.
While interpreting lateral variations of seismic velocity in terms of composition and temperature is a major goal [MCD+12, TDRY04], to determine the bulk composition the current challenge appears to be concurrently fitting absolute P and S wave velocities and incorporate the significant uncertainties in mineral physical parameters).
The BurnMan Tutorial¶
This tutorial introduces users to the main classes and some of the important functions implemented in BurnMan. The tutorial is split into several parts:
The BurnMan Tutorial
Part 1: Material Classes¶
This file is part of BurnMan  a thermoelastic and thermodynamic toolkit for the Earth and Planetary Sciences
Copyright (C) 2012  2021 by the BurnMan team, released under the GNU GPL v2 or later.
Introduction¶
This ipython notebook is the first in a series designed to introduce new users to the code structure and functionalities present in BurnMan.
Demonstrates
burnman.Mineral: Equations of state, property modification schemes, initialization, interrogating properties at given pressure and temperature.
burnman.CombinedMineral: Initialization (otherwise similar to mineral).
burnman.Solution: Formulations (ideal, asymmetric, subregular), initialization, interrogating properties at given pressure, temperature and composition.
burnman.Composite: Initialization, interrogating properties at given pressure, temperature, phase proportions and using different seismic averaging schemes.
Everything in BurnMan and in this tutorial is defined in SI units.
Acknowledgments¶
The authors are grateful to M. Ghiorso for useful discussions. R.M. would like to thank B. Watterson for the concept of Planet Zog (see Part 3: Layers and Planets).
This project was initiated at, and followup research support was received through, CIDER (NSF FESD grant 1135452). The development of BurnMan has been supported by the Computational Infrastructure for Geodynamics initiative (CIG), through the Science and Technologies Funding Council (U.K.) under Award No. ST/R001332/1 and through the Natural Environment Research Council (U.K.) Large Grant MCsquared (Award No. NE/T012633/1). The authors have also received support from the University of California  Davis.
Getting started with BurnMan¶
Our first task is to import BurnMan. If you haven’t yet installed the current version, you can do this by typing pip install e .
from the toplevel directory of the repository. Alternatively, you could just uncomment (remove the leading #
from) the first line in the following code block. The second line in the code block imports the BurnMan module.
[1]:
#!pip install q e ..
import burnman
Warning: No module named 'cdd'. For full functionality of BurnMan, please install pycddlib.
Types of BurnMan Material objects¶
The BurnMan package allows users to define and use objects that represent different kinds of Materials. The most important classes of Material are named Mineral, Solution and Composite. In the following subsections, we show how users can create objects of each type, set their state (pressure and temperature) and composition, and interrogate them for their material properties.
Mineral objects¶
Mineral objects are the building blocks for more complex objects in BurnMan. These objects are intended to represent minerals (or melts, or fluids) of fixed composition, with a welldefined equation of state that defines the relationship between the current state (pressure and temperature) of the mineral and its thermodynamic potentials and derivatives (such as Gibbs energy, volume and entropy).
Mineral objects are initialized with a “params” dictionary containing all of the parameters required by the desired equation of state and an optional “property modifiers” argument. Here we initialize a generic Mineral, just to show the general structure:
[2]:
mineral_object = burnman.Mineral(params={}, property_modifiers=[[]])
The required keys in the parameters dictionary depends on the equation of state, described in the following section.
Equations of state¶
BurnMan identifies the desired equation of state by checking the value of the string parameter “equation_of_state” in the parameters dictionary that is passed as an argument to the Mineral class. BurnMan currently contains implementations of the following static equations of state (i.e., equations of state with no temperature dependence):  “bm2”, “bm3” and “bm4”: BirchMurnaghan (2nd, 3rd and 4th order)  “mt”: Modified Tait  “morse”: Morse potential  “vinet”: Vinet (originally the Rydberg equation of state)  “rkprime”: Reciprocal Kprime
And also the following thermal equations of state:  “mgd2” and “mgd3”: MieDebyeGrueneisen equation of state (second and third order expansions for the shear modulus)  “slb2” and “slb3”: Stixrude and LithgowBertelloni (2011; second and third order expansions for the shear modulus)  “hp_tmt”: Thermal Modified Tait (Holland and Powell; 2011)  “dks_s”: de Koker and Stixrude (2013; solids)  “dks_l”: de Koker and Stixrude (2013; liquids)  “cork”: Compensated RedlichKwong equation of state  “aa”: Anderson and Ahrens (1998)  “brosh_calphad”: Brosh et al. (2007)
Each equation of state assumes the presence of a different set of keys in the “params” dictionary. These keys are checked and validated on initialization. Two important parameters for most mineral objects are the chemical formula of the mineral and its molar mass, which can be calculated using the functions “dictionarize_formula” and “formula_mass”.
[3]:
from burnman.tools.chemistry import dictionarize_formula, formula_mass
forsterite_formula = dictionarize_formula('Mg2SiO4')
molar_mass = formula_mass(forsterite_formula)
print(f'Formula in dictionary form: {forsterite_formula}')
print(f'Molar mass: {molar_mass:.5f} kg')
Formula in dictionary form: {'Mg': 2.0, 'Si': 1.0, 'O': 4.0}
Molar mass: 0.14069 kg
Below, we demonstrate the creation of a forsterite object for the Stixrude and LithgowBertelloni (2011) equation of state which uses a 3rd order expansion for the shear modulus (equation_of_state = ‘slb3’).
[4]:
forsterite_params = {'name': 'Forsterite',
'formula': forsterite_formula,
'equation_of_state': 'slb3',
'F_0': 2055403.0,
'V_0': 4.3603e05,
'K_0': 1.279555e+11,
'Kprime_0': 4.21796,
'Debye_0': 809.1703,
'grueneisen_0': 0.99282,
'q_0': 2.10672,
'G_0': 81.6e9,
'Gprime_0': 1.46257,
'eta_s_0': 2.29972,
'n': sum(forsterite_formula.values()),
'molar_mass': molar_mass}
forsterite = burnman.Mineral(params=forsterite_params)
print(forsterite.params)
{'name': 'Forsterite', 'formula': {'Mg': 2.0, 'Si': 1.0, 'O': 4.0}, 'equation_of_state': 'slb3', 'F_0': 2055403.0, 'V_0': 4.3603e05, 'K_0': 127955500000.0, 'Kprime_0': 4.21796, 'Debye_0': 809.1703, 'grueneisen_0': 0.99282, 'q_0': 2.10672, 'G_0': 81600000000.0, 'Gprime_0': 1.46257, 'eta_s_0': 2.29972, 'n': 7.0, 'molar_mass': 0.14069310000000002, 'T_0': 300.0, 'E_0': 0.0, 'P_0': 0.0}
Property modifiers¶
Thermodynamic models of minerals can include modifications to the properties predicted by the underlying equation of state. Often, this is done to approximate physical processes which were neglected during development of the equation of state.
BurnMan allows users to apply modifications to the Gibbs energy \(\mathcal{G}(P, T)\) via a “property_modifiers” attribute. This attribute takes the form of a list of different modifiers, which are themselves each composed of a list consisting of an identifying string and a dictionary containing the required parameters. The following modifiers are currently implemented in BurnMan:  “linear”: Linear in pressure and temperature (i.e. \(\Delta \mathcal{G} = \Delta \mathcal{E}  T \Delta S + P \Delta V\))  “landau”: A Landau orderdisorder transition (Putnis, 1992)  “landau_hp”: A Landau orderdisorder transition (Holland and Powell, 2011)  “bragg_williams”: A BraggWilliams orderdisorder transition  “magnetic_chs”: A magnetic orderdisorder transition
The following code block implements the Landau transition for hematite as given in the Holland et al. (2018) database (ds6.33).
[5]:
# Dictionary with parameters for equation of state.
hematite_params = {'name': 'hematite',
'formula': {'Fe': 2.0, 'O': 3.0},
'equation_of_state': 'hp_tmt',
'H_0': 825420.0,
'S_0': 87.4,
'V_0': 3.027e05,
'Cp': [163.9, 0.0, 2257200.0, 657.6],
'a_0': 2.79e05,
'K_0': 223000e6,
'Kprime_0': 4.04,
'Kdprime_0': 1.8e11,
'n': 5.0,
'molar_mass': 0.1596882}
# Dictionary for Landau transition which modifieds the Gibbs energy
hematite_property_modifiers = [['landau_hp', {'P_0': 1.e5,
'T_0': 298.15,
'Tc_0': 955.0,
'S_D': 15.6,
'V_D': 0.0}]]
# Initialise mineral with a property modifier
hematite = burnman.Mineral(params=hematite_params,
property_modifiers=hematite_property_modifiers)
The CombinedMineral class¶
In petrology, we are often interested in phases for which we have little experimental or theoretical information. One common example is when we want to approximate the properties of an ordered phase relative to its disordered counterparts. In many cases, a reasonable procedure is to make a mechanical mixture of the known phases, such that they are the correct composition of the unknown phase, and then apply a linear correction to the Gibbs energy of the phase (i.e. \(\Delta \mathcal{G} = \Delta \mathcal{E}  T \Delta S + P \Delta V\)). In BurnMan, we do this using the CombinedMineral class. In the lines below, we create an ordered orthopyroxene with composition MgFeSi\(_2\)O\(_6\) from a 50:50 mixture of enstatite and ferrosilite. We make this compound 6 kJ/mol more stable than the mechanical mixture.
[6]:
from burnman import CombinedMineral
from burnman.minerals import HP_2011_ds62
fe_mg_orthopyroxene = CombinedMineral(name='ordered ferroenstatite',
mineral_list=[HP_2011_ds62.en(),
HP_2011_ds62.fs()],
molar_amounts=[0.5, 0.5],
free_energy_adjustment=[6.e3, 0., 0.])
print(fe_mg_orthopyroxene.formula)
Counter({'O': 6.0, 'Si': 2.0, 'Mg': 1.0, 'Fe': 1.0})
Implemented datasets¶
BurnMan already includes several published mineral datasets. These can be found in the “minerals” subdirectory of BurnMan, and include the Stixrude and LithgowBertelloni (2011) dataset, Holland et al. (2013) and Jennings et al. (2015) datasets. The script misc/table_mineral_library.py will list all minerals included for each equation of state.
Mineral objects can be initialised from these datasets using commands like the following:
[7]:
from burnman import minerals
fo_SLB = minerals.SLB_2011.forsterite()
fo_HP = minerals.HP_2011_ds62.fo()
The “()” at the end of each call indicate that these commands initialize an object from a class. This is done so that users can create multiple distinct objects with the same material parameters but potentially different states.
Interrogating the properties of a Mineral object at a given pressure and temperature¶
Once an object has been initialised, the user can set its state (pressure in Pa and temperature in K):
[8]:
fo_SLB.set_state(1.e9, 1000)
The thermodynamic and seismic properties of the object at those conditions can then by returned as attributes of the object (returned in SI units):
[9]:
from burnman.tools.chemistry import formula_to_string
print(f'{fo_SLB.name} (SLB2011)')
print(f'Formula: {formula_to_string(fo_SLB.formula)}')
print(f'Gibbs: {fo_SLB.gibbs:.3e} J/mol')
print(f'S: {fo_SLB.molar_entropy:.3e} J/K/mol')
print(f'V_p: {fo_SLB.v_p:.3e} m/s')
Forsterite (SLB2011)
Formula: Mg2SiO4
Gibbs: 2.148e+06 J/mol
S: 2.747e+02 J/K/mol
V_p: 8.308e+03 m/s
It is common for users to want to return properties over a range of states. BurnMan makes this easy through the “evaluate” method of the Material class. In the following lines, we investigate the properties of quartz as defined by Stixrude and Lithgow Bertelloni (2011):
[10]:
# Necessary imports for creating the P, T arrays and plotting
import numpy as np
import matplotlib.pyplot as plt
# Temperatures and pressures at which to evaluate properties
temperatures = np.linspace(300., 1300., 1001)
pressures = 1.e5*np.ones_like(temperatures)
# Desired properties
properties = ['molar_heat_capacity_p', 'v_phi']
# Initialize a quartz object and print the parameters dictionary
# and property modifiers
qtz = minerals.SLB_2011.quartz()
print(qtz.params)
print(qtz.property_modifiers)
# Return arrays containing the desired properties
isobaric_heat_capacities, bulk_sound_velocities = qtz.evaluate(properties, pressures, temperatures)
# Plot the properties
fig = plt.figure(figsize=(8, 4))
fig.suptitle('Heat capacity and bulk sound velocity across alpha>beta transition in quartz')
ax = [fig.add_subplot(1, 2, i) for i in range(1, 3)]
ax[0].plot(temperatures, isobaric_heat_capacities)
ax[1].plot(temperatures, bulk_sound_velocities/1000.)
for i in range(2):
ax[i].set_xlabel('Temperature (K)')
ax[0].set_ylabel('Isobaric heat capacity (J/K/mol)')
ax[1].set_ylabel('Bulk sound velocity (km/s)')
ax[0].set_ylim(50., 150.)
ax[1].set_ylim(1.5, 5.)
fig.tight_layout()
fig.subplots_adjust(top=0.88)
{'name': 'Quartz', 'formula': {'Si': 1.0, 'O': 2.0}, 'equation_of_state': 'slb3', 'F_0': 858853.4, 'V_0': 2.367003e05, 'K_0': 49547430000.0, 'Kprime_0': 4.33155, 'Debye_0': 816.3307, 'grueneisen_0': 0.00296, 'q_0': 1.0, 'G_0': 44856170000.0, 'Gprime_0': 0.95315, 'eta_s_0': 2.36469, 'n': 3.0, 'molar_mass': 0.0600843, 'T_0': 300.0, 'E_0': 0.0, 'P_0': 0.0}
[['landau', {'Tc_0': 847.0, 'S_D': 5.164, 'V_D': 1.222e06}]]
In the figure produced by the code above, one can see the effect of the alpha>beta transition in quartz (also known as the quartz inversion). This is a displacive transition from trigonal to hexagonal symmetry that can be modelled (as done here) by a Landautype transition.
To give users confidence that BurnMan is outputting accurate derivative properties, we include a tool that checks the thermodynamic selfconsistency of the equation of state:
[11]:
from burnman.tools.eos import check_eos_consistency
assert(check_eos_consistency(qtz, 1.e9, 1000., verbose=False))
Solution objects¶
The Solution class in BurnMan allows users to define materials that have crystal chemical sites that can be occupied by multiple species. The species occupancies on the sites can be varied by changing the proportions of a number of ‘’independent’’ endmembers.
Currently, all the models implemented in BurnMan are of “BraggWilliams” type (Bragg and Williams, 1930s). This means that the thermodynamic properties of solutions are uniquely defined by the average occupancies of species on each site (Myhill and Connolly, 2021). These models are therefore unable to accurately account for local interactions that give rise to shortrange order, but BraggWilliams models with multiple sites can provide a firstorder approximation to the energetic effects of longrange order. From a microscopic standpoint, the distinction is an artificial one, since longrange order arises from local interactions (e.g. Bethe 1934, 1935). However, a detailed treatment of orderdisorder requires either complex models (Kikuchi, 1951) or abinitio simulations that have their own dedicated software packages (e.g. VASP).
Model formalisms¶
The simplest solution model implemented in BurnMan is the ideal solution model. In this solution, the excess Gibbs energy of mixing is purely configurational in nature. Mathematically:
\(\Delta \mathcal{G}^{\text{mix}} = RT \sum_i m_i \sum_j p_{ij} \ln p_{ij}\)
where \(R\) is the gas constant (in J/K/mol), \(T\) is the temperature (in K), \(m_i\) is the multiplicity of site \(i\) and \(p_{ij}\) is the proportion of species \(j\) on site \(i\).
The following code initializes a binary ideal solution model between pyrope and almandine. The solution_type argument defines the other arguments that are reuired by the model. The endmembers argument must be a list of lists containing all the endmembers of the solution. The first item in each inner list should contain a Mineral object corresponding to the endmember of interest. The second item should be a chemical formula. The exchange sites in this formula are denoted by square brackets, followed by the multiplicities \(m_i\) (optional if \(m_i = 1\)). Everything not contained within an “[…]m” block is ignored.
[12]:
from burnman import Solution, minerals
ideal_garnet = Solution(name = 'Ideal pyropealmandine garnet',
solution_type = 'ideal',
endmembers = [[minerals.HP_2011_ds62.py(),
'[Mg]3[Al]2Si3O12'],
[minerals.HP_2011_ds62.alm(),
'[Fe]3[Al]2Si3O12']],
molar_fractions = [0.5, 0.5])
Each species on each site must begin with a capital letter, but the string does not have to correspond to an element. Partial occupancies are allowed. For example, [Al0.5Fef0.5] would be valid for an endmember with both \(\text{Al}^{3+}\) and \(\text{Fe}^{3+}\) on the same site. [Vac] would be valid for a vacancy, whereas [v] would not. The molar_fractions argument in the initialization above is optional.
As of BurnMan 1.1, ideality does not require that the multiplicity of each site is fixed between endmembers. This allows the implementation of Temkintype models (Temkin, 1945), which are used to model melts. For example, here is the ideal part of the melt model proposed by Holland et al. (2018), restricted to the CMAS system (nonideality and changes in endmember energies have been ignored):
Another solution model formalism implemented in BurnMan is the asymmetric model (Holland and Powell, 2003). This model is an extension of the regular solution model (Wohl, 1946), with the extension following van Laar (1906).
The general idea is that there are interaction energies \(W_{ij}\) associated with each binary in the solution model. Those binary interactions are modified by a set of \(\alpha_k\) parameters, such that:
\(\Delta \mathcal{G}^{\text{mix}} = \Delta \mathcal{G}^{\text{mix,ideal}} + \Delta \mathcal{G}^{\text{mix,nonideal}}\)
and
\(\Delta \mathcal{G}^{\text{mix,nonideal}} = (\sum_l \alpha_l p_l) (\sum_i^{n1} \sum_{j=i+1}^n \phi_i \phi_j \frac{2 w_{ij}}{\alpha_i + \alpha_j}\))
where
\(\phi_i = \frac{\alpha_i p_i}{\sum_k \alpha_k p_k}\)
\(\Delta \mathcal{G}^{\text{mix,ideal}}\) is caalculated using the same expressions as for the ideal model. The following code block initialises an asymmetric garnet model in the system CFMASO (taken from the Holland and Powell, 2011 dataset; ds6.2).
[13]:
g2 = Solution(name='asymmetric garnet (ThermoCalc ds6.2)',
solution_type='asymmetric',
endmembers = [[minerals.HP_2011_ds62.py(), '[Mg]3[Al]2Si3O12'],
[minerals.HP_2011_ds62.alm(), '[Fe]3[Al]2Si3O12'],
[minerals.HP_2011_ds62.gr(), '[Ca]3[Al]2Si3O12'],
[minerals.HP_2011_ds62.andr(), '[Ca]3[Fe]2Si3O12']],
alphas = [1.0, 1.0, 2.7, 2.7],
energy_interaction = [[2.5e3, 31.e3, 53.2e3],
[5.e3, 37.24e3],
[2.e3]])
In the case that all \(\alpha_i\) are equal to each other, the asymmetric model becomes a symmetric model. BurnMan allows users to specify this type of model by setting solution_type=’symmetric’. In this case, the alphas argument does not need to be specified.
BurnMan also contains an implementation of the subregular model (Helffrich and Wood, 1989). In this model, the excess nonideal Gibbs energy is expressed as a cubic polynomial in endmember proportions:
\(G^* = p_i G^{\text{mbr}}_i + p_i p_j W_{ij}^{\text{binary}} (1 + p_j  p_i)/2 + p_i p_j p_k W_{ijk}^{\text{ternary}}\)
where the binary terms must be equal to zero when \(i=j\) and the ternary terms can only be nonzero for \(i<j<k\).
[14]:
# Parameters from Ganguly et al. (1996), converted to SI units
Wh_1bar = [[[2117., 695.], [9834., 21627.], [12083., 12083.]],
[[6773., 873.], [539., 539.]],
[[0., 0.]]]
Wv = [[[0.07e5, 0.], [0.058e5, 0.012e5], [0.04e5, 0.03e5]],
[[0.03e5, 0.], [0.04e5, 0.01e5]],
[[0., 0.]]]
Ws = [[[0., 0.], [5.78, 5.78], [7.67, 7.67]],
[[1.69, 1.69], [0., 0.]],
[[0., 0.]]]
# We now convert from 1 bar enthalpy, entropy and volume (1 cation basis)
# to energy, entropy and volume (3 cation basis)
mult = lambda x, n: [[[v*n for v in i] for i in j] for j in x]
add = lambda x, y: [[[x[i][j][k] + y[i][j][k] for k in range(len(x[i][j]))]
for j in range(len(x[i]))] for i in range(len(x))]
energy_interaction = add(mult(Wh_1bar, 3.), mult(Wv, 3.e5))
volume_interaction = mult(Wv, 3.)
entropy_interaction = mult(Ws, 3.)
g3 = Solution(name='CFMnMAS garnet (Ganguly et al., 1996)',
solution_type='subregular',
endmembers=[[minerals.HP_2011_ds62.py(), '[Mg]3[Al]2Si3O12'],
[minerals.HP_2011_ds62.alm(), '[Fe]3[Al]2Si3O12'],
[minerals.HP_2011_ds62.gr(), '[Ca]3[Al]2Si3O12'],
[minerals.HP_2011_ds62.spss(), '[Mn]3[Al]2Si3O12']],
energy_interaction=energy_interaction,
volume_interaction=volume_interaction,
entropy_interaction=entropy_interaction,
energy_ternary_terms = [[0, 1, 2, 0.e3]],
volume_ternary_terms = [[0, 1, 2, 0.e6]],
entropy_ternary_terms = [[0, 1, 2, 0.]])
For the model above, the ternary terms are all equal to zero, but in general, the structure of the ternary_terms input is a list of lists, where the inner list corresponds to the three endmember indices \(i<j<k\) followed by the value of the interaction term.
Implemented datasets¶
As for the endmembers, BurnMan already contains solution models from multiple datasets. One of these is the Stixrude and LithgowBertelloni (2011) dataset. In the codeblock below, we initialize an object using the threeendmember (FMAS) bridgmanite solution from that publication.
[15]:
bdg = minerals.SLB_2011.mg_fe_bridgmanite()
bdg.endmembers
[15]:
[[<burnman.minerals.SLB_2011.mg_perovskite at 0x7fc22e861350>, '[Mg][Si]O3'],
[<burnman.minerals.SLB_2011.fe_perovskite at 0x7fc230c16ed0>, '[Fe][Si]O3'],
[<burnman.minerals.SLB_2011.al_perovskite at 0x7fc26c2fc9d0>, '[Al][Al]O3']]
Interrogating the properties of a Solution object at a given composition, pressure and temperature¶
Unlike Minerals, Solutions can vary in composition. Before interrogating properties, it is therefore necessary to set both the state (using the method “set_state”) and the composition (using the method “set_composition”).
[16]:
bdg.set_composition([0.8, 0.1, 0.1])
bdg.set_state(30.e9, 2000.)
print(f'Bridgmanite composition: {dict(bdg.formula)}')
print(f'Bridgmanite volume: {bdg.V:.3e} m')
print(f'Bridgmanite endmember partial Gibbs energies:')
for pG in bdg.partial_gibbs:
print(f' {pG/1e3:.3e} kJ/mol')
Bridgmanite composition: {'Mg': 0.8, 'Si': 0.9, 'O': 3.0, 'Fe': 0.1, 'Al': 0.2}
Bridgmanite volume: 2.315e05 m
Bridgmanite endmember partial Gibbs energies:
9.690e+02 kJ/mol
6.868e+02 kJ/mol
1.143e+03 kJ/mol
The evaluate method can be used to output properties of solutions at constant composition:
[17]:
pressures = np.linspace(30.e9, 130.e9, 101)
temperatures = 2000.*np.ones_like(pressures)
densities, p_wave_velocities = bdg.evaluate(['rho', 'v_p'], pressures, temperatures)
# The following lines do the plotting
fig = plt.figure(figsize=(8, 4))
ax = [fig.add_subplot(1, 2, i) for i in range(1, 3)]
ax[0].plot(pressures/1.e9, densities)
ax[1].plot(pressures/1.e9, p_wave_velocities/1000.)
for i in range(2):
ax[i].set_xlabel('Pressure (GPa)')
ax[0].set_ylabel('Densities (kg/m^3)')
ax[1].set_ylabel('P wave velocity (km/s)')
fig.tight_layout()
The “set_composition” method must be used every time the Solution composition is to be changed:
[18]:
import numpy as np
import matplotlib.pyplot as plt
comp = np.linspace(1e5, 1.1e5, 101)
bdg_excess_gibbs_400 = np.empty_like(comp)
bdg_excess_gibbs_800 = np.empty_like(comp)
bdg_excess_gibbs_1200 = np.empty_like(comp)
bdg_activities_400 = np.empty((101, 3))
bdg_activities_800 = np.empty((101, 3))
bdg_activities_1200 = np.empty((101, 3))
pressure = 1.e9
for i, c in enumerate(comp):
molar_fractions = [1.0  c, c, 0.]
bdg.set_composition(molar_fractions)
bdg.set_state(pressure, 400.)
bdg_excess_gibbs_400[i] = bdg.excess_gibbs
bdg_activities_400[i] = bdg.activities
bdg.set_state(pressure, 800.)
bdg_excess_gibbs_800[i] = bdg.excess_gibbs
bdg_activities_800[i] = bdg.activities
bdg.set_state(pressure, 1200.)
bdg_excess_gibbs_1200[i] = bdg.excess_gibbs
bdg_activities_1200[i] = bdg.activities
fig = plt.figure(figsize=(8, 4))
ax = [fig.add_subplot(1, 2, i) for i in range(1, 3)]
ax[0].plot(comp, bdg_excess_gibbs_400/1000., 'r', linewidth=1., label='400 K')
ax[0].plot(comp, bdg_excess_gibbs_800/1000., 'g', linewidth=1., label='800 K')
ax[0].plot(comp, bdg_excess_gibbs_1200/1000., 'b', linewidth=1., label='1200 K')
ax[1].plot(comp, bdg_activities_1200[:,0], 'b', linewidth=1., label='MgSiO$_3$')
ax[1].plot(comp, bdg_activities_1200[:,1], 'b:', linewidth=1., label='FeSiO$_3$')
ax[0].set_title("MgSiO$_3$FeSiO$_3$ bridgmanite join")
ax[0].set_ylabel("Excess gibbs free energy of solution (kJ/mol)")
ax[1].set_ylabel("Endmember activities")
for i in range(2):
ax[i].set_xlabel("Molar FeSiO$_3$ fraction")
ax[i].legend(loc='lower left')
fig.tight_layout()
plt.show()
Composite objects¶
The third major class derived from Material is the Composite class. This class is designed to represent collections of other Materials, which are typically objects of type Mineral and Solution, but can be any Material, including other Composites. A list of Materials is passed to the argument “phases” on initialization of a Composite object.
The properties of Composite materials can be interrogated in a similar manner to Solutions. To set the molar fractions of the material, either pass an array of fractions as the argument “fractions” on initialization, or use the “set_fractions” method at any time after initialization. As before, “set_state” is used to set the pressure and temperature of the object.
[19]:
from burnman import Composite
bdg = minerals.SLB_2011.mg_fe_bridgmanite()
fper = minerals.SLB_2011.ferropericlase()
bdg.set_composition([0.9, 0.1, 0.0])
fper.set_composition([0.8, 0.2])
assemblage = Composite(phases=[bdg, fper],
fractions=[0.5, 0.5],
fraction_type='molar',
name='rock')
pressure = 30.e9
temperature = 2000.
assemblage.set_fractions([0.5, 0.5])
assemblage.set_state(pressure, temperature)
print(f'Assemblage density at {pressure/1.e9} GPa and {temperature} K: {assemblage.density:.1f} kg/m^3')
Assemblage density at 30.0 GPa and 2000.0 K: 4491.0 kg/m^3
As for the Solution class, the properties of Composites at fixed phase fraction and phase compositions can be obtined using the evaluate method of the class. In the following code block, we evaluate the seismic properties of our 50:50 molar mixture of bridgmanite and ferropericlase from 30 to 130 GPa at 2000 K.
In this code block, we also demonstrate the ability for BurnMan to switch between different schemes for averaging seismic properties. Available schemes are “Voigt”, “Reuss”, “VoigtReussHill” (the standard scheme, which averages the Voigt and Reuss averages), “HashinShtrikmanLower” and “HashinShtrikmanUpper”. The figure shows the difference between the Reuss (lower) and Voigt (upper) bounds on the Pwave and Swave velocities.
[20]:
pressures = np.linspace(30.e9, 130.e9, 101)
temperatures = np.ones(101)*2000.
assemblage.set_averaging_scheme('VoigtReussHill')
densities, Vp_VRH, Vs_VRH = assemblage.evaluate(['rho', 'v_p', 'v_s'],
pressures, temperatures)
assemblage.set_averaging_scheme('Reuss')
Vp_R, Vs_R = assemblage.evaluate(['v_p', 'v_s'], pressures, temperatures)
assemblage.set_averaging_scheme('Voigt')
Vp_V, Vs_V = assemblage.evaluate(['v_p', 'v_s'], pressures, temperatures)
fig = plt.figure(figsize=(8, 4))
ax = [fig.add_subplot(1, 2, i) for i in range(1, 3)]
ax[0].plot(pressures/1.e9, densities)
ax[1].fill_between(pressures/1.e9, Vp_R/1.e3, Vp_V/1.e3, alpha=0.2)
ax[1].fill_between(pressures/1.e9, Vs_R/1.e3, Vs_V/1.e3, alpha=0.2)
ax[1].plot(pressures/1.e9, Vs_R/1.e3, color='orange', linewidth=0.5)
ax[1].plot(pressures/1.e9, Vs_V/1.e3, color='orange', linewidth=0.5, label='$V_S$')
ax[1].plot(pressures/1.e9, Vp_R/1.e3, color='blue', linewidth=0.5)
ax[1].plot(pressures/1.e9, Vp_V/1.e3, color='blue', linewidth=0.5, label='$V_P$')
for i in range(2):
ax[i].set_xlabel('Pressure (GPa)')
ax[0].set_ylabel('Density (kg/m$3$)')
ax[1].set_ylabel('Velocities (km/s)')
fig.tight_layout()
plt.show()
In the next part of the tutorial, we will look at BurnMan’s Composition class.
The BurnMan Tutorial
Part 2: The Composition Class¶
This file is part of BurnMan  a thermoelastic and thermodynamic toolkit for the Earth and Planetary Sciences
Copyright (C) 2012  2021 by the BurnMan team, released under the GNU GPL v2 or later.
Introduction¶
This ipython notebook is the second in a series designed to introduce new users to the code structure and functionalities present in BurnMan.
Demonstrates
burnman.Composition: Defining Composition objects, converting between molar, weight and atomic amounts, changing component bases. and modifying compositions.
Everything in BurnMan and in this tutorial is defined in SI units.
The Composition class¶
It is quite common in petrology to want to perform simple manipulations on chemical compositions. These manipulations might include:  converting between molar and weight percent of oxides or elements  changing from one compositional basis to another (e.g. ‘FeO’ and ‘Fe2O3’ to ‘Fe’ and ‘O’)  adding new chemical components to an existing composition in specific proportions with existing components.
These operations are easy to perform in Excel (for example), but errors are surprisingly common, and are even present in published literature. BurnMan’s Composition class is designed to make some of these common tasks easy and hopefully less error prone. Composition objects are initialised with a dictionary of component amounts (in any format), followed by a string that indicates whether that composition is given in “molar” amounts or “weight” (more technically mass, but weight is a more commonly used word in chemistry).
[1]:
from burnman import Composition
olivine_composition = Composition({'MgO': 1.8,
'FeO': 0.2,
'SiO2': 1.}, 'molar')
Warning: No module named 'cdd'. For full functionality of BurnMan, please install pycddlib.
After initialization, the “print” method can be used to directly print molar, weight or atomic amounts. Optional variables control the print precision and normalization of amounts.
[2]:
olivine_composition.print('molar', significant_figures=4,
normalization_component='SiO2', normalization_amount=1.)
olivine_composition.print('weight', significant_figures=4,
normalization_component='total', normalization_amount=1.)
olivine_composition.print('atomic', significant_figures=4,
normalization_component='total', normalization_amount=7.)
Molar composition
FeO: 0.2000
MgO: 1.8000
SiO2: 1.0000
Weight composition
FeO: 0.0977
MgO: 0.4935
SiO2: 0.4087
Atomic composition
Fe: 0.2000
Mg: 1.8000
O: 4.0000
Si: 1.0000
Let’s do something a little more complicated. When we’re making a starting mix for petrological experiments, we often have to add additional components. For example, we add iron as Fe2O3 even if we want a reduced oxide starting mix, because FeO is not a stable stoichiometric compound.
Here we show how to use BurnMan to create such mixes. In this case, let’s say we want to create a KLB1 starting mix (Takahashi, 1986). We know the weight proportions of the various oxides (including only components in the NCFMAS system):
[3]:
KLB1 = Composition({'SiO2': 44.48,
'Al2O3': 3.59,
'FeO': 8.10,
'MgO': 39.22,
'CaO': 3.44,
'Na2O': 0.30}, 'weight')
However, this composition is not the composition we wish to make in the lab. We need to make the following changes:  \(\text{CaO}\) and \(\text{Na}_2\text{O}\) should be added as \(\text{CaCO}_3\) and \(\text{Na}_2\text{CO}_3\).  \(\text{FeO}\) should be added as \(\text{Fe}_2\text{O}_3\)
First, we change the bulk composition to satisfy these requirements. The molar amounts of the existing components are stored in a dictionary “molar_composition”, and can be used to determine the amounts of CO2 and O to add to the bulk composition:
[4]:
CO2_molar = KLB1.molar_composition['CaO'] + KLB1.molar_composition['Na2O']
O_molar = KLB1.molar_composition['FeO']*0.5
KLB1.add_components(composition_dictionary = {'CO2': CO2_molar,
'O': O_molar},
unit_type = 'molar')
Then we can change the component set to the oxidised, carbonated compounds and print the desired starting compositions, for 2 g total mass:
[5]:
KLB1.change_component_set(['Na2CO3', 'CaCO3', 'Fe2O3', 'MgO', 'Al2O3', 'SiO2'])
KLB1.print('weight', significant_figures=4, normalization_amount=2.)
Weight composition
Al2O3: 0.0697
CaCO3: 0.1193
Fe2O3: 0.1749
MgO: 0.7620
Na2CO3: 0.0100
SiO2: 0.8642
And that’s it! The next tutorial will be on making Layer and Planet objects for planetary science applications.
The BurnMan Tutorial
Part 3: Layers and Planets¶
This file is part of BurnMan  a thermoelastic and thermodynamic toolkit for the Earth and Planetary Sciences
Copyright (C) 2012  2021 by the BurnMan team, released under the GNU GPL v2 or later.
Introduction¶
This ipython notebook is the third in a series designed to introduce new users to the code structure and functionalities present in BurnMan.
Demonstrates
burnman.Layer
burnman.Planet
Everything in BurnMan and in this tutorial is defined in SI units.
Building a planet¶
Planets are, to a good first approximation, layered like an onion. They typically have a core, a mantle and a crust, stratified according to density. Because temperature and pressure are both strongly correlated with depth, major phase transitions are (again, to first order) reasonably treated as being continuous and horizontal.
On Earth, these approximations have led to 1D models of properties, such as PREM (Dziewonski and Anderson, 1981) and AK135 (Kennett, Engdahl and Buland, 1995). These models can be used as a starting point for studies investigating the possible composition and temperature of Earth’s deep interior.
On other planets, we do not yet have data anywhere near as good as for Earth, and so interior structure is less wellknown. What we do have is gravity at the surface, mass of the planet, and moment of inertia. So the question is then, what might those things tell us about interior mineralogy?
BurnMan includes Layer and Planet classes to help investigate these questions.
The Layer class¶
The Layer class in BurnMan is designed to represent a spherical shell of a planet. That shell is made of a BurnMan Material object. The user can specify how the pressure and temperature evolve within the Layer.
The following code block creates a lower_mantle layer of mg_bridgmanite and periclase in 80:20 molar proportions, and assigns an adiabatic profile with a temperature of 1900 K at the top of the layer. The pressure is set to be selfconsistent; i.e. the gravity and pressure are adjusted to ensure that pressure is hydrostatic.
[1]:
import numpy as np
import matplotlib.pyplot as plt
import burnman
from burnman import Mineral, PerplexMaterial, Composite, Layer, Planet
from burnman import minerals
depths = np.linspace(2890e3, 670e3, 20)
rock = Composite([minerals.SLB_2011.mg_bridgmanite(),
minerals.SLB_2011.periclase()],
[0.8, 0.2])
lower_mantle = Layer(name='Lower Mantle', radii=6371.e3depths)
lower_mantle.set_material(rock)
lower_mantle.set_temperature_mode(temperature_mode='adiabatic',
temperature_top=1900.)
lower_mantle.set_pressure_mode(pressure_mode='selfconsistent',
pressure_top=23.8e9,
gravity_bottom=10.7)
# The "make" method does the calculations to make the pressure and gravity selfconsistent.
lower_mantle.make()
fig = plt.figure(figsize=(8, 8))
ax = [fig.add_subplot(2, 2, i) for i in range(1, 5)]
ax[0].plot(lower_mantle.pressure/1.e9, 6371.lower_mantle.radii/1.e3)
ax[1].plot(lower_mantle.temperature, 6371.lower_mantle.radii/1.e3)
ax[2].plot(lower_mantle.gravity, 6371.lower_mantle.radii/1.e3)
ax[3].plot(lower_mantle.bullen, 6371.lower_mantle.radii/1.e3)
for i in range(3):
ax[i].set_ylim(6371.lower_mantle.radii[0]/1.e3,
6371.lower_mantle.radii[1]/1.e3)
ax[i].set_ylabel('Depth (km)')
ax[0].set_xlabel('Pressure (GPa)')
ax[1].set_xlabel('Temperature (K)')
ax[2].set_xlabel('Gravity (m/s$^2$)')
ax[3].set_xlabel('Bullen parameter')
fig.set_tight_layout(True)
Warning: No module named 'cdd'. For full functionality of BurnMan, please install pycddlib.
The Planet class¶
In a 1D Planet, the pressure, gravity, temperature and temperature gradient at the interfaces between layers must be continuous. In BurnMan, it is possible to collect layers together into a Planet, and have the “make” method of Planet work out how to ensure continuity (at least for pressure, gravity and temperature; for the sake of flexibility, temperature gradient is allowed to be discontinuous).
In the following example, we build Planet Zog, a planet similar to Earth but a little simpler in mineralogical makeup. First, we create an adiabatic inner core. The inner core probably isn’t adiabatic, but this is largely unimportant for the innermost layer.
[2]:
from burnman import Composition
from burnman.tools.chemistry import formula_mass
# Compositions from midpoints of Hirose et al. (2021), ignoring carbon and hydrogen
inner_core_composition = Composition({'Fe': 94.4, 'Ni': 5., 'Si': 0.55, 'O': 0.05}, 'weight')
outer_core_composition = Composition({'Fe': 90., 'Ni': 5., 'Si': 2., 'O': 3.}, 'weight')
for c in [inner_core_composition, outer_core_composition]:
c.renormalize('atomic', 'total', 1.)
inner_core_elemental_composition = dict(inner_core_composition.atomic_composition)
outer_core_elemental_composition = dict(outer_core_composition.atomic_composition)
inner_core_molar_mass = formula_mass(inner_core_elemental_composition)
outer_core_molar_mass = formula_mass(outer_core_elemental_composition)
[3]:
icb_radius = 1220.e3
inner_core = Layer('inner core', radii=np.linspace(0., icb_radius, 21))
hcp_iron = minerals.SE_2015.hcp_iron()
params = hcp_iron.params
params['name'] = 'modified solid iron'
params['formula'] = inner_core_elemental_composition
params['molar_mass'] = inner_core_molar_mass
delta_V = 2.0e7
inner_core_material = Mineral(params=params,
property_modifiers=[['linear',
{'delta_E': 0.,
'delta_S': 0.,
'delta_V': delta_V}]])
# check that the new inner core material does what we expect:
hcp_iron.set_state(200.e9, 4000.)
inner_core_material.set_state(200.e9, 4000.)
assert np.abs(delta_V  (inner_core_material.V  hcp_iron.V)) < 1.e12
inner_core.set_material(inner_core_material)
inner_core.set_temperature_mode('adiabatic')
Now, we create an adiabatic outer core.
[4]:
cmb_radius = 3480.e3
outer_core = Layer('outer core', radii=np.linspace(icb_radius, cmb_radius, 21))
liq_iron = minerals.SE_2015.liquid_iron()
params = liq_iron.params
params['name'] = 'modified liquid iron'
params['formula'] = outer_core_elemental_composition
params['molar_mass'] = outer_core_molar_mass
delta_V = 2.3e7
outer_core_material = Mineral(params=params,
property_modifiers=[['linear',
{'delta_E': 0.,
'delta_S': 0.,
'delta_V': delta_V}]])
# check that the new inner core material does what we expect:
liq_iron.set_state(200.e9, 4000.)
outer_core_material.set_state(200.e9, 4000.)
assert np.abs(delta_V  (outer_core_material.V  liq_iron.V)) < 1.e12
outer_core.set_material(outer_core_material)
outer_core.set_temperature_mode('adiabatic')
Now, we assume that there is a single mantle layer that is convecting. We import a PerpleX input table that contains the material properties of pyrolite for this layer. We’ll use this for the convecting mantle layer. We apply a perturbed adiabatic temperature profile.
[5]:
from burnman import BoundaryLayerPerturbation
lab_radius = 6171.e3 # 200 km thick lithosphere
lab_temperature = 1550.
convecting_mantle_radii = np.linspace(cmb_radius, lab_radius, 101)
convecting_mantle = Layer('convecting mantle', radii=convecting_mantle_radii)
# Import a low resolution PerpleX data table.
fname = '../tutorial/data/pyrolite_perplex_table_lo_res.dat'
pyrolite = PerplexMaterial(fname, name='pyrolite')
convecting_mantle.set_material(pyrolite)
# Here we add a thermal boundary layer perturbation, assuming that the
# lower mantle has a Rayleigh number of 1.e7, and that the basal thermal
# boundary layer has a temperature jump of 840 K and the top
# boundary layer has a temperature jump of 60 K.
tbl_perturbation = BoundaryLayerPerturbation(radius_bottom=cmb_radius,
radius_top=lab_radius,
rayleigh_number=1.e7,
temperature_change=900.,
boundary_layer_ratio=60./900.)
# Onto this perturbation, we add a linear superadiabaticity term according
# to Anderson (he settled on 200 K over the lower mantle)
dT_superadiabatic = 300.*(convecting_mantle_radii  convecting_mantle_radii[1])/(convecting_mantle_radii[0]  convecting_mantle_radii[1])
convecting_mantle_tbl = (tbl_perturbation.temperature(convecting_mantle_radii)
+ dT_superadiabatic)
convecting_mantle.set_temperature_mode('perturbedadiabatic',
temperatures=convecting_mantle_tbl)
And the lithosphere has a userdefined conductive gradient.
[6]:
moho_radius = 6341.e3
moho_temperature = 620.
dunite = minerals.SLB_2011.mg_fe_olivine(molar_fractions=[0.92, 0.08])
lithospheric_mantle = Layer('lithospheric mantle',
radii=np.linspace(lab_radius, moho_radius, 31))
lithospheric_mantle.set_material(dunite)
lithospheric_mantle.set_temperature_mode('userdefined',
np.linspace(lab_temperature,
moho_temperature, 31))
Finally, we assume the crust has the density of andesine ~ 40% anorthite
[7]:
planet_radius = 6371.e3
surface_temperature = 300.
andesine = minerals.SLB_2011.plagioclase(molar_fractions=[0.4, 0.6])
crust = Layer('crust', radii=np.linspace(moho_radius, planet_radius, 11))
crust.set_material(andesine)
crust.set_temperature_mode('userdefined',
np.linspace(moho_temperature,
surface_temperature, 11))
Everything is ready! Let’s make our planet from its consistuent layers.
[8]:
planet_zog = Planet('Planet Zog',
[inner_core, outer_core,
convecting_mantle, lithospheric_mantle,
crust], verbose=True)
planet_zog.make()
Iteration 1 maximum relative pressure error: 9.5e01
Iteration 2 maximum relative pressure error: 4.8e01
Iteration 3 maximum relative pressure error: 2.1e01
Iteration 4 maximum relative pressure error: 8.5e02
Iteration 5 maximum relative pressure error: 3.4e02
Iteration 6 maximum relative pressure error: 1.3e02
Iteration 7 maximum relative pressure error: 5.2e03
Iteration 8 maximum relative pressure error: 2.0e03
Iteration 9 maximum relative pressure error: 7.9e04
Iteration 10 maximum relative pressure error: 3.1e04
Iteration 11 maximum relative pressure error: 1.2e04
Iteration 12 maximum relative pressure error: 4.7e05
Iteration 13 maximum relative pressure error: 1.8e05
Iteration 14 maximum relative pressure error: 7.1e06
Now we output the mass of the planet and moment of inertia and the mass of the individual layers:
[9]:
earth_mass = 5.972e24
earth_moment_of_inertia_factor = 0.3307
print(f'mass = {planet_zog.mass:.3e} (Earth = {earth_mass:.3e})')
print(f'moment of inertia factor= {planet_zog.moment_of_inertia_factor:.4f} '
f'(Earth = {earth_moment_of_inertia_factor:.4f})')
print('Layer mass fractions:')
for layer in planet_zog.layers:
print(f'{layer.name}: {layer.mass / planet_zog.mass:.3f}')
mass = 5.962e+24 (Earth = 5.972e+24)
moment of inertia factor= 0.3307 (Earth = 0.3307)
Layer mass fractions:
inner core: 0.016
outer core: 0.309
convecting mantle: 0.621
lithospheric mantle: 0.047
crust: 0.007
BurnMan contains some utility functions that output data in a format readable by seismic codes such as Axisem and Mineos.
[10]:
from burnman.tools.output_seismo import write_axisem_input
from burnman.tools.output_seismo import write_mineos_input
write_axisem_input([convecting_mantle, lithospheric_mantle, crust], modelname='zog_silicates', plotting=True)
write_mineos_input([convecting_mantle, lithospheric_mantle, crust], modelname='zog_silicates', plotting=True)
# Now we delete the newlycreated files. If you want them, comment out these lines.
import os
os.remove('axisem_zog_silicates.txt')
os.remove('mineos_zog_silicates.txt')
Writing axisem_zog_silicates.txt ...
Writing mineos_zog_silicates.txt ...
Let’s compare the properties of this planet to PREM
[11]:
import warnings
prem = burnman.seismic.PREM()
premdepth = prem.internal_depth_list()
premradii = 6371.e3  premdepth
with warnings.catch_warnings(record=True) as w:
eval = prem.evaluate(['density', 'pressure', 'gravity', 'v_s', 'v_p'])
premdensity, prempressure, premgravity, premvs, premvp = eval
print(w[1].message)
Gravity is not given in PREM and is now being computed. This will only work when density is defined for the entire planet. Use at your own risk.
Also create the Anzellini et al. (2013) geotherm:
[12]:
from scipy.interpolate import interp1d
d = np.loadtxt('../tutorial/data/Anzellini_2013_geotherm.dat')
Anz_interp = interp1d(d[:,0]*1.e9, d[:,1])
Finally, plot the 1D structure of the planet
[13]:
fig = plt.figure(figsize=(8, 5))
ax = [fig.add_subplot(2, 2, i) for i in range(1, 5)]
bounds = np.array([[layer.radii[0]/1.e3, layer.radii[1]/1.e3]
for layer in planet_zog.layers])
maxy = [15, 400, 12, 7000]
for bound in bounds:
for i in range(4):
ax[i].fill_betweenx([0., maxy[i]],
[bound[0], bound[0]],
[bound[1], bound[1]], alpha=0.2)
ax[0].plot(planet_zog.radii / 1.e3, planet_zog.density / 1.e3,
label=planet_zog.name)
ax[0].plot(premradii / 1.e3, premdensity / 1.e3, linestyle=':', label='PREM')
ax[0].set_ylabel('Density ($10^3$ kg/m$^3$)')
ax[0].legend()
# Make a subplot showing the calculated pressure profile
ax[1].plot(planet_zog.radii / 1.e3, planet_zog.pressure / 1.e9)
ax[1].plot(premradii / 1.e3, prempressure / 1.e9, linestyle=':')
ax[1].set_ylabel('Pressure (GPa)')
# Make a subplot showing the calculated gravity profile
ax[2].plot(planet_zog.radii / 1.e3, planet_zog.gravity)
ax[2].plot(premradii / 1.e3, premgravity, linestyle=':')
ax[2].set_ylabel('Gravity (m/s$^2)$')
ax[2].set_xlabel('Radius (km)')
# Make a subplot showing the calculated temperature profile
ax[3].plot(planet_zog.radii / 1.e3, planet_zog.temperature)
ax[3].set_ylabel('Temperature (K)')
ax[3].set_xlabel('Radius (km)')
ax[3].set_ylim(0.,)
# Finally, let's overlay some geotherms onto our model
# geotherm
labels = ['Stacey (1977)',
'Brown and Shankland (1981)',
'Anderson (1982)',
'Alfe et al. (2007)',
'Anzellini et al. (2013)']
short_labels = ['S1977',
'BS1981',
'A1982',
'A2007',
'A2013']
ax[3].plot(planet_zog.radii / 1.e3,
burnman.geotherm.stacey_continental(planet_zog.depth),
linestyle=':', label=short_labels[0])
mask = planet_zog.depth > 269999.
ax[3].plot(planet_zog.radii[mask] / 1.e3,
burnman.geotherm.brown_shankland(planet_zog.depth[mask]),
linestyle=':', label=short_labels[1])
ax[3].plot(planet_zog.radii / 1.e3,
burnman.geotherm.anderson(planet_zog.depth),
linestyle=':', label=short_labels[2])
ax[3].scatter([planet_zog.layers[0].radii[1] / 1.e3,
planet_zog.layers[1].radii[1] / 1.e3],
[5400., 4000.],
linestyle=':', label=short_labels[3])
mask = planet_zog.pressure < 330.e9
temperatures = Anz_interp(planet_zog.pressure[mask])
ax[3].plot(planet_zog.radii[mask] / 1.e3, temperatures,
linestyle=':', label=short_labels[4])
ax[3].legend()
for i in range(2):
ax[i].set_xticklabels([])
for i in range(4):
ax[i].set_xlim(0., max(planet_zog.radii) / 1.e3)
ax[i].set_ylim(0., maxy[i])
fig.set_tight_layout(True)
plt.show()
And that’s it! Next time, we’ll look at some of BurnMan’s fitting routines.
The BurnMan Tutorial
Part 4: Fitting¶
This file is part of BurnMan  a thermoelastic and thermodynamic toolkit for the Earth and Planetary Sciences
Copyright (C) 2012  2021 by the BurnMan team, released under the GNU GPL v2 or later.
Introduction¶
This ipython notebook is the fourth in a series designed to introduce new users to the code structure and functionalities present in BurnMan.
Demonstrates
burnman.optimize.eos_fitting.fit_PTV_data
burnman.optimize.composition_fitting.fit_composition_to_solution
burnman.optimize.composition_fitting.fit_phase_proportions_to_bulk_composition
Everything in BurnMan and in this tutorial is defined in SI units.
[1]:
import burnman
import numpy as np
import matplotlib.pyplot as plt
Warning: No module named 'cdd'. For full functionality of BurnMan, please install pycddlib.
Fitting parameters for an equation of state to experimental data¶
BurnMan contains leastsquares optimization functions that fit model parameters to data. There are two helper functions especially for use in fitting Mineral parameters to experimental data; these are burnman.optimize.eos_fitting.fit_PTp_data
(which can fit multiple kinds of data at the same time), and burnman.optimize.eos_fitting.fit_PTV_data
, which specifically fits only pressuretemperaturevolume data.
An extended example of fitting various kinds of data, outlier removal and detailed analysis can be found in examples/example_fit_eos.py
. In this tutorial, we shall focus solely on fitting room temperature pressuretemperaturevolume data. Specifically, the data we will fit is experimental volumes of stishovite, taken from Andrault et al. (2003). This data is provided in the form [P (GPa), V (Angstrom^3) and sigma_V (Angstrom^3)].
[2]:
PV = np.array([[0.0001, 46.5126, 0.0061],
[1.168, 46.3429, 0.0053],
[2.299, 46.1756, 0.0043],
[3.137, 46.0550, 0.0051],
[4.252, 45.8969, 0.0045],
[5.037, 45.7902, 0.0053],
[5.851, 45.6721, 0.0038],
[6.613, 45.5715, 0.0050],
[7.504, 45.4536, 0.0041],
[8.264, 45.3609, 0.0056],
[9.635, 45.1885, 0.0042],
[11.69, 44.947, 0.002],
[17.67, 44.264, 0.002],
[22.38, 43.776, 0.003],
[29.38, 43.073, 0.009],
[37.71, 42.278, 0.008],
[46.03, 41.544, 0.017],
[52.73, 40.999, 0.009],
[26.32, 43.164, 0.006],
[30.98, 42.772, 0.005],
[34.21, 42.407, 0.003],
[38.45, 42.093, 0.004],
[43.37, 41.610, 0.004],
[47.49, 41.280, 0.007]])
print(f'{len(PV)} data points loaded successfully.')
24 data points loaded successfully.
BurnMan works exclusively in SI units, so we have to convert from GPa to Pa, and Angstrom per cell into molar volume in m^3. The fitting function also takes covariance matrices as input, so we have to build those matrices.
[3]:
from burnman.tools.unitcell import molar_volume_from_unit_cell_volume
Z = 2. # number of formula units per unit cell in stishovite
PTV = np.array([PV[:,0]*1.e9,
298.15 * np.ones_like(PV[:,0]),
molar_volume_from_unit_cell_volume(PV[:,1], Z)]).T
# Here, we assume that the pressure uncertainties are equal to 3% of the total pressure,
# that the temperature uncertainties are negligible, and take the unit cell volume
# uncertainties from the paper.
# We also assume that the uncertainties in pressure and volume are uncorrelated.
nul = np.zeros_like(PTV[:,0])
PTV_covariances = np.array([[0.03*PTV[:,0], nul, nul],
[nul, nul, nul],
[nul, nul, molar_volume_from_unit_cell_volume(PV[:,2], Z)]]).T
PTV_covariances = np.power(PTV_covariances, 2.)
The next code block creates a Mineral object (stv
), providing starting guesses for the parameters. The user selects which parameters they wish to fit, and which they wish to keep fixed. The parameters of the Mineral object are automatically updated during fitting.
Finally, the optimized parameter values and their variances are printed to screen.
[4]:
stv = burnman.minerals.HP_2011_ds62.stv()
params = ['V_0', 'K_0', 'Kprime_0']
fitted_eos = burnman.optimize.eos_fitting.fit_PTV_data(stv, params, PTV, PTV_covariances, verbose=False)
print('Optimized equation of state for stishovite:')
burnman.utils.misc.pretty_print_values(fitted_eos.popt, fitted_eos.pcov, fitted_eos.fit_params)
print('')
Optimized equation of state for stishovite:
V_0: (1.4005 +/ 0.0001) x 1e05
K_0: (3.12 +/ 0.03) x 1e+11
Kprime_0: (4.4 +/ 0.2) x 1e+00
The fitted_eos
object contains a lot of useful information about the fit. In the next code block, we fit the corner plot of the covariances, showing the tradeoffs in different parameters.
[5]:
import matplotlib
matplotlib.rc('axes.formatter', useoffset=False) # turns offset off, makes for a more readable plot
fig = burnman.nonlinear_fitting.corner_plot(fitted_eos.popt, fitted_eos.pcov,
params)
axes = fig[1]
axes[1][0].set_xlabel('$V_0$ ($\\times 10^{5}$ m$^3$)')
axes[1][1].set_xlabel('$K_0$ ($\\times 10^{11}$ Pa)')
axes[0][0].set_ylabel('$K_0$ ($\\times 10^{11}$ Pa)')
axes[1][0].set_ylabel('$K\'_0$')
plt.show()
We now plot our optimized equation of state against the original data. BurnMan also includes a useful function burnman.optimize.nonlinear_fitting.confidence_prediction_bands
that can be used to calculate the nth percentile confidence and prediction bounds on a function given a model using the delta method.
[6]:
from burnman.utils.misc import attribute_function
from burnman.optimize.nonlinear_fitting import confidence_prediction_bands
T = 298.15
pressures = np.linspace(1.e5, 60.e9, 101)
temperatures = T*np.ones_like(pressures)
volumes = stv.evaluate(['V'], pressures, temperatures)[0]
PTVs = np.array([pressures, temperatures, volumes]).T
# Calculate the 95% confidence and prediction bands
cp_bands = confidence_prediction_bands(model=fitted_eos,
x_array=PTVs,
confidence_interval=0.95,
f=attribute_function(stv, 'V'),
flag='V')
plt.fill_between(pressures/1.e9, cp_bands[2] * 1.e6, cp_bands[3] * 1.e6,
color=[0.75, 0.25, 0.55], label='95% prediction bands')
plt.fill_between(pressures/1.e9, cp_bands[0] * 1.e6, cp_bands[1] * 1.e6,
color=[0.75, 0.95, 0.95], label='95% confidence bands')
plt.plot(PTVs[:,0] / 1.e9, PTVs[:,2] * 1.e6, label='Optimized fit for stishovite')
plt.errorbar(PTV[:,0] / 1.e9, PTV[:,2] * 1.e6,
xerr=np.sqrt(PTV_covariances[:,0,0]) / 1.e9,
yerr=np.sqrt(PTV_covariances[:,2,2]) * 1.e6,
linestyle='None', marker='o', label='Andrault et al. (2003)')
plt.ylabel("Volume (cm$^3$/mol)")
plt.xlabel("Pressure (GPa)")
plt.legend(loc="upper right")
plt.title("Stishovite EoS (room temperature)")
plt.show()
We can also calculate the confidence and prediction bands for any other property of the mineral. In the code block below, we calculate and plot the optimized isothermal bulk modulus and its uncertainties.
[7]:
cp_bands = confidence_prediction_bands(model=fitted_eos,
x_array=PTVs,
confidence_interval=0.95,
f=attribute_function(stv, 'K_T'),
flag='V')
plt.fill_between(pressures/1.e9, (cp_bands[0])/1.e9, (cp_bands[1])/1.e9, color=[0.75, 0.95, 0.95], label='95% confidence band')
plt.plot(pressures/1.e9, (cp_bands[0] + cp_bands[1])/2.e9, color='b', label='Best fit')
plt.ylabel("Bulk modulus (GPa)")
plt.xlabel("Pressure (GPa)")
plt.legend(loc="upper right")
plt.title("Stishovite EoS; uncertainty in bulk modulus (room temperature)")
plt.show()
Finding the best fit endmember proportions of a solution given a bulk composition¶
Let’s now turn our focus to a different kind of fitting. It is common in petrology to have a bulk composition of a phase (provided, for example, by electron probe microanalysis), and want to turn this composition into a formula that satisfies stoichiometric constraints. This can be formulated as a constrained, weighted least squares problem, and BurnMan can be used to solve these problems using the function burnman.optimize.composition_fitting.fit_composition_to_solution
.
In the following example, we shall create a model garnet composition, and then fit that to the Jennings and Holland (2015) garnet solution model. First, let’s look at the solution model endmembers (pyrope, almandine, grossular, andradite and knorringite):
[8]:
from burnman import minerals
gt = minerals.JH_2015.garnet()
print(f'Endmembers: {gt.endmember_names}')
print(f'Elements: {gt.elements}')
print('Stoichiometric matrix:')
print(gt.stoichiometric_matrix)
Endmembers: ['py', 'alm', 'gr', 'andr', 'knor']
Elements: ['Ca', 'Mg', 'Cr', 'Fe', 'Al', 'Si', 'O']
Stoichiometric matrix:
Matrix([[0, 3, 0, 0, 2, 3, 12], [0, 0, 0, 3, 2, 3, 12], [3, 0, 0, 0, 2, 3, 12], [3, 0, 0, 2, 0, 3, 12], [0, 3, 2, 0, 0, 3, 12]])
Now, let’s create a model garnet composition. A unique composition can be determined with the species Fe (total), Ca, Mg, Cr, Al, Si and Fe3+, all given in mole amounts. On top of this, we add some random noise (using a fixed seed so that the composition is reproducible).
[9]:
fitted_variables = ['Fe', 'Ca', 'Mg', 'Cr', 'Al', 'Si', 'Fe3+']
variable_values = np.array([1.1, 2., 0., 0, 1.9, 3., 0.1])
variable_covariances = np.eye(7)*0.01*0.01
# Add some noise.
v_err = np.random.rand(7)
np.random.seed(100)
variable_values = np.random.multivariate_normal(variable_values,
variable_covariances)
Importantly, Fe3+ isn’t an element or a sitespecies of the solution model, so we need to provide the linear conversion from Fe3+ to elements and/or site species. In this case, Fe3+ resides only on the second site (Site B), and the JH_2015.gt model has labelled Fe3+ on that site as Fef. Therefore, the conversion is simply Fe3+ = Fef_B.
[10]:
variable_conversions = {'Fe3+': {'Fef_B': 1.}}
Now we’re ready to do the fitting. The following line is all that is required, and yields as output the optimized parameters, the corresponding covariance matrix and the residual.
[11]:
from burnman.optimize.composition_fitting import fit_composition_to_solution
popt, pcov, res = fit_composition_to_solution(gt,
fitted_variables,
variable_values,
variable_covariances,
variable_conversions)
Finally, the optimized parameters can be used to set the composition of the solution model, and the optimized parameters printed to stdout.
[12]:
# We can set the composition of gt using the optimized parameters
gt.set_composition(popt)
# Print the optimized parameters and principal uncertainties
print('Molar fractions:')
for i in range(len(popt)):
print(f'{gt.endmember_names[i]}: '
f'{gt.molar_fractions[i]:.3f} +/ '
f'{np.sqrt(pcov[i][i]):.3f}')
print(f'Weighted residual: {res:.3f}')
Molar fractions:
py: 0.004 +/ 0.005
alm: 0.328 +/ 0.003
gr: 0.619 +/ 0.004
andr: 0.048 +/ 0.004
knor: 0.000 +/ 0.004
Weighted residual: 0.519
As in the equation of state fitting, a corner plot of the covariances can also be plotted.
[13]:
fig = burnman.nonlinear_fitting.corner_plot(popt, pcov, gt.endmember_names)
Fitting phase proportions to a bulk composition¶
Another common constrained weighted least squares problem involves fitting phase proportions, given their individual compositions and the overall bulk composition. This is particularly important in experimental petrology, where the bulk composition is known from a starting composition. In these cases, the residual after fitting is often used to assess whether the sample remained a closed system during the experiment.
In the following example, we take phase compositions and the bulk composition reported from high pressure experiments on a Martian mantle composition by Bertka and Fei (1997), and use these to calculate phase proportions in Mars mantle, and the quality of the experiments.
First, some tedious data preparation…
[14]:
import itertools
# Load and transpose input data
filename = '../burnman/data/input_fitting/Bertka_Fei_1997_mars_mantle.dat'
with open(filename) as f:
column_names = f.readline().strip().split()[1:]
data = np.genfromtxt(filename, dtype=None, encoding='utf8')
data = list(map(list, itertools.zip_longest(*data, fillvalue=None)))
# The first six columns are compositions given in weight % oxides
compositions = np.array(data[:6])
# The first row is the bulk composition
bulk_composition = compositions[:, 0]
# Load all the data into a dictionary
data = {column_names[i]: np.array(data[i])
for i in range(len(column_names))}
# Make ordered lists of samples (i.e. experiment ID) and phases
samples = []
phases = []
for i in range(len(data['sample'])):
if data['sample'][i] not in samples:
samples.append(data['sample'][i])
if data['phase'][i] not in phases:
phases.append(data['phase'][i])
samples.remove("bulk_composition")
phases.remove("bulk")
# Get the indices of all the phases present in each sample
sample_indices = [[i for i in range(len(data['sample']))
if data['sample'][i] == sample]
for sample in samples]
# Get the run pressures of each experiment
pressures = np.array([data['pressure'][indices[0]] for indices in sample_indices])
The following code block loops over each of the compositions, and finds the best weight proportions and uncertainties on those proportions.
[15]:
from burnman.optimize.composition_fitting import fit_phase_proportions_to_bulk_composition
# Create empty arrays to store the weight proportions of each phase,
# and the principal uncertainties (we do not use the covariances here,
# although they are calculated)
weight_proportions = np.zeros((len(samples), len(phases)))*np.NaN
weight_proportion_uncertainties = np.zeros((len(samples),
len(phases)))*np.NaN
residuals = []
# Loop over the samples, fitting phase proportions
# to the provided bulk composition
for i, sample in enumerate(samples):
# This line does the heavy lifting
popt, pcov, res = fit_phase_proportions_to_bulk_composition(compositions[:, sample_indices[i]],
bulk_composition)
residuals.append(res)
# Fill the correct elements of the weight_proportions
# and weight_proportion_uncertainties arrays
sample_phases = [data['phase'][i] for i in sample_indices[i]]
for j, phase in enumerate(sample_phases):
weight_proportions[i, phases.index(phase)] = popt[j]
weight_proportion_uncertainties[i, phases.index(phase)] = np.sqrt(pcov[j][j])
Finally, we plot the data.
[16]:
fig = plt.figure(figsize=(6, 6))
ax = [fig.add_subplot(4, 1, 1)]
ax.append(fig.add_subplot(4, 1, (2, 4)))
for i, phase in enumerate(phases):
ebar = plt.errorbar(pressures, weight_proportions[:, i],
yerr=weight_proportion_uncertainties[:, i],
fmt="none", zorder=2)
ax[1].scatter(pressures, weight_proportions[:, i], label=phase, zorder=3)
ax[0].set_title('Phase proportions in the Martian Mantle (Bertka and Fei, 1997)')
ax[0].scatter(pressures, residuals)
for i in range(2):
ax[i].set_xlim(0., 40.)
ax[1].set_ylim(0., 1.)
ax[0].set_ylabel('Residual')
ax[1].set_xlabel('Pressure (GPa)')
ax[1].set_ylabel('Phase fraction (wt %)')
ax[1].legend()
fig.set_tight_layout(True)
plt.show()
We can see from this plot that most of the residuals are below one, indicating that the probe analyses consistent with the bulk composition. Three analyses have higher residuals, which may indicate a problem with the experiments, or with the analyses around the wadsleyite field.
The phase proportions also show some nice trends; clinopyroxene weight percentage increases with pressure at the expense of orthopyroxene. Garnet / majorite percentage increases sharply as clinopyroxene is exhausted at 1416 GPa.
And we’re done! Next time, we’ll look at how to determine equilibrium assemblages in BurnMan.
The BurnMan Tutorial
Part 5: Equilibrium problems¶
This file is part of BurnMan  a thermoelastic and thermodynamic toolkit for the Earth and Planetary Sciences
Copyright (C) 2012  2021 by the BurnMan team, released under the GNU GPL v2 or later.
Introduction¶
This ipython notebook is the fifth in a series designed to introduce new users to the code structure and functionalities present in BurnMan.
Demonstrates
burnman.equilibrate, an experimental function that determines the bulk elemental composition, pressure, temperature, phase proportions and compositions of an assemblage subject to userdefined constraints.
Everything in BurnMan and in this tutorial is defined in SI units.
Phase equilibria¶
What BurnMan does and doesn’t do¶
Members of the BurnMan Team are often asked whether BurnMan does Gibbs energy minimization. The short answer to that is no, for three reasons: 1) Python is illsuited to such computationally intensive problems. 2) There are many pieces of software already in the community that do Gibbs energy minimization, including but not limited to: PerpleX, HeFESTo, Theriak Domino, MELTS, ENKI, FactSAGE (proprietary), and MMAEoS. 3) Gibbs minimization is a hard problem. The bruteforce pseudocompound/simplex technique employed by Perple_X is the only globally robust method, but clever techniques have to be used to make the computations tractable, and the solution found is generally only a (very close) approximation to the true minimum assemblage. More refined Newton / higher order schemes (e.g. HeFESTo, MELTS, ENKI) provide an exact solution, but can get stuck in local minima or even fail to find a solution.
So, with those things in mind, what does BurnMan do? Well, because BurnMan can compute the Gibbs energy and analytical derivatives of composite materials, it is well suited to solving the equilibrium relations for fixed assemblages. This is done using the burnman.equilibrate
function, which acts in a similar (but slightly more general) way to the THERMOCALC software developed by Tim Holland, Roger Powell and coworkers. Essentially, one chooses an assemblage (e.g. olivine + garnet +
orthopyroxene) and some equality constraints (typically related to bulk composition, pressure, temperature, entropy, volume, phase proportions or phase compositions) and the equilibrate
function attempts to find the remaining unknowns that satisfy those constraints.
In a sense, then, the equilibrate
function is simultaneously more powerful and more limited than Gibbs minimization techniques. It allows the user to investigate and plot metastable reactions, and quickly obtain answers to questions like “at what pressure does wadsleyite first become stable along a given isentrope?”. However, it is not designed to create PT tables of equilibrium assemblages. If a user wishes to do this for a complex problem, we refer them to other existing codes. BurnMan
also contains a useful utility material called burnman.PerplexMaterial
that is specifically designed to read in and interrogate PT data from PerpleX.
There are a couple more caveats to bear in mind. Firstly, the equilibrate
function is experimental and can certainly be improved. Equilibrium problems are highly nonlinear, and sometimes solvers struggle to find a solution. If you have a better, more robust way of solving these problems, we would love to hear from you! Secondly, the equilibrate
function is not completely free from the curse of multiple roots  sometimes there is more than one solution to the equilibrium problem, and
BurnMan (and indeed any equilibrium software) may find one a metastable root.
Equilibrating at fixed bulk composition¶
Fixed bulk composition problems are most similar to those asked by Gibbs minimization software like HeFESTo. Essentially, the only difference is that rather than allowing the assemblage to change to minimize the Gibbs energy, the assemblage is instead fixed.
In the following code block, we calculate the equilibrium assemblage of olivine, orthopyroxene and garnet for a mantle composition in the system NCFMAS at 10 GPa and 1500 K.
[1]:
import numpy as np
import matplotlib.pyplot as plt
import burnman
from burnman import equilibrate
from burnman.minerals import SLB_2011
# Set the pressure, temperature and composition
pressure = 3.e9
temperature = 1500.
composition = {'Na': 0.02, 'Fe': 0.2, 'Mg': 2.0, 'Si': 1.9,
'Ca': 0.2, 'Al': 0.4, 'O': 6.81}
# Create the assemblage
gt = SLB_2011.garnet()
ol = SLB_2011.mg_fe_olivine()
opx = SLB_2011.orthopyroxene()
assemblage = burnman.Composite(phases=[ol, opx, gt],
fractions=[0.7, 0.1, 0.2],
name='NCFMAS olopxgt assemblage')
# The solver uses the current compositions of each solution as a starting guess,
# so we have to set them here
ol.set_composition([0.93, 0.07])
opx.set_composition([0.8, 0.1, 0.05, 0.05])
gt.set_composition([0.8, 0.1, 0.05, 0.03, 0.02])
equality_constraints = [('P', 10.e9), ('T', 1500.)]
sol, prm = equilibrate(composition, assemblage, equality_constraints)
print(f'It is {sol.success} that equilibrate was successful')
print(sol.assemblage)
# The total entropy of the assemblage is the molar entropy
# multiplied by the number of moles in the assemblage
entropy = sol.assemblage.S*sol.assemblage.n_moles
Warning: No module named 'cdd'. For full functionality of BurnMan, please install pycddlib.
It is True that equilibrate was successful
Composite: NCFMAS olopxgt assemblage
P, T: 1e+10 Pa, 1500 K
Phase and endmember fractions:
olivine: 0.4971
Forsterite: 0.9339
Fayalite: 0.0661
orthopyroxene: 0.2925
Enstatite: 0.8640
Ferrosilite: 0.0687
Mg_Tschermaks: 0.0005
Ortho_Diopside: 0.0668
garnet: 0.2104
Pyrope: 0.4458
Almandine: 0.1239
Grossular: 0.2607
Mg_Majorite: 0.1258
Jd_Majorite: 0.0437
Each equality constraint can be a list of constraints, in which case equilibrate will loop over them. In the next code block we change the equality constraints to be a series of pressures which correspond to the total entropy obtained from the previous solve.
[2]:
equality_constraints = [('P', np.linspace(3.e9, 13.e9, 21)),
('S', entropy)]
sols, prm = equilibrate(composition, assemblage, equality_constraints)
The object sols
is now a 1D list of solution objects. Each one of these contains an equilibrium assemblage object that can be interrogated for any properties:
[3]:
data = np.array([[sol.assemblage.pressure,
sol.assemblage.temperature,
sol.assemblage.p_wave_velocity,
sol.assemblage.shear_wave_velocity,
sol.assemblage.molar_fractions[0],
sol.assemblage.molar_fractions[1],
sol.assemblage.molar_fractions[2]]
for sol in sols if sol.success])
The next code block plots these properties.
[4]:
fig = plt.figure(figsize=(12, 4))
ax = [fig.add_subplot(1, 3, i) for i in range(1, 4)]
P, T, V_p, V_s = data.T[:4]
phase_proportions = data.T[4:]
ax[0].plot(P/1.e9, T)
ax[1].plot(P/1.e9, V_p/1.e3)
ax[1].plot(P/1.e9, V_s/1.e3)
for i in range(3):
ax[2].plot(P/1.e9, phase_proportions[i], label=sol.assemblage.phases[i].name)
for i in range(3):
ax[i].set_xlabel('Pressure (GPa)')
ax[0].set_ylabel('Temperature (K)')
ax[1].set_ylabel('Seismic velocities (km/s)')
ax[2].set_ylabel('Molar phase proportions')
ax[2].legend()
plt.show()
From the above figure, we can see that the proportion of orthopyroxene is decreasing rapidly and is exhausted near 13 GPa. In the next code block, we determine the exact pressure at which orthopyroxene is exhausted.
[5]:
equality_constraints = [('phase_fraction', [opx, 0.]),
('S', entropy)]
sol, prm = equilibrate(composition, assemblage, equality_constraints)
print(f'Orthopyroxene is exhausted from the assemblage at {sol.assemblage.pressure/1.e9:.2f} GPa, {sol.assemblage.temperature:.2f} K.')
Orthopyroxene is exhausted from the assemblage at 13.04 GPa, 1511.64 K.
Equilibrating while allowing bulk composition to vary¶
[6]:
# Initialize the minerals we will use in this example.
ol = SLB_2011.mg_fe_olivine()
wad = SLB_2011.mg_fe_wadsleyite()
rw = SLB_2011.mg_fe_ringwoodite()
# Set the starting guess compositions for each of the solutions
ol.set_composition([0.90, 0.10])
wad.set_composition([0.90, 0.10])
rw.set_composition([0.80, 0.20])
First, we find the compositions of the three phases at the univariant.
[7]:
T = 1600.
composition = {'Fe': 0.2, 'Mg': 1.8, 'Si': 1.0, 'O': 4.0}
assemblage = burnman.Composite([ol, wad, rw], [1., 0., 0.])
equality_constraints = [('T', T),
('phase_fraction', (ol, 0.0)),
('phase_fraction', (rw, 0.0))]
free_compositional_vectors = [{'Mg': 1., 'Fe': 1.}]
sol, prm = equilibrate(composition, assemblage, equality_constraints,
free_compositional_vectors,
verbose=False)
if not sol.success:
raise Exception('Could not find solution for the univariant using '
'provided starting guesses.')
P_univariant = sol.assemblage.pressure
phase_names = [sol.assemblage.phases[i].name for i in range(3)]
x_fe_mbr = [sol.assemblage.phases[i].molar_fractions[1] for i in range(3)]
print(f'Univariant pressure at {T:.0f} K: {P_univariant/1.e9:.3f} GPa')
print('Fe2SiO4 concentrations at the univariant:')
for i in range(3):
print(f'{phase_names[i]}: {x_fe_mbr[i]:.2f}')
Univariant pressure at 1600 K: 12.002 GPa
Fe2SiO4 concentrations at the univariant:
olivine: 0.22
wadsleyite: 0.37
ringwoodite: 0.50
Now we solve for the stable sections of the three binary loops
[8]:
output = []
for (m1, m2, x_fe_m1) in [[ol, wad, np.linspace(x_fe_mbr[0], 0.001, 20)],
[ol, rw, np.linspace(x_fe_mbr[0], 0.999, 20)],
[wad, rw, np.linspace(x_fe_mbr[1], 0.001, 20)]]:
assemblage = burnman.Composite([m1, m2], [1., 0.])
# Reset the compositions of the two phases to have compositions
# close to those at the univariant point
m1.set_composition([1.x_fe_mbr[1], x_fe_mbr[1]])
m2.set_composition([1.x_fe_mbr[1], x_fe_mbr[1]])
# Also set the pressure and temperature
assemblage.set_state(P_univariant, T)
# Here our equality constraints are temperature,
# the phase fraction of the second phase,
# and we loop over the composition of the first phase.
equality_constraints = [('T', T),
('phase_composition',
(m1, [['Mg_A', 'Fe_A'],
[0., 1.], [1., 1.], x_fe_m1])),
('phase_fraction', (m2, 0.0))]
sols, prm = equilibrate(composition, assemblage,
equality_constraints,
free_compositional_vectors,
verbose=False)
# Process the solutions
out = np.array([[sol.assemblage.pressure,
sol.assemblage.phases[0].molar_fractions[1],
sol.assemblage.phases[1].molar_fractions[1]]
for sol in sols if sol.success])
output.append(out)
output = np.array(output)
Finally, we do some plotting
[9]:
fig = plt.figure()
ax = [fig.add_subplot(1, 1, 1)]
color='purple'
# Plot the line connecting the three phases
ax[0].plot([x_fe_mbr[0], x_fe_mbr[2]],
[P_univariant/1.e9, P_univariant/1.e9], color=color)
for i in range(3):
if i == 0:
ax[0].plot(output[i,:,1], output[i,:,0]/1.e9, color=color, label=f'{T} K')
else:
ax[0].plot(output[i,:,1], output[i,:,0]/1.e9, color=color)
ax[0].plot(output[i,:,2], output[i,:,0]/1.e9, color=color)
ax[0].fill_betweenx(output[i,:,0]/1.e9, output[i,:,1], output[i,:,2],
color=color, alpha=0.2)
ax[0].text(0.1, 6., 'olivine', horizontalalignment='left')
ax[0].text(0.015, 14.2, 'wadsleyite', horizontalalignment='left',
bbox=dict(facecolor='white',
edgecolor='white',
boxstyle='round,pad=0.2'))
ax[0].text(0.9, 15., 'ringwoodite', horizontalalignment='right')
ax[0].set_xlim(0., 1.)
ax[0].set_ylim(0.,20.)
ax[0].set_xlabel('p(Fe$_2$SiO$_4$)')
ax[0].set_ylabel('Pressure (GPa)')
ax[0].legend()
plt.show()
And we’re done!
CIDER Tutorial 2014¶
 The tutorial for BurnMan presented at CIDER 2014 consists of three separate units:
CIDER 2014 BurnMan Tutorial — step 1¶
In this first part of the tutorial we will acquaint ourselves with a basic script for calculating the elastic properties of a mantle mineralogical model.
In general, there are three portions of this script:
1) Define a set of pressures and temperatures at which we want to calculate elastic properties
2) Setup a composite of minerals (or “rock”) and calculate its elastic properties at those pressures and temperatures.
3) Plot those elastic properties, and compare them to a seismic model, in this case PREM
The script is basically already written, and should run as is by typing:
python step_1.py
on the command line. However, the mineral model for the rock is not very realistic, and you will want to change it to one that is more in accordance with what we think the bulk composition of Earth’s lower mantle is.
When run (without putting in a more realistic composition), the program produces the following image:
Your goal in this tutorial is to improve this awful fit…
CIDER 2014 BurnMan Tutorial — step 2¶
In this second part of the tutorial we try to get a closer fit to our 1D seismic reference model. In the simple Mg, Si, and O model that we used in step 1 there was one free parameter, namely phase_1_fraction, which goes between zero and one.
In this script we want to explore how good of a fit to PREM we can get by varying this fraction. We create a simple function that calculates a misfit between PREM and our mineral model as a function of phase_1_fraction, and then plot this misfit function to try to find a best model.
This script may be run by typing
python step_2.py
Without changing any input, the program should produce the following image showing the misfit as a function of perovskite content:
CIDER 2014 BurnMan Tutorial — step 3¶
In the previous two steps of the tutorial we tried to find a very simple mineralogical model that best fit the 1D seismic model PREM. But we know that there is consideral uncertainty in many of the mineral physical parameters that control how the elastic properties of minerals change with pressure and temperature. In this step we explore how uncertainties in these parameters might affect the conclusions you draw.
The strategy here is to make many different “realizations” of the rock that you determined was the closest fit to PREM, where each realization has its mineral physical parameters perturbed by a small amount, hopefully related to the uncertainty in that parameter. In particular, we will look at how perturbations to \(K_{0}^{'}\) and \(G_{0}^{'}\) (the pressure derivatives of the bulk and shear modulus, respectively) change the calculated 1D seismic profiles.
This script may be run by typing
python step_3.py
After changing the standard deviations for \(K_{0}^{'}\) and \(G_{0}^{'}\) to 0.2, the following figure of velocities for 1000 realizations is produced:
Examples¶
BurnMan comes with a large collection of example programs under examples/. Below you can find a summary of the different examples. They are grouped into three categories: Class examples, Simple Examples and More Advanced Examples. The Class examples introduce the main class types in BurnMan, and covers similar ground to The BurnMan Tutorial but in a little more detail. Simple Examples provides introductions to some of the seismic and chemical functions. More Advanced Examples covers functionality which is more suited to research projects.
Finally, we also include the scripts that were used for all computations and figures in the 2014 BurnMan paper in the misc/ folder, see Reproducing Cottaar, Heister, Rose and Unterborn (2014).
Class examples¶
The following is a list of examples that introduce the main classes of BurnMan objects:
example_geotherms
, andexample_mineral¶
This example shows how to create mineral objects in BurnMan, and how to output their thermodynamic and thermoelastic quantities.
Mineral objects are the building blocks for more complex objects in BurnMan. These objects are intended to represent minerals (or melts, or fluids) of fixed composition, with a well defined equation of state that defines the relationship between the current state (pressure and temperature) of the mineral and its thermodynamic potentials and derivatives (such as volume and entropy).
Mineral objects are initialized with a dictionary containing all of the parameters required by the desired equation of state. BurnMan contains implementations of may equations of state (Equations of state).
Uses:
Demonstrates:
Different ways to define an endmember
How to set state
How to output thermodynamic and thermoelastic properties
Resulting figure:
example_gibbs_modifiers¶
This example script demonstrates the modifications to the gibbs free energy (and derivatives) that can be applied as masks over the results from the equations of state.
These modifications currently take the forms:
Landau corrections (implementations of Putnis (1992) and Holland and Powell (2011)
BraggWilliams corrections (implementation of Holland and Powell (1996))
Linear (a simple delta_E + delta_V*P  delta_S*T
Magnetic (Chin, Hertzman and Sundman (1987))
Uses:
Demonstrates:
creating a mineral with excess contributions
calculating thermodynamic properties
Resulting figures:
example_solution¶
This example shows how to create different solution models and output thermodynamic and thermoelastic quantities.
There are four main types of solution currently implemented in BurnMan:
Ideal solutions
Symmmetric solutions
Asymmetric solutions
Subregular solutions
These solutions can potentially deal with:
Disordered endmembers (more than one element on a crystallographic site)
Site vacancies
More than one valence/spin state of the same element on a site
Uses:
Demonstrates:
Different ways to define a solution
How to set composition and state
How to output thermodynamic and thermoelastic properties
Resulting figures:
example_composite¶
This example demonstrates the functionalities of the burnman.Composite class.
Uses:
Demonstrates:
How to initialize a composite object containing minerals and solutions
How to set state and composition of composite objects
How to interrogate composite objects for their compositional, thermodynamic and thermoelastic properties.
How to use the stoichiometric and reaction affinity methods to solve simple thermodynamic equilibrium problems.
Resulting figures:
example_mineral¶
This example shows how to create mineral objects in BurnMan, and how to output their thermodynamic and thermoelastic quantities.
Mineral objects are the building blocks for more complex objects in BurnMan. These objects are intended to represent minerals (or melts, or fluids) of fixed composition, with a well defined equation of state that defines the relationship between the current state (pressure and temperature) of the mineral and its thermodynamic potentials and derivatives (such as volume and entropy).
Mineral objects are initialized with a dictionary containing all of the parameters required by the desired equation of state. BurnMan contains implementations of may equations of state (Equations of state).
Uses:
Demonstrates:
Different ways to define an endmember
How to set state
How to output thermodynamic and thermoelastic properties
example_calibrants¶
This example demonstrates the use of BurnMan’s library of pressure calibrants. These calibrants are strippeddown versions of BurnMan’s minerals, in that they are only designed to return pressure as a function of volume and temperature or volume as a function of pressure and temperature.
Uses:
burnman.calibrants.tools.pressure_to_pressure()
Decker (1971) calibration for NaCl
Demonstrates:
Use of the Calibrant class
Conversion between pressures given by two different calibrations.
example_anisotropy¶
This example illustrates the basic functions required to convert an elastic stiffness tensor into elastic properties.
Specifically uses:
Demonstrates:
anisotropic functions
Resulting figure:
example_anisotropic_mineral¶
This example illustrates how to create and interrogate an AnisotropicMineral object.
Specifically uses:
Demonstrates:
anisotropic functions
Resulting figure:
example_geotherms¶
This example shows each of the geotherms currently possible with BurnMan. These are:
Brown and Shankland, 1981 [BS81]
Anderson, 1982 [And82]
Watson and Baxter, 2007 [WB07]
linear extrapolation
Read in from file from user
Adiabatic from potential temperature and choice of mineral
Uses:
burnman.geotherm.brown_shankland()
burnman.geotherm.anderson()
input geotherm file input_geotherm/example_geotherm.txt (optional)
burnman.Composite
for adiabatDemonstrates:
the available geotherms
Resulting figure:
example_composition¶
This example script demonstrates the use of BurnMan’s Composition class.
Uses:
Demonstrates:
Creating an instance of the Composition class with a molar or weight composition
Printing weight, molar, atomic compositions
Renormalizing compositions
Modifying the independent set of components
Modifying compositions by adding and removing components
Simple Examples¶
 The following is a list of simple examples:
example_beginner¶
This example script is intended for absolute beginners to BurnMan. We cover importing BurnMan modules, creating a composite material, and calculating its seismic properties at lower mantle pressures and temperatures. Afterwards, we plot it against a 1D seismic model for visual comparison.
Uses:
burnman.seismic.PREM
burnman.geotherm.brown_shankland()
Demonstrates:
creating basic composites
calculating thermoelastic properties
seismic comparison
Resulting figure:
example_seismic¶
Shows the various ways to input seismic models (\(V_s, V_p, V_{\phi}, \rho\)) as a function of depth (or pressure) as well as different velocity model libraries available within Burnman:
This example will first calculate or read in a seismic model and plot the model along the defined pressure range. The example also illustrates how to import a seismic model of your choice, here shown by importing AK135 [KEB95].
Uses:
Demonstrates:
Utilization of library seismic models within BurnMan
Input of userdefined seismic models
Resulting figures:
example_composite_seismic_velocities¶
This example shows how to create different minerals, how to compute seismic velocities, and how to compare them to a seismic reference model.
There are many different ways in BurnMan to combine minerals into a composition. Here we present a couple of examples:
Two minerals mixed in simple mole fractions. Can be chosen from the BurnMan libraries or from user defined minerals (see example_user_input_material)
Example with three minerals
Using preset solutions
Defining your own solution
To turn a method of mineral creation “on” the first if statement above the method must be set to True, with all others set to False.
Note: These minerals can include a spin transition in (Mg,Fe)O, see example_spintransition.py for explanation of how to implement this
Uses:
Demonstrates:
Different ways to define a composite
Using minerals and solutions
Compare computations to seismic models
Resulting figure:
example_averaging¶
This example shows the effect of different averaging schemes. Currently four averaging schemes are available:
VoightReussHill
Voight averaging
Reuss averaging
HashinShtrikman averaging
See [WDOConnell76] Journal of Geophysics and Space Physics for explanations of each averaging scheme.
Specifically uses:
Demonstrates:
implemented averaging schemes
Resulting figure:
example_chemical_potentials¶
This example shows how to obtain chemical potentials and associated properties from an assemblage.
Demonstrates:
How to calculate chemical potentials of an assemblage.
How to compute fugacities and relative fugacities.
Resulting figure:
More Advanced Examples¶
 Advanced examples:
example_spintransition¶
This example shows the different minerals that are implemented with a spin transition. Minerals with spin transition are implemented by defining two separate minerals (one for the low and one for the high spin state). Then a third dynamic mineral is created that switches between the two previously defined minerals by comparing the current pressure to the transition pressure.
Specifically uses:
burnman.mineral_helpers.HelperSpinTransition()
Demonstrates:
implementation of spin transition in (Mg,Fe)O at user defined pressure
Resulting figure:
example_spintransition_thermal¶
This example illustrates how to create a nonideal solution model for (Mg,Fe^{HS},Fe^{LS})O ferropericlase that has a gradual spin transition at finite temperature. First, we define the MgO endmember and two endmembers for the low and high spin states of FeO. Then we create a regular/symmetric solution that incorporates all three endmembers. The modified solution class contains a method called set_equilibrium_composition, which calculates the equilibrium proportions of the low and high spin phases at the desired bulk composition, pressure and temperature.
In this example, we neglect the elastic component of mixing. We also implicitly apply the BraggWilliams approximation (i.e., we assume that there is no shortrange order by only incorporating interactions that are a function of the average occupancy of species on each distinct site). Furthermore, the one site model [Mg,Fe^{HS},Fe^{LS}]O explicitly precludes long range order.
Specifically uses:
Demonstrates:
implementation of gradual spin transition in (Mg,Fe)O at a userdefined pressure and temperature
Resulting figure:
example_user_input_material¶
Shows user how to input a mineral of his/her choice without usint the library and which physical values need to be input for BurnMan to calculate \(V_P, V_\Phi, V_S\) and density at depth.
Specifically uses:
Demonstrates:
how to create your own minerals
example_optimize_pv¶
Vary the amount perovskite vs. ferropericlase and compute the error in the seismic data against PREM. For more extensive comments on this setup, see tutorial/step_2.py
Uses:
burnman.seismic.PREM
burnman.geotherm.brown_shankland()
Demonstrates:
compare errors between models
loops over models
Resulting figure:
example_compare_all_methods¶
This example demonstrates how to call each of the individual calculation methodologies that exist within BurnMan. See below for current options. This example calculates seismic velocity profiles for the same set of minerals and a plot of \(V_s, V_\phi\) and \(\rho\) is produce for the user to compare each of the different methods.
Specifically uses:
Demonstrates:
Each method for calculating velocity profiles currently included within BurnMan
Resulting figure:
example_build_planet¶
For Earth we have wellconstrained onedimensional density models. This allows us to calculate pressure as a function of depth. Furthermore, petrologic data and assumptions regarding the convective state of the planet allow us to estimate the temperature.
For planets other than Earth we have much less information, and in particular we know almost nothing about the pressure and temperature in the interior. Instead, we tend to have measurements of things like mass, radius, and momentofinertia. We would like to be able to make a model of the planet’s interior that is consistent with those measurements.
However, there is a difficulty with this. In order to know the density of the planetary material, we need to know the pressure and temperature. In order to know the pressure, we need to know the gravity profile. And in order to the the gravity profile, we need to know the density. This is a nonlinear problem which requires us to iterate to find a selfconsistent solution.
This example allows the user to define layers of planets of known outer radius and self consistently solve for the density, pressure and gravity profiles. The calculation will iterate until the difference between central pressure calculations are less than 1e5. The planet class in BurnMan (../burnman/planet.py) allows users to call multiple properties of the model planet after calculations, such as the mass of an individual layer, the total mass of the planet and the moment if inertia. See planets.py for information on each of the parameters which can be called.
Uses:
Demonstrates:
setting up a planet
computing its selfconsistent state
computing various parameters for the planet
seismic comparison
Resulting figure:
example_fit_composition¶
This example shows how to fit compositional data to a solution model, how to partition a bulk composition between phases of known composition, and how to assess goodness of fit.
Uses:
burnman.optimize.composition_fitting.fit_phase_proportions_to_bulk_composition()
burnman.optimize.composition_fitting.fit_composition_to_solution()
Demonstrates:
Fitting compositional data to a solution model
Partitioning of a bulk composition between phases of known composition
Assessing goodness of fit.
Resulting figure:
example_fit_data¶
This example demonstrates BurnMan’s functionality to fit various mineral physics data to an EoS of the user’s choice.
Please note also the separate file example_fit_eos.py, which can be viewed as a more advanced example in the same general field.
teaches:  least squares fitting
Resulting figures:
example_fit_eos¶
This example demonstrates BurnMan’s functionality to fit data to an EoS of the user’s choice.
The first example deals with simple PVT fitting. The second example illustrates how powerful it can be to provide nonPVT constraints to the same fitting problem.
teaches:  least squares fitting
Last seven resulting figures:
example_fit_solution¶
This example demonstrates how to fit parameters for solution models using a range of compositionallyvariable experimental data.
The example in this file deals with finding optimized parameters for the forsteritefayalite binary using a mixture of volume and seismic velocity data.
teaches:  least squares fitting for solution data
Resulting figures:
example_equilibrate¶
This example demonstrates how BurnMan may be used to calculate the equilibrium state for an assemblage of a fixed bulk composition given two constraints. Each constraint has the form [<constraint type>, <constraint>], where <constraint type> is one of the strings: P, T, S, V, X, PT_ellipse, phase_fraction, or phase_composition. The <constraint> object should either be a float or an array of floats for P, T, S, V (representing the desired pressure, temperature, entropy or volume of the material). If the constraint type is X (a generic constraint on the solution vector) then the constraint c is represented by the following equality: np.dot(c[0], x)  c[1]. If the constraint type is PT_ellipse, the equality is given by norm(([P, T]  c[0])/c[1])  1. The constraint_type phase_fraction assumes a tuple of the phase object (which must be one of the phases in the burnman.Composite) and a float or vector corresponding to the phase fractions. Finally, a phase_composition constraint has the format (site_names, n, d, v), where site names dictates the sites involved in the equality constraint. The equality constraint is given by n*x/d*x = v, where x are the site occupancies and n and d are fixed vectors of site coefficients. So, one could for example choose a constraint ([Mg_A, Fe_A], [1., 0.], [1., 1.], [0.5]) which would correspond to equal amounts Mg and Fe on the A site.
This script provides a number of examples, which can be turned on and off with a series of boolean variables. In order of complexity:
run_aluminosilicates: Creates the classic aluminosilicate diagram involving univariate reactions between andalusite, sillimanite and kyanite.
run_ordering: Calculates the state of order of Jennings and Holland (2015) orthopyroxene in the simple enfs binary at 1 bar.
run_gt_solvus: Demonstrates the shape of the pyropegrossular solvus.
run_fper_ol: Calculates the equilibrium MgFe partitioning between ferropericlase and olivine.
run_fixed_ol_composition: Calculates the composition of wadsleyite in equilibrium with olivine of a fixed composition at a fixed pressure.
run_upper_mantle: Calculates the equilibrium compositions and phase proportions for an olopxgt composite in an NCFMAS bulk composition.
run_lower_mantle: Calculates temperatures and assemblage properties along an isentrope in the lower mantle. Includes calculations of the postperovskitein and bridgmaniteout lines.
run_olivine_polymorphs: Produces a PT pseudosection for a fo90 composition.
Uses:
Resulting figures:
example_olivine_binary¶
This example demonstrates how BurnMan may be used to calculate the equilibrium binary phase diagram for the three olivine polymorphs (olivine, wadsleyite and ringwoodite).
The calculations use the equilibrate function. Unlike the examples in example_equilibrate.py, which are constrained to a fixed bulk composition, the bulk composition is allowed to vary along the vector [n_Mg  n_Fe].
Uses:
Resulting figures:
Reproducing Cottaar, Heister, Rose and Unterborn (2014)¶
In this section we include the scripts that were used for all computations and figures in the 2014 BurnMan paper: Cottaar, Heister, Rose & Unterborn (2014) [CHRU14]
paper_averaging¶
This script reproduces [CHRU14], Figure 2.
This example shows the effect of different averaging schemes. Currently four averaging schemes are available: 1. VoightReussHill 2. Voight averaging 3. Reuss averaging 4. HashinShtrikman averaging
See [WDOConnell76] for explanations of each averaging scheme.
requires:  geotherms  compute seismic velocities
teaches:  averaging
paper_benchmark¶
This script reproduces the benchmark in [CHRU14], Figure 3.
paper_fit_data¶
This script reproduces [CHRU14] Figure 4.
This example demonstrates BurnMan’s functionality to fit thermoelastic data to both 2nd and 3rd orders using the EoS of the user’s choice at 300 K. User’s must create a file with \(P, T\) and \(V_s\). See input_minphys/ for example input files.
requires:  compute seismic velocities
teaches:  averaging
paper_incorrect_averaging¶
This script reproduces [CHRU14], Figure 5. Attempt to reproduce Figure 6.12 from [Mur13]
paper_opt_pv¶
This script reproduces [CHRU14], Figure 6. Vary the amount perovskite vs. ferropericlase and compute the error in the seismic data against PREM.
requires:  creating minerals  compute seismic velocities  geotherms  seismic models  seismic comparison
teaches:  compare errors between models  loops over models
paper_onefit¶
This script reproduces [CHRU14], Figure 7. It shows an example for a best fit for a pyrolitic model within mineralogical error bars.
paper_uncertain¶
This script reproduces [CHRU14], Figure 8. It shows the sensitivity of the velocities to various mineralogical parameters.
Misc or work in progress¶
example_grid¶
This example shows how to evaluate seismic quantities on a \(P,T\) grid.
example_woutput¶
This example explains how to perform the basic i/o of BurnMan. A method of calculation is chosen, a composite mineral/material (see example_composition.py for explanation of this process) is created in the class “rock,” finally a geotherm is created and seismic velocities calculated.
Postcalculation, the results are written to a simple text file to plot/manipulate at the user’s whim.
requires:  creating minerals  compute seismic velocities  geotherms
teaches:  output computed seismic data to file
Autogenerated Full API¶
Materials¶
Burnman operates on materials (type Material
)
most prominently in the form of minerals
(Mineral
) and composites (Composite
).
Material Base Class¶
 class burnman.Material[source]¶
Bases:
object
Base class for all materials. The main functionality is unroll() which returns a list of objects of type
Mineral
and their molar fractions. This class is available asburnman.Material
.The user needs to call set_method() (once in the beginning) and set_state() before querying the material with unroll() or density().
 property name¶
Humanreadable name of this material.
By default this will return the name of the class, but it can be set to an arbitrary string. Overriden in Mineral.
 set_method(method)[source]¶
Set the averaging method. See Averaging Schemes for details.
Notes
Needs to be implemented in derived classes.
 to_string()[source]¶
Returns a humanreadable name of this material. The default implementation will return the name of the class, which is a reasonable default.
 Returns
 namestring
Name of this material.
 print_minerals_of_current_state()[source]¶
Print a humanreadable representation of this Material at the current P, T as a list of minerals. This requires set_state() has been called before.
 set_state(pressure, temperature)[source]¶
Set the material to the given pressure and temperature.
 Parameters
 pressurefloat
The desired pressure in [Pa].
 temperaturefloat
The desired temperature in [K].
 set_state_with_volume(volume, temperature, pressure_guesses=[0.0, 10000000000.0])[source]¶
This function acts similarly to set_state, but takes volume and temperature as input to find the pressure. In order to ensure selfconsistency, this function does not use any pressure functions from the material classes, but instead finds the pressure using the brentq rootfinding method.
 Parameters
 volumefloat
The desired molar volume of the mineral [m^3].
 temperaturefloat
The desired temperature of the mineral [K].
 pressure_guesseslist of floats (default: [5.e9, 10.e9])
The initial low and high guesses for bracketing of the pressure [Pa]. These guesses should preferably bound the correct pressure, but do not need to do so. More importantly, they should not lie outside the valid region of the equation of state.
 reset()[source]¶
Resets all cached material properties.
It is typically not required for the user to call this function.
 unroll()[source]¶
Unroll this material into a list of
burnman.Mineral
and their molar fractions. All averaging schemes then operate on this list of minerals. Note that the return value of this function may depend on the current state (temperature, pressure). Returns
 fractionslist of float
List of molar fractions, should sum to 1.0.
 mineralslist of
burnman.Mineral
List of minerals.
Notes
Needs to be implemented in derived classes.
 evaluate(vars_list, pressures, temperatures)[source]¶
Returns an array of material properties requested through a list of strings at given pressure and temperature conditions. At the end it resets the set_state to the original values. The user needs to call set_method() before.
 Parameters
 vars_listlist of strings
Variables to be returned for given conditions
 pressuresndlist or ndarray of float
ndimensional array of pressures in [Pa].
 temperaturesndlist or ndarray of float
ndimensional array of temperatures in [K].
 Returns
 outputarray of array of float
Array returning all variables at given pressure/temperature values. output[i][j] is property vars_list[j] and temperatures[i] and pressures[i].
 property pressure¶
Returns current pressure that was set with
set_state()
. Returns
 pressurefloat
Pressure in [Pa].
Notes
Aliased with
P()
.
 property temperature¶
Returns current temperature that was set with
set_state()
. Returns
 temperaturefloat
Temperature in [K].
Notes
Aliased with
T()
.
 property molar_internal_energy¶
Returns the molar internal energy of the mineral.
 Returns
 molar_internal_energyfloat
The internal energy in [J/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
energy()
.
 property molar_gibbs¶
Returns the molar Gibbs free energy of the mineral.
 Returns
 molar_gibbsfloat
Gibbs free energy in [J/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
gibbs()
.
 property molar_helmholtz¶
Returns the molar Helmholtz free energy of the mineral.
 Returns
 molar_helmholtzfloat
Helmholtz free energy in [J/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
helmholtz()
.
 property molar_mass¶
Returns molar mass of the mineral.
 Returns
 molar_massfloat
Molar mass in [kg/mol].
Notes
Needs to be implemented in derived classes.
 property molar_volume¶
Returns molar volume of the mineral.
 Returns
 molar_volumefloat
Molar volume in [m^3/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
V()
.
 property density¶
Returns the density of this material.
 Returns
 densityfloat
The density of this material in [kg/m^3].
Notes
Needs to be implemented in derived classes.
Aliased with
rho()
.
 property molar_entropy¶
Returns molar entropy of the mineral.
 Returns
 molar_entropyfloat
Entropy in [J/K/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
S()
.
 property molar_enthalpy¶
Returns molar enthalpy of the mineral.
 Returns
 molar_enthalpyfloat
Enthalpy in [J/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
H()
.
 property isothermal_bulk_modulus¶
Returns isothermal bulk modulus of the material.
 Returns
 isothermal_bulk_modulusfloat
Bulk modulus in [Pa].
Notes
Needs to be implemented in derived classes.
Aliased with
K_T()
.
 property adiabatic_bulk_modulus¶
Returns the adiabatic bulk modulus of the mineral.
 Returns
 adiabatic_bulk_modulusfloat
Adiabatic bulk modulus in [Pa].
Notes
Needs to be implemented in derived classes.
Aliased with
K_S()
.
 property isothermal_compressibility¶
Returns isothermal compressibility of the mineral (or inverse isothermal bulk modulus).
 Returns
 (K_T)^1float
Compressibility in [1/Pa].
Notes
Needs to be implemented in derived classes.
Aliased with
beta_T()
.
 property adiabatic_compressibility¶
Returns adiabatic compressibility of the mineral (or inverse adiabatic bulk modulus).
 Returns
 adiabatic_compressibilityfloat
adiabatic compressibility in [1/Pa].
Notes
Needs to be implemented in derived classes.
Aliased with
beta_S()
.
 property shear_modulus¶
Returns shear modulus of the mineral.
 Returns
 shear_modulusfloat
Shear modulus in [Pa].
Notes
Needs to be implemented in derived classes.
Aliased with
beta_G()
.
 property p_wave_velocity¶
Returns P wave speed of the mineral.
 Returns
 p_wave_velocityfloat
P wave speed in [m/s].
Notes
Needs to be implemented in derived classes.
Aliased with
v_p()
.
 property bulk_sound_velocity¶
Returns bulk sound speed of the mineral.
 Returns
 bulk sound velocity: float
Sound velocity in [m/s].
Notes
Needs to be implemented in derived classes.
Aliased with
v_phi()
.
 property shear_wave_velocity¶
Returns shear wave speed of the mineral.
 Returns
 shear_wave_velocityfloat
Wave speed in [m/s].
Notes
Needs to be implemented in derived classes.
Aliased with
v_s()
.
 property grueneisen_parameter¶
Returns the grueneisen parameter of the mineral.
 Returns
 grfloat
Grueneisen parameters [unitless].
Notes
Needs to be implemented in derived classes.
Aliased with
gr()
.
 property thermal_expansivity¶
Returns thermal expansion coefficient of the mineral.
 Returns
 alphafloat
Thermal expansivity in [1/K].
Notes
Needs to be implemented in derived classes.
Aliased with
alpha()
.
 property molar_heat_capacity_v¶
Returns molar heat capacity at constant volume of the mineral.
 Returns
 molar_heat_capacity_vfloat
Heat capacity in [J/K/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
C_v()
.
 property molar_heat_capacity_p¶
Returns molar heat capacity at constant pressure of the mineral.
 Returns
 molar_heat_capacity_pfloat
Heat capacity in [J/K/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
C_p()
.
 property P¶
Alias for
pressure()
 property T¶
Alias for
temperature()
 property energy¶
Alias for
molar_internal_energy()
 property helmholtz¶
Alias for
molar_helmholtz()
 property gibbs¶
Alias for
molar_gibbs()
 property V¶
Alias for
molar_volume()
 property S¶
Alias for
molar_entropy()
 property H¶
Alias for
molar_enthalpy()
 property K_T¶
Alias for
isothermal_bulk_modulus()
 property K_S¶
Alias for
adiabatic_bulk_modulus()
 property beta_T¶
Alias for
isothermal_compressibility()
 property beta_S¶
Alias for
adiabatic_compressibility()
 property G¶
Alias for
shear_modulus()
 property v_p¶
Alias for
p_wave_velocity()
 property v_phi¶
Alias for
bulk_sound_velocity()
 property v_s¶
Alias for
shear_wave_velocity()
 property gr¶
Alias for
grueneisen_parameter()
 property alpha¶
Alias for
thermal_expansivity()
 property C_v¶
Alias for
molar_heat_capacity_v()
 property C_p¶
Alias for
molar_heat_capacity_p()
Perple_X Class¶
 class burnman.PerplexMaterial(tab_file, name='Perple_X material')[source]¶
Bases:
Material
This is the base class for a PerpleX material. States of the material can only be queried after setting the pressure and temperature using set_state().
Instances of this class are initialised with a 2D PerpleX tab file. This file should be in the standard format (as output by werami), and should have columns with the following names: ‘rho,kg/m3’, ‘alpha,1/K’, ‘beta,1/bar’, ‘Ks,bar’, ‘Gs,bar’, ‘v0,km/s’, ‘vp,km/s’, ‘vs,km/s’, ‘s,J/K/kg’, ‘h,J/kg’, ‘cp,J/K/kg’, ‘V,J/bar/mol’. The order of these names is not important.
Properties of the material are determined by linear interpolation from the PerpleX grid. They are all returned in SI units on a molar basis, even though the PerpleX tab file is not in these units.
This class is available as
burnman.PerplexMaterial
. property name¶
Humanreadable name of this material.
By default this will return the name of the class, but it can be set to an arbitrary string. Overriden in Mineral.
 set_state()¶
(copied from set_state):
Set the material to the given pressure and temperature.
 Parameters
 pressurefloat
The desired pressure in [Pa].
 temperaturefloat
The desired temperature in [K].
 property molar_volume¶
Returns molar volume of the mineral.
 Returns
 molar_volumefloat
Molar volume in [m^3/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
V()
.
 property molar_enthalpy¶
Returns molar enthalpy of the mineral.
 Returns
 molar_enthalpyfloat
Enthalpy in [J/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
H()
.
 property molar_entropy¶
Returns molar entropy of the mineral.
 Returns
 molar_entropyfloat
Entropy in [J/K/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
S()
.
 property isothermal_bulk_modulus¶
Returns isothermal bulk modulus of the material.
 Returns
 isothermal_bulk_modulusfloat
Bulk modulus in [Pa].
Notes
Needs to be implemented in derived classes.
Aliased with
K_T()
.
 property adiabatic_bulk_modulus¶
Returns the adiabatic bulk modulus of the mineral.
 Returns
 adiabatic_bulk_modulusfloat
Adiabatic bulk modulus in [Pa].
Notes
Needs to be implemented in derived classes.
Aliased with
K_S()
.
 property molar_heat_capacity_p¶
Returns molar heat capacity at constant pressure of the mineral.
 Returns
 molar_heat_capacity_pfloat
Heat capacity in [J/K/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
C_p()
.
 property thermal_expansivity¶
Returns thermal expansion coefficient of the mineral.
 Returns
 alphafloat
Thermal expansivity in [1/K].
Notes
Needs to be implemented in derived classes.
Aliased with
alpha()
.
 property shear_modulus¶
Returns shear modulus of the mineral.
 Returns
 shear_modulusfloat
Shear modulus in [Pa].
Notes
Needs to be implemented in derived classes.
Aliased with
beta_G()
.
 property p_wave_velocity¶
Returns P wave speed of the mineral.
 Returns
 p_wave_velocityfloat
P wave speed in [m/s].
Notes
Needs to be implemented in derived classes.
Aliased with
v_p()
.
 property bulk_sound_velocity¶
Returns bulk sound speed of the mineral.
 Returns
 bulk sound velocity: float
Sound velocity in [m/s].
Notes
Needs to be implemented in derived classes.
Aliased with
v_phi()
.
 property shear_wave_velocity¶
Returns shear wave speed of the mineral.
 Returns
 shear_wave_velocityfloat
Wave speed in [m/s].
Notes
Needs to be implemented in derived classes.
Aliased with
v_s()
.
 property molar_gibbs¶
Returns the molar Gibbs free energy of the mineral.
 Returns
 molar_gibbsfloat
Gibbs free energy in [J/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
gibbs()
.
 property molar_mass¶
Returns molar mass of the mineral.
 Returns
 molar_massfloat
Molar mass in [kg/mol].
Notes
Needs to be implemented in derived classes.
 property density¶
Returns the density of this material.
 Returns
 densityfloat
The density of this material in [kg/m^3].
Notes
Needs to be implemented in derived classes.
Aliased with
rho()
.
 property molar_internal_energy¶
Returns the molar internal energy of the mineral.
 Returns
 molar_internal_energyfloat
The internal energy in [J/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
energy()
.
 property molar_helmholtz¶
Returns the molar Helmholtz free energy of the mineral.
 Returns
 molar_helmholtzfloat
Helmholtz free energy in [J/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
helmholtz()
.
 property isothermal_compressibility¶
Returns isothermal compressibility of the mineral (or inverse isothermal bulk modulus).
 Returns
 (K_T)^1float
Compressibility in [1/Pa].
Notes
Needs to be implemented in derived classes.
Aliased with
beta_T()
.
 property adiabatic_compressibility¶
Returns adiabatic compressibility of the mineral (or inverse adiabatic bulk modulus).
 Returns
 adiabatic_compressibilityfloat
adiabatic compressibility in [1/Pa].
Notes
Needs to be implemented in derived classes.
Aliased with
beta_S()
.
 property molar_heat_capacity_v¶
Returns molar heat capacity at constant volume of the mineral.
 Returns
 molar_heat_capacity_vfloat
Heat capacity in [J/K/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
C_v()
.
 property grueneisen_parameter¶
Returns the grueneisen parameter of the mineral.
 Returns
 grfloat
Grueneisen parameters [unitless].
Notes
Needs to be implemented in derived classes.
Aliased with
gr()
.
 property C_p¶
Alias for
molar_heat_capacity_p()
 property C_v¶
Alias for
molar_heat_capacity_v()
 property G¶
Alias for
shear_modulus()
 property H¶
Alias for
molar_enthalpy()
 property K_S¶
Alias for
adiabatic_bulk_modulus()
 property K_T¶
Alias for
isothermal_bulk_modulus()
 property P¶
Alias for
pressure()
 property S¶
Alias for
molar_entropy()
 property T¶
Alias for
temperature()
 property V¶
Alias for
molar_volume()
 property alpha¶
Alias for
thermal_expansivity()
 property beta_S¶
Alias for
adiabatic_compressibility()
 property beta_T¶
Alias for
isothermal_compressibility()
 copy()¶
 debug_print(indent='')¶
Print a humanreadable representation of this Material.
 property energy¶
Alias for
molar_internal_energy()
 evaluate(vars_list, pressures, temperatures)¶
Returns an array of material properties requested through a list of strings at given pressure and temperature conditions. At the end it resets the set_state to the original values. The user needs to call set_method() before.
 Parameters
 vars_listlist of strings
Variables to be returned for given conditions
 pressuresndlist or ndarray of float
ndimensional array of pressures in [Pa].
 temperaturesndlist or ndarray of float
ndimensional array of temperatures in [K].
 Returns
 outputarray of array of float
Array returning all variables at given pressure/temperature values. output[i][j] is property vars_list[j] and temperatures[i] and pressures[i].
 property gibbs¶
Alias for
molar_gibbs()
 property gr¶
Alias for
grueneisen_parameter()
 property helmholtz¶
Alias for
molar_helmholtz()
 property pressure¶
Returns current pressure that was set with
set_state()
. Returns
 pressurefloat
Pressure in [Pa].
Notes
Aliased with
P()
.
 print_minerals_of_current_state()¶
Print a humanreadable representation of this Material at the current P, T as a list of minerals. This requires set_state() has been called before.
 reset()¶
Resets all cached material properties.
It is typically not required for the user to call this function.
 set_method(method)¶
Set the averaging method. See Averaging Schemes for details.
Notes
Needs to be implemented in derived classes.
 set_state_with_volume(volume, temperature, pressure_guesses=[0.0, 10000000000.0])¶
This function acts similarly to set_state, but takes volume and temperature as input to find the pressure. In order to ensure selfconsistency, this function does not use any pressure functions from the material classes, but instead finds the pressure using the brentq rootfinding method.
 Parameters
 volumefloat
The desired molar volume of the mineral [m^3].
 temperaturefloat
The desired temperature of the mineral [K].
 pressure_guesseslist of floats (default: [5.e9, 10.e9])
The initial low and high guesses for bracketing of the pressure [Pa]. These guesses should preferably bound the correct pressure, but do not need to do so. More importantly, they should not lie outside the valid region of the equation of state.
 property temperature¶
Returns current temperature that was set with
set_state()
. Returns
 temperaturefloat
Temperature in [K].
Notes
Aliased with
T()
.
 to_string()¶
Returns a humanreadable name of this material. The default implementation will return the name of the class, which is a reasonable default.
 Returns
 namestring
Name of this material.
 unroll()¶
Unroll this material into a list of
burnman.Mineral
and their molar fractions. All averaging schemes then operate on this list of minerals. Note that the return value of this function may depend on the current state (temperature, pressure). Returns
 fractionslist of float
List of molar fractions, should sum to 1.0.
 mineralslist of
burnman.Mineral
List of minerals.
Notes
Needs to be implemented in derived classes.
 property v_p¶
Alias for
p_wave_velocity()
 property v_phi¶
Alias for
bulk_sound_velocity()
 property v_s¶
Alias for
shear_wave_velocity()
Minerals¶
Endmembers¶
 class burnman.Mineral(params=None, property_modifiers=None)[source]¶
Bases:
Material
This is the base class for all minerals. States of the mineral can only be queried after setting the pressure and temperature using set_state(). The method for computing properties of the material is set using set_method(). This is done during initialisation if the param ‘equation_of_state’ has been defined. The method can be overridden later by the user.
This class is available as
burnman.Mineral
.If deriving from this class, set the properties in self.params to the desired values. For more complicated materials you can overwrite set_state(), change the params and then call set_state() from this class.
All the material parameters are expected to be in plain SI units. This means that the elastic moduli should be in Pascals and NOT Gigapascals, and the Debye temperature should be in K not C. Additionally, the reference volume should be in m^3/(mol molecule) and not in unit cell volume and ‘n’ should be the number of atoms per molecule. Frequently in the literature the reference volume is given in Angstrom^3 per unit cell. To convert this to m^3/(mol of molecule) you should multiply by 10^(30) * N_a / Z, where N_a is Avogadro’s number and Z is the number of formula units per unit cell. You can look up Z in many places, including www.mindat.org
 property name¶
Humanreadable name of this material.
By default this will return the name of the class, but it can be set to an arbitrary string. Overriden in Mineral.
 set_method(equation_of_state)[source]¶
Set the equation of state to be used for this mineral. Takes a string corresponding to any of the predefined equations of state: ‘bm2’, ‘bm3’, ‘mgd2’, ‘mgd3’, ‘slb2’, ‘slb3’, ‘mt’, ‘hp_tmt’, or ‘cork’. Alternatively, you can pass a user defined class which derives from the equation_of_state base class. After calling set_method(), any existing derived properties (e.g., elastic parameters or thermodynamic potentials) will be out of date, so set_state() will need to be called again.
 unroll()[source]¶
Unroll this material into a list of
burnman.Mineral
and their molar fractions. All averaging schemes then operate on this list of minerals. Note that the return value of this function may depend on the current state (temperature, pressure). Returns
 fractionslist of float
List of molar fractions, should sum to 1.0.
 mineralslist of
burnman.Mineral
List of minerals.
Notes
Needs to be implemented in derived classes.
 set_state()¶
(copied from set_state):
Set the material to the given pressure and temperature.
 Parameters
 pressurefloat
The desired pressure in [Pa].
 temperaturefloat
The desired temperature in [K].
 property molar_gibbs¶
Returns the molar Gibbs free energy of the mineral.
 Returns
 molar_gibbsfloat
Gibbs free energy in [J/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
gibbs()
.
 property molar_volume¶
Returns molar volume of the mineral.
 Returns
 molar_volumefloat
Molar volume in [m^3/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
V()
.
 property molar_entropy¶
Returns molar entropy of the mineral.
 Returns
 molar_entropyfloat
Entropy in [J/K/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
S()
.
 property isothermal_bulk_modulus¶
Returns isothermal bulk modulus of the material.
 Returns
 isothermal_bulk_modulusfloat
Bulk modulus in [Pa].
Notes
Needs to be implemented in derived classes.
Aliased with
K_T()
.
 property molar_heat_capacity_p¶
Returns molar heat capacity at constant pressure of the mineral.
 Returns
 molar_heat_capacity_pfloat
Heat capacity in [J/K/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
C_p()
.
 property thermal_expansivity¶
Returns thermal expansion coefficient of the mineral.
 Returns
 alphafloat
Thermal expansivity in [1/K].
Notes
Needs to be implemented in derived classes.
Aliased with
alpha()
.
 property shear_modulus¶
Returns shear modulus of the mineral.
 Returns
 shear_modulusfloat
Shear modulus in [Pa].
Notes
Needs to be implemented in derived classes.
Aliased with
beta_G()
.
 property formula¶
Returns the chemical formula of the Mineral class
 property molar_mass¶
Returns molar mass of the mineral.
 Returns
 molar_massfloat
Molar mass in [kg/mol].
Notes
Needs to be implemented in derived classes.
 property density¶
Returns the density of this material.
 Returns
 densityfloat
The density of this material in [kg/m^3].
Notes
Needs to be implemented in derived classes.
Aliased with
rho()
.
 property molar_internal_energy¶
Returns the molar internal energy of the mineral.
 Returns
 molar_internal_energyfloat
The internal energy in [J/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
energy()
.
 property molar_helmholtz¶
Returns the molar Helmholtz free energy of the mineral.
 Returns
 molar_helmholtzfloat
Helmholtz free energy in [J/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
helmholtz()
.
 property molar_enthalpy¶
Returns molar enthalpy of the mineral.
 Returns
 molar_enthalpyfloat
Enthalpy in [J/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
H()
.
 property adiabatic_bulk_modulus¶
Returns the adiabatic bulk modulus of the mineral.
 Returns
 adiabatic_bulk_modulusfloat
Adiabatic bulk modulus in [Pa].
Notes
Needs to be implemented in derived classes.
Aliased with
K_S()
.
 property isothermal_compressibility¶
Returns isothermal compressibility of the mineral (or inverse isothermal bulk modulus).
 Returns
 (K_T)^1float
Compressibility in [1/Pa].
Notes
Needs to be implemented in derived classes.
Aliased with
beta_T()
.
 property adiabatic_compressibility¶
Returns adiabatic compressibility of the mineral (or inverse adiabatic bulk modulus).
 Returns
 adiabatic_compressibilityfloat
adiabatic compressibility in [1/Pa].
Notes
Needs to be implemented in derived classes.
Aliased with
beta_S()
.
 property p_wave_velocity¶
Returns P wave speed of the mineral.
 Returns
 p_wave_velocityfloat
P wave speed in [m/s].
Notes
Needs to be implemented in derived classes.
Aliased with
v_p()
.
 property bulk_sound_velocity¶
Returns bulk sound speed of the mineral.
 Returns
 bulk sound velocity: float
Sound velocity in [m/s].
Notes
Needs to be implemented in derived classes.
Aliased with
v_phi()
.
 property shear_wave_velocity¶
Returns shear wave speed of the mineral.
 Returns
 shear_wave_velocityfloat
Wave speed in [m/s].
Notes
Needs to be implemented in derived classes.
Aliased with
v_s()
.
 property C_p¶
Alias for
molar_heat_capacity_p()
 property C_v¶
Alias for
molar_heat_capacity_v()
 property G¶
Alias for
shear_modulus()
 property H¶
Alias for
molar_enthalpy()
 property K_S¶
Alias for
adiabatic_bulk_modulus()
 property K_T¶
Alias for
isothermal_bulk_modulus()
 property P¶
Alias for
pressure()
 property S¶
Alias for
molar_entropy()
 property T¶
Alias for
temperature()
 property V¶
Alias for
molar_volume()
 property alpha¶
Alias for
thermal_expansivity()
 property beta_S¶
Alias for
adiabatic_compressibility()
 property beta_T¶
Alias for
isothermal_compressibility()
 copy()¶
 property energy¶
Alias for
molar_internal_energy()
 evaluate(vars_list, pressures, temperatures)¶
Returns an array of material properties requested through a list of strings at given pressure and temperature conditions. At the end it resets the set_state to the original values. The user needs to call set_method() before.
 Parameters
 vars_listlist of strings
Variables to be returned for given conditions
 pressuresndlist or ndarray of float
ndimensional array of pressures in [Pa].
 temperaturesndlist or ndarray of float
ndimensional array of temperatures in [K].
 Returns
 outputarray of array of float
Array returning all variables at given pressure/temperature values. output[i][j] is property vars_list[j] and temperatures[i] and pressures[i].
 property gibbs¶
Alias for
molar_gibbs()
 property gr¶
Alias for
grueneisen_parameter()
 property grueneisen_parameter¶
Returns the grueneisen parameter of the mineral.
 Returns
 grfloat
Grueneisen parameters [unitless].
Notes
Needs to be implemented in derived classes.
Aliased with
gr()
.
 property helmholtz¶
Alias for
molar_helmholtz()
 property pressure¶
Returns current pressure that was set with
set_state()
. Returns
 pressurefloat
Pressure in [Pa].
Notes
Aliased with
P()
.
 print_minerals_of_current_state()¶
Print a humanreadable representation of this Material at the current P, T as a list of minerals. This requires set_state() has been called before.
 reset()¶
Resets all cached material properties.
It is typically not required for the user to call this function.
 set_state_with_volume(volume, temperature, pressure_guesses=[0.0, 10000000000.0])¶
This function acts similarly to set_state, but takes volume and temperature as input to find the pressure. In order to ensure selfconsistency, this function does not use any pressure functions from the material classes, but instead finds the pressure using the brentq rootfinding method.
 Parameters
 volumefloat
The desired molar volume of the mineral [m^3].
 temperaturefloat
The desired temperature of the mineral [K].
 pressure_guesseslist of floats (default: [5.e9, 10.e9])
The initial low and high guesses for bracketing of the pressure [Pa]. These guesses should preferably bound the correct pressure, but do not need to do so. More importantly, they should not lie outside the valid region of the equation of state.
 property temperature¶
Returns current temperature that was set with
set_state()
. Returns
 temperaturefloat
Temperature in [K].
Notes
Aliased with
T()
.
 property v_p¶
Alias for
p_wave_velocity()
 property v_phi¶
Alias for
bulk_sound_velocity()
 property v_s¶
Alias for
shear_wave_velocity()
Solutions¶
 class burnman.Solution(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)[source]¶
Bases:
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
burnman.Solution
. It uses an instance ofburnman.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
 namestring
Name of the solution
 solution_typestring
String determining which SolutionModel to use. One of ‘mechanical’, ‘ideal’, ‘symmetric’, ‘asymmetric’ or ‘subregular’.
 endmemberslist of lists
List of endmembers in this solution. The first item of each list should be a
burnman.Mineral
object. The second item should be a string with the site formula of the endmember. molar_fractionsnumpy array (optional)
The molar fractions of each endmember in the solution. Can be reset using the set_composition() method.
 property name¶
Humanreadable name of this material.
By default this will return the name of the class, but it can be set to an arbitrary string. Overriden in Mineral.
 set_composition(molar_fractions)[source]¶
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.
 set_method(method)[source]¶
Set the equation of state to be used for this mineral. Takes a string corresponding to any of the predefined equations of state: ‘bm2’, ‘bm3’, ‘mgd2’, ‘mgd3’, ‘slb2’, ‘slb3’, ‘mt’, ‘hp_tmt’, or ‘cork’. Alternatively, you can pass a user defined class which derives from the equation_of_state base class. After calling set_method(), any existing derived properties (e.g., elastic parameters or thermodynamic potentials) will be out of date, so set_state() will need to be called again.
 set_state(pressure, temperature)[source]¶
(copied from set_state):
Set the material to the given pressure and temperature.
 Parameters
 pressurefloat
The desired pressure in [Pa].
 temperaturefloat
The desired temperature in [K].
 property formula¶
Returns molar chemical formula of the solution.
 property activities¶
Returns a list of endmember activities [unitless].
 property activity_coefficients¶
Returns a list of endmember activity coefficients (gamma = activity / ideal activity) [unitless].
 property molar_internal_energy¶
Returns molar internal energy of the mineral [J/mol]. Aliased with self.energy
 property excess_partial_gibbs¶
Returns excess partial molar gibbs free energy [J/mol]. Property specific to solutions.
 property excess_partial_volumes¶
Returns excess partial volumes [m^3]. Property specific to solutions.
 property excess_partial_entropies¶
Returns excess partial entropies [J/K]. Property specific to solutions.
 property partial_gibbs¶
Returns endmember partial molar gibbs free energy [J/mol]. Property specific to solutions.
 property partial_volumes¶
Returns endmember partial volumes [m^3]. Property specific to solutions.
 property partial_entropies¶
Returns endmember partial entropies [J/K]. Property specific to solutions.
 property excess_gibbs¶
Returns molar excess gibbs free energy [J/mol]. Property specific to solutions.
 property gibbs_hessian¶
Returns an array containing the second compositional derivative of the Gibbs free energy [J]. Property specific to solutions.
 property entropy_hessian¶
Returns an array containing the second compositional derivative of the entropy [J/K]. Property specific to solutions.
 property volume_hessian¶
Returns an array containing the second compositional derivative of the volume [m^3]. Property specific to solutions.
 property molar_gibbs¶
Returns molar Gibbs free energy of the solution [J/mol]. Aliased with self.gibbs.
 property molar_helmholtz¶
Returns molar Helmholtz free energy of the solution [J/mol]. Aliased with self.helmholtz.
 property molar_mass¶
Returns molar mass of the solution [kg/mol].
 property excess_volume¶
Returns excess molar volume of the solution [m^3/mol]. Specific property for solutions.
 property molar_volume¶
Returns molar volume of the solution [m^3/mol]. Aliased with self.V.
 property density¶
Returns density of the solution [kg/m^3]. Aliased with self.rho.
 property excess_entropy¶
Returns excess molar entropy [J/K/mol]. Property specific to solutions.
 property molar_entropy¶
Returns molar entropy of the solution [J/K/mol]. Aliased with self.S.
 property excess_enthalpy¶
Returns excess molar enthalpy [J/mol]. Property specific to solutions.
 property molar_enthalpy¶
Returns molar enthalpy of the solution [J/mol]. Aliased with self.H.
 property isothermal_bulk_modulus¶
Returns isothermal bulk modulus of the solution [Pa]. Aliased with self.K_T.
 property adiabatic_bulk_modulus¶
Returns adiabatic bulk modulus of the solution [Pa]. Aliased with self.K_S.
 property isothermal_compressibility¶
Returns isothermal compressibility of the solution. (or inverse isothermal bulk modulus) [1/Pa]. Aliased with self.K_T.
 property adiabatic_compressibility¶
Returns adiabatic compressibility of the solution. (or inverse adiabatic bulk modulus) [1/Pa]. Aliased with self.K_S.
 property shear_modulus¶
Returns shear modulus of the solution [Pa]. Aliased with self.G.
 property p_wave_velocity¶
Returns P wave speed of the solution [m/s]. Aliased with self.v_p.
 property bulk_sound_velocity¶
Returns bulk sound speed of the solution [m/s]. Aliased with self.v_phi.
 property shear_wave_velocity¶
Returns shear wave speed of the solution [m/s]. Aliased with self.v_s.
 property grueneisen_parameter¶
Returns grueneisen parameter of the solution [unitless]. Aliased with self.gr.
 property thermal_expansivity¶
Returns thermal expansion coefficient (alpha) of the solution [1/K]. Aliased with self.alpha.
 property molar_heat_capacity_v¶
Returns molar heat capacity at constant volume of the solution [J/K/mol]. Aliased with self.C_v.
 property C_p¶
Alias for
molar_heat_capacity_p()
 property C_v¶
Alias for
molar_heat_capacity_v()
 property G¶
Alias for
shear_modulus()
 property H¶
Alias for
molar_enthalpy()
 property K_S¶
Alias for
adiabatic_bulk_modulus()
 property K_T¶
Alias for
isothermal_bulk_modulus()
 property P¶
Alias for
pressure()
 property S¶
Alias for
molar_entropy()
 property T¶
Alias for
temperature()
 property V¶
Alias for
molar_volume()
 property alpha¶
Alias for
thermal_expansivity()
 property beta_S¶
Alias for
adiabatic_compressibility()
 property beta_T¶
Alias for
isothermal_compressibility()
 copy()¶
 debug_print(indent='')¶
Print a humanreadable representation of this Material.
 property energy¶
Alias for
molar_internal_energy()
 evaluate(vars_list, pressures, temperatures)¶
Returns an array of material properties requested through a list of strings at given pressure and temperature conditions. At the end it resets the set_state to the original values. The user needs to call set_method() before.
 Parameters
 vars_listlist of strings
Variables to be returned for given conditions
 pressuresndlist or ndarray of float
ndimensional array of pressures in [Pa].
 temperaturesndlist or ndarray of float
ndimensional array of temperatures in [K].
 Returns
 outputarray of array of float
Array returning all variables at given pressure/temperature values. output[i][j] is property vars_list[j] and temperatures[i] and pressures[i].
 property gibbs¶
Alias for
molar_gibbs()
 property gr¶
Alias for
grueneisen_parameter()
 property helmholtz¶
Alias for
molar_helmholtz()
 property molar_heat_capacity_p¶
Returns molar heat capacity at constant pressure of the solution [J/K/mol]. Aliased with self.C_p.
 property pressure¶
Returns current pressure that was set with
set_state()
. Returns
 pressurefloat
Pressure in [Pa].
Notes
Aliased with
P()
.
 print_minerals_of_current_state()¶
Print a humanreadable representation of this Material at the current P, T as a list of minerals. This requires set_state() has been called before.
 reset()¶
Resets all cached material properties.
It is typically not required for the user to call this function.
 set_state_with_volume(volume, temperature, pressure_guesses=[0.0, 10000000000.0])¶
This function acts similarly to set_state, but takes volume and temperature as input to find the pressure. In order to ensure selfconsistency, this function does not use any pressure functions from the material classes, but instead finds the pressure using the brentq rootfinding method.
 Parameters
 volumefloat
The desired molar volume of the mineral [m^3].
 temperaturefloat
The desired temperature of the mineral [K].
 pressure_guesseslist of floats (default: [5.e9, 10.e9])
The initial low and high guesses for bracketing of the pressure [Pa]. These guesses should preferably bound the correct pressure, but do not need to do so. More importantly, they should not lie outside the valid region of the equation of state.
 property temperature¶
Returns current temperature that was set with
set_state()
. Returns
 temperaturefloat
Temperature in [K].
Notes
Aliased with
T()
.
 to_string()¶
Returns the name of the mineral class
 unroll()¶
Unroll this material into a list of
burnman.Mineral
and their molar fractions. All averaging schemes then operate on this list of minerals. Note that the return value of this function may depend on the current state (temperature, pressure). Returns
 fractionslist of float
List of molar fractions, should sum to 1.0.
 mineralslist of
burnman.Mineral
List of minerals.
Notes
Needs to be implemented in derived classes.
 property v_p¶
Alias for
p_wave_velocity()
 property v_phi¶
Alias for
bulk_sound_velocity()
 property v_s¶
Alias for
shear_wave_velocity()
 property stoichiometric_matrix[source]¶
A sympy Matrix where each element M[i,j] corresponds to the number of atoms of element[j] in endmember[i].
 property stoichiometric_array[source]¶
An array where each element arr[i,j] corresponds to the number of atoms of element[j] in endmember[i].
 property reaction_basis[source]¶
An array where each element arr[i,j] corresponds to the number of moles of endmember[j] involved in reaction[i].
 property independent_element_indices[source]¶
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).
 property dependent_element_indices[source]¶
The element indices not included in the independent list.
 class burnman.ElasticSolution(name=None, solution_type=None, endmembers=None, energy_interaction=None, pressure_interaction=None, entropy_interaction=None, energy_ternary_terms=None, pressure_ternary_terms=None, entropy_ternary_terms=None, alphas=None, excess_helmholtz_function=None, molar_fractions=None)[source]¶
Bases:
Mineral
This is the base class for all Elastic solutions. Site occupancies, endmember activities and the constant and volume 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() and set_composition.
This class is available as
burnman.ElasticSolution
. It uses an instance ofburnman.ElasticSolutionModel
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 V derivatives in J/K/mol and Pa/mol.
The parameters are relevant to all Elastic solution models. Please see the documentation for individual models for details about other parameters.
 Parameters
 namestring
Name of the solution
 solution_typestring
String determining which SolutionModel to use. One of ‘mechanical’, ‘ideal’, ‘symmetric’, ‘asymmetric’ or ‘subregular’.
 endmemberslist of lists
List of endmembers in this solution. The first item of each list should be a
burnman.Mineral
object. The second item should be a string with the site formula of the endmember. molar_fractionsnumpy array (optional)
The molar fractions of each endmember in the solution. Can be reset using the set_composition() method.
 property name¶
Humanreadable name of this material.
By default this will return the name of the class, but it can be set to an arbitrary string. Overriden in Mineral.
 set_composition(molar_fractions)[source]¶
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.
 set_method(method)[source]¶
Set the equation of state to be used for this mineral. Takes a string corresponding to any of the predefined equations of state: ‘bm2’, ‘bm3’, ‘mgd2’, ‘mgd3’, ‘slb2’, ‘slb3’, ‘mt’, ‘hp_tmt’, or ‘cork’. Alternatively, you can pass a user defined class which derives from the equation_of_state base class. After calling set_method(), any existing derived properties (e.g., elastic parameters or thermodynamic potentials) will be out of date, so set_state() will need to be called again.
 set_state(pressure, temperature)[source]¶
(copied from set_state):
Set the material to the given pressure and temperature.
 Parameters
 pressurefloat
The desired pressure in [Pa].
 temperaturefloat
The desired temperature in [K].
 property formula¶
Returns molar chemical formula of the solution.
 property activities¶
Returns a list of endmember activities [unitless].
 property activity_coefficients¶
Returns a list of endmember activity coefficients (gamma = activity / ideal activity) [unitless].
 property molar_internal_energy¶
Returns molar internal energy of the mineral [J/mol]. Aliased with self.energy
 property partial_gibbs¶
Returns endmember partial molar Gibbs energy at constant pressure [J/mol]. Property specific to solutions.
 property partial_volumes¶
Returns endmember partial molar volumes [m^3/mol]. Property specific to solutions.
 property partial_entropies¶
Returns endmember partial molar entropies [J/K/mol]. Property specific to solutions.
 property gibbs_hessian¶
Returns an array containing the second compositional derivative of the Gibbs energy at constant pressure [J/mol]. Property specific to solutions.
 property molar_helmholtz¶
Returns molar Helmholtz energy of the solution [J/mol]. Aliased with self.helmholtz.
 property molar_gibbs¶
Returns molar Gibbs free energy of the solution [J/mol]. Aliased with self.gibbs.
 property molar_mass¶
Returns molar mass of the solution [kg/mol].
 property excess_pressure¶
Returns excess pressure of the solution [Pa]. Specific property for solutions.
 property molar_volume¶
Returns molar volume of the solution [m^3/mol]. Aliased with self.V.
 property density¶
Returns density of the solution [kg/m^3]. Aliased with self.rho.
 property excess_entropy¶
Returns excess molar entropy [J/K/mol]. Property specific to solutions.
 property molar_entropy¶
Returns molar entropy of the solution [J/K/mol]. Aliased with self.S.
 property excess_enthalpy¶
Returns excess molar enthalpy [J/mol]. Property specific to solutions.
 property molar_enthalpy¶
Returns molar enthalpy of the solution [J/mol]. Aliased with self.H.
 property isothermal_bulk_modulus¶
Returns isothermal bulk modulus of the solution [Pa]. Aliased with self.K_T.
 property adiabatic_bulk_modulus¶
Returns adiabatic bulk modulus of the solution [Pa]. Aliased with self.K_S.
 property isothermal_compressibility¶
Returns isothermal compressibility of the solution. (or inverse isothermal bulk modulus) [1/Pa]. Aliased with self.K_T.
 property adiabatic_compressibility¶
Returns adiabatic compressibility of the solution. (or inverse adiabatic bulk modulus) [1/Pa]. Aliased with self.K_S.
 property shear_modulus¶
Returns shear modulus of the solution [Pa]. Aliased with self.G.
 property p_wave_velocity¶
Returns P wave speed of the solution [m/s]. Aliased with self.v_p.
 property bulk_sound_velocity¶
Returns bulk sound speed of the solution [m/s]. Aliased with self.v_phi.
 property shear_wave_velocity¶
Returns shear wave speed of the solution [m/s]. Aliased with self.v_s.
 property grueneisen_parameter¶
Returns grueneisen parameter of the solution [unitless]. Aliased with self.gr.
 property thermal_expansivity¶
Returns thermal expansion coefficient (alpha) of the solution [1/K]. Aliased with self.alpha.
 property molar_heat_capacity_v¶
Returns molar heat capacity at constant volume of the solution [J/K/mol]. Aliased with self.C_v.
 property C_p¶
Alias for
molar_heat_capacity_p()
 property C_v¶
Alias for
molar_heat_capacity_v()
 property G¶
Alias for
shear_modulus()
 property H¶
Alias for
molar_enthalpy()
 property K_S¶
Alias for
adiabatic_bulk_modulus()
 property K_T¶
Alias for
isothermal_bulk_modulus()
 property P¶
Alias for
pressure()
 property S¶
Alias for
molar_entropy()
 property T¶
Alias for
temperature()
 property V¶
Alias for
molar_volume()
 property alpha¶
Alias for
thermal_expansivity()
 property beta_S¶
Alias for
adiabatic_compressibility()
 property beta_T¶
Alias for
isothermal_compressibility()
 copy()¶
 debug_print(indent='')¶
Print a humanreadable representation of this Material.
 property energy¶
Alias for
molar_internal_energy()
 evaluate(vars_list, pressures, temperatures)¶
Returns an array of material properties requested through a list of strings at given pressure and temperature conditions. At the end it resets the set_state to the original values. The user needs to call set_method() before.
 Parameters
 vars_listlist of strings
Variables to be returned for given conditions
 pressuresndlist or ndarray of float
ndimensional array of pressures in [Pa].
 temperaturesndlist or ndarray of float
ndimensional array of temperatures in [K].
 Returns
 outputarray of array of float
Array returning all variables at given pressure/temperature values. output[i][j] is property vars_list[j] and temperatures[i] and pressures[i].
 property gibbs¶
Alias for
molar_gibbs()
 property gr¶
Alias for
grueneisen_parameter()
 property helmholtz¶
Alias for
molar_helmholtz()
 property molar_heat_capacity_p¶
Returns molar heat capacity at constant pressure of the solution [J/K/mol]. Aliased with self.C_p.
 property pressure¶
Returns current pressure that was set with
set_state()
. Returns
 pressurefloat
Pressure in [Pa].
Notes
Aliased with
P()
.
 print_minerals_of_current_state()¶
Print a humanreadable representation of this Material at the current P, T as a list of minerals. This requires set_state() has been called before.
 reset()¶
Resets all cached material properties.
It is typically not required for the user to call this function.
 set_state_with_volume(volume, temperature, pressure_guesses=[0.0, 10000000000.0])¶
This function acts similarly to set_state, but takes volume and temperature as input to find the pressure. In order to ensure selfconsistency, this function does not use any pressure functions from the material classes, but instead finds the pressure using the brentq rootfinding method.
 Parameters
 volumefloat
The desired molar volume of the mineral [m^3].
 temperaturefloat
The desired temperature of the mineral [K].
 pressure_guesseslist of floats (default: [5.e9, 10.e9])
The initial low and high guesses for bracketing of the pressure [Pa]. These guesses should preferably bound the correct pressure, but do not need to do so. More importantly, they should not lie outside the valid region of the equation of state.
 property temperature¶
Returns current temperature that was set with
set_state()
. Returns
 temperaturefloat
Temperature in [K].
Notes
Aliased with
T()
.
 to_string()¶
Returns the name of the mineral class
 unroll()¶
Unroll this material into a list of
burnman.Mineral
and their molar fractions. All averaging schemes then operate on this list of minerals. Note that the return value of this function may depend on the current state (temperature, pressure). Returns
 fractionslist of float
List of molar fractions, should sum to 1.0.
 mineralslist of
burnman.Mineral
List of minerals.
Notes
Needs to be implemented in derived classes.
 property v_p¶
Alias for
p_wave_velocity()
 property v_phi¶
Alias for
bulk_sound_velocity()
 property v_s¶
Alias for
shear_wave_velocity()
 property stoichiometric_matrix[source]¶
A sympy Matrix where each element M[i,j] corresponds to the number of atoms of element[j] in endmember[i].
 property stoichiometric_array[source]¶
An array where each element arr[i,j] corresponds to the number of atoms of element[j] in endmember[i].
 property reaction_basis[source]¶
An array where each element arr[i,j] corresponds to the number of moles of endmember[j] involved in reaction[i].
 property independent_element_indices[source]¶
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).
 property dependent_element_indices[source]¶
The element indices not included in the independent list.
 burnman.ElasticSolidSolution¶
alias of
ElasticSolution
Mineral helpers¶
 class burnman.classes.mineral_helpers.HelperSpinTransition(transition_pressure, ls_mat, hs_mat)[source]¶
Bases:
Composite
Helper class that makes a mineral that switches between two materials (for low and high spin) based on some transition pressure [Pa]
 set_state(pressure, temperature)[source]¶
Update the material to the given pressure [Pa] and temperature [K].
 property C_p¶
Alias for
molar_heat_capacity_p()
 property C_v¶
Alias for
molar_heat_capacity_v()
 property G¶
Alias for
shear_modulus()
 property H¶
Alias for
molar_enthalpy()
 property K_S¶
Alias for
adiabatic_bulk_modulus()
 property K_T¶
Alias for
isothermal_bulk_modulus()
 property P¶
Alias for
pressure()
 property S¶
Alias for
molar_entropy()
 property T¶
Alias for
temperature()
 property V¶
Alias for
molar_volume()
 property adiabatic_bulk_modulus¶
Returns adiabatic bulk modulus of the mineral [Pa] Aliased with self.K_S
 property adiabatic_compressibility¶
Returns isothermal compressibility of the composite (or inverse isothermal bulk modulus) [1/Pa] Aliased with self.beta_S
 property alpha¶
Alias for
thermal_expansivity()
 property beta_S¶
Alias for
adiabatic_compressibility()
 property beta_T¶
Alias for
isothermal_compressibility()
 property bulk_sound_velocity¶
Returns bulk sound speed of the composite [m/s] Aliased with self.v_phi
 chemical_potential(components=None)¶
Returns the chemical potentials of the currently defined components in the composite. Raises an exception if the assemblage is not equilibrated.
 Parameters
 components: list of dictionaries (optional)
List of formulae of the desired components. If not specified, the method uses the components specified by a previous call to set_components.
 Returns
 chemical_potential: numpy array of floats
The chemical potentials of the desired components in the equilibrium composite.
 property compositional_null_basis¶
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 composite.
 copy()¶
 property density¶
Compute the density of the composite based on the molar volumes and masses Aliased with self.rho
 property dependent_element_indices¶
The element indices not included in the independent list.
 property elements¶
A list of the elements which could be contained in the composite, returned in the IUPAC element order.
 property endmember_formulae¶
A list of the formulae in the composite.
 property endmember_names¶
A list of the endmember names contained in the composite. Mineral names are returned as given in Mineral.name. Solution endmember names are given in the format Mineral.name in Solution.name.
 property endmember_partial_gibbs¶
Returns the partial Gibbs energies for all the endmember minerals in the Composite
 property endmembers_per_phase¶
A list of integers corresponding to the number of endmembers stored within each phase.
 property energy¶
Alias for
molar_internal_energy()
 property equilibrated¶
Returns True if the reaction affinities are all zero within a given tolerance given by self.equilibrium_tolerance.
 evaluate(vars_list, pressures, temperatures)¶
Returns an array of material properties requested through a list of strings at given pressure and temperature conditions. At the end it resets the set_state to the original values. The user needs to call set_method() before.
 Parameters
 vars_listlist of strings
Variables to be returned for given conditions
 pressuresndlist or ndarray of float
ndimensional array of pressures in [Pa].
 temperaturesndlist or ndarray of float
ndimensional array of temperatures in [K].
 Returns
 outputarray of array of float
Array returning all variables at given pressure/temperature values. output[i][j] is property vars_list[j] and temperatures[i] and pressures[i].
 property formula¶
Returns molar chemical formula of the composite
 property gibbs¶
Alias for
molar_gibbs()
 property gr¶
Alias for
grueneisen_parameter()
 property grueneisen_parameter¶
Returns grueneisen parameter of the composite [unitless] Aliased with self.gr
 property helmholtz¶
Alias for
molar_helmholtz()
 property independent_element_indices¶
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)
 property isothermal_bulk_modulus¶
Returns isothermal bulk modulus of the composite [Pa] Aliased with self.K_T
 property isothermal_compressibility¶
Returns isothermal compressibility of the composite (or inverse isothermal bulk modulus) [1/Pa] Aliased with self.beta_T
 property molar_enthalpy¶
Returns enthalpy of the mineral [J] Aliased with self.H
 property molar_entropy¶
Returns enthalpy of the mineral [J] Aliased with self.S
 property molar_gibbs¶
Returns molar Gibbs free energy of the composite [J/mol] Aliased with self.gibbs
 property molar_heat_capacity_p¶
Returns molar heat capacity at constant pressure of the composite [J/K/mol] Aliased with self.C_p
 property molar_heat_capacity_v¶
Returns molar_heat capacity at constant volume of the composite [J/K/mol] Aliased with self.C_v
 property molar_helmholtz¶
Returns molar Helmholtz free energy of the mineral [J/mol] Aliased with self.helmholtz
 property molar_internal_energy¶
Returns molar internal energy of the mineral [J/mol] Aliased with self.energy
 property molar_mass¶
Returns molar mass of the composite [kg/mol]
 property molar_volume¶
Returns molar volume of the composite [m^3/mol] Aliased with self.V
 property n_elements¶
Returns the total number of distinct elements which might be in the composite.
 property n_endmembers¶
Returns the number of endmembers in the composite.
 property n_reactions¶
The number of reactions in reaction_basis.
 property name¶
Humanreadable name of this material.
By default this will return the name of the class, but it can be set to an arbitrary string. Overriden in Mineral.
 property p_wave_velocity¶
Returns P wave speed of the composite [m/s] Aliased with self.v_p
 property pressure¶
Returns current pressure that was set with
set_state()
. Returns
 pressurefloat
Pressure in [Pa].
Notes
Aliased with
P()
.
 print_minerals_of_current_state()¶
Print a humanreadable representation of this Material at the current P, T as a list of minerals. This requires set_state() has been called before.
 property reaction_affinities¶
Returns the affinities corresponding to each reaction in reaction_basis
 property reaction_basis¶
An array where each element arr[i,j] corresponds to the number of moles of endmember[j] involved in reaction[i].
 property reaction_basis_as_strings¶
Returns a list of string representations of all the reactions in reaction_basis.
 property reduced_stoichiometric_array¶
The stoichiometric array including only the independent elements
 reset()¶
Resets all cached material properties.
It is typically not required for the user to call this function.
 set_averaging_scheme(averaging_scheme)¶
Set the averaging scheme for the moduli in the composite. Default is set to VoigtReussHill, when Composite is initialized.
 set_components(components)¶
Sets the components and components_array attributes of the composite material. The components attribute is a list of dictionaries containing the chemical formulae of the components. The components_array attribute is a 2D numpy array describing the linear transformation between endmember amounts and component amounts.
The components do not need to be linearly independent, not do they need to form a complete basis set for the composite. However, it must be possible to obtain the composition of each component from a linear sum of the endmember compositions of the composite. For example, if the composite was composed of MgSiO3 and Mg2SiO4, SiO2 would be a valid component, but Si would not. The method raises an exception if any of the chemical potentials are not defined by the assemblage.
 Parameters
 components: list of dictionaries
List of formulae of the components.
 set_fractions(fractions, fraction_type='molar')¶
Change the fractions of the phases of this Composite. Resets cached properties
 Parameters
 fractions: list or numpy array of floats
molar or mass fraction for each phase.
 fraction_type: ‘molar’ or ‘mass’
specify whether molar or mass fractions are specified.
 set_method(method)¶
set the same equation of state method for all the phases in the composite
 set_state_with_volume(volume, temperature, pressure_guesses=[0.0, 10000000000.0])¶
This function acts similarly to set_state, but takes volume and temperature as input to find the pressure. In order to ensure selfconsistency, this function does not use any pressure functions from the material classes, but instead finds the pressure using the brentq rootfinding method.
 Parameters
 volumefloat
The desired molar volume of the mineral [m^3].
 temperaturefloat
The desired temperature of the mineral [K].
 pressure_guesseslist of floats (default: [5.e9, 10.e9])
The initial low and high guesses for bracketing of the pressure [Pa]. These guesses should preferably bound the correct pressure, but do not need to do so. More importantly, they should not lie outside the valid region of the equation of state.
 property shear_modulus¶
Returns shear modulus of the mineral [Pa] Aliased with self.G
 property shear_wave_velocity¶
Returns shear wave speed of the composite [m/s] Aliased with self.v_s
 property stoichiometric_array¶
An array where each element arr[i,j] corresponds to the number of atoms of element[j] in endmember[i].
 property stoichiometric_matrix¶
An sympy Matrix where each element M[i,j] corresponds to the number of atoms of element[j] in endmember[i].
 property temperature¶
Returns current temperature that was set with
set_state()
. Returns
 temperaturefloat
Temperature in [K].
Notes
Aliased with
T()
.
 property thermal_expansivity¶
Returns thermal expansion coefficient of the composite [1/K] Aliased with self.alpha
 to_string()¶
return the name of the composite
 unroll()¶
Unroll this material into a list of
burnman.Mineral
and their molar fractions. All averaging schemes then operate on this list of minerals. Note that the return value of this function may depend on the current state (temperature, pressure). Returns
 fractionslist of float
List of molar fractions, should sum to 1.0.
 mineralslist of
burnman.Mineral
List of minerals.
Notes
Needs to be implemented in derived classes.
 property v_p¶
Alias for
p_wave_velocity()
 property v_phi¶
Alias for
bulk_sound_velocity()
 property v_s¶
Alias for
shear_wave_velocity()
Anisotropic materials¶
 class burnman.AnisotropicMaterial(rho, cijs)[source]¶
Bases:
Material
A base class for anisotropic elastic materials. The base class is initialised with a density and a full isentropic stiffness tensor in Voigt notation. It can then be interrogated to find the values of different properties, such as bounds on seismic velocities. There are also several functions which can be called to calculate properties along directions oriented with respect to the isentropic elastic tensor.
See [MHS11] and https://docs.materialsproject.org/methodology/elasticity/ for mathematical descriptions of each function.
 property isentropic_stiffness_tensor¶
 property full_isentropic_stiffness_tensor¶
 property isentropic_compliance_tensor¶
 property full_isentropic_compliance_tensor¶
 property density¶
Returns the density of this material.
 Returns
 densityfloat
The density of this material in [kg/m^3].
Notes
Needs to be implemented in derived classes.
Aliased with
rho()
.
 property isentropic_bulk_modulus_voigt¶
Computes the isentropic bulk modulus (Voigt bound)
 property isentropic_bulk_modulus_reuss¶
Computes the isentropic bulk modulus (Reuss bound)
 property isentropic_bulk_modulus_vrh¶
Computes the isentropic bulk modulus (VoigtReussHill average)
 property isentropic_shear_modulus_voigt¶
Computes the isentropic shear modulus (Voigt bound)
 property isentropic_shear_modulus_reuss¶
Computes the isentropic shear modulus (Reuss bound)
 property isentropic_shear_modulus_vrh¶
Computes the shear modulus (VoigtReussHill average)
 property isentropic_universal_elastic_anisotropy¶
Compute the universal elastic anisotropy
 property isentropic_isotropic_poisson_ratio¶
Compute mu, the isotropic Poisson ratio (a description of the laterial response to loading)
 christoffel_tensor(propagation_direction)[source]¶
Computes the Christoffel tensor from an elastic stiffness tensor and a propagation direction for a seismic wave relative to the stiffness tensor
T_ik = C_ijkl n_j n_l
 isentropic_linear_compressibility(direction)[source]¶
Computes the linear isentropic compressibility in a given direction relative to the stiffness tensor
 isentropic_youngs_modulus(direction)[source]¶
Computes the isentropic Youngs modulus in a given direction relative to the stiffness tensor
 isentropic_shear_modulus(plane_normal, shear_direction)[source]¶
Computes the isentropic shear modulus on a plane in a given shear direction relative to the stiffness tensor
 isentropic_poissons_ratio(axial_direction, lateral_direction)[source]¶
Computes the isentropic poisson ratio given loading and response directions relative to the stiffness tensor
 wave_velocities(propagation_direction)[source]¶
Computes the compressional wave velocity, and two shear wave velocities in a given propagation direction
Returns two lists, containing the wave speeds and directions of particle motion relative to the stiffness tensor
 property C_p¶
Alias for
molar_heat_capacity_p()
 property C_v¶
Alias for
molar_heat_capacity_v()
 property G¶
Alias for
shear_modulus()
 property H¶
Alias for
molar_enthalpy()
 property K_S¶
Alias for
adiabatic_bulk_modulus()
 property K_T¶
Alias for
isothermal_bulk_modulus()
 property P¶
Alias for
pressure()
 property S¶
Alias for
molar_entropy()
 property T¶
Alias for
temperature()
 property V¶
Alias for
molar_volume()
 property adiabatic_bulk_modulus¶
Returns the adiabatic bulk modulus of the mineral.
 Returns
 adiabatic_bulk_modulusfloat
Adiabatic bulk modulus in [Pa].
Notes
Needs to be implemented in derived classes.
Aliased with
K_S()
.
 property adiabatic_compressibility¶
Returns adiabatic compressibility of the mineral (or inverse adiabatic bulk modulus).
 Returns
 adiabatic_compressibilityfloat
adiabatic compressibility in [1/Pa].
Notes
Needs to be implemented in derived classes.
Aliased with
beta_S()
.
 property alpha¶
Alias for
thermal_expansivity()
 property beta_S¶
Alias for
adiabatic_compressibility()
 property beta_T¶
Alias for
isothermal_compressibility()
 property bulk_sound_velocity¶
Returns bulk sound speed of the mineral.
 Returns
 bulk sound velocity: float
Sound velocity in [m/s].
Notes
Needs to be implemented in derived classes.
Aliased with
v_phi()
.
 copy()¶
 debug_print(indent='')¶
Print a humanreadable representation of this Material.
 property energy¶
Alias for
molar_internal_energy()
 evaluate(vars_list, pressures, temperatures)¶
Returns an array of material properties requested through a list of strings at given pressure and temperature conditions. At the end it resets the set_state to the original values. The user needs to call set_method() before.
 Parameters
 vars_listlist of strings
Variables to be returned for given conditions
 pressuresndlist or ndarray of float
ndimensional array of pressures in [Pa].
 temperaturesndlist or ndarray of float
ndimensional array of temperatures in [K].
 Returns
 outputarray of array of float
Array returning all variables at given pressure/temperature values. output[i][j] is property vars_list[j] and temperatures[i] and pressures[i].
 property gibbs¶
Alias for
molar_gibbs()
 property gr¶
Alias for
grueneisen_parameter()
 property grueneisen_parameter¶
Returns the grueneisen parameter of the mineral.
 Returns
 grfloat
Grueneisen parameters [unitless].
Notes
Needs to be implemented in derived classes.
Aliased with
gr()
.
 property helmholtz¶
Alias for
molar_helmholtz()
 property isothermal_bulk_modulus¶
Returns isothermal bulk modulus of the material.
 Returns
 isothermal_bulk_modulusfloat
Bulk modulus in [Pa].
Notes
Needs to be implemented in derived classes.
Aliased with
K_T()
.
 property isothermal_compressibility¶
Returns isothermal compressibility of the mineral (or inverse isothermal bulk modulus).
 Returns
 (K_T)^1float
Compressibility in [1/Pa].
Notes
Needs to be implemented in derived classes.
Aliased with
beta_T()
.
 property molar_enthalpy¶
Returns molar enthalpy of the mineral.
 Returns
 molar_enthalpyfloat
Enthalpy in [J/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
H()
.
 property molar_entropy¶
Returns molar entropy of the mineral.
 Returns
 molar_entropyfloat
Entropy in [J/K/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
S()
.
 property molar_gibbs¶
Returns the molar Gibbs free energy of the mineral.
 Returns
 molar_gibbsfloat
Gibbs free energy in [J/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
gibbs()
.
 property molar_heat_capacity_p¶
Returns molar heat capacity at constant pressure of the mineral.
 Returns
 molar_heat_capacity_pfloat
Heat capacity in [J/K/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
C_p()
.
 property molar_heat_capacity_v¶
Returns molar heat capacity at constant volume of the mineral.
 Returns
 molar_heat_capacity_vfloat
Heat capacity in [J/K/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
C_v()
.
 property molar_helmholtz¶
Returns the molar Helmholtz free energy of the mineral.
 Returns
 molar_helmholtzfloat
Helmholtz free energy in [J/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
helmholtz()
.
 property molar_internal_energy¶
Returns the molar internal energy of the mineral.
 Returns
 molar_internal_energyfloat
The internal energy in [J/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
energy()
.
 property molar_mass¶
Returns molar mass of the mineral.
 Returns
 molar_massfloat
Molar mass in [kg/mol].
Notes
Needs to be implemented in derived classes.
 property molar_volume¶
Returns molar volume of the mineral.
 Returns
 molar_volumefloat
Molar volume in [m^3/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
V()
.
 property name¶
Humanreadable name of this material.
By default this will return the name of the class, but it can be set to an arbitrary string. Overriden in Mineral.
 property p_wave_velocity¶
Returns P wave speed of the mineral.
 Returns
 p_wave_velocityfloat
P wave speed in [m/s].
Notes
Needs to be implemented in derived classes.
Aliased with
v_p()
.
 property pressure¶
Returns current pressure that was set with
set_state()
. Returns
 pressurefloat
Pressure in [Pa].
Notes
Aliased with
P()
.
 print_minerals_of_current_state()¶
Print a humanreadable representation of this Material at the current P, T as a list of minerals. This requires set_state() has been called before.
 reset()¶
Resets all cached material properties.
It is typically not required for the user to call this function.
 set_method(method)¶
Set the averaging method. See Averaging Schemes for details.
Notes
Needs to be implemented in derived classes.
 set_state(pressure, temperature)¶
Set the material to the given pressure and temperature.
 Parameters
 pressurefloat
The desired pressure in [Pa].
 temperaturefloat
The desired temperature in [K].
 set_state_with_volume(volume, temperature, pressure_guesses=[0.0, 10000000000.0])¶
This function acts similarly to set_state, but takes volume and temperature as input to find the pressure. In order to ensure selfconsistency, this function does not use any pressure functions from the material classes, but instead finds the pressure using the brentq rootfinding method.
 Parameters
 volumefloat
The desired molar volume of the mineral [m^3].
 temperaturefloat
The desired temperature of the mineral [K].
 pressure_guesseslist of floats (default: [5.e9, 10.e9])
The initial low and high guesses for bracketing of the pressure [Pa]. These guesses should preferably bound the correct pressure, but do not need to do so. More importantly, they should not lie outside the valid region of the equation of state.
 property shear_modulus¶
Returns shear modulus of the mineral.
 Returns
 shear_modulusfloat
Shear modulus in [Pa].
Notes
Needs to be implemented in derived classes.
Aliased with
beta_G()
.
 property shear_wave_velocity¶
Returns shear wave speed of the mineral.
 Returns
 shear_wave_velocityfloat
Wave speed in [m/s].
Notes
Needs to be implemented in derived classes.
Aliased with
v_s()
.
 property temperature¶
Returns current temperature that was set with
set_state()
. Returns
 temperaturefloat
Temperature in [K].
Notes
Aliased with
T()
.
 property thermal_expansivity¶
Returns thermal expansion coefficient of the mineral.
 Returns
 alphafloat
Thermal expansivity in [1/K].
Notes
Needs to be implemented in derived classes.
Aliased with
alpha()
.
 to_string()¶
Returns a humanreadable name of this material. The default implementation will return the name of the class, which is a reasonable default.
 Returns
 namestring
Name of this material.
 unroll()¶
Unroll this material into a list of
burnman.Mineral
and their molar fractions. All averaging schemes then operate on this list of minerals. Note that the return value of this function may depend on the current state (temperature, pressure). Returns
 fractionslist of float
List of molar fractions, should sum to 1.0.
 mineralslist of
burnman.Mineral
List of minerals.
Notes
Needs to be implemented in derived classes.
 property v_p¶
Alias for
p_wave_velocity()
 property v_phi¶
Alias for
bulk_sound_velocity()
 property v_s¶
Alias for
shear_wave_velocity()
 class burnman.AnisotropicMineral(isotropic_mineral, cell_parameters, anisotropic_parameters, psi_function=None, orthotropic=None)[source]¶
Bases:
Mineral
,AnisotropicMaterial
A class implementing the anisotropic mineral equation of state described in [Myhill22]. This class is derived from both Mineral and AnisotropicMaterial, and inherits most of the methods from these classes.
Instantiation of an AnisotropicMineral takes three required arguments; a reference Mineral (i.e. a standard isotropic mineral which provides volume as a function of pressure and temperature), cell_parameters, which give the lengths of the molar cell vectors and the angles between them (see
cell_parameters_to_vectors()
), and an anisotropic parameters object, which should be either a 4D array of anisotropic parameters or a dictionary of parameters which describe the anisotropic behaviour of the mineral. For a description of the physical meaning of the parameters in the 4D array, please refer to the code or to the original paper.If the user chooses to define their parameters as a dictionary, they must also provide a function to the psi_function argument that describes how to compute the tensors Psi, dPsidf and dPsidPth (in Voigt form). The function arguments should be f, Pth and params, in that order. The output variables Psi, dPsidf and dPsidth must be returned in that order in a tuple. The user should also explicitly state whether the material is orthotropic or not by supplying a boolean to the orthotropic argument.
States of the mineral can only be queried after setting the pressure and temperature using set_state().
This class is available as
burnman.AnisotropicMineral
.All the material parameters are expected to be in plain SI units. This means that the elastic moduli should be in Pascals and NOT Gigapascals. Additionally, the cell parameters should be in m/(mol formula unit) and not in unit cell lengths. To convert unit cell lengths given in Angstrom to molar cell parameters you should multiply by 10^(10) * (N_a / Z)^1/3, where N_a is Avogadro’s number and Z is the number of formula units per unit cell. You can look up Z in many places, including www.mindat.org.
Finally, it is assumed that the unit cell of the anisotropic material is aligned in a particular way relative to the coordinate axes (the anisotropic_parameters are defined relative to the coordinate axes). The crystallographic aaxis is assumed to be parallel to the first spatial coordinate axis, and the crystallographic baxis is assumed to be perpendicular to the third spatial coordinate axis.
 property name¶
Humanreadable name of this material.
By default this will return the name of the class, but it can be set to an arbitrary string. Overriden in Mineral.
 set_state()¶
(copied from set_state):
Set the material to the given pressure and temperature.
 Parameters
 pressurefloat
The desired pressure in [Pa].
 temperaturefloat
The desired temperature in [K].
 property deformation_gradient_tensor¶
 Returns
 deformation_gradient_tensor2D numpy array
The deformation gradient tensor describing the deformation of the mineral from its undeformed state (i.e. the state at the reference pressure and temperature).
 property unrotated_cell_vectors¶
 Returns
 unrotated_cell_vectors2D numpy array
The vectors of the cell constructed from one mole of formula units after deformation of the mineral from its undeformed state (i.e. the state at the reference pressure and temperature). Each vector is given in [m]. See the documentation for the function
cell_parameters_to_vectors()
for the assumed relationships between the cell vectors and spatial coordinate axes.
 property deformed_coordinate_frame¶
 Returns
 deformed_coordinate_frame2D numpy array
The orientations of the three spatial coordinate axes after deformation of the mineral. For orthotropic minerals, this is equal to the identity matrix, as hydrostatic stresses only induce rotations in monoclinic and triclinic crystals.
 property rotation_matrix¶
 Returns
 rotation_matrix2D numpy array
The matrix required to rotate the properties of the deformed mineral into the deformed coordinate frame. For orthotropic minerals, this is equal to the identity matrix.
 property cell_vectors¶
 Returns
 cell_vectors2D numpy array
The vectors of the cell constructed from one mole of formula units. Each vector is given in [m]. See the documentation for the function
cell_parameters_to_vectors()
for the assumed relationships between the cell vectors and spatial coordinate axes.
 property cell_parameters¶
 Returns
 cell_parameters1D numpy array
The molar cell parameters of the mineral, given in standard form: [\(a\), \(b\), \(c\), \(\alpha\), \(\beta\), \(\gamma\)], where the first three floats are the lengths of the vectors in [m] defining the cell constructed from one mole of formula units. The last three floats are angles between vectors (given in radians). See the documentation for the function
cell_parameters_to_vectors()
for the assumed relationships between the cell vectors and spatial coordinate axes.
 property shear_modulus¶
Anisotropic minerals do not (in general) have a single shear modulus. This function returns a NotImplementedError. Users should instead consider directly querying the elements in the isothermal_stiffness_tensor or isentropic_stiffness_tensor.
 property isothermal_bulk_modulus¶
Anisotropic minerals do not have a single isothermal bulk modulus. This function returns a NotImplementedError. Users should instead consider either using isothermal_bulk_modulus_reuss, isothermal_bulk_modulus_voigt, or directly querying the elements in the isothermal_stiffness_tensor.
 property isentropic_bulk_modulus¶
Anisotropic minerals do not have a single isentropic bulk modulus. This function returns a NotImplementedError. Users should instead consider either using isentropic_bulk_modulus_reuss, isentropic_bulk_modulus_voigt (both derived from Anisotropicmineral), or directly querying the elements in the isentropic_stiffness_tensor.
 property isothermal_bulk_modulus_reuss¶
Returns isothermal bulk modulus of the material.
 Returns
 isothermal_bulk_modulusfloat
Bulk modulus in [Pa].
Notes
Needs to be implemented in derived classes.
Aliased with
K_T()
.
 property isothermal_bulk_modulus_voigt¶
 Returns
 isothermal_bulk_modulus_voigtfloat
The Voigt bound on the isothermal bulk modulus in [Pa].
 property isothermal_compressibility_reuss¶
 Returns
 isothermal_compressibility_reussfloat
The Reuss bound on the isothermal compressibility in [1/Pa].
 property isothermal_compliance_tensor¶
 Returns
 isothermal_compliance_tensor2D numpy array
The isothermal compliance tensor [1/Pa] in Voigt form (\(\mathbb{S}_{\text{T} pq}\)).
 property thermal_expansivity_tensor¶
 Returns
 thermal_expansivity_tensor2D numpy array
The tensor of thermal expansivities [1/K].
 property isothermal_stiffness_tensor¶
 Returns
 isothermal_stiffness_tensor2D numpy array
The isothermal stiffness tensor [Pa] in Voigt form (\(\mathbb{C}_{\text{T} pq}\)).
 property full_isothermal_compliance_tensor¶
 Returns
 full_isothermal_stiffness_tensor4D numpy array
The isothermal compliance tensor [1/Pa] in standard form (\(\mathbb{S}_{\text{T} ijkl}\)).
 property full_isothermal_stiffness_tensor¶
 Returns
 full_isothermal_stiffness_tensor4D numpy array
The isothermal stiffness tensor [Pa] in standard form (\(\mathbb{C}_{\text{T} ijkl}\)).
 property full_isentropic_compliance_tensor¶
 Returns
 full_isentropic_stiffness_tensor4D numpy array
The isentropic compliance tensor [1/Pa] in standard form (\(\mathbb{S}_{\text{N} ijkl}\)).
 property isentropic_compliance_tensor¶
 Returns
 isentropic_compliance_tensor2D numpy array
The isentropic compliance tensor [1/Pa] in Voigt form (\(\mathbb{S}_{\text{N} pq}\)).
 property isentropic_stiffness_tensor¶
 Returns
 isentropic_stiffness_tensor2D numpy array
The isentropic stiffness tensor [Pa] in Voigt form (\(\mathbb{C}_{\text{N} pq}\)).
 property full_isentropic_stiffness_tensor¶
 Returns
 full_isentropic_stiffness_tensor4D numpy array
The isentropic stiffness tensor [Pa] in standard form (\(\mathbb{C}_{\text{N} ijkl}\)).
 property grueneisen_tensor¶
 Returns
 grueneisen_tensor2D numpy array
The grueneisen tensor. This is defined by [BM67] as \(\mathbb{C}_{\text{N} ijkl} \alpha_{kl} V/C_{P}\).
 property grueneisen_parameter¶
Anisotropic minerals do not (in general) have a single grueneisen parameter. This function returns a NotImplementedError. Users should instead consider directly querying the elements in the grueneisen_tensor.
 property isothermal_compressibility_tensor¶
 Returns
 isothermal_compressibility_tensor2D numpy array
The isothermal compressibility tensor.
 property isentropic_compressibility_tensor¶
 Returns
 isentropic_compressibility_tensor2D numpy array
The isentropic compressibility tensor.
 property thermal_stress_tensor¶
 Returns
 thermal stress2D numpy array
The change in stress with temperature at constant strain.
 property C_p¶
Alias for
molar_heat_capacity_p()
 property C_v¶
Alias for
molar_heat_capacity_v()
 property G¶
Alias for
shear_modulus()
 property H¶
Alias for
molar_enthalpy()
 property K_S¶
Alias for
adiabatic_bulk_modulus()
 property K_T¶
Alias for
isothermal_bulk_modulus()
 property P¶
Alias for
pressure()
 property S¶
Alias for
molar_entropy()
 property T¶
Alias for
temperature()
 property V¶
Alias for
molar_volume()
 property adiabatic_bulk_modulus¶
Returns the adiabatic bulk modulus of the mineral.
 Returns
 adiabatic_bulk_modulusfloat
Adiabatic bulk modulus in [Pa].
Notes
Needs to be implemented in derived classes.
Aliased with
K_S()
.
 property adiabatic_compressibility¶
Returns adiabatic compressibility of the mineral (or inverse adiabatic bulk modulus).
 Returns
 adiabatic_compressibilityfloat
adiabatic compressibility in [1/Pa].
Notes
Needs to be implemented in derived classes.
Aliased with
beta_S()
.
 property alpha¶
Alias for
thermal_expansivity()
 property beta_S¶
Alias for
adiabatic_compressibility()
 property beta_T¶
Alias for
isothermal_compressibility()
 property bulk_sound_velocity¶
Returns bulk sound speed of the mineral.
 Returns
 bulk sound velocity: float
Sound velocity in [m/s].
Notes
Needs to be implemented in derived classes.
Aliased with
v_phi()
.
 christoffel_tensor(propagation_direction)¶
Computes the Christoffel tensor from an elastic stiffness tensor and a propagation direction for a seismic wave relative to the stiffness tensor
T_ik = C_ijkl n_j n_l
 copy()¶
 debug_print(indent='')¶
Print a humanreadable representation of this Material.
 property density¶
Returns the density of this material.
 Returns
 densityfloat
The density of this material in [kg/m^3].
Notes
Needs to be implemented in derived classes.
Aliased with
rho()
.
 property energy¶
Alias for
molar_internal_energy()
 evaluate(vars_list, pressures, temperatures)¶
Returns an array of material properties requested through a list of strings at given pressure and temperature conditions. At the end it resets the set_state to the original values. The user needs to call set_method() before.
 Parameters
 vars_listlist of strings
Variables to be returned for given conditions
 pressuresndlist or ndarray of float
ndimensional array of pressures in [Pa].
 temperaturesndlist or ndarray of float
ndimensional array of temperatures in [K].
 Returns
 outputarray of array of float
Array returning all variables at given pressure/temperature values. output[i][j] is property vars_list[j] and temperatures[i] and pressures[i].
 property formula¶
Returns the chemical formula of the Mineral class
 property gibbs¶
Alias for
molar_gibbs()
 property gr¶
Alias for
grueneisen_parameter()
 property helmholtz¶
Alias for
molar_helmholtz()
 property isentropic_bulk_modulus_reuss¶
Computes the isentropic bulk modulus (Reuss bound)
 property isentropic_bulk_modulus_voigt¶
Computes the isentropic bulk modulus (Voigt bound)
 property isentropic_bulk_modulus_vrh¶
Computes the isentropic bulk modulus (VoigtReussHill average)
 property isentropic_isotropic_poisson_ratio¶
Compute mu, the isotropic Poisson ratio (a description of the laterial response to loading)
 isentropic_linear_compressibility(direction)¶
Computes the linear isentropic compressibility in a given direction relative to the stiffness tensor
 isentropic_poissons_ratio(axial_direction, lateral_direction)¶
Computes the isentropic poisson ratio given loading and response directions relative to the stiffness tensor
 isentropic_shear_modulus(plane_normal, shear_direction)¶
Computes the isentropic shear modulus on a plane in a given shear direction relative to the stiffness tensor
 property isentropic_shear_modulus_reuss¶
Computes the isentropic shear modulus (Reuss bound)
 property isentropic_shear_modulus_voigt¶
Computes the isentropic shear modulus (Voigt bound)
 property isentropic_shear_modulus_vrh¶
Computes the shear modulus (VoigtReussHill average)
 property isentropic_universal_elastic_anisotropy¶
Compute the universal elastic anisotropy
 isentropic_youngs_modulus(direction)¶
Computes the isentropic Youngs modulus in a given direction relative to the stiffness tensor
 property isothermal_compressibility¶
Returns isothermal compressibility of the mineral (or inverse isothermal bulk modulus).
 Returns
 (K_T)^1float
Compressibility in [1/Pa].
Notes
Needs to be implemented in derived classes.
Aliased with
beta_T()
.
 property molar_enthalpy¶
Returns molar enthalpy of the mineral.
 Returns
 molar_enthalpyfloat
Enthalpy in [J/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
H()
.
 property molar_entropy¶
Returns molar entropy of the mineral.
 Returns
 molar_entropyfloat
Entropy in [J/K/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
S()
.
 property molar_gibbs¶
Returns the molar Gibbs free energy of the mineral.
 Returns
 molar_gibbsfloat
Gibbs free energy in [J/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
gibbs()
.
 property molar_heat_capacity_p¶
Returns molar heat capacity at constant pressure of the mineral.
 Returns
 molar_heat_capacity_pfloat
Heat capacity in [J/K/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
C_p()
.
 property molar_heat_capacity_v¶
Returns molar heat capacity at constant volume of the mineral.
 Returns
 molar_heat_capacity_vfloat
Heat capacity in [J/K/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
C_v()
.
 property molar_helmholtz¶
Returns the molar Helmholtz free energy of the mineral.
 Returns
 molar_helmholtzfloat
Helmholtz free energy in [J/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
helmholtz()
.
 property molar_internal_energy¶
Returns the molar internal energy of the mineral.
 Returns
 molar_internal_energyfloat
The internal energy in [J/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
energy()
.
 property molar_isometric_heat_capacity¶
 Returns
 molar_isometric_heat_capacityfloat
The molar heat capacity at constant strain.
 property molar_mass¶
Returns molar mass of the mineral.
 Returns
 molar_massfloat
Molar mass in [kg/mol].
Notes
Needs to be implemented in derived classes.
 property molar_volume¶
Returns molar volume of the mineral.
 Returns
 molar_volumefloat
Molar volume in [m^3/mol].
Notes
Needs to be implemented in derived classes.
Aliased with
V()
.
 property p_wave_velocity¶
Returns P wave speed of the mineral.
 Returns
 p_wave_velocityfloat
P wave speed in [m/s].
Notes
Needs to be implemented in derived classes.
Aliased with
v_p()
.
 property pressure¶
Returns current pressure that was set with
set_state()
. Returns
 pressurefloat
Pressure in [Pa].
Notes
Aliased with
P()
.
 print_minerals_of_current_state()¶
Print a humanreadable representation of this Material at the current P, T as a list of minerals. This requires set_state() has been called before.
 reset()¶
Resets all cached material properties.
It is typically not required for the user to call this function.
 set_method(equation_of_state)¶
Set the equation of state to be used for this mineral. Takes a string corresponding to any of the predefined equations of state: ‘bm2’, ‘bm3’, ‘mgd2’, ‘mgd3’, ‘slb2’, ‘slb3’, ‘mt’, ‘hp_tmt’, or ‘cork’. Alternatively, you can pass a user defined class which derives from the equation_of_state base class. After calling set_method(), any existing derived properties (e.g., elastic parameters or thermodynamic potentials) will be out of date, so set_state() will need to be called again.
 set_state_with_volume(volume, temperature, pressure_guesses=[0.0, 10000000000.0])¶
This function acts similarly to set_state, but takes volume and temperature as input to find the pressure. In order to ensure selfconsistency, this function does not use any pressure functions from the material classes, but instead finds the pressure using the brentq rootfinding method.
 Parameters
 volumefloat
The desired molar volume of the mineral [m^3].
 temperaturefloat
The desired temperature of the mineral [K].
 pressure_guesseslist of floats (default: [5.e9, 10.e9])
The initial low and high guesses for bracketing of the pressure [Pa]. These guesses should preferably bound the correct pressure, but do not need to do so. More importantly, they should not lie outside the valid region of the equation of state.
 property shear_wave_velocity¶
Returns shear wave speed of the mineral.
 Returns
 shear_wave_velocityfloat
Wave speed in [m/s].
Notes
Needs to be implemented in derived classes.
Aliased with
v_s()
.
 property temperature¶
Returns current temperature that was set with
set_state()
. Returns
 temperaturefloat
Temperature in [K].
Notes
Aliased with
T()
.
 property thermal_expansivity¶
Returns thermal expansion coefficient of the mineral.
 Returns
 alphafloat
Thermal expansivity in [1/K].
Notes
Needs to be implemented in derived classes.
Aliased with
alpha()
.
 to_string()¶
Returns the name of the mineral class
 unroll()¶
Unroll this material into a list of
burnman.Mineral
and their molar fractions. All averaging schemes then operate on this list of minerals. Note that the return value of this function may depend on the current state (temperature, pressure). Returns
 fractionslist of float
List of molar fractions, should sum to 1.0.
 mineralslist of
burnman.Mineral
List of minerals.
Notes
Needs to be implemented in derived classes.
 property v_p¶
Alias for
p_wave_velocity()
 property v_phi¶
Alias for
bulk_sound_velocity()
 property v_s¶
Alias for
shear_wave_velocity()
 wave_velocities(propagation_direction)¶
Computes the compressional wave velocity, and two shear wave velocities in a given propagation direction
Returns two lists, containing the wave speeds and directions of particle motion relative to the stiffness tensor
 burnman.cell_parameters_to_vectors(cell_parameters)[source]¶
Converts cell parameters to unit cell vectors.
 Parameters
 cell_parameters1D numpy array
An array containing the three lengths of the unit cell vectors [m], and the three angles [degrees]. The first angle (\(\alpha\)) corresponds to the angle between the second and the third cell vectors, the second (\(\beta\)) to the angle between the first and third cell vectors, and the third (\(\gamma\)) to the angle between the first and second vectors.
 Returns
 M2D numpy array
The three vectors defining the parallelopiped cell [m]. This function assumes that the first cell vector is colinear with the xaxis, and the second is perpendicular to the zaxis, and the third is defined in a righthanded sense.
 burnman.cell_vectors_to_parameters(M)[source]¶
Converts unit cell vectors to cell parameters.
 Parameters
 M2D numpy array
The three vectors defining the parallelopiped cell [m]. This function assumes that the first cell vector is colinear with the xaxis, the second is perpendicular to the zaxis, and the third is defined in a righthanded sense.
 Returns
 cell_parameters1D numpy array
An array containing the three lengths of the unit cell vectors [m], and the three angles [degrees]. The first angle (\(\alpha\)) corresponds to the angle between the second and the third cell vectors, the second (\(\beta\)) to the angle between the first and third cell vectors, and the third (\(\gamma\)) to the angle between the first and second vectors.
Composites¶
 class burnman.Composite(phases, fractions=None, fraction_type='molar', name='Unnamed composite')[source]¶
Bases:
Material
Base class for a composite material. The static phases can be minerals or materials, meaning composite can be nested arbitrarily.
The fractions of the phases can be input as either ‘molar’ or ‘mass’ during instantiation, and modified (or initialised) after this point by using set_fractions.
This class is available as
burnman.Composite
. property name¶
Humanreadable name of this material.
By default this will return the name of the class, but it can be set to an arbitrary string. Overriden in Mineral.
 set_fractions(fractions, fraction_type='molar')[source]¶
Change the fractions of the phases of this Composite. Resets cached properties
 Parameters
 fractions: list or numpy array of floats
molar or mass fraction for each phase.
 fraction_type: ‘molar’ or ‘mass’
specify whether molar or mass fractions are specified.
 set_method(method)[source]¶
set the same equation of state method for all the phases in the composite
 set_averaging_scheme(averaging_scheme)[source]¶
Set the averaging scheme for the moduli in the composite. Default is set to VoigtReussHill, when Composite is initialized.
 set_state(pressure, temperature)[source]¶
Update the material to the given pressure [Pa] and temperature [K].
 unroll()[source]¶
Unroll this material into a list of
burnman.Mineral
and their molar fractions. All averaging schemes then operate on this list of minerals. Note that the return value of this function may depend on the current state (temperature, pressure). Returns
 fractionslist of float
List of molar fractions, should sum to 1.0.
 mineralslist of
burnman.Mineral
List of minerals.
Notes
Needs to be implemented in derived classes.
 property formula¶
Returns molar chemical formula of the composite
 property molar_internal_energy¶
Returns molar internal energy of the mineral [J/mol] Aliased with self.energy
 property molar_gibbs¶
Returns molar Gibbs free energy of the composite [J/mol] Aliased with self.gibbs
 property molar_helmholtz¶
Returns molar Helmholtz free energy of the mineral [J/mol] Aliased with self.helmholtz
 property molar_volume¶
Returns molar volume of the composite [m^3/mol] Aliased with self.V
 property molar_mass¶
Returns molar mass of the composite [kg/mol]
 property density¶
Compute the density of the composite based on the molar volumes and masses Aliased with self.rho
 property molar_entropy¶
Returns enthalpy of the mineral [J] Aliased with self.S
 property molar_enthalpy¶
Returns enthalpy of the mineral [J] Aliased with self.H
 property isothermal_bulk_modulus¶
Returns isothermal bulk modulus of the composite [Pa] Aliased with self.K_T
 property adiabatic_bulk_modulus¶
Returns adiabatic bulk modulus of the mineral [Pa] Aliased with self.K_S
 property isothermal_compressibility¶
Returns isothermal compressibility of the composite (or inverse isothermal bulk modulus) [1/Pa] Aliased with self.beta_T
 property adiabatic_compressibility¶
Returns isothermal compressibility of the composite (or inverse isothermal bulk modulus) [1/Pa] Aliased with self.beta_S
 property shear_modulus¶
Returns shear modulus of the mineral [Pa] Aliased with self.G
 property p_wave_velocity¶
Returns P wave speed of the composite [m/s] Aliased with self.v_p
 property bulk_sound_velocity¶
Returns bulk sound speed of the composite [m/s] Aliased with self.v_phi
 property shear_wave_velocity¶
Returns shear wave speed of the composite [m/s] Aliased with self.v_s
 property grueneisen_parameter¶
Returns grueneisen parameter of the composite [unitless] Aliased with self.gr
 property thermal_expansivity¶
Returns thermal expansion coefficient of the composite [1/K] Aliased with self.alpha
 property molar_heat_capacity_v¶
Returns molar_heat capacity at constant volume of the composite [J/K/mol] Aliased with self.C_v
 property molar_heat_capacity_p¶
Returns molar heat capacity at constant pressure of the composite [J/K/mol] Aliased with self.C_p
 property endmember_partial_gibbs¶
Returns the partial Gibbs energies for all the endmember minerals in the Composite
 property reaction_affinities¶
Returns the affinities corresponding to each reaction in reaction_basis
 property equilibrated¶
Returns True if the reaction affinities are all zero within a given tolerance given by self.equilibrium_tolerance.
 set_components(components)[source]¶
Sets the components and components_array attributes of the composite material. The components attribute is a list of dictionaries containing the chemical formulae of the components. The components_array attribute is a 2D numpy array describing the linear transformation between endmember amounts and component amounts.
The components do not need to be linearly independent, not do they need to form a complete basis set for the composite. However, it must be possible to obtain the composition of each component from a linear sum of the endmember compositions of the composite. For example, if the composite was composed of MgSiO3 and Mg2SiO4, SiO2 would be a valid component, but Si would not. The method raises an exception if any of the chemical potentials are not defined by the assemblage.
 Parameters
 components: list of dictionaries
List of formulae of the components.
 chemical_potential(components=None)[source]¶
Returns the chemical potentials of the currently defined components in the composite. Raises an exception if the assemblage is not equilibrated.
 Parameters
 components: list of dictionaries (optional)
List of formulae of the desired components. If not specified, the method uses the components specified by a previous call to set_components.
 Returns
 chemical_potential: numpy array of floats
The chemical potentials of the desired components in the equilibrium composite.
 property stoichiometric_matrix[source]¶
An sympy Matrix where each element M[i,j] corresponds to the number of atoms of element[j] in endmember[i].
 property stoichiometric_array[source]¶
An array where each element arr[i,j] corresponds to the number of atoms of element[j] in endmember[i].
 property reaction_basis[source]¶
An array where each element arr[i,j] corresponds to the number of moles of endmember[j] involved in reaction[i].
 property reaction_basis_as_strings[source]¶
Returns a list of string representations of all the reactions in reaction_basis.
 property independent_element_indices[source]¶
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)
 property dependent_element_indices[source]¶
The element indices not included in the independent list.
 property reduced_stoichiometric_array[source]¶
The stoichiometric array including only the independent elements
 property compositional_null_basis[source]¶
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 composite.
 property endmember_names[source]¶
A list of the endmember names contained in the composite. Mineral names are returned as given in Mineral.name. Solution endmember names are given in the format Mineral.name in Solution.name.
 property endmembers_per_phase[source]¶
A list of integers corresponding to the number of endmembers stored within each phase.
 property elements[source]¶
A list of the elements which could be contained in the composite, returned in the IUPAC element order.
 property n_elements[source]¶
Returns the total number of distinct elements which might be in the composite.
 property C_p¶
Alias for
molar_heat_capacity_p()
 property C_v¶
Alias for
molar_heat_capacity_v()
 property G¶
Alias for
shear_modulus()
 property H¶
Alias for
molar_enthalpy()
 property K_S¶
Alias for
adiabatic_bulk_modulus()
 property K_T¶
Alias for
isothermal_bulk_modulus()
 property P¶
Alias for
pressure()
 property S¶
Alias for
molar_entropy()
 property T¶
Alias for
temperature()
 property V¶
Alias for
molar_volume()
 property alpha¶
Alias for
thermal_expansivity()
 property beta_S¶
Alias for
adiabatic_compressibility()
 property beta_T¶
Alias for
isothermal_compressibility()
 copy()¶
 property energy¶
Alias for
molar_internal_energy()
 evaluate(vars_list, pressures, temperatures)¶
Returns an array of material properties requested through a list of strings at given pressure and temperature conditions. At the end it resets the set_state to the original values. The user needs to call set_method() before.
 Parameters
 vars_listlist of strings
Variables to be returned for given conditions
 pressuresndlist or ndarray of float
ndimensional array of pressures in [Pa].
 temperaturesndlist or ndarray of float
ndimensional array of temperatures in [K].
 Returns
 outputarray of array of float
Array returning all variables at given pressure/temperature values. output[i][j] is property vars_list[j] and temperatures[i] and pressures[i].
 property gibbs¶
Alias for
molar_gibbs()