Tools

Burnman has a number of high-level tools to help achieve common goals. Several of these have already been described in previous sections.

Unit cell

Equilibration

burnman.tools.equilibration.calculate_constraints(assemblage, n_free_compositional_vectors)

This function calculates the linear inequality constraints bounding the valid parameter space for a given assemblage.

The constraints are as follows:

• Pressure and temperature must be positive

• All phase fractions must be positive

• All site-species occupancies must be positive

The constraints are stored in a vector (b) and matrix (A). The sign convention is chosen such that the constraint is satisfied if A.x + b < eps.

Parameters:

assemblage (burnman.Composite) – The assemblage for which the constraints are calculated.

Returns:

The constraints vector and matrix.

Return type:

tuple

burnman.tools.equilibration.get_parameters(assemblage, n_free_compositional_vectors=0)

Gets the starting parameters vector (x) for the current equilibrium problem. These are:

• pressure

• temperature

• absolute amount of each phase. if a phase is a solution with >1 endmember, the following parameters are the mole fractions of the independent endmembers in the solution, except for the first endmember (as the mole fractions must sum to one).

Parameters:

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

Returns:

The current values of all the parameters.

Return type:

numpy.array

burnman.tools.equilibration.get_endmember_amounts(assemblage)

Gets the absolute amounts of all the endmembers in the solution.

Parameters:

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

Returns:

The current amounts of all the endmembers.

Return type:

numpy.array

burnman.tools.equilibration.set_compositions_and_state_from_parameters(assemblage, parameters)

Sets the phase compositions, amounts and state of the assemblage from a list of parameter values.

Parameters:

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

• parameters (numpy.array) – The current parameter values.

burnman.tools.equilibration.default_F_tolerances(assemblage, equality_constraints, n_atoms)

Returns the default tolerances for each component of F(x).

burnman.tools.equilibration.F(x, assemblage, equality_constraints, reduced_composition_vector, reduced_free_composition_vectors)

The vector-valued function for which the root is sought. The first two vector values depend on the equality_constraints chosen. For example, if

• eq[i][0] = ‘P’, F[i] = P - eq[i][1]

• eq[i][0] = ‘T’, F[i] = T - eq[i][1]

• eq[i][0] = ‘S’, F[i] = entropy - eq[i][1]

• eq[i][0] = ‘V’, F[i] = volume - eq[i][1]

• eq[i][0] = ‘X’, np.dot(eq[i][1][0], x) - eq[i][1][1]

The next set of vector values correspond to the reaction affinities. The final set of vector values correspond to the bulk composition constraints.

Parameters:

• x (numpy array) – Parameter values for the equilibrium problem to be solved.

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

• equality_constraints (list of lists) – A list of the equality constraints (see above for valid formats).

• reduced_composition_vector (numpy.array) – The amounts of the independent elements.

• reduced_free_composition_vectors (2D numpy.array) – The amounts of the independent elements in each of the free_compositional_vectors.

Returns:

The vector corresponding to F(x).

Return type:

numpy.array

burnman.tools.equilibration.jacobian(x, assemblage, equality_constraints, reduced_free_composition_vectors)

The Jacobian of the vector-valued function \(F\) for which the root is sought (\(\partial F / \partial x\)). See documentation for F() and get_parameters() (which return \(F\) and \(x\) respectively) for more details.

Parameters:

• x (numpy.array) – Parameter values for the equilibrium problem to be solved.

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

• equality_constraints (list of lists) – A list of the equality constraints (see documentation for burnman.tools.equilbration.F()).

• reduced_free_composition_vectors (2D numpy array) – The amounts of the independent elements in each of the free_compositional_vectors.

Returns:

The Jacobian for the equilibrium problem.

Return type:

2D numpy.array

burnman.tools.equilibration.lambda_bounds(dx, x, endmembers_per_phase)

Returns the lambda bounds for the damped affine invariant modification to Newton’s method for nonlinear problems (Deuflhard, 1974;1975;2004).

Parameters:

• dx (numpy.array) – The proposed newton step.

• x (numpy.array) – Parameter values for the equilibrium problem to be solved.

