Phase Example: Feldspar (Solid Solution)

Demonstrates phase properties that can be called by the ThermoEngine package.

This code was originally written by Mark Ghiorso and modified by Suzanne Birner in 2025.

Open this code in an executable MyBinder instance (MyBinder links may be slow to load– please be patient!):

Initialization

Import necessary packages.

import thermoengine
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

Import a feldspar phase object.

phase_calculator = thermoengine.model.get_phase_calculator('MELTS_Feldspar_EG90')
Feldspar = thermoengine.phases.SolutionPhase('Fsp', phase_calculator)

Phase Information

Return basic properties of the feldspar phase:

print (f"Phase Name: {Feldspar.props['phase_name']}")
print (f"Abbreviaton: {Feldspar.props['abbrev']}")
print (f"Endmember Names: {Feldspar.props['endmember_name']}")
print (f"Endmember Formulas: {Feldspar.props['formula']}")
print (f"Molecular Weights: {Feldspar.props['molwt']}")

Phase Name: Feldspar
Abbreviaton: Fsp
Endmember Names: ['Albite' 'Potassium_Feldspar' 'Anorthite']
Endmember Formulas: ['NaAlSi3O8' 'KAlSi3O8' 'Al2CaSi2O8']
Molecular Weights: [262.22301 278.33524 278.20928]

Composition

Specify the composition of this instantiation of the feldspar phase: (Sandine composition in wt% oxides, taken from Deer, Howie and Zussman, Table 30, entry 9: Sanidine from a rhyolite in Mitchell Mesa, Texas)

mol_oxides = thermoengine.core.chem.format_mol_oxide_comp({'SiO2':67.27,
                                                           'Al2O3':18.35,
                                                           'FeO':0.92,
                                                           'CaO':0.15,
                                                           'Na2O':6.45,
                                                           'K2O':7.05,
                                                           'H2O':0.16},
                                                           convert_grams_to_moles=True)

Calculations

Convert analytical composition to moles of endmember components:

moles_end,oxide_res = Feldspar.calc_endmember_comp(mol_oxide_comp=mol_oxides, method='intrinsic', output_residual=True)
for i in range(0,Feldspar.props['endmember_num']):
    print (f"Mole number of {Feldspar.props['endmember_name'][i]} = {moles_end[i]:.5f}")
if not Feldspar.test_endmember_comp(moles_end):
    print ("Calculated composition is infeasible!")

print (f"Total moles of endmembers: {Feldspar.convert_endmember_comp(moles_end,output='total_moles'):.5f}")
mol_elm = Feldspar.convert_endmember_comp(moles_end,output='moles_elements')
print (f"Mole fractions of endmembers: {Feldspar.convert_endmember_comp(moles_end,output='mole_fraction')}")
print ('Total grams of phase: ', Feldspar.convert_elements(mol_elm, output='total_grams'))
print("\n")

Mole number of Albite = 0.20814
Mole number of Potassium_Feldspar = 0.14969
Mole number of Anorthite = 0.00267
Total moles of endmembers: 0.36050
Mole fractions of endmembers: [0.57735337 0.41522697 0.00741966]
Total grams of phase:  96.98566962270826

Compute chemical potentials of endmember components:

t = 1000 # in K
p = 1000 # in bar

names = Feldspar.props['endmember_name']
forms = Feldspar.props['formula']

chem_pots = Feldspar.chem_potential(t, p, mol=moles_end)

for i in range(0,Feldspar.props['endmember_num']):
    print (f"Chemical potential of {names[i]} ({forms[i]}) = {chem_pots[i]:.2f} J")

Chemical potential of Albite (NaAlSi3O8) = -4265563.32 J
Chemical potential of Potassium_Feldspar (KAlSi3O8) = -4305773.96 J
Chemical potential of Anorthite (Al2CaSi2O8) = -4566600.22 J

Compute Gibbs Free Energy and compositional derivatives:

print ('Gibbs free energy = {0:12.2f} J'.format(Feldspar.gibbs_energy(t, p, mol=moles_end)))
print ("")
dgdm = Feldspar.gibbs_energy(t, p, mol=moles_end, deriv={'dmol':1})[0]
d2gdm2 = Feldspar.gibbs_energy(t, p, mol=moles_end, deriv={'dmol':2})[0]
for i in range (0, Feldspar.props['endmember_num']):
    print ('dg/dm[{0:>2d}] = {1:13.6e}     d2gdm2[{0:>2d}][*]:  '.format(i, dgdm[i], i), end=' ')
    for j in range (0, Feldspar.props['endmember_num']):
        print ('{0:13.6e}'.format(d2gdm2[i, j]), end=' ')
    print()
