Korzhinskii potential minimization (T, P, \(\mu\)SiO2 constrained)¶
import numpy as np
import scipy.optimize as opt
import scipy.linalg as lin
import sys
from thermoengine import core, phases, model, equilibrate
np.set_printoptions(linewidth=200, precision=1)
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¶
Note that quartz saturation will be imposed, so adding quartz as a system phase is redundant and merely moditors the stauration condition.
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 SiO2¶
def muSiO2(t, p, state):
return Quartz.gibbs_energy(t, p)
Instantiate class instance and run calculation¶
equil = equilibrate.Equilibrate(elm_sys, phs_sys, lagrange_l=[({'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: Water
Quad (000) norm: 2.0919339116256e-02 Lin (022) step: 1.1404184587577e+00 func: -6.3463025951454e+05
Quad (001) norm: 3.3488088202158e-03 Lin (034) step: 1.1687446083504e+00 func: -6.3463005405689e+05
Quad (002) norm: 5.5404363799077e-04 Lin (024) step: 1.1640376709056e+00 func: -6.3463004850062e+05
Quad (003) norm: 9.1181655016418e-05 Lin (014) step: 1.1642086859434e+00 func: -6.3463004832265e+05
Quad (004) norm: 1.4962742159588e-05 Lin (033) step: 1.1664159924080e+00 func: -6.3463004833999e+05
Quad (005) norm: 2.4903355411261e-06 Lin (035) step: 1.1721712654619e+00 func: -6.3463004834073e+05
Minimal energy termination of quadratic loop.
Add: Quartz
Quad (000) norm: 4.2865017842861e-07 Lin (037) step: 1.0227907354696e+00 func: -6.3463989322211e+05
Quad (001) norm: 9.6509648395378e-09 Lin (038) step: -3.8716905600240e-01 func: -6.3463989322079e+05
Quad (002) norm: 1.3387516320064e-08 Lin (035) step: 1.6393202493172e+00 func: -6.3463989322079e+05
Quad (003) norm: 8.5590033366465e-09 Lin (040) step: -1.0509103027852e+00 func: -6.3463989322079e+05
Quad (004) norm: 1.7554178136231e-08 Lin (038) step: 8.9452579470777e-01 func: -6.3463989322079e+05
Quad (005) norm: 1.8506933462712e-09 Lin (050) step: 9.0631518900942e-03 func: -6.3463989322079e+05
Minimal energy termination of quadratic loop.
T = 776.85 °C, P = 175.0 MPa
Liquid moles: 1.711813 grams: 108.664
SiO2 form: SiO2 X: 0.6808 wt% SiO2 74.34
TiO2 form: TiO2 X: 0.0006 wt% TiO2 0.07
Al2O3 form: Al2O3 X: 0.0414 wt% Al2O3 11.50
Fe2O3 form: Fe2O3 X: 0.0008 wt% Fe2O3 0.19
MgCr2O4 form: MgCr2O4 X: 0.0000 wt% FeO 0.44
Fe2SiO4 form: Fe2SiO4 X: 0.0019 wt% MgO 0.03
MnSi0.5O2 form: MnSi0.5O2 X: 0.0000 wt% CaO 0.40
Mg2SiO4 form: Mg2SiO4 X: 0.0002 wt% Na2O 3.66
NiSi0.5O2 form: NiSi0.5O2 X: 0.0000 wt% K2O 4.49
CoSi0.5O2 form: CoSi0.5O2 X: 0.0000 wt% H2O 4.88
CaSiO3 form: CaSiO3 X: 0.0045
Na2SiO3 form: Na2SiO3 X: 0.0375
KAlSiO4 form: KAlSiO4 X: 0.0605
Ca3(PO4)2 form: Ca3(PO4)2 X: 0.0000
H2O form: H2O X: 0.1719
Feldspar affn: 264.01
albite form: NaAlSi3O8 X: 0.7411
anorthite form: CaAl2Si2O8 X: 0.1888
sanidine form: KAlSi3O8 X: 0.0701
Water moles: 0.011074 grams: 0.200
Quartz moles: 0.000010 grams: 0.001
Pickup runs use previously computed state
state = equil.execute(t-5.0, p, state=state, debug=0, stats=True)
state.print_state()
Add: Feldspar
Quad (000) norm: 4.1530329504509e-01 Lin (037) step: -5.5625677045860e-02 func: -6.5238951381008e+05
Quad (001) norm: 7.0429498604093e-01 Lin (038) step: 4.4370623407854e-03 func: -6.5447193926378e+05
Quad (002) norm: 5.6770585823860e-01 Lin (037) step: -9.2767951349095e-03 func: -6.5736313225043e+05
Quad (003) norm: 6.3586523505206e-01 Lin (040) step: -1.0150575486173e-03 func: -6.5737880325226e+05
Quad (004) norm: 5.7201759044269e-01 Lin (037) step: 3.9933561616004e-03 func: -6.5757808862571e+05
Quad (005) norm: 9.5799555960644e-01 Lin (038) step: -6.0865053398050e-03 func: -6.5894764733790e+05
Quad (006) norm: 2.1868223831377e+00 Lin (043) step: 3.5392868187752e-04 func: -6.5902061774685e+05
Quad (007) norm: 2.1471153681086e+00 Lin (048) step: 3.8726518890751e-05 func: -6.5902800631958e+05
Quad (008) norm: 2.1369697719763e+00 Lin (058) step: 1.7962899250578e-07 func: -6.5902803994307e+05
Remove: Water
Quad (000) norm: 1.8711143615073e+00 Lin (039) step: 3.9933561294381e-03 func: -6.5973127077201e+05
Quad (001) norm: 7.3595725016869e-01 Lin (037) step: -2.4315330541724e-01 func: -7.7684298032846e+05
Quad (002) norm: 1.4478457722711e+00 Lin (037) step: 1.4139300658851e-02 func: -7.8112221305756e+05
Quad (003) norm: 1.0430285089303e+00 Lin (037) step: 1.2725370572772e-02 func: -7.9113032478277e+05
Quad (004) norm: 6.5323615637559e-01 Lin (037) step: -1.0466952493557e-01 func: -8.2691962540789e+05
Quad (005) norm: 6.6461155267123e-01 Lin (040) step: -2.6200409720147e-03 func: -8.2796047317863e+05
Quad (006) norm: 6.8415457564177e-01 Lin (040) step: -2.6200409720147e-03 func: -8.2924105627869e+05
Quad (007) norm: 6.9378571255856e-01 Lin (040) step: -2.6200409720147e-03 func: -8.3051766289862e+05
Quad (008) norm: 7.0495145670417e-01 Lin (040) step: -2.6200409720147e-03 func: -8.3175832252417e+05
Quad (009) norm: 7.2036448957476e-01 Lin (040) step: -2.6200409720147e-03 func: -8.3295174507994e+05
Quad (010) norm: 7.4171934579853e-01 Lin (040) step: -2.6200409720147e-03 func: -8.3408697267458e+05
Quad (011) norm: 7.7081352312931e-01 Lin (040) step: -2.6200409647262e-03 func: -8.3515203750756e+05
Quad (012) norm: 8.0939421471903e-01 Lin (040) step: -2.6200409647262e-03 func: -8.3613426427993e+05
Quad (013) norm: 8.5880776705343e-01 Lin (040) step: -2.6200409647262e-03 func: -8.3702087385821e+05
Quad (014) norm: 9.1962373385636e-01 Lin (040) step: -2.6200409566279e-03 func: -8.3779976099576e+05
Quad (015) norm: 9.9134479158062e-01 Lin (040) step: -2.6200409566279e-03 func: -8.3846000656735e+05
Quad (016) norm: 1.0724205980903e+00 Lin (040) step: -2.6200409502138e-03 func: -8.3900026475082e+05
Quad (017) norm: 8.4147166109310e-16
Add: Water
Quad (000) norm: 5.9873014135237e-01 Lin (042) step: -4.3029465263668e-05 func: -8.3902618254186e+05
Quad (001) norm: 5.5438854091657e-01 Lin (060) step: -4.2374165324941e-06 func: -8.3902846091424e+05
Quad (002) norm: 5.7445840878890e-01 Lin (053) step: -1.3094953623177e-07 func: -8.3902853755442e+05
Remove: Water
Quad (000) norm: 8.2723861329952e-16
Add: Water
Quad (000) norm: 5.7256978130451e-01 Lin (055) step: -4.7810516978195e-05 func: -8.3905503858508e+05
Quad (001) norm: 5.7502051822156e-01 Lin (065) step: -2.4640465361673e-07 func: -8.3905517349952e+05
Remove: Water
Quad (000) norm: 8.1456766203904e-16
Add: Water
Quad (000) norm: 5.9557552355154e-01 Lin (042) step: -4.7810516682650e-05 func: -8.3908321531022e+05
Quad (001) norm: 5.7081950517278e-01 Lin (055) step: -6.2632755853170e-08 func: -8.3908325171622e+05
Remove: Water
Quad (000) norm: 8.0245280288361e-16
Add: Water
Quad (000) norm: 7.6453909807180e-16
T = 771.85 °C, P = 175.0 MPa
Liquid moles: 1.355664 grams: 81.544
SiO2 form: SiO2 X: 0.6603 wt% SiO2 75.37
TiO2 form: TiO2 X: 0.0007 wt% TiO2 0.10
Al2O3 form: Al2O3 X: 0.0186 wt% Al2O3 9.27
Fe2O3 form: Fe2O3 X: 0.0010 wt% Fe2O3 0.25
MgCr2O4 form: MgCr2O4 X: 0.0000 wt% FeO 0.58
Fe2SiO4 form: Fe2SiO4 X: 0.0024 wt% MgO 0.04
MnSi0.5O2 form: MnSi0.5O2 X: 0.0000 wt% CaO 0.00
Mg2SiO4 form: Mg2SiO4 X: 0.0003 wt% Na2O 2.00
NiSi0.5O2 form: NiSi0.5O2 X: 0.0000 wt% K2O 5.65
CoSi0.5O2 form: CoSi0.5O2 X: 0.0000 wt% H2O 6.74
CaSiO3 form: CaSiO3 X: 0.0000
Na2SiO3 form: Na2SiO3 X: 0.0194
KAlSiO4 form: KAlSiO4 X: 0.0722
Ca3(PO4)2 form: Ca3(PO4)2 X: 0.0000
H2O form: H2O X: 0.2252
Feldspar moles: 0.089326 grams: 23.638
albite form: NaAlSi3O8 X: 0.8499 wt% SiO2 66.17
anorthite form: CaAl2Si2O8 X: 0.0858 wt% Al2O3 20.92
sanidine form: KAlSi3O8 X: 0.0643 wt% CaO 1.82
wt% Na2O 9.95
wt% K2O 1.14
Water moles: 0.000010 grams: 0.000
Quartz moles: 0.000012 grams: 0.001
state = equil.execute(t-10.0, p, state=state, debug=0, stats=True)
state.print_state()
Quad (000) norm: 5.6688245635369e-01 Lin (043) step: -4.7810516655517e-05 func: -8.3831498931258e+05
Quad (001) norm: 7.5407306663823e-16
T = 766.85 °C, P = 175.0 MPa
Liquid moles: 1.342725 grams: 80.765
SiO2 form: SiO2 X: 0.6570 wt% SiO2 75.14
TiO2 form: TiO2 X: 0.0007 wt% TiO2 0.10
Al2O3 form: Al2O3 X: 0.0187 wt% Al2O3 9.35
Fe2O3 form: Fe2O3 X: 0.0010 wt% Fe2O3 0.26
MgCr2O4 form: MgCr2O4 X: 0.0000 wt% FeO 0.59
Fe2SiO4 form: Fe2SiO4 X: 0.0025 wt% MgO 0.04
MnSi0.5O2 form: MnSi0.5O2 X: 0.0000 wt% CaO 0.00
Mg2SiO4 form: Mg2SiO4 X: 0.0003 wt% Na2O 2.01
NiSi0.5O2 form: NiSi0.5O2 X: 0.0000 wt% K2O 5.71
CoSi0.5O2 form: CoSi0.5O2 X: 0.0000 wt% H2O 6.81
CaSiO3 form: CaSiO3 X: 0.0000
Na2SiO3 form: Na2SiO3 X: 0.0196
KAlSiO4 form: KAlSiO4 X: 0.0729
Ca3(PO4)2 form: Ca3(PO4)2 X: 0.0000
H2O form: H2O X: 0.2274
Feldspar moles: 0.089337 grams: 23.641
albite form: NaAlSi3O8 X: 0.8499 wt% SiO2 66.17
anorthite form: CaAl2Si2O8 X: 0.0858 wt% Al2O3 20.92
sanidine form: KAlSi3O8 X: 0.0643 wt% CaO 1.82
wt% Na2O 9.95
wt% K2O 1.14
Water moles: 0.000001 grams: 0.000
Quartz moles: 0.000012 grams: 0.001
state = equil.execute(t-15.0, p, state=state, debug=0, stats=True)
state.print_state()
Quad (000) norm: 7.5572812753852e-16
T = 761.85 °C, P = 175.0 MPa
Liquid moles: 1.329912 grams: 79.995
SiO2 form: SiO2 X: 0.6537 wt% SiO2 74.90
TiO2 form: TiO2 X: 0.0008 wt% TiO2 0.10
Al2O3 form: Al2O3 X: 0.0189 wt% Al2O3 9.44
Fe2O3 form: Fe2O3 X: 0.0010 wt% Fe2O3 0.26
MgCr2O4 form: MgCr2O4 X: 0.0000 wt% FeO 0.59
Fe2SiO4 form: Fe2SiO4 X: 0.0025 wt% MgO 0.04
MnSi0.5O2 form: MnSi0.5O2 X: 0.0000 wt% CaO 0.00
Mg2SiO4 form: Mg2SiO4 X: 0.0003 wt% Na2O 2.03
NiSi0.5O2 form: NiSi0.5O2 X: 0.0000 wt% K2O 5.76
CoSi0.5O2 form: CoSi0.5O2 X: 0.0000 wt% H2O 6.88
CaSiO3 form: CaSiO3 X: 0.0000
Na2SiO3 form: Na2SiO3 X: 0.0197
KAlSiO4 form: KAlSiO4 X: 0.0736
Ca3(PO4)2 form: Ca3(PO4)2 X: 0.0000
H2O form: H2O X: 0.2296
Feldspar moles: 0.089337 grams: 23.641
albite form: NaAlSi3O8 X: 0.8499 wt% SiO2 66.17
anorthite form: CaAl2Si2O8 X: 0.0858 wt% Al2O3 20.92
sanidine form: KAlSi3O8 X: 0.0643 wt% CaO 1.82
wt% Na2O 9.95
wt% K2O 1.14
Water moles: 0.000001 grams: 0.000
Quartz moles: 0.000012 grams: 0.001
state = equil.execute(t-20.0, p, state=state, debug=0, stats=True)
state.print_state()
Quad (000) norm: 7.5733316417540e-16
T = 756.85 °C, P = 175.0 MPa
Liquid moles: 1.317201 grams: 79.232
SiO2 form: SiO2 X: 0.6504 wt% SiO2 74.66
TiO2 form: TiO2 X: 0.0008 wt% TiO2 0.10
Al2O3 form: Al2O3 X: 0.0191 wt% Al2O3 9.53
Fe2O3 form: Fe2O3 X: 0.0010 wt% Fe2O3 0.26
MgCr2O4 form: MgCr2O4 X: 0.0000 wt% FeO 0.60
Fe2SiO4 form: Fe2SiO4 X: 0.0025 wt% MgO 0.04
MnSi0.5O2 form: MnSi0.5O2 X: 0.0000 wt% CaO 0.00
Mg2SiO4 form: Mg2SiO4 X: 0.0003 wt% Na2O 2.05
NiSi0.5O2 form: NiSi0.5O2 X: 0.0000 wt% K2O 5.82
CoSi0.5O2 form: CoSi0.5O2 X: 0.0000 wt% H2O 6.94
CaSiO3 form: CaSiO3 X: 0.0000
Na2SiO3 form: Na2SiO3 X: 0.0199
KAlSiO4 form: KAlSiO4 X: 0.0743
Ca3(PO4)2 form: Ca3(PO4)2 X: 0.0000
H2O form: H2O X: 0.2318
Feldspar moles: 0.089337 grams: 23.641
albite form: NaAlSi3O8 X: 0.8499 wt% SiO2 66.17
anorthite form: CaAl2Si2O8 X: 0.0858 wt% Al2O3 20.92
sanidine form: KAlSi3O8 X: 0.0643 wt% CaO 1.82
wt% Na2O 9.95
wt% K2O 1.14
Water moles: 0.000001 grams: 0.000
Quartz moles: 0.000012 grams: 0.001