• endmembers_per_phase (list of int) – A list of the number of endmembers in each phase.

Returns:

Minimum and maximum allowed fractions of the full newton step (dx).

Return type:

tuple of floats

burnman.tools.equilibration.phase_fraction_constraints(phase, assemblage, fractions, prm)

Converts a phase fraction constraint into standard linear form that can be processed by the root finding problem.

We start with a single fraction or an array of fractions for a particular phase (\(n_p / \sum n = f\)). These are then converted into the “X” form of constraint by multiplying by \(\sum n\) and moving all terms to the LHS of the equation:

\(-f n_0 - f n_1 - \ldots + (1-f) n_p - \ldots = 0\)

This form is less readable, but easier to use as a constraint in a nonlinear solve.

Parameters:

• phase (burnman.Solution or burnman.Mineral) – The phase for which the fraction is to be constrained

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

• fractions (numpy.array) – The phase fractions to be satified at equilibrium.

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

Returns:

The phase fraction constraints.

Return type:

list

burnman.tools.equilibration.phase_composition_constraints(phase, assemblage, constraints, prm)

Converts a phase composition constraint into standard linear form that can be processed by the root finding problem.

We start with constraints in the form (site_names, n, d, v), where \((n x)/ (d x) = v\) 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 function converts the user-defined vectors of site constraints \(n\) and \(d\) into vectors of endmember proportion constraints \(n'\) and \(d'\), such that \((n' x)/ (d' x) = v\). This is done via linear transformation using the site occupancy matrix provided by burnman.Solution. By multiplying by the denominator, we have the following scalar comparison: \((n' x) = v (d' x)\)

The equilibration function does not use the proportion of the first endmember (as the endmember proportions must sum to one), and so we split \(x\), \(n'\) and \(d'\) into the first element and following elements: \((n'_0 x_0 + n'_i x_i) = v (d'_0 x_0 + d'_i x_i)\) where \(i\) is taken over all elements apart from the first.

With some more rearranging we can express the constraint in standard linear form: \((n'_0 (1 - \sum_j x_j) + n'_i x_i) = v (d'_0 (1 - \sum_j x_j) + d'_i x_i)\)

\((n'_0 + (n'_i - 1_i n'_0) x_i) = v (d'_0 + (d'_i - 1_i d'_0) x_i)\)

\((((n'_i - 1_i n'_0) - v(d'_i - 1_i d'_0)) x_i) = (v d'_0 - n'_0)\)

This form is less readable, but easier to use as a constraint in a nonlinear solve.

Parameters:

• phase (burnman.Solution) – The phase for which the composition is to be constrained.

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

• constraints (tuple) – The desired constraints in the form: site_names (list of strings), numerator (numpy.array), denominator (numpy.array), values (numpy.array).

Returns:

The phase composition constraints in standard form.

Return type:

list

burnman.tools.equilibration.get_equilibration_parameters(assemblage, composition, free_compositional_vectors)

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

Parameters:

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

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

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

Returns:

A tuple with attributes:

• parameter_names: list of strings, the names of the parameters

  for the current equilibrium problem

• default_tolerances: numpy.array, the default tolerances for each parameter

• n_parameters: int, the number of parameters for the current equilibrium problem

• phase_amount_indices: list of int, the indices of the parameters that

  correspond to phase amounts

• bulk_composition_vector: numpy.array, the bulk composition vector

• free_compositional_vectors: 2D numpy.array, the free compositional vectors

• reduced_composition_vector: numpy.array, the bulk composition vector

  reduced to the independent elements

• reduced_free_composition_vectors: 2D numpy.array, the free compositional vectors

  reduced to the independent elements

• constraint_vector: numpy.array, vector for the linear constraints bounding the

  valid parameter space (b in Ax + b < eps)

• constraint_matrix: 2D numpy.array, matrix for the linear constraints bounding the

  valid parameter space (A in Ax + b < eps)

Return type:

namedtuple

burnman.tools.equilibration.process_eq_constraints(equality_constraints, assemblage, prm)

A function that processes the equality constraints into a form that can be processed by the F and jacobian functions.

