PhaseObjC - DEW Quartz Solubility¶
Required python code to load the PhaseObjC library. The library libphaseobjc.dylib (see build instructions in README.md) must be locatable to teh system in a standard location (by default /usr/local/lib)
from ctypes import cdll
from ctypes import util
from rubicon.objc import ObjCClass, objc_method
cdll.LoadLibrary(util.find_library('phaseobjc'))
/Users/ghiorso/anaconda3/lib/python3.7/site-packages/rubicon/objc/ctypes_patch.py:24: UserWarning: rubicon.objc.ctypes_patch has only been tested with Python 3.4 through 3.6. The current version is sys.version_info(major=3, minor=7, micro=6, releaselevel='final', serial=0). Most likely things will work properly, but you may experience crashes if Python's internals have changed significantly.
.format(sys.version_info)
<CDLL '/usr/local/lib/libphaseobjc.dylib', handle 7fe0a04dcfb0 at 0x7fe0785a8650>
Create a python reference to the DEWFluid class and instantiate an instance of that class.¶
In Objective-C the code has the form:
obj = [[CpxBerman alloc] init]
and in python:
DEWFluid = ObjCClass('DEWFluid')
obj = DEWFluid.alloc().init()
Solution Protocol Functions¶
–> Retrieves the number of endmember components in the system
(NSUInteger)numberOfSolutionComponents
–> Retrieves the number of species (dependent endmembers with positive mole fractions) in the soluton
(NSUInteger)numberOfSolutionSpecies
Note that the number of components (nc) may be the same as the number of endmember species (ns) if the solution does not involve speciation (complexing) or if the solid solution is not a reciprocal solution.
nc = obj.numberOfSolutionComponents()
print ('Number of components = ', nc)
ns = obj.numberOfSolutionSpecies()
print ('Number of species = ', ns)
Number of components = 17
Number of species = 92
Return name, formula and molecular weight of each end-member component in the solution. Note teh use of the PhaseBase class (for further info see the Stoichiometric Phase notebook examples)
–> Retrieves superclass instance of PhaseBase object for component at specified index
(id)componentAtIndex:(NSUInteger)index
PhaseBase = ObjCClass('PhaseBase')
print ('Component name'.ljust(20), '\t', 'Formula'.ljust(20), '\t', 'Molecular weight (g/mol)')
for i in range(0, nc):
component = obj.componentAtIndex_(i)
print (component.phaseName.ljust(20), '\t', component.phaseFormula.ljust(20), '\t', component.mw)
Component name Formula Molecular weight (g/mol)
Water H2O 18.0152
CO2,aq CO2 44.0098
O2,aq O2 31.9988
HF,aq HF 20.006303
NaOH,aq NaOH 39.99707
Mg(OH)2,aq Mg(OH)2 58.319599999999994
HAlO2,aq HAlO2 59.98824
SiO2,aq SiO2 60.0843
H3PO4,aq H3PO4 97.99506
SO2,aq SO2 64.0588
HCl,aq HCl 36.4609
KOH,aq KOH 56.1093
Ca(OH)2,aq Ca(OH)2 74.0946
H2CrO4,aq H2CrO4 118.0094
Mn(OH)2,aq Mn(OH)2 88.9526
Fe(OH)2,aq Fe(OH)2 89.8616
Co(OH)2,aq Co(OH)2 92.9478
Allocate a “c”-type pointer to a double precision one-dimensional array, and initialize the array to hold the moles of each component in the solution
import ctypes
m = (ctypes.c_double*nc)()
ctypes.cast(m, ctypes.POINTER(ctypes.c_double))
m[ 0] = 55.55
m[ 1] = 0.0
m[ 2] = 0.0
m[ 3] = 0.0
m[ 4] = 0.0
m[ 5] = 0.0
m[ 6] = 0.0
m[ 7] = 0.04223 # 0.03 at 400 C and 2 kb, 0.3 at 600 C and 10 kb
m[ 8] = 0.0
m[ 9] = 0.0
m[10] = 0.0
m[11] = 0.0
m[12] = 0.0
m[13] = 0.0
m[14] = 0.0
m[15] = 0.0
m[16] = 0.0
for i in range (0, nc):
component = obj.componentAtIndex_(i)
if m[i] > 0.0:
print ('moles of (', component.phaseName.ljust(20), ') = ', m[i])
moles of ( Water ) = 55.55
moles of ( SiO2,aq ) = 0.04223
if (obj.testPermissibleValuesOfComponents_(m) == 1):
print ('mole input is feasible')
else:
print ('mole input is infeasible')
print ('Total moles = ', obj.totalMolesFromMolesOfComponents_(m))
mole_frac_pointer = obj.convertMolesToMoleFractions_(m)
print ('component name'.ljust(20), 'component formula'.ljust(20), 'mole fraction')
for i in range (0, nc):
if m[i] > 0.0:
print (obj.componentAtIndex_(i).phaseName.ljust(20), obj.componentAtIndex_(i).phaseFormula.ljust(20), mole_frac_pointer.valueAtIndex_(i))
moles_pointer = obj.convertMolesToElements_(m)
ne = moles_pointer.size
formula = ''
for i in range(0, ne):
value = moles_pointer.valueAtIndex_(i)
if value != 0.0:
name = PhaseBase.elementNameFromAtomicNumber_(i)
formula = formula + name + '(' + str(value) + ')'
print ('Solution formula = ', formula)
mole input is feasible
Total moles = 55.59223
component name component formula mole fraction
Water H2O 0.9992403614677806
SiO2,aq SiO2 0.0007596385322193406
Solution formula = H(111.1)O(55.63446)Si(0.04223)
Change convention for Berman properties¶
Quartz.enableGibbsFreeEnergyReferenceStateUsed()
Quartz.disableQuartzCorrectionUsed()
BermanGTrPr = -856288.0
t = 673.15
p = 1000.0
activity = obj.getActivityFromMolesOfComponents_andT_andP_(m, t, p)
potential = obj.getChemicalPotentialFromMolesOfComponents_andT_andP_(m, t, p)
print ('Gibbs free energy (J) = ', obj.getGibbsFreeEnergyFromMolesOfComponents_andT_andP_(m, t, p))
print ('component, activity, chemical potential, mu0')
mu0_quartz = qtz.getGibbsFreeEnergyFromT_andP_(t, p)
for i in range (0, nc):
if m[i] > 0.0:
component = obj.componentAtIndex_(i)
g = component.getGibbsFreeEnergyFromT_andP_(t, p)
print ("{0:>20s} {1:10.6f} {2:15.2f} {3:15.2f}".format(component.phaseName, activity.valueAtIndex_(i), potential.valueAtIndex_(i), g))
print ("{0:>20s} {1:15.2f}".format('Quartz', mu0_quartz))
import numpy as np
muSpecies = obj.chemicalPotentialsOfSpeciesFromMolesOfComponents_andT_andP_(m, t, p)
print ('species, activity, mu, mu0')
for i in range (0, ns):
if muSpecies.valueAtIndex_(i) != 0.0:
pure = obj.componentAtIndex_(i)
g = pure.getGibbsFreeEnergyFromT_andP_(t, p)
print("{0:>20s} {1:10.6e} {2:15.2f} {3:15.2f}".format(obj.nameOfSolutionSpeciesAtIndex_(i), np.exp((muSpecies.valueAtIndex_(i)-g)/8.3143/t), muSpecies.valueAtIndex_(i), g))
print ('deltaG quartz -> SiO2 (aq)', muSpecies.valueAtIndex_(7)-mu0_quartz)
print ("at log (molality) = {0:10.6e}".format(np.log10(m[7])))
Gibbs free energy (J) = -15269944.699001541
component, activity, chemical potential, mu0
Water 0.999394 -274219.43 -274216.04
SiO2,aq 0.025098 -877459.59 -856835.60
Quartz -878750.13
species, activity, mu, mu0
Water 9.993938e-01 -274219.43 -274216.04
SiO2,aq 2.509769e-02 -877459.59 -856835.60
H+ 7.652831e-06 -65932.40 0.00
OH- 4.452470e-06 -208287.04 -139323.34
HSiO3- 3.200361e-06 -1085746.63 -1014934.90
Si2O4,aq 8.535953e-03 -1754919.18 -1728259.14
deltaG quartz -> SiO2 (aq) 1290.539609355852
at log (molality) = -1.374379e+00
print ("{0:>10s} {1:10.6f} {2:10.6f} {3:10.6f}".format('G TP TPr TrPr', qtz.getGibbsFreeEnergyFromT_andP_(t, p), qtz.getGibbsFreeEnergyFromT_andP_(t, 1.0), qtz.getGibbsFreeEnergyFromT_andP_(298.15, 1.0)))
print ("{0:>10s} {1:10.6f} {2:10.6f} {3:10.6f}".format('H TP TPr TrPr', qtz.getEnthalpyFromT_andP_(t, p), qtz.getEnthalpyFromT_andP_(t, 1.0), qtz.getEnthalpyFromT_andP_(298.15, 1.0)))
print ("{0:>10s} {1:10.6f} {2:10.6f} {3:10.6f}".format('S TP TPr TrPr', qtz.getEntropyFromT_andP_(t, p), qtz.getEntropyFromT_andP_(t, 1.0), qtz.getEntropyFromT_andP_(298.15, 1.0)))
print ("{0:>10s} {1:10.6f} {2:10.6f} {3:10.6f}".format('C TP TPr TrPr', qtz.getHeatCapacityFromT_andP_(t, p), qtz.getHeatCapacityFromT_andP_(t, 1.0), qtz.getHeatCapacityFromT_andP_(298.15, 1.0)))
HTP = qtz.getEnthalpyFromT_andP_(t, p)
HTrPr = qtz.getEnthalpyFromT_andP_(298.15, 1.0)
STP = qtz.getEntropyFromT_andP_(t,p)
STrPr = qtz.getEntropyFromT_andP_(298.15, 1.0)
GTP = qtz.getGibbsFreeEnergyFromT_andP_(t, p)
GTrPr = qtz.getGibbsFreeEnergyFromT_andP_(298.15, 1.0)
print ('G(T,P)-G(Tr,Pr) = ', GTP - GTrPr)
BermanG = -856288
BermanH = -910700
print ('Berman deltaG-deltaH formation at 298 K and 1 bar', BermanG - BermanH)
print ('BermanG + integral from PhaseObjC', BermanG+GTP-GTrPr)
print ('Dimitri at 400C and 1KB deltaG quartz = ', -210026.5*4.184)
print ('Dimitri at 400C and 1KB deltaG SiO2,aq = ', -204755.13*4.184)
print ('Dimitri at 400C and 1KB deltaG Si2O4 = ', -413054.54*4.184)
G TP TPr TrPr -878750.130116 -881047.101142 -856287.786356
H TP TPr TrPr -820003.153422 -822172.288621 -843929.046342
S TP TPr TrPr 87.271747 87.461654 41.451417
C TP TPr TrPr 67.531020 68.196114 44.844355
G(T,P)-G(Tr,Pr) = -22462.343760243733
Berman deltaG-deltaH formation at 298 K and 1 bar 54412
BermanG + integral from PhaseObjC -878750.3437602438
Dimitri at 400C and 1KB deltaG quartz = -878750.876
Dimitri at 400C and 1KB deltaG SiO2,aq = -856695.46392
Dimitri at 400C and 1KB deltaG Si2O4 = -1728220.19536