ENKI

Enthalpy potential minimization (S, P, constrained)

Closed system; crystallization of a rhyolitic liquid using rhyolite-MELTS

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

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 entropy

Note that the entropy is equivalent to $ - :raw-latex:`frac{{partial G}}{{partial T}}`$ - Run an equilibration step at fixed T,P - Calculate the entropy of the system - Define a function to set the entropy for subsequent equilibration steps

equil = equilibrate.Equilibrate(elm_sys, phs_sys)
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.8609503260090e-02 Lin (019) step:  9.4395802431652e-01 func: -1.7280794403704e+06
Quad (001) norm:  1.7020967377192e-08 Lin (026) step: -3.9941348412705e-01 func: -1.7280794403704e+06
Quad (002) norm:  2.3819371335756e-08 Lin (037) step:  9.4323348836889e-01 func: -1.7280794403704e+06
Quad (003) norm:  1.3521403531091e-09 Lin (032) step: -4.7213659552033e-01 func: -1.7280794403704e+06
Quad (004) norm:  1.9905400629586e-09 Lin (039) step: -1.2872026626417e+00 func: -1.7280794403704e+06
Quad (005) norm:  4.5527674310389e-09 Lin (039) step: -1.7924442664609e+00 func: -1.7280794403704e+06
Minimal energy termination of quadratic loop.

Add: Feldspar
Quad (000) norm:  5.8343755255746e-03 Lin (020) step:  9.7809006932692e-01 func: -1.7280795838328e+06
Quad (001) norm:  2.2098427302856e-04 Lin (012) step:  1.0471013041669e+00 func: -1.7280795887688e+06
Quad (002) norm:  9.8888798416875e-05 Lin (025) step:  1.0023917655135e+00 func: -1.7280795888147e+06
Quad (003) norm:  1.4361892154363e-07 Lin (029) step:  9.9860621867021e-01 func: -1.7280795888147e+06
Quad (004) norm:  1.6471976415762e-10 Lin (037) step: -1.0717369638414e+00 func: -1.7280795888147e+06
Quad (005) norm:  3.4122477608691e-10 Lin (034) step: -1.0565415854491e+00 func: -1.7280795888147e+06
Minimal energy termination of quadratic loop.


T =     776.85 °C, P =      175.0 MPa
Liquid          moles:   1.639690 grams: 104.628
           SiO2 form:  SiO2           X:  0.6745  wt%    SiO2   73.73
           TiO2 form:  TiO2           X:  0.0006  wt%    TiO2    0.08
          Al2O3 form:  Al2O3          X:  0.0424  wt%   Al2O3   11.82
          Fe2O3 form:  Fe2O3          X:  0.0008  wt%   Fe2O3    0.20
        MgCr2O4 form:  MgCr2O4        X:  0.0000  wt%     FeO    0.45
        Fe2SiO4 form:  Fe2SiO4        X:  0.0020  wt%     MgO    0.03
      MnSi0.5O2 form:  MnSi0.5O2      X:  0.0000  wt%     CaO    0.39
        Mg2SiO4 form:  Mg2SiO4        X:  0.0002  wt%    Na2O    3.76
      NiSi0.5O2 form:  NiSi0.5O2      X:  0.0000  wt%     K2O    4.66
      CoSi0.5O2 form:  CoSi0.5O2      X:  0.0000  wt%     H2O    4.89
         CaSiO3 form:  CaSiO3         X:  0.0044
        Na2SiO3 form:  Na2SiO3        X:  0.0387
        KAlSiO4 form:  KAlSiO4        X:  0.0631
      Ca3(PO4)2 form:  Ca3(PO4)2      X:  0.0000
            H2O form:  H2O            X:  0.1732
Feldspar        moles:   0.002127 grams:   0.567
         albite form:  NaAlSi3O8      X:  0.7392  wt%    SiO2   63.41
      anorthite form:  CaAl2Si2O8     X:  0.1886  wt%   Al2O3   22.75
       sanidine form:  KAlSi3O8       X:  0.0723  wt%     CaO    3.97
                                                  wt%    Na2O    8.60
                                                  wt%     K2O    1.28
Water           moles:   0.021370 grams:   0.385
Quartz          affn:     134.38
delta_dGdT = 0.0
dGdT = state.dGdT(t,p)
def con(t, p, state):
    return dGdT + delta_dGdT

Instantiate class instance and run calculation

equil = equilibrate.Equilibrate(elm_sys, phs_sys, lagrange_l=[('T',con)])
state = equil.execute(t, p, bulk_comp=blk_cmp, debug=0, stats=True)
state.print_state()
Add: Water
Quad (000) norm:  8.4109760121255e-02 Lin (020) step: -7.6247405768250e-01 func: -1.4419214302114e+06
Quad (001) norm:  5.5409228636560e-07 Lin (040) step:  3.8706352769574e-01 func: -1.4419214302114e+06
Quad (002) norm:  3.3961756355139e-07 Lin (033) step: -7.4383200643679e-01 func: -1.4419214302114e+06
Quad (003) norm:  5.9224454039414e-07 Lin (039) step:  5.0200500853460e-01 func: -1.4419214302114e+06
Quad (004) norm:  2.9492666607495e-07 Lin (045) step: -2.9039028052303e-02 func: -1.4419214302114e+06
Quad (005) norm:  3.0349128388339e-07 Lin (039) step: -5.0155281000757e-01 func: -1.4419214302114e+06
Minimal energy termination of quadratic loop.