This function has two main tasks: it turns phase_fraction and phase_composition constraints into standard linear constraints in the solution parameters. It also turns vector-valued constraints into a list of scalar-valued constraints.

Parameters:

• equality_constraints (list) – A list of equality constraints. For valid types of constraints, please see the documentation for burnman.equilibrate().

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

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

Returns:

Equality constraints in a form that can be processed by the F and jacobian functions.

Return type:

list of lists

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

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

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

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

• P: float or numpy.array of

  pressures [Pa]

• T: float or numpy.array of

  temperatures [K]

• S: float or numpy.array of

  entropies [J/K]

• V: float or numpy.array of

  volumes [m:math:^3]

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

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

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

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

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

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

Parameters:

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

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

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

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

• tol (float or numpy.array) – The tolerance for the nonlinear solver. This can be a single float, which is then applied to all parameters, or a numpy array with the same length as the number of parameters in the solve (2 + number of phases + number of free compositional vectors). The default is None, which sets the tolerances to 1 Pa for pressure, 1e-3 K for temperature, 1e-6 for phase amounts and 1e-6 for free compositional vectors. One of the termination criteria for the nonlinear solver is that the absolute change in each parameter is less than the corresponding tolerance.

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

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

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

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

Returns:

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

Return type:

tuple

Optimal thermobarometry

burnman.tools.thermobarometry.get_reaction_matrix(assemblage, small_fraction_tol=0.0)

Set up a matrix of independent reactions for an assemblage, possibly reducing the number of endmembers by excluding components with small molar fractions.

Parameters:

• assemblage (Assemblage) – The mineral assemblage for which to build the reaction matrix.

• small_fraction_tol (float, optional, default 0.0) – If > 0.0, reduces the number of endmembers in solution phases by transforming to a smaller set of independent endmembers using a greedy algorithm and excluding those with molar fractions smaller than this value during the inversion.

Returns:

The reaction matrix for the assemblage, returned in the original endmember basis.

Return type:

2D np.array

burnman.tools.thermobarometry.assemblage_set_state_from_params(assemblage, params)

Set the state of the assemblage (P, T, compositions) from the given list of parameters.

Parameters:

• assemblage (Assemblage) – The mineral assemblage for which to set the state. If there are free compositional vectors for any solution phases (e.g., to account for an unknown state of order at fixed composition), the corresponding phases must have the attributes baseline_composition and free_compositional_vectors set, and the length of params must match the number of free compositional vectors plus two (for P and T).

• params (list or np.array) – List of parameters, where the first two are pressure (Pa) and temperature (K), and any additional parameters correspond to compositional degrees of freedom. Each compositional degree of freedom corresponds to a free compositional vector in one of the solution phases in the assemblage.

Returns:

None

Return type:

None

burnman.tools.thermobarometry.assemblage_affinity_covariance_matrix(assemblage, reaction_matrix, dataset_covariances=None, include_state_uncertainties=False)

Compute the covariance matrix of the affinities of the independent reactions for the given assemblage. This can include contributions from compositional uncertainties, dataset uncertainties, and from uncertainties in pressure and temperature.

Parameters:

• assemblage (Assemblage) – The mineral assemblage for which to compute the covariance matrix. Should have its state (P, T, compositions) set prior to calling this function.

• reaction_matrix (2D np.array) – The reaction matrix for the assemblage.

• dataset_covariances (dict, with keys 'endmember_names' and 'covariance_matrix'. Default is None, in which case only compositional uncertainties are considered.) – The covariance data from the thermodynamic dataset.

• include_state_uncertainties (bool) – If True, includes the contribution from uncertainties in pressure and temperature to the covariance matrix. If True, the assemblage must have the attribute state_covariances.

Returns:

Covariance matrix of the affinities of the independent reactions.

Return type:

2D np.array

burnman.tools.thermobarometry.assemblage_affinity_misfit(assemblage, reaction_matrix, dataset_covariances=None, include_state_uncertainties=False)

Compute the objective misfit function (chi-squared) for given P and T.

Parameters:

• assemblage (Assemblage) – The mineral assemblage for which to compute the misfit. Should have its state (P, T, compositions) set prior to calling this function.

