Korzhinskii potential minimization (T, P, \(\mu\)SiO2, \(\mu\)Albite, \(\mu\)Anorthite, \(\mu\)Sanidine constrained)¶
import numpy as np
import scipy.optimize as opt
import scipy.linalg as lin
import scipy as sci
import sys
from thermoengine import core, phases, model, equilibrate
#np.set_printoptions(linewidth=200, precision=1)
np.set_printoptions(linewidth=200, precision=2)
Create phases for equilibrium assemblages¶
modelDB = model.Database(liq_mod='v1.0')
Liquid = modelDB.get_phase('Liq')
Feldspar = modelDB.get_phase('Fsp')
Quartz = modelDB.get_phase('Qz')
The Berman model database provides the SWIM water model by default. Instead, override that choice by instantiating the MELTS 1.0.2 water model directly.
Water = phases.PurePhase('WaterMelts', 'H2O', calib=False)
Define elements in system and phases in system¶
elm_sys = ['H','O','Na','Mg','Al','Si','P','K','Ca','Ti','Cr','Mn','Fe','Co','Ni']
phs_sys = [Liquid, Feldspar, Water, Quartz]
Composition of the system¶
This is a high-silica rhyolite
grm_oxides = {
'SiO2': 77.5,
'TiO2': 0.08,
'Al2O3': 12.5,
'Fe2O3': 0.207,
'Cr2O3': 0.0,
'FeO': 0.473,
'MnO': 0.0,
'MgO': 0.03,
'NiO': 0.0,
'CoO': 0.0,
'CaO': 0.43,
'Na2O': 3.98,
'K2O': 4.88,
'P2O5': 0.0,
'H2O': 5.5
}
Cast this composition as moles of elements for input to the Equilibrate class
mol_oxides = core.chem.format_mol_oxide_comp(grm_oxides, convert_grams_to_moles=True)
moles_end,oxide_res = Liquid.calc_endmember_comp(
mol_oxide_comp=mol_oxides, method='intrinsic', output_residual=True)
if not Liquid.test_endmember_comp(moles_end):
print ("Calculated composition is infeasible!")
mol_elm = Liquid.covert_endmember_comp(moles_end,output='moles_elements')
blk_cmp = []
for elm in elm_sys:
index = core.chem.PERIODIC_ORDER.tolist().index(elm)
blk_cmp.append(mol_elm[index])
blk_cmp = np.array(blk_cmp)
Function to constrain the chemical potential of silica¶
def muSiO2(t, p, state):
return Quartz.gibbs_energy(t, p)
Function to constrain the chemical potential of feldspar components¶
t = 1050.0
p = 1750.0
moles_fld = np.array([0.4388, 0.0104, 0.5508])
muFld = Feldspar.gibbs_energy(t,p,mol=moles_fld,deriv={"dmol":1})
muFld[0][0],muFld[0][1],muFld[0][2]
(-4286486.8312951205, -4587946.630682982, -4325071.949624277)
def muAb(t, p, state):
return muFld[0][0]
def muAn(t, p, state):
return muFld[0][1]
def muSn(t, p, state):
return muFld[0][2]
Instantiate class instance and run calculation¶
equil = equilibrate.Equilibrate(elm_sys, phs_sys,
lagrange_l=[({'Na':1.0,'Al':1.0,'Si':3.0,'O':8.0},muAb),
({'Ca':1.0,'Al':2.0,'Si':2.0,'O':8.0},muAn),
({ 'K':1.0,'Al':1.0,'Si':3.0,'O':8.0},muSn),
({'Si':1.0, 'O':2.0},muSiO2)])
t = 1050.0
p = 1750.0
state = equil.execute(t, p, bulk_comp=blk_cmp, debug=0, stats=True)
state.print_state()
Add: Quartz
Quad (000) norm: 0.0000000000000e+00
T = 776.85 °C, P = 175.0 MPa
Liquid moles: 1.932513 grams: 124.994
SiO2 form: SiO2 X: 0.6898 wt% SiO2 74.39
TiO2 form: TiO2 X: 0.0005 wt% TiO2 0.06
Al2O3 form: Al2O3 X: 0.0400 wt% Al2O3 11.71
Fe2O3 form: Fe2O3 X: 0.0007 wt% Fe2O3 0.17
MgCr2O4 form: MgCr2O4 X: 0.0000 wt% FeO 0.38
Fe2SiO4 form: Fe2SiO4 X: 0.0017 wt% MgO 0.02
MnSi0.5O2 form: MnSi0.5O2 X: 0.0000 wt% CaO 0.17
Mg2SiO4 form: Mg2SiO4 X: 0.0002 wt% Na2O 3.70
NiSi0.5O2 form: NiSi0.5O2 X: 0.0000 wt% K2O 4.99
CoSi0.5O2 form: CoSi0.5O2 X: 0.0000 wt% H2O 4.40
CaSiO3 form: CaSiO3 X: 0.0019
Na2SiO3 form: Na2SiO3 X: 0.0386
KAlSiO4 form: KAlSiO4 X: 0.0686
Ca3(PO4)2 form: Ca3(PO4)2 X: 0.0000
H2O form: H2O X: 0.1580
Feldspar affn: 0.12
albite form: NaAlSi3O8 X: 0.4429
anorthite form: CaAl2Si2O8 X: 0.0107
sanidine form: KAlSi3O8 X: 0.5464
Water affn: 756.37
Quartz moles: 0.000010 grams: 0.001
Pickup runs use previously computed state