Note
Go to the end to download the full example code.
Calculator: Solid Phase Properties#
This notebook calculates the activity of endmembers given the composition of a solid phase. It replaces the functionality previously implemented by the MELTS Supplemental Calculator (https://melts.ofm-research.org/CalcForms/)
This calculator is a work in progress – please let us know if you find any errors!
Open this code in an executable MyBinder instance (MyBinder links may be slow to load– please be patient!):
Information#
The following solution phases are implemented in this calculator:
‘Ol’ (olivine)
‘Fsp’ (feldspar)
‘Opx’ (orthopyroxene)
‘Cpx’ (clinopyroxene)
‘SplS’ (spinel solid solution)
‘Gt’ (garnet)
This notebook is experimental and has not yet been extensively tested! Proceed with caution!
User Input#
Define phase composition in wt% oxides, temperature, and pressure below. Components not included in the phase model will be ignored.
Once you have input your data, click “Run” –> “Run All Cells” at the top of the notebook. This will run the entire notebook and produce an output table at the bottom of the file.
phase = 'Opx'
# oxide_wt_dict = {
# 'SiO2': 0.000,
# 'TiO2': 0.118399,
# 'Al2O3': 31.4741,
# 'Fe2O3': 3.623535369,
# 'Cr2O3': 33.6659,
# 'FeO': 15.77978837,
# 'MnO': 0.248043,
# 'MgO': 13.0798,
# 'NiO': 0.140058,
# 'CoO': 0.000,
# 'CaO': 0.019049,
# 'Na2O': 0.0039,
# 'K2O': 0.000,
# 'P2O5': 0.000,
# }
oxide_wt_dict = {
'SiO2': 55.564,
'TiO2': 0.045,
'Al2O3': 3.436,
'Fe2O3': 0.000,
'Cr2O3': 1.011,
'FeO': 5.453,
'MnO': 0.116,
'MgO': 32.341,
'NiO': 0.105,
'CoO': 0.000,
'CaO': 1.569,
'Na2O': 0.077,
'K2O': 0.000,
'P2O5': 0.000,
}
T_K = 794.12 + 273.15
P_bar = 6000
# Advanced options
database='MELTS_v1_2'
Calculations#
The user does not need to modify this code– it will run automatically when “Run All Cells” is triggered.
Initialization#
###-----------------------------------------------------------------------------------------
### 0. Import modules
###-----------------------------------------------------------------------------------------
import numpy as np
import pandas as pd
from thermoengine import model, rockychem
###-----------------------------------------------------------------------------------------
### 1. Create Phase Instance and Calculate Endmember Fractions
###-----------------------------------------------------------------------------------------
phs_obj = model.Database(database_name=database).get_phase(phase)
phs_end_frac = rockychem.convert_phase_comp(phase, database=database, wt_oxides=oxide_wt_dict, to='mol_endmems')
if phase == 'SplS':
# Custom method for calculating endmembers; more accurate for magnetite activity than default endmember calculation
mol_oxides = rockychem.convert_bulk_comp(wt_oxides=oxide_wt_dict, to='mol_oxides')
mol_cat = pd.Series([mol_oxides['TiO2'],
2*mol_oxides['Al2O3'],
2*mol_oxides['Cr2O3'],
2*mol_oxides['Fe2O3'],
mol_oxides['FeO'],
mol_oxides['MgO']],
['Ti', 'Al', 'Cr', 'Fe3', 'Fe2', 'Mg'])
cat_O = pd.Series([2, 3/2, 3/2, 3/2, 1, 1],
['Ti', 'Al', 'Cr', 'Fe3', 'Fe2', 'Mg'])
mol_O = sum(mol_cat*cat_O)
norm_cat = 4*mol_cat/mol_O
magnesioaluminate = norm_cat['Mg']
chromite = norm_cat['Cr']/2
ulvospinel = norm_cat['Ti']
magnetite = norm_cat['Fe3']/2
hercynite = norm_cat['Fe2'] - (chromite + 2*ulvospinel + magnetite)
phs_end_frac = np.array([hercynite, magnesioaluminate, chromite, ulvospinel, magnetite])
phs_end_frac = phs_end_frac/np.sum(phs_end_frac)
Phase Properties#
#mass_scaling = 100 / phs_endmems.dot(phs_endmems.endmember_molwts)
info = pd.Series({
'Molar Mass (g/mol)': phs_end_frac.dot(phs_obj.endmember_molwts),
'Molar Volume (cm^3/mol)': phs_obj.volume(T_K, P_bar, mol=phs_end_frac) * 10,
'Density (g/cm^3)': phs_end_frac.dot(phs_obj.endmember_molwts)/(phs_obj.volume(T_K, P_bar, mol=phs_end_frac) * 10),
'Molar Heat Capacity (J/K/mol)':(phs_obj.heat_capacity(T_K, P_bar, mol=phs_end_frac)),
'Molar Gibbs Energy (kJ/mol)': (phs_obj.gibbs_energy(T_K, P_bar, mol=phs_end_frac) / 1000)
})
pd.DataFrame(info, columns=[f'{phs_obj.phase_name} Properties'])
Endmember Properties#
###-----------------------------------------------------------------------------------------
### 2. Calculate Gibbs Free Energies of Formation (Composition-Independent)
###-----------------------------------------------------------------------------------------
end_array = np.identity(len(phs_end_frac))
G_array = np.zeros_like(phs_end_frac)
for i in range(len(phs_end_frac)):
G = phs_obj.gibbs_energy(T_K,P_bar,mol=end_array[i,:])
G_array[i] = G
###-----------------------------------------------------------------------------------------
### 3. Calculate Activities (Composition-Dependent)
###-----------------------------------------------------------------------------------------
activities = phs_obj.activity(T_K, P_bar, mol=phs_end_frac)
###-----------------------------------------------------------------------------------------
### 4. Combine and Display Data
###-----------------------------------------------------------------------------------------
data = {'Endmember': phs_obj.endmember_names,
'Formula': phs_obj.endmember_formulas,
'Mole fraction': pd.Series(phs_end_frac),
'Gibbs Energy (kJ)': G_array/1000,
'Activity': activities}
endmember_df = pd.DataFrame(data)
pd.set_option('display.max_columns', None)
pd.set_option('display.max_rows', None)
endmember_df
print(np.sum(G_array*phs_end_frac)/1000)
-3272.0930543122986
Total running time of the script: (0 minutes 0.125 seconds)