• reaction_matrix (2D np.array) – The reaction matrix for the assemblage.

• dataset_covariances (dict, with keys 'endmember_names' and 'covariance_matrix'. Default is None, in which case only compositional uncertainties are considered.) – The covariance data from the thermodynamic dataset.

• include_state_uncertainties (bool) – If True, includes the contribution from uncertainties in pressure and temperature to the covariance matrix. If True, the assemblage must have the attribute state_covariances.

Returns:

Chi-squared misfit value.

Return type:

float

burnman.tools.thermobarometry.assemblage_state_misfit(assemblage)

Compute the objective misfit function (chi-squared) for given P and T based on prior expectations for P and T and their covariance. The assemblage must have attributes state_priors and state_inverse_covariances.

Parameters:

assemblage (Assemblage) – The mineral assemblage for which to compute the misfit.

Returns:

Chi-squared misfit value.

Return type:

float

burnman.tools.thermobarometry.estimate_conditions(assemblage, dataset_covariances=None, include_state_misfit=False, guessed_conditions=[1000000000.0, 1000.0], pressure_bounds=[100000.0, 400000000000.0], temperature_bounds=[300.0, 4000.0], P_scaling=1000000.0, small_fraction_tol=0.0, max_it=100)

Perform a least-squares inversion to find the optimal pressure and temperature for a given mineral assemblage. Algorithm modified from [PH94].

Parameters:

• assemblage (Assemblage) – The mineral assemblage for which to perform the inversion. Each solution phase in the assemblage must have its composition set along with its compositional covariance matrix (called compositional_covariances). If there are compositional degrees of freedom, they can be added by setting the attribute free_compositional_vectors on the relevant solution phases, where each row of the array corresponds to a free compositional vector, and the columns correspond to the amounts of endmembers of that phase in each vector.

• dataset_covariances (dict, with keys 'endmember_names' and 'covariance_matrix'. Default is None, in which case only compositional uncertainties are considered.) – The covariance data describing the uncertainties in endmember energies, taken from the thermodynamic dataset. For example, the covariance for the Holland et al. (2018) dataset can be obtained using burnman.minerals.HGP_2018_ds633.cov().

• include_state_misfit (bool) – If True, includes the misfit from prior expectations on P and T. The assemblage must also have attributes state_priors and state_inverse_covariances.

• guessed_conditions (np.array of shape (2,), optional, default None) – Initial guess for pressure (Pa) and temperature (K). If not provided, the initial guess will be taken from the current state of the assemblage.

• pressure_bounds (list of two floats) – Bounds for pressure (Pa) during optimization.

• temperature_bounds (list of two floats) – Bounds for temperature (K) during optimization.

• P_scaling (float) – Scaling factor for pressure to improve numerical stability.

• small_fraction_tol (float, optional, default 0.0) – If > 0.0, reduces the number of endmembers in solution phases by transforming to a smaller set of independent endmembers using a greedy algorithm and excluding those with molar fractions smaller than this value during the inversion.

• max_it (int, optional, default 100) – Maximum number of iterations for the optimization algorithm.

Returns:

Result object from the optimization containing the optimal conditions (x, which includes P, T, and any free compositional vectors) and other properties of the solution. These include the covariance matrix of the estimated parameters (xcov), the standard deviations of the estimated parameters (var), the correlation matrix (xcorr), the affinities of the independent reactions at the optimal conditions (affinities), the covariance matrix of the affinities (acov), the affinities weighted by the inverse square root of their covariance matrix (weighted_affinities), the partial derivatives of the affinities with respect to P and T (dadx), the number of independent reactions (n_reactions), the number of fitted parameters (n_params), the degrees of freedom (degrees_of_freedom), the chi-squared value (fun), the reduced chi-squared value (reduced_chisqr, given by the reduced chi-squared divided by the number of degrees of freedom), and the fit quality (fit, given by the square root of the reduced chi-squared value).

Return type:

OptimizeResult

Equations of state

burnman.tools.eos.check_eos_consistency(m, P=1000000000.0, T=300.0, tol=0.0001, verbose=False, including_shear_properties=True, equilibration_function=None)