Add: Feldspar
Quad (000) norm:  8.2433571075966e-01 Lin (024) step:  9.8350376266276e-01 func: -1.4419216218052e+06
Quad (001) norm:  1.8528134113517e-02 Lin (010) step:  1.0311516961913e+00 func: -1.4419216272040e+06
Quad (002) norm:  1.1294075179570e-02 Lin (026) step:  1.0033290064394e+00 func: -1.4419216272453e+06
Quad (003) norm:  2.5808106760253e-05 Lin (038) step:  9.7213570108658e-01 func: -1.4419216272453e+06
Quad (004) norm:  7.1851279402016e-07 Lin (039) step:  1.0557287173513e+00 func: -1.4419216272453e+06
Quad (005) norm:  4.0253001638508e-08 Lin (041) step:  5.5728093004448e-02 func: -1.4419216272453e+06
Minimal energy termination of quadratic loop.


T =     776.85 °C, P =      175.0 MPa
Liquid          moles:   1.639690 grams: 104.628
           SiO2 form:  SiO2           X:  0.6745  wt%    SiO2   73.73
           TiO2 form:  TiO2           X:  0.0006  wt%    TiO2    0.08
          Al2O3 form:  Al2O3          X:  0.0424  wt%   Al2O3   11.82
          Fe2O3 form:  Fe2O3          X:  0.0008  wt%   Fe2O3    0.20
        MgCr2O4 form:  MgCr2O4        X:  0.0000  wt%     FeO    0.45
        Fe2SiO4 form:  Fe2SiO4        X:  0.0020  wt%     MgO    0.03
      MnSi0.5O2 form:  MnSi0.5O2      X:  0.0000  wt%     CaO    0.39
        Mg2SiO4 form:  Mg2SiO4        X:  0.0002  wt%    Na2O    3.76
      NiSi0.5O2 form:  NiSi0.5O2      X:  0.0000  wt%     K2O    4.66
      CoSi0.5O2 form:  CoSi0.5O2      X:  0.0000  wt%     H2O    4.89
         CaSiO3 form:  CaSiO3         X:  0.0044
        Na2SiO3 form:  Na2SiO3        X:  0.0387
        KAlSiO4 form:  KAlSiO4        X:  0.0631
      Ca3(PO4)2 form:  Ca3(PO4)2      X:  0.0000
            H2O form:  H2O            X:  0.1732
Feldspar        moles:   0.002127 grams:   0.567
         albite form:  NaAlSi3O8      X:  0.7392  wt%    SiO2   63.41
      anorthite form:  CaAl2Si2O8     X:  0.1886  wt%   Al2O3   22.75
       sanidine form:  KAlSi3O8       X:  0.0723  wt%     CaO    3.97
                                                  wt%    Na2O    8.60
                                                  wt%     K2O    1.28
Water           moles:   0.021370 grams:   0.385
Quartz          affn:     134.38

Pickup runs use previously computed state

delta_dGdT = 5.0
state = equil.execute(t, p, state=state, debug=0)
state.print_state()
Minimal energy termination of quadratic loop.
Minimal energy termination of quadratic loop.

T =     762.82 °C, P =      175.0 MPa
Liquid          moles:   1.325625 grams:  84.947
           SiO2 form:  SiO2           X:  0.6770  wt%    SiO2   73.80
           TiO2 form:  TiO2           X:  0.0008  wt%    TiO2    0.09
          Al2O3 form:  Al2O3          X:  0.0399  wt%   Al2O3   11.65
          Fe2O3 form:  Fe2O3          X:  0.0010  wt%   Fe2O3    0.24
        MgCr2O4 form:  MgCr2O4        X:  0.0000  wt%     FeO    0.56
        Fe2SiO4 form:  Fe2SiO4        X:  0.0025  wt%     MgO    0.04
      MnSi0.5O2 form:  MnSi0.5O2      X:  0.0000  wt%     CaO    0.27
        Mg2SiO4 form:  Mg2SiO4        X:  0.0003  wt%    Na2O    3.64
      NiSi0.5O2 form:  NiSi0.5O2      X:  0.0000  wt%     K2O    4.90
      CoSi0.5O2 form:  CoSi0.5O2      X:  0.0000  wt%     H2O    4.82
         CaSiO3 form:  CaSiO3         X:  0.0031
        Na2SiO3 form:  Na2SiO3        X:  0.0376
        KAlSiO4 form:  KAlSiO4        X:  0.0666
      Ca3(PO4)2 form:  Ca3(PO4)2      X:  0.0000
            H2O form:  H2O            X:  0.1713
Feldspar        moles:   0.025831 grams:   6.875
         albite form:  NaAlSi3O8      X:  0.7556  wt%    SiO2   64.95
      anorthite form:  CaAl2Si2O8     X:  0.1231  wt%   Al2O3   21.51
       sanidine form:  KAlSi3O8       X:  0.1213  wt%     CaO    2.59
                                                  wt%    Na2O    8.80
                                                  wt%     K2O    2.15
Water           moles:   0.078248 grams:   1.410
Quartz          moles:   0.107219 grams:   6.442
Feldspar        moles:   0.021754 grams:   5.906
         albite form:  NaAlSi3O8      X:  0.4246  wt%    SiO2   66.03
      anorthite form:  CaAl2Si2O8     X:  0.0165  wt%   Al2O3   19.09
       sanidine form:  KAlSi3O8       X:  0.5589  wt%     CaO    0.34
                                                  wt%    Na2O    4.85
                                                  wt%     K2O    9.70