print ("")
d3gdm3 = Feldspar.gibbs_energy(t, p, mol=moles_end, deriv={'dmol':3})[0]
for i in range (0, Feldspar.props['endmember_num']):
    print('d3gdm3[{0:>2d}][*][*]: '.format(i), end=' ')
    for j in range (0, Feldspar.props['endmember_num']):
        for k in range (0, Feldspar.props['endmember_num']):
            print ('{0:13.6e}'.format(d3gdm3[i, j, k]), end=' ')
        print ('  ', end=' ')
    print ()

Gibbs free energy =  -1544554.89 J

dg/dm[ 0] = -4.265563e+06     d2gdm2[ 0][*]:    2.751184e+03 -2.538435e+03 -7.202198e+04
dg/dm[ 1] = -4.305774e+06     d2gdm2[ 1][*]:   -2.538435e+03  2.698688e+03  4.649893e+04
dg/dm[ 2] = -4.566600e+06     d2gdm2[ 2][*]:   -7.202198e+04  4.649893e+04  3.002095e+06

d3gdm3[ 0][*][*]:  -3.780052e+04  2.817484e+04  3.360941e+05     2.817484e+04 -2.226127e+04  2.438818e+03     3.360941e+05  2.438818e+03  6.370287e+05
d3gdm3[ 1][*][*]:   2.817484e+04 -2.226127e+04  2.438818e+03    -2.226127e+04  1.847594e+04 -3.106709e+05     2.438818e+03 -3.106709e+05 -1.878906e+05
d3gdm3[ 2][*][*]:   3.360941e+05  2.438818e+03  6.370287e+05     2.438818e+03 -3.106709e+05 -1.878906e+05     6.370287e+05 -1.878906e+05 -1.161426e+09

Plot Gibbs Free Energy versus temperature:

T_array = np.linspace(250.0, 1200.0, 100, endpoint=True)
G_array = Feldspar.gibbs_energy(T_array, 1000.0, mol=moles_end)
plt.plot(T_array, G_array)
plt.ylabel('GFE (J)')
plt.xlabel('T (K)')
plt.show()

Compute Entropy, Enthalpy, and molar derivatives:

print ('Entropy = {0:12.2f} J/K'.format(Feldspar.entropy(t, p, mol=moles_end)))
print ("")
dsdm = Feldspar.entropy(t, p, mol=moles_end, deriv={'dmol':1})[0]
d2sdm2 = Feldspar.entropy(t, p, mol=moles_end, deriv={'dmol':2})[0]
for i in range (0, Feldspar.props['endmember_num']):
    print ('ds/dm[{0:>2d}] = {1:13.6e}     d2sdm2[{0:>2d}][*]:  '.format(i, dsdm[i], i), end=' ')
    for j in range (0, Feldspar.props['endmember_num']):
        print ('{0:13.6e}'.format(d2sdm2[i, j]), end=' ')
    print()

Entropy =       199.58 J/K

ds/dm[ 0] =  5.506183e+02     d2sdm2[ 0][*]:   -2.691158e+01  3.697437e+01  2.489871e+01
ds/dm[ 1] =  5.573605e+02     d2sdm2[ 1][*]:    3.697437e+01 -5.177324e+01  2.026651e+01
ds/dm[ 2] =  5.768298e+02     d2sdm2[ 2][*]:    2.489871e+01  2.026651e+01 -3.071644e+03

Compute Heat Capacity and its derivatives:

print ('{0:<10s}{1:13.6f} {2:<15s}'.format('Cp', Feldspar.heat_capacity(t, p, mol=moles_end), 'J/K'))
print ("")
dcpdm = Feldspar.heat_capacity(t, p, mol=moles_end, deriv={'dmol':1})[0]
for i in range (0, Feldspar.props['endmember_num']):
    print ('dCp/dm[{0:>2d}]    = {1:13.6e} J/K-m'.format(i, dcpdm[i]))

Cp           117.868897 J/K

dCp/dm[ 0]    =  3.313120e+02 J/K-m
dCp/dm[ 1]    =  3.210186e+02 J/K-m
dCp/dm[ 2]    =  3.208837e+02 J/K-m

Compute Volume:

print ('{0:<10s}{1:13.6f} {2:<15s}'.format('Volume', Feldspar.volume(t, p, mol=moles_end), 'J/bar'))

Volume         3.843450 J/bar