Checks that numerical derivatives of the Gibbs energy of a mineral or composite under given conditions are equal to those provided analytically by the equation of state.

Parameters:

• m (burnman.Mineral or burnman.Composite) – The material for which the equation of state is to be checked for consistency.

• P (float) – The pressure at which to check consistency.

• T (float) – The temperature at which to check consistency.

• tol (float) – The fractional tolerance for each of the checks.

• verbose (bool) – Decide whether to print information about each check.

• including_shear_properties (bool) – Decide whether to check shear information, which is pointless for liquids and equations of state without shear modulus parameterizations.

• equilibration_function (function) – Function to internally equilibrate object. Called after every set_state. Takes the mineral object as its only argument.

Returns:

Boolean stating whether all checks have passed.

Return type:

bool

burnman.tools.eos.check_anisotropic_eos_consistency(m, P=1000000000.0, T=300.0, tol=0.0001, verbose=False, equilibration_function=None)

Checks that numerical derivatives of the Gibbs energy of an anisotropic mineral under given conditions are equal to those provided analytically by the equation of state.

Parameters:

• m (burnman.AnisotropicMineral) – The anisotropic mineral for which the equation of state is to be checked for consistency.

• P (float) – The pressure at which to check consistency.

• T (float) – The temperature at which to check consistency.

• tol (float) – The fractional tolerance for each of the checks.

• verbose (bool) – Decide whether to print information about each check.

• equilibration_function (function) – Function to internally equilibrate object. Called after every set_state. Takes the mineral object as its only argument.

Returns:

Boolean stating whether all checks have passed.

Return type:

bool

Solution

burnman.tools.solution.transform_solution_to_new_basis(solution, new_basis, n_mbrs=None, solution_name=None, endmember_names=None, molar_fractions=None)

Transforms a solution model from one endmember basis to another. Returns a new Solution object.

Parameters:

• solution (burnman.Solution object) – The original solution object.

• new_basis (2D numpy array) – The new endmember basis, given as amounts of the old endmembers.

• n_mbrs (float, optional) – The number of endmembers in the new solution (defaults to the length of new_basis).

• solution_name (str, optional) – A name corresponding to the new solution.

• endmember_names (list of str, optional) – A list corresponding to the names of the new endmembers.

• molar_fractions (numpy.array, optional) – Fractions of the new endmembers in the new solution.

Returns:

The transformed solution.

Return type:

burnman.Solution object

Chemistry

burnman.tools.chemistry.fugacity(standard_material, assemblage)

Calculates the fugacity of a standard material in another assemblage.

Note

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

Parameters:

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

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

Returns:

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

Return type:

float

burnman.tools.chemistry.relative_fugacity(component_formula, assemblage, reference_assemblage)

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

Note

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

Parameters:

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

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

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

Returns:

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

Return type:

float

burnman.tools.chemistry.equilibrium_pressure(minerals, stoichiometry, temperature, pressure_initial_guess=100000.0)

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

Parameters:

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

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

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

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

Returns:

The equilibrium pressure of the reaction [Pa].

Return type:

float

burnman.tools.chemistry.equilibrium_temperature(minerals, stoichiometry, pressure, temperature_initial_guess=1000.0)

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

Parameters:

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

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

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

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

Returns:

The equilibrium temperature of the reaction [K].

Return type:

float

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

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

Parameters:

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

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

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

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

Returns:

The equilibrium temperature of the reaction [K].

Return type:

float

burnman.tools.chemistry.hugoniot(mineral, P_ref, T_ref, pressures, reference_mineral=None)

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

Parameters:

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

• P_ref (float) – Reference pressure [Pa]

• T_ref (float) – Reference temperature [K]

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

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

Returns:

The Hugoniot temperatures and volumes at the given pressures.

Return type:

tuple of numpy.arrays

burnman.tools.chemistry.reactions_from_stoichiometric_matrix(stoichiometric_matrix)

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

Parameters:

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

Returns:

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

Return type:

2D numpy array

burnman.tools.chemistry.reactions_from_formulae(formulae, compound_names, return_strings=True)

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

Parameters:

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

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

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

