Calculator: Liquid Properties

This notebook calculates the activity of endmembers within a liquid. 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!):

User Input

1. Input liquid composition in wt% oxides, temperature, pressure, and log(fO2) below

2. Once you have input your data in the cell above, click “Run” –> “Run All Cells” at the top of the notebook. This will run the entire notebook and produce the output table at the bottom of the file.

oxide_wt_dict = {
    'SiO2':    48.4012,
    'TiO2':     0.9208,
    'Al2O3':   14.9349,
    'FeO':      9.0825,
    'MnO':      0.1587,
    'MgO':      8.1285,
    'CaO':      14.1368,
    'Na2O':      1.3124,
    'K2O':       0.1462,
    'P2O5':      0.001
    }

TK = 1500 + 273.15 # temperature value in K
Pbar = 15000 # pressure value in bar
logfO2 = None # can replace this with a value for absolute logfO2 (not relative to a buffer). Otherwise, calculates with given Fe2O3 value.

Calculations

The user does not need to modify this code– it will run automatically when “Run All Cells” is triggered.

###-----------------------------------------------------------------------------------------
### 0. Import modules
###-----------------------------------------------------------------------------------------

import numpy as np
import pandas as pd
from thermoengine import model, equilibrate, chemistry
from thermoengine.core import chem
from thermoengine import redox

###-----------------------------------------------------------------------------------------
### 1. Calculate oxide moles
###-----------------------------------------------------------------------------------------

oxide_moles = chem.format_mol_oxide_comp(oxide_wt_dict, convert_grams_to_moles=True)
oxide_moles_ser = pd.Series(oxide_moles, index=chemistry.OXIDES)
oxide_dict_for_fO2 = {'Liq': oxide_moles}

###-----------------------------------------------------------------------------------------
### 2. Adjust composition to be consistent with log fO2 input
###-----------------------------------------------------------------------------------------

if logfO2 is not None:
    ln_Fe_oxide_ratio = redox.redox_state(TK, Pbar, oxide_comp=oxide_dict_for_fO2, logfO2=logfO2,
                phase_of_interest='Liq', method='Kress91')
    Fe2O3_FeO = np.exp(ln_Fe_oxide_ratio)
    oxide_moles_ser = redox.impose_redox_ratio(Fe2O3_FeO, oxide_moles_ser)

###-----------------------------------------------------------------------------------------
### 3. Create liquid phase instance
###-----------------------------------------------------------------------------------------

liquid = model.Database(database_name="MELTS_v1_2").get_phase('Liq')
liq_end = liquid.calc_endmember_comp(oxide_moles_ser.values, method='least_squares',
                            output_residual=False, normalize=False,
                            decimals=10)

###-----------------------------------------------------------------------------------------
### 4. Calculate Gibbs Free Energies of Formation
###-----------------------------------------------------------------------------------------

end_array = np.identity(liquid.endmember_num)
G_array = np.zeros(liquid.endmember_num)
for i in range(liquid.endmember_num):
    G = liquid.gibbs_energy(TK,Pbar,mol=end_array[i,:])
    G_array[i] = G

###-----------------------------------------------------------------------------------------
### 5. Calculate Activities
###-----------------------------------------------------------------------------------------

activities = liquid.activity(TK, Pbar, mol=liq_end)

###-----------------------------------------------------------------------------------------
### 6. Combine and display data
###-----------------------------------------------------------------------------------------

data = {'Endmember': liquid.endmember_names, 'Gibbs Energy (kJ)': G_array/1000, 'Activity': activities}
endmember_df = pd.DataFrame(data)
pd.set_option('display.max_rows', None)
endmember_df