Returns:

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

Return type:

2D numpy array or list of strings

Partitioning

burnman.tools.partitioning.calculate_nakajima_fp_pv_partition_coefficient(pressure, temperature, bulk_composition_mol, initial_distribution_coefficient)

Calculate the partitioning of iron between periclase and bridgmanite as given by Nakajima et al., 2012.

Parameters:

• pressure (float) – Equilibrium pressure [Pa]

• temperature (float) – Equilibrium temperature [K]

• bulk_composition_mol (dict) – Bulk composition [mol]. Only Mg, Fe, and Si are assumed to explicitly affect the partitioning with Al playing an implicit role.

• initial_distribution_coefficient (float) – The distribution coefficient (Kd_0) at 25 GPa and 0 K.

Returns:

The proportion of Fe in ferropericlase and perovskite, respectively

Return type:

tuple

Plotting

burnman.tools.plot.pretty_plot()

Makes pretty plots. Overwrites the matplotlib default settings to allow for better fonts. Slows down plotting

burnman.tools.plot.plot_projected_elastic_properties(mineral, plot_types, axes, n_zenith=31, n_azimuth=91, n_divs=100, n_rots=181, plot_vectors=True)

Parameters:

• mineral (burnman.Mineral) – Mineral object on which calculations should be done

• plot_types (list of str) – Plot types must be one of the following * ‘vp’ - V$_{P}$ (km/s) * ‘vs1’ - ‘V$_{S1}$ (km/s) * ‘vs2’ - V$_{S2}$ (km/s) * ‘vp/vs1’ - V$_{P}$/V$_{S1}$ * ‘vp/vs2’ - V$_{P}$/V$_{S2}$ * ‘s anisotropy’ - S-wave anisotropy (%), 200(vs1s - vs2s)/(vs1s + vs2s)) * ‘linear compressibility’ - Linear compressibility (GPa$^{-1}$) * ‘youngs modulus’ - Youngs Modulus (GPa) * ‘minimum poisson ratio’ - Minimum Poisson ratio * ‘maximum poisson ratio’ - Maximum Poisson ratio

• axes (matplotlib.pyplot.axes objects) – axes objects to be modified. Must be initialised with projection=’polar’.

• n_zenith (int) – Number of zeniths (plot resolution).

• n_azimuth (int) – Number of azimuths (plot resolution).

• n_divs (int) – Number of divisions for the color scale.

• n_rots (int) – Number of along-vector rotations to find minimum and maximum Poisson ratios.

• plot_vectors (bool) – Whether or not to plot vectors for $V_S$ and Poisson axial directions

Seismic

burnman.tools.seismic.attenuation_correction(v_p, v_s, v_phi, Qs, Qphi)

Applies the attenuation correction following [MBR+07], page 4. This is simplified, and there is also currently no 1D Q model implemented. The correction, however, only slightly reduces the velocities, and can be ignored for our current applications. Arguably, it might not be as relevant when comparing computations to PREM for periods of 1s as is implemented here. Called from burnman.tools.seismic.attenuation_correction()

Parameters:

• v_p (float) – P wave velocity in [m/s].

• v_s (float) – S wave velocitiy in [m/s].

• v_phi (float) – Bulk sound velocity in [m/s].

• Qs (float) – Shear quality factor [dimensionless].

• Qphi (float) – Bulk quality factor [dimensionless].

Returns:

Corrected P wave, S wave and bulk sound velocities in [m/s].

Return type:

tuple

Output for seismology

burnman.tools.output_seismo.write_tvel_file(planet_or_layer, modelname='burnmanmodel', background_model=None)

Function to write input file for obspy travel time calculations. Helper function which uses tvel_formatted_data_and_header()

Parameters:

• planet_or_layer (burnman.Planet or burnman.Layer.) – Planet or layer to write out to tvel file

• modelname (str) – Basename of file to write to (suffix .tvel added by function).

• background_model (burnman.seismic.Seismic1DModel) – 1D seismic model to fill in parts of planet (likely to be an earth model) that aren’t defined by layer (only need when using burnman.Layer)

burnman.tools.output_seismo.write_axisem_input(layers, modelname='burnmanmodel_foraxisem', axisem_ref='axisem_prem_ani_noocean.txt', plotting=False, verbose=True)

Function to write input file for AXISEM (www.axisem.info). Helper function which uses axisem_formatted_data_and_reference()

The input can be a single layer, or a list of layers taken from a planet (planet.layers). Currently this function will implement explicit discontinuities between layers in the seismic model. Currently this function is only set for Earth.

Parameters:

• layers (list of one or more burnman.Layer) – List of layers to put in AXISEM file.

• modelname (str) – Name of model, appears in name of output file.

• axisem_ref (str) – Reference file, used to copy the header and for the rest of the planet, in the case of a burnman.Layer.

• plotting (bool) – Choose whether to show plot of the old model and replaced model.

burnman.tools.output_seismo.write_mineos_input(layers, modelname='burnmanmodel_formineos', mineos_ref='mineos_prem_noocean.txt', plotting=False, verbose=True)

Function to write input file for Mineos (https://geodynamics.org/cig/software/mineos/). Helper function which uses mineos_formatted_data_and_reference()

Note: currently, this function only honors the discontinuities already in the synthetic input file, so it is best to only replace certain layers with burnman values

Parameters:

• layers (list of one or more burnman.Layer) – List of layers to put in Mineos file.

• modelname (str) – Name of model, appears in name of output file.

• mineos_ref (str) – Reference file, used to copy the header and for the rest of the planet, in the case of a burnman.Layer.

• plotting (bool) – Choose whether to show plot of the old model and replaced model.

burnman.tools.output_seismo.tvel_formatted_data_and_header(planet_or_layer, background_model=None)

Formats data and creates a header for obspy travel time calculations. Note: Because density isn’t defined for most 1D seismic models, densities are output as zeroes. The tvel format has a column for density, but this column is not used by obspy for travel time calculations.

Parameters:

• planet_or_layer (burnman.Planet or burnman.Layer.) – Planet or layer to write out to tvel file

• background_model (burnman.seismic.Seismic1DModel) – 1D seismic model to fill in parts of planet (likely to be an earth model) that aren’t defined by layer (only need when using burnman.Layer)

burnman.tools.output_seismo.axisem_formatted_data_and_reference(layers, axisem_ref='axisem_prem_ani_noocean.txt', plotting=False)

Formats data for AXISEM (www.axisem.info). The input can be a single layer, or a list of layers taken from a planet (planet.layers). Currently this function will implement explicit discontinuities between layers in the seismic model. Currently this function is only set for Earth.

Parameters:

• layers (list of one or more burnman.Layer) – List of layers to put in AXISEM file.

• axisem_ref (str) – Reference file, used to copy the header and for the rest of the planet, in the case of a burnman.Layer.

• plotting (bool) – Choose whether to show plot of the old model and replaced model.

burnman.tools.output_seismo.mineos_formatted_data_and_reference(layers, mineos_ref='mineos_prem_noocean.txt', plotting=False)

Formats data and reference data for Mineos (https://geodynamics.org/cig/software/mineos/).

Note: currently, this function only honors the discontinuities already in the synthetic input file, so it is best to only replace certain layers with burnman values

Parameters:

• layers (list of one or more burnman.Layer) – List of layers to put in Mineos file.

• mineos_ref (str) – Reference file, used to copy the header and for the rest of the planet, in the case of a burnman.Layer.

• plotting (bool) – Choose whether to show plot of the old model and replaced model.

Polytopes

burnman.tools.polytope.solution_polytope_from_charge_balance(charges, charge_total, return_fractions=False)

Creates a polytope object from a list of the charges for each species on each site and the total charge for all site-species.

Parameters:

• charges (2D list of floats) – 2D list containing the total charge for species j on site i, including the site multiplicity. So, for example, a solution with the site formula [Mg,Fe]3[Mg,Al,Si]2Si3O12 would have the following list: [[6., 6.], [4., 6., 8.]].

• charge_total (float) – The total charge for all site-species per formula unit. The example given above would have charge_total = 12.

• return_fractions (bool) – Determines whether the created polytope object returns its attributes (such as endmember occupancies) as fractions or as floats. Default is False.

Returns:

A polytope object corresponding to the parameters provided.

Return type:

burnman.polytope.MaterialPolytope object

burnman.tools.polytope.solution_polytope_from_endmember_occupancies(endmember_occupancies, return_fractions=False)

Creates a polytope object from a list of independent endmember occupancies.

Parameters:

• endmember_occupancies (2D numpy array) – 2D list containing the site-species occupancies j for endmember i. So, for example, a solution with independent endmembers [Mg]3[Al]2Si3O12, [Mg]3[Mg0.5Si0.5]2Si3O12, [Fe]3[Al]2Si3O12 might have the following array: [[1., 0., 1., 0., 0.], [1., 0., 0., 0.5, 0.5], [0., 1., 1., 0., 0.]], where the order of site-species is Mg_A, Fe_A, Al_B, Mg_B, Si_B.

• return_fractions (bool) – Determines whether the created polytope object returns its attributes (such as endmember occupancies) as fractions or as floats.

Returns:

A polytope object corresponding to the parameters provided.

Return type:

burnman.polytope.MaterialPolytope object

burnman.tools.polytope.composite_polytope_at_constrained_composition(composite, composition, return_fractions=False)

Creates a polytope object from a Composite object and a composition. This polytope describes the complete set of valid composite endmember amounts that satisfy the compositional constraints.

Parameters:

• composite (burnman.Composite object) – A composite containing one or more Solution and Mineral objects.

• composition (dict) – A dictionary containing the amounts of each element.

• return_fractions (bool) – Determines whether the created polytope object returns its attributes (such as endmember occupancies) as fractions or as floats.

Returns:

A polytope object corresponding to the parameters provided.

Return type:

burnman.polytope.MaterialPolytope object

burnman.tools.polytope.reorder_to_echelon(A)

Reorders the rows of a matrix A so that it is in echelon form. :param A: The input matrix. :type A: 2D numpy array :returns: The reordered matrix. :rtype: 2D numpy array

burnman.tools.polytope.simplify_composite_with_composition(composite, composition)

Takes a composite and a composition, and returns the simplest composite object that spans the solution space at the given composition.

For example, if the composition is given as {‘Mg’: 2., ‘Si’: 1.5, ‘O’: 5.}, and the composite is given as a mix of Mg,Fe olivine and pyroxene solutions, this function will return a composite that only contains the Mg-bearing endmembers.

If the solutions have endmember proportions that are consistent with the bulk composition, the site occupancies of the new solution models are set to the values in the original models.

Parameters:

• composite (burnman.Composite object) – The initial Composite object.

• composition (dict) – A dictionary containing the amounts of each element.

Returns:

The simplified Composite object

Return type:

burnman.Composite object

burnman.tools.polytope.greedy_independent_endmember_selection(endmember_site_occupancies, site_occupancies, small_fraction_tol=0.0, norm_tol=1e-12)

Greedy algorithm to select independent endmembers from a solution to approximate given site occupancies through a non-negative linear combination of endmember site occupancies.

This function starts with the full site occupancies and then iteratively selects endmembers to approximate those site occupancies. It loops through all possible endmembers in a number of steps. At each step the algorithm selects the endmember that can be subtracted in the largest amount from the current residual site occupancies without making any site occupancy negative. The process continues until either no endmember can be subtracted in an amount greater than fraction_tol, or the norm of the residual site occupancies is less than tol.

Parameters:

• endmember_site_occupancies (np.ndarray) – A 2D array of shape (m, n), where m is the number of endmembers and n is the number of sites. Each row corresponds to the site occupancies of an endmember.

• site_occupancies (np.ndarray) – A 1D array of length n, representing the target site occupancies to approximate.

• small_fraction_tol (float, optional, default 0.0) – Algorithm stops if no endmember can be added with a molar fraction larger than this value.

• norm_tol (float, optional, default 1e-12) – Algorithm stops if the norm of the residual site occupancies is less than this value.

Returns:

indices of selected endmembers, their fractions, and the final residual site occupancies.

Return type:

tuple(list[int], list[float], np.ndarray)