General Usage
Import pygcc
from importlib.metadata import version
print(version("pygcc"))
from pygcc.pygcc_utils import *
---------------------------------------------------------------------------
StopIteration Traceback (most recent call last)
File ~/.asdf/installs/python/3.11.14/lib/python3.11/importlib/metadata/__init__.py:563, in Distribution.from_name(cls, name)
562 try:
--> 563 return next(cls.discover(name=name))
564 except StopIteration:
StopIteration:
During handling of the above exception, another exception occurred:
PackageNotFoundError Traceback (most recent call last)
Cell In[1], line 2
1 from importlib.metadata import version
----> 2 print(version("pygcc"))
3
4 from pygcc.pygcc_utils import *
File ~/.asdf/installs/python/3.11.14/lib/python3.11/importlib/metadata/__init__.py:1009, in version(distribution_name)
1002 def version(distribution_name):
1003 """Get the version string for the named package.
1004
1005 :param distribution_name: The name of the distribution package to query.
1006 :return: The version string for the package as defined in the package's
1007 "Version" metadata key.
1008 """
-> 1009 return distribution(distribution_name).version
File ~/.asdf/installs/python/3.11.14/lib/python3.11/importlib/metadata/__init__.py:982, in distribution(distribution_name)
976 def distribution(distribution_name):
977 """Get the ``Distribution`` instance for the named package.
978
979 :param distribution_name: The name of the distribution package as a string.
980 :return: A ``Distribution`` instance (or subclass thereof).
981 """
--> 982 return Distribution.from_name(distribution_name)
File ~/.asdf/installs/python/3.11.14/lib/python3.11/importlib/metadata/__init__.py:565, in Distribution.from_name(cls, name)
563 return next(cls.discover(name=name))
564 except StopIteration:
--> 565 raise PackageNotFoundError(name)
PackageNotFoundError: No package metadata was found for pygcc
Read database by specifying the direct- or sequential- access and source database
Using the default sequential-access database - speq21.dat
ps = db_reader(sourcedb = 'thermo.2021', sourceformat = 'gwb')
# ps.dbaccessdic, ps.sourcedic, ps.specielist
Using the user-specified sequential-access database
ps = db_reader(dbaccess = './database/slop07.dat', sourcedb = 'thermo.2021', sourceformat = 'gwb')
# ps.dbaccessdic, ps.sourcedic, ps.specielist
Duplicate found for species "AlOH++" in slop07.dat
Duplicate found for species "MgCO3(aq)" in slop07.dat
Duplicate found for species "CaCO3(aq)" in slop07.dat
Duplicate found for species "SrCO3(aq)" in slop07.dat
Duplicate found for species "BaCO3(aq)" in slop07.dat
Duplicate found for species "BaF+" in slop07.dat
Duplicate found for species "Sr(Succ)(aq)" in slop07.dat
Duplicate found for species "Sc(Glut)+(aq)" in slop07.dat
Duplicate found for species "AMP2-" in slop07.dat
Duplicate found for species "HAMP-" in slop07.dat
Duplicate found for species "+H2AMP-" in slop07.dat
Example: Calculate water properties
With IAPWS95
water = iapws95(T = np.array([ 0.01, 25, 60, 100, 150, 200, 250, 300]), P = 200)
print('Density:', water.rho)
print('Gibbs Energy:', water.G)
print('Enthalpy:', water.H)
print('Entropy:', water.S)
Density: [1009.73580931 1005.83998999 991.7058825 967.43838406 927.69054221
877.9652082 816.08878805 734.71208475]
Gibbs Energy: [-56194.26867312 -56592.33889817 -57211.87993461 -58000.04190115
-59093.00868798 -60294.40876454 -61596.26593847 -62995.11966287]
Enthalpy: [-68680.33499723 -68236.29564253 -67613.18053963 -66897.34619192
-65991.9279271 -65062.66635227 -64087.83220987 -63021.29280159]
Entropy: [15.14005087 16.69550855 18.67153539 20.7005717 22.97720284 25.05223509
27.00976077 28.9549847 ]
With ZhangDuan
water = ZhangDuan(T= np.array([1000, 1050]), P = np.array([1000, 2000]))
print(water.rho, water.G)
[175.8286937 314.74510733] [-88767.03099 -89082.94223142]
Example: Calculate water dielectric constants
dielect = water_dielec(T= np.array([1000, 1050]), P = np.array([1000, 2000]), Dielec_method = 'DEW')
dielect.E, dielect.rhohat, dielect.Ah, dielect.Bh
(array([1.19637079, 2.52660344]),
array([0.17582869, 0.31474511]),
array([12.8720893 , 5.29642074]),
array([0.5403405 , 0.48797969]))
dielect = water_dielec(T= np.array([100, 150]), P = np.array([100, 200]), Dielec_method = 'JN91')
dielect.E, dielect.rhohat, dielect.Ah, dielect.Bh
(array([55.83017195, 44.79388023]),
array([0.96293375, 0.92769054]),
array([0.59551378, 0.67352582]),
array([0.34191435, 0.35183609]))
Example: Calculate quartz properties
With Maier-Kelly powerlaw
Temp = np.array([100, 200, 400, 600, 800, 900])
ps = db_reader()
sup = heatcap( T = Temp, P = 1, Species = 'Quartz', Species_ppt = ps.dbaccessdic['Quartz'],
method = 'SUPCRT')
sup.dG, sup.dCp
(array([-205319.43164294, -206728.65136032, -210399.56629393,
-215044.47020067, -220518.04896025, -223492.57970966]),
array([12.34074489, 13.89377788, 16.14397572, 16.103911 , 16.491911 ,
16.685911 ]))
With Holland and Power’s formulation
Temp = np.array([100, 200, 400, 600, 800, 900])
pshp = db_reader(dbHP_dir = './database/supcrtbl.dat')
hp = heatcap( T = Temp, P = 10, Species = 'Quartz', Species_ppt = pshp.dbaccessdic['Quartz'],
method = 'HP11')
hp.dG, hp.dCp
Duplicate found for species "Fluorphlogopite" in supcrtbl.dat
Duplicate found for species "Arsenic" in supcrtbl.dat
Duplicate found for species "Arsenolite" in supcrtbl.dat
(array([-205491.44716714, -206900.38310661, -210575.79120014,
-215228.32290096, -220682.44526921, -223649.25812766]),
array([12.4075447 , 13.99685509, 16.16426364, 16.05341748, 16.6660167 ,
16.90252028]))
With Berman’s formulation
Temp = np.array([100, 200, 400, 600, 800, 900])
ps = db_reader(dbBerman_dir = './database/berman.dat')
bm = heatcap( T = Temp, P = 1, Species = 'Quartz', Species_ppt = ps.dbaccessdic['Quartz'],
method = 'berman88')
bm.dG, bm.dCp
Duplicate found for species "Cordierite" in berman.dat
(array([-205498.95550508, -206907.11144614, -210575.29239247,
-215221.44900321, -220654.77893829, -223612.45803826]),
array([12.323813 , 13.87414095, 16.29997608, 16.24449507, 16.72930639,
16.90352372]))
import matplotlib.pyplot as plt
fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.plot(Temp, sup.dG*4.182/1e3, 'r', Temp, hp.dG*4.182/1e3, 'b', Temp, bm.dG*4.182/1e3, 'k')
ax1.set_xlabel('Temperature [°C]'); ax1.set_ylabel('Gibbs Energy [kJ/mol] ')
plt.legend(['SUPCRT', 'HP11', 'Berman88'])
ax2.plot(Temp, hp.dG*4.182/1e3 - sup.dG*4.182/1e3, 'r', Temp, bm.dG*4.182/1e3 - sup.dG*4.182/1e3, 'b')
ax2.set_xlabel('Temperature [°C]'); ax2.set_ylabel('Delta Gibbs Energy [kJ/mol] ')
plt.legend(['SUPCRT - HP11', 'SUPCRT - Berman88'])
fig.tight_layout()
Example: Calculate olivine solid solutions
calc = calcRxnlogK( X = 0.85,T = np.array([300, 400, 450]), P = np.array([200, 200, 200]),
Specie = 'olivine', dbaccessdic = ps.dbaccessdic, densityextrap = True)
calc.logK, calc.Rxn
(array([ 9.2646654 , 23.54844166, 19.7311314 ]),
{'type': 'ol',
'name': 'Fo85',
'formula': 'Mg1.70Fe0.30Si1O4',
'MW': 150.1557,
'min': ['Mg1.70Fe0.30Si1O4',
' R&H95, Stef2001',
np.float64(-466862.12264868076),
nan,
np.float64(26.21037219328013),
44.049,
24.058078393881452,
0.017393236137667304,
-890893.8814531548,
171.38145315487571,
-3.6586998087954105e-06],
'spec': ['H+', 'Mg++', 'Fe++', 'SiO2(aq)', 'H2O'],
'coeff': [-4, 1.7, 0.30000000000000004, 1, 2],
'nSpec': 5,
'V': 44.049,
'source': ' R&H95, Stef2001',
'elements': ['1.7000', 'Mg', '0.3000', 'Fe', '1.0000', 'Si', '4.0000', 'O']})
Example: Calculate glauconite mineral properties
Temp = np.array([0.01, 20, 30, 40, 60, 80, 100, 150])
calc = calcRxnlogK(T = Temp, P = 'T', Specie = 'clay', dbaccessdic = ps.dbaccessdic, densityextrap = True,
elem = ['Glauconite', '3.654', '.687', '1.175', '0.079', '0.337', '0.679', '0.212', '0.0359', '0'])
calc.logK, calc.Rxn
(array([ 7.33889672, 6.08136856, 5.34485756, 4.59351803, 3.11784876,
1.73061272, 0.45134995, -2.30826164]),
{'type': 'Smectite',
'name': 'Glauconite',
'formula': 'K0.68Na0.21Fe0.08Ca0.04Mg0.34Al0.69FeIII1.18Si3.654O10(OH)2',
'MW': 426.25271512,
'min': ['K0.68Na0.21Fe0.08Ca0.04Mg0.34Al0.69FeIII1.18Si3.654O10(OH)2',
'B2015, B2021',
np.float64(-1173001.1408482206),
np.float64(-1256479.2241755025),
np.float64(90.78782850627151),
np.float64(141.34069593680297),
np.float64(77.70879420984171),
np.float64(83.36446524554509),
np.float64(-17.841879278113826)],
'struct_formula': '(K0.679Na0.212Ca0.036)(Mg0.337Fe0.079Al0.341FeIII1.175)(Si3.654Al0.346)O10(OH)2',
'spec': ['H+',
'K+',
'Na+',
'Fe++',
'Ca++',
'Mg++',
'Al+++',
'Fe+++',
'SiO2(aq)',
'H2O'],
'coeff': [np.float64(-7.381),
0.679,
0.212,
0.079,
0.0359,
0.337,
0.687,
1.175,
3.654,
4.692],
'nSpec': 10,
'V': np.float64(141.34069593680297),
'dG': np.float64(-4907.836773308955),
'dHf': np.float64(-5257.109073950302),
'S': np.float64(379.85627447024),
'Cp': 345.17439398884756,
'source': 'B2015, B2021',
'elements': ['0.0359',
'Ca',
'0.6790',
'K',
'0.2120',
'Na',
'1.2540',
'Fe',
'0.3370',
'Mg',
'0.6870',
'Al',
'3.6540',
'Si',
'2.0000',
'H',
'12.0000',
'O']})
Example: Create new reactions and calculate equilibrium constants
An example with pyrite, pyrrhotite, magnetite (PPM) reaction
Pyrite + 4 H2O + 2 Pyrrhotite → 4 H2S(aq) + Magnetite
Then we can include the reaction in sourcedic, one of the output of db_reader, which is a dictionary of list of reaction coefficients and species. An example with the format for sourcedic is as follows:
ps.sourcedic[‘Name’] = [‘formula’, number of reactants in the reaction, ‘coefficient of specie 1’, ‘specie 1’, ‘coefficient of specie 2’, ‘specie 2’]
AB → 0.5 A2(aq) + B
ps.sourcedic[‘AB’] = [‘AB’, 2, ‘0.5’, ‘A2(aq), ‘1’, ‘B’]
ps = db_reader(sourcedb = './database/thermo.2021.dat', sourceformat = 'gwb')
ps.sourcedic['Pyrite'] = ['', 4, '4', 'H2S(aq)', '1', 'Magnetite', '-2', 'Pyrrhotite', '-4', 'H2O']
Temp = np.array([300.0000, 325, 350.0000, 400.0000, 415, 425.0000, 435, 450.0000])
Press = 500*np.ones(np.size(Temp))
log_K_PPM = calcRxnlogK( T = Temp, P = Press, Specie = 'Pyrite', dbaccessdic = ps.dbaccessdic,
sourcedic = ps.sourcedic, specielist = ps.specielist).logK
log_K_PPM
array([-9.61842122, -8.38742387, -7.21969753, -4.97772897, -4.29319566,
-3.81890541, -3.32234015, -2.52289322])
Example: Generate GWB thermodynamic database
Note: GWB requires exactly 8 temperature and pressure pairs
# Vectors for Temperature (C) and Pressure (bar) inputs
T = np.array([ 0.010, 25.0000 , 60.0000, 100.0000, 120.0000, 150.0000, 250.0000, 300.0000])
P = 350*np.ones(np.size(T))
nCa = 1
write GWB using default sourced database, with inclusion of solid_solution and clay thermo properties
%%time
write_database(T = T, P = P, cpx_Ca = nCa, solid_solution = 'Yes', clay_thermo = 'Yes',
dataset = 'GWB')
Success, your new GWB database is ready for download
CPU times: user 2.95 s, sys: 49.2 ms, total: 2.99 s
Wall time: 3.08 s
<pygcc.pygcc_utils.write_database at 0x10f3a8ec0>
write GWB using user-specified sourced database, with inclusion of solid_solution and clay thermo properties
%%time
src_db = os.path.join(os.getcwd(), 'database', 'thermo.29Sep15.dat')
write_database(T = T, P = 175, cpx_Ca = nCa, solid_solution = 'Yes', clay_thermo = 'Yes',
sourcedb = src_db, dataset = 'GWB')
Success, your new GWB database is ready for download
CPU times: user 2.49 s, sys: 38.6 ms, total: 2.53 s
Wall time: 2.67 s
<pygcc.pygcc_utils.write_database at 0x10f3c5310>
write GWB using Jan2020/Apr20 formatted sourced database
%%time
src_db = os.path.join(os.getcwd(), 'database', 'thermo.com.tdat')
write_database(T = T, P = 125, cpx_Ca = nCa, solid_solution = True, clay_thermo = True,
sourcedb = src_db, dataset = 'GWB')
Success, your new GWB database is ready for download
CPU times: user 2.84 s, sys: 31.2 ms, total: 2.87 s
Wall time: 2.9 s
<pygcc.pygcc_utils.write_database at 0x118fd7ed0>
write GWB using Jan2020/Apr20 formatted sourced database with logK as polynomial coefficients, using Tmax and Tmin
%%time
src_db = os.path.join(os.getcwd(), 'database', 'thermo.com.tdat')
write_database(T = [0, 300], P = 150, clay_thermo = True, logK_form = 'polycoeffs',
sourcedb = src_db, dataset = 'GWB')
Success, your new GWB database is ready for download
CPU times: user 2.63 s, sys: 28.2 ms, total: 2.65 s
Wall time: 2.69 s
<pygcc.pygcc_utils.write_database at 0x1192196e0>
write GWB using Mar21 formatted sourced database with logK as polynomial coefficients
%%time
src_db = os.path.join(os.getcwd(), 'database', 'thermo_latest.tdat')
write_database(T = [0, 300], P = 145, clay_thermo = True, logK_form = 'polycoeffs',
sourcedb = src_db, dataset = 'GWB')
Success, your new GWB database is ready for download
CPU times: user 2.67 s, sys: 33.6 ms, total: 2.7 s
Wall time: 2.76 s
<pygcc.pygcc_utils.write_database at 0x119219f30>
write GWB using user-specified sourced database and direct-access database (slop07) and FGL97 dielectric constant
%%time
db_db = os.path.join(os.getcwd(), 'database', 'slop07.dat')
src_db = os.path.join(os.getcwd(), 'database', 'thermo.29Sep15.dat')
write_database(T = [0, 400], P = 300, cpx_Ca = 0.5, solid_solution = 'Yes', Dielec_method = 'FGL97',
dbaccess = db_db, sourcedb = src_db, dataset = 'GWB')
Duplicate found for species "AlOH++" in slop07.dat
Duplicate found for species "MgCO3(aq)" in slop07.dat
Duplicate found for species "CaCO3(aq)" in slop07.dat
Duplicate found for species "SrCO3(aq)" in slop07.dat
Duplicate found for species "BaCO3(aq)" in slop07.dat
Duplicate found for species "BaF+" in slop07.dat
Duplicate found for species "Sr(Succ)(aq)" in slop07.dat
Duplicate found for species "Sc(Glut)+(aq)" in slop07.dat
Duplicate found for species "AMP2-" in slop07.dat
Duplicate found for species "HAMP-" in slop07.dat
Duplicate found for species "+H2AMP-" in slop07.dat
Success, your new GWB database is ready for download
CPU times: user 2.41 s, sys: 50 ms, total: 2.46 s
Wall time: 2.47 s
<pygcc.pygcc_utils.write_database at 0x10f3ef0b0>
write GWB using default sourced database and direct-access database (slop07 with Berman mineral data) and FGL97 dielectric constant
%%time
db_berman = os.path.join(os.getcwd(), 'database', 'berman.dat')
db_db = os.path.join(os.getcwd(), 'database', 'slop07.dat')
write_database(T = [0, 340], P = 150, Dielec_method = 'FGL97', dbaccess = db_db,
dbBerman_dir = db_berman, dataset = 'GWB')
Duplicate found for species "AlOH++" in slop07.dat
Duplicate found for species "MgCO3(aq)" in slop07.dat
Duplicate found for species "CaCO3(aq)" in slop07.dat
Duplicate found for species "SrCO3(aq)" in slop07.dat
Duplicate found for species "BaCO3(aq)" in slop07.dat
Duplicate found for species "BaF+" in slop07.dat
Duplicate found for species "Sr(Succ)(aq)" in slop07.dat
Duplicate found for species "Sc(Glut)+(aq)" in slop07.dat
Duplicate found for species "AMP2-" in slop07.dat
Duplicate found for species "HAMP-" in slop07.dat
Duplicate found for species "+H2AMP-" in slop07.dat
Duplicate found for species "Cordierite" in berman.dat
Success, your new GWB database is ready for download
CPU times: user 1.9 s, sys: 55.3 ms, total: 1.95 s
Wall time: 2.01 s
<pygcc.pygcc_utils.write_database at 0x10fc34490>
write GWB using Mar21 formatted sourced database and direct-access database (speq21 with HP mineral data) and DEW dielectric constant
%%time
write_database(T = [0, 275], P = 1100, Dielec_method = 'DEW', dbaccess = os.path.join(os.getcwd(), 'database', 'speq21.dat'),
dbHP_dir = os.path.join(os.getcwd(), 'database', 'supcrtbl.dat'),
sourcedb = os.path.join(os.getcwd(), 'database', 'thermo_latest.tdat'), dataset = 'GWB')
Duplicate found for species "Fluorphlogopite" in supcrtbl.dat
Duplicate found for species "Arsenic" in supcrtbl.dat
Duplicate found for species "Arsenolite" in supcrtbl.dat
/Users/noahawolayo/Documents/pygcc_build/src/pygcc/pygcc_utils.py:1924: UserWarning: Warning: Using SiO2(aq) data from Sverjensky et al. (2014) GCA. This thermodynamic model requires that you also add data for Si2O4(aq) to achieve accurate solubility calculations.Please ensure thermodynamic data for Si2O4(aq) are added to both your source GWB or EQ3/6 database AND your source direct-access database.
warnings.warn(
Success, your new GWB database is ready for download
CPU times: user 1.48 s, sys: 74.2 ms, total: 1.55 s
Wall time: 1.63 s
<pygcc.pygcc_utils.write_database at 0x10fc345a0>
write GWB using user-specified sourced Pitzer database and default direct-access database along the saturation curve
%%time
write_database(T = [0, 350], P = 'T', dataset = 'GWB', sourcedb = os.path.join(os.getcwd(), 'database', 'thermo_hmw.tdat') )
Success, your new GWB database is ready for download
CPU times: user 869 ms, sys: 22.2 ms, total: 892 ms
Wall time: 931 ms
<pygcc.pygcc_utils.write_database at 0x118eb5ae0>
write GWB using user-specified sourced EQ3/6 Pitzer database and default direct-access database
%%time
write_database(T = [0, 350], P = 250, dataset = 'GWB', sourcedb = os.path.join(os.getcwd(), 'database', 'data0.hmw'),
sourceformat = 'EQ36')
Success, your new GWB database is ready for download
CPU times: user 749 ms, sys: 10.5 ms, total: 759 ms
Wall time: 773 ms
<pygcc.pygcc_utils.write_database at 0x119216580>
Example: Generate EQ3/6 thermodynamic database
write EQ3/6 using default sourced database
%%time
write_database(T = T, P = P, cpx_Ca = 1, solid_solution = 'Yes', clay_thermo = 'Yes',
dataset = 'EQ36')
Success, your new EQ3/6 database is ready for download
CPU times: user 3.03 s, sys: 128 ms, total: 3.16 s
Wall time: 5.09 s
<pygcc.pygcc_utils.write_database at 0x119217bd0>
write EQ3/6 user-specified sourced database
%%time
write_database(T = [0, 400], P = 350, cpx_Ca = 0.5, solid_solution = 'Yes', clay_thermo = 'Yes',
sourcedb = os.path.join(os.getcwd(), 'database', 'data0.geo'), dataset = 'EQ36', sourcedb_codecs = 'latin-1')
Success, your new EQ3/6 database is ready for download
CPU times: user 2.84 s, sys: 90.1 ms, total: 2.93 s
Wall time: 3.09 s
<pygcc.pygcc_utils.write_database at 0x119217ac0>
write EQ3/6 using user-specified sourced Pitzer database
%%time
write_database(T = [0, 350], P = 200, sourcedb = os.path.join(os.getcwd(), 'database', 'data0.hmw'), dataset = 'EQ36')
Success, your new EQ3/6 database is ready for download
CPU times: user 717 ms, sys: 22 ms, total: 739 ms
Wall time: 746 ms
<pygcc.pygcc_utils.write_database at 0x119217240>
write EQ3/6 user-specified sourced database using FGL97 dielectric constant
%%time
write_database(T = [0, 400], P = 300, cpx_Ca = 0.1, sourcedb = os.path.join(os.getcwd(), 'database', 'data0.geo'),
dataset = 'EQ36', solid_solution = 'Yes', Dielec_method = 'FGL97', clay_thermo = 'Yes')
Success, your new EQ3/6 database is ready for download
CPU times: user 2.76 s, sys: 31.9 ms, total: 2.8 s
Wall time: 2.88 s
<pygcc.pygcc_utils.write_database at 0x104237350>
write EQ3/6 user-specified sourced database using DEW model
%%time
Temp = np.array([50, 100, 150, 300, 450, 500, 600, 700])
write_database(T = Temp, P = 1500, sourcedb = os.path.join(os.getcwd(), 'database', 'data0.geo'), dataset = 'EQ36',
Dielec_method = 'DEW')
Success, your new EQ3/6 database is ready for download
CPU times: user 1.39 s, sys: 29.5 ms, total: 1.42 s
Wall time: 1.45 s
<pygcc.pygcc_utils.write_database at 0x11a784d10>
write EQ3/6 using default sourced GWB database and direct-access database
%%time
write_database(T = np.array([0.010, 25, 60, 100, 150, 200, 250, 300]), P = 250, dataset = 'EQ36',
sourcedb = 'thermo.2021', sourceformat = 'gwb', solid_solution = True, clay_thermo = True,
print_msg = True)
Success, your new EQ3/6 database is ready for download
Database for EQ36 generated successfully using JN91 dielectric constant, thermo.2021.dat gwb source database, solid solution included with cpx excluded, clay thermodynamics included
CPU times: user 11.3 s, sys: 172 ms, total: 11.5 s
Wall time: 11.8 s
<pygcc.pygcc_utils.write_database at 0x11a7858c0>
Example: Generate PHREEQC thermodynamic database
write PHREEQC using default sourced database
%%time
write_database(T = T, P = P, cpx_Ca = 1, solid_solution = 'Yes', clay_thermo = 'Yes', dataset = 'PHREEQC')
Success, your new PHREEQC database is ready for download
CPU times: user 1.9 s, sys: 17.7 ms, total: 1.92 s
Wall time: 1.92 s
<pygcc.pygcc_utils.write_database at 0x11a7857b0>
write PHREEQC user-specified sourced database
%%time
write_database(T = [0, 600], P = 350, cpx_Ca = 0.5, solid_solution = True, clay_thermo = True, sourceformat = 'PHREEQC',
sourcedb = os.path.join(os.getcwd(), 'database', 'PHREEQC_ThermoddemV1.10_15Dec2020.dat'), dataset = 'PHREEQC')
/Users/noahawolayo/Documents/pygcc_build/src/pygcc/pygcc_utils.py:4891: UserWarning: Some temperature and pressure points are out of aqueous species HKF eqns regions of applicability, hence, density extrapolation has been applied
warnings.warn('Some temperature and pressure points are out of aqueous species HKF eqns regions of applicability, hence, density extrapolation has been applied')
Success, your new PHREEQC database is ready for download
CPU times: user 1min 9s, sys: 725 ms, total: 1min 10s
Wall time: 1min 13s
<pygcc.pygcc_utils.write_database at 0x11a786030>
write PHREEQC using user-specified sourced Pitzer database
%%time
write_database(T = [0, 350], P = 200, sourcedb = os.path.join(os.getcwd(), 'database', 'pitzer.dat'), sourceformat = 'PHREEQC', dataset = 'PHREEQC')
Success, your new PHREEQC database is ready for download
CPU times: user 693 ms, sys: 16.7 ms, total: 710 ms
Wall time: 811 ms
<pygcc.pygcc_utils.write_database at 0x11a787ac0>
write PHREEQC user-specified sourced database using FGL97 dielectric constant
%%time
write_database(T = [0, 600], P = 300, cpx_Ca = 0.1, sourcedb = os.path.join(os.getcwd(), 'database', 'PHREEQC_ThermoddemV1.10_15Dec2020.dat'), dataset = 'PHREEQC',
solid_solution = True, Dielec_method = 'FGL97', clay_thermo = True, sourceformat = 'PHREEQC')
Success, your new PHREEQC database is ready for download
CPU times: user 1min 4s, sys: 417 ms, total: 1min 5s
Wall time: 1min 6s
<pygcc.pygcc_utils.write_database at 0x10f5046b0>
write PHREEQC user-specified sourced database using DEW model
%%time
T = np.array([50, 100, 150, 300, 450, 500, 600, 700])
write_database(T = T, P = 1500, sourcedb = os.path.join(os.getcwd(), 'database', 'PHREEQC_ThermoddemV1.10_15Dec2020.dat'),
sourceformat = 'PHREEQC', dataset = 'PHREEQC', Dielec_method = 'DEW')
Success, your new PHREEQC database is ready for download
CPU times: user 3.3 s, sys: 50 ms, total: 3.35 s
Wall time: 3.64 s
<pygcc.pygcc_utils.write_database at 0x1185a1e10>
Example: Generate ToughReact thermodynamic database
write ToughReact using user-specified EQ3/6 database
%%time
write_database(T = T, P = P, cpx_Ca = nCa, solid_solution = 'Yes', sourcedb = os.path.join(os.getcwd(), 'database', 'data0.dat'),
dataset = 'ToughReact', sourceformat = 'EQ36')
/Users/noahawolayo/Documents/pygcc_build/src/pygcc/pygcc_utils.py:5750: UserWarning: Some temperature and pressure points are out of aqueous species HKF eqns regions of applicability, hence, density extrapolation has been applied
warnings.warn('Some temperature and pressure points are out of aqueous species HKF eqns regions of applicability, hence, density extrapolation has been applied')
Success, your new ToughReact database is ready for download
CPU times: user 43.4 s, sys: 437 ms, total: 43.8 s
Wall time: 44.6 s
<pygcc.pygcc_utils.write_database at 0x11d36ae00>
write ToughReact using user-specified GWB database
%%time
write_database(T = [0, 350], P = 250, cpx_Ca = 0.25, clay_thermo = 'Yes', dataset = 'ToughReact',
sourcedb = os.path.join(os.getcwd(), 'database', 'thermo.com.tdat'), sourceformat = 'GWB')
Success, your new ToughReact database is ready for download
CPU times: user 1.62 s, sys: 20.2 ms, total: 1.64 s
Wall time: 1.67 s
<pygcc.pygcc_utils.write_database at 0x1185a0490>
write ToughReact using EQ3/6 user-specified Ptizer sourced database and JN91 dielectric constant
%%time
write_database(T = [0, 300], P = 200, sourceformat = 'EQ36', sourcedb = os.path.join(os.getcwd(), 'database', 'data0.fmt'),
dataset = 'ToughReact', Dielec_method = 'JN91')
Success, your new ToughReact database is ready for download
CPU times: user 865 ms, sys: 5.8 ms, total: 871 ms
Wall time: 876 ms
<pygcc.pygcc_utils.write_database at 0x11d36b8a0>
Example: Generate Pflotran thermodynamic database
write Pflotran using user-specified EQ3/6 database
%%time
write_database(T = T, P = P, clay_thermo = 'Yes', sourcedb = os.path.join(os.getcwd(), 'database', 'data0.dat'),
dataset = 'Pflotran', sourceformat = 'EQ36')
/Users/noahawolayo/Documents/pygcc_build/src/pygcc/pygcc_utils.py:5313: UserWarning: Some temperature and pressure points are out of aqueous species HKF eqns regions of applicability, hence, density extrapolation has been applied
warnings.warn('Some temperature and pressure points are out of aqueous species HKF eqns regions of applicability, hence, density extrapolation has been applied')
Success, your new Pflotran database is ready for download
CPU times: user 37.4 s, sys: 455 ms, total: 37.8 s
Wall time: 38.4 s
<pygcc.pygcc_utils.write_database at 0x11e358160>
write Pflotran using user-specified GWB database
%%time
write_database(T = [0, 350], P = 250, cpx_Ca = 0.1, solid_solution = True, clay_thermo = True,
sourcedb = os.path.join(os.getcwd(), 'database', 'thermo.com.tdat'), dataset = 'Pflotran', sourceformat = 'GWB')
Success, your new Pflotran database is ready for download
CPU times: user 2.64 s, sys: 20 ms, total: 2.66 s
Wall time: 2.69 s
<pygcc.pygcc_utils.write_database at 0x10f507680>
Example: Calculate clay mineral thermodynamics
#%% specify the direct access thermodynamic database
db_dic = db_reader(dbaccess = os.path.join(os.getcwd(), 'database', 'speq21.dat'))
db_dic = db_dic.dbaccessdic
folder_to_save = 'output'
if os.path.exists(os.path.join(os.getcwd(), folder_to_save)) == False:
os.makedirs(os.path.join(os.getcwd(), folder_to_save))
fid = open('./output/logK_test.txt', 'w')
T = np.array([ 0.01, 25, 60, 100, 150, 200, 250, 300])
TK_range = T + 273.15
logKRxn = calcRxnlogK(T = T, P = 'T', Specie = 'Clay', dbaccessdic = db_dic,
elem = ['Clinochlore', '3', '2', '0', '0', '5', '0', '0', '0', '0', '0'],
densityextrap = True)
logK, Rxn = logKRxn.logK, logKRxn.Rxn
# output in EQ36 format
outputfmt(fid, logK, Rxn, dataset = 'EQ36')
# output in GWB format
outputfmt(fid, logK, Rxn, dataset = 'GWB')
# output in PHREEQC format
outputfmt(fid, logK, Rxn, *TK_range, dataset = 'PHREEQC')
# output in Pflotran format
outputfmt(fid, logK, Rxn, dataset = 'Pflotran')
# output in ToughReact format
outputfmt(fid, logK, Rxn, dataset = 'ToughReact')
fid.close()
Example: Calculate plagioclase solid-solution thermodynamics
#%% specify the direct access thermodynamic database
db_dic = db_reader(dbaccess = os.path.join(os.getcwd(), 'database', 'speq21.dat'))
db_dic = db_dic.dbaccessdic
folder_to_save = 'output'
if os.path.exists(os.path.join(os.getcwd(), folder_to_save)) == False:
os.makedirs(os.path.join(os.getcwd(), folder_to_save))
fid = open('./output/logK_test.txt', 'w')
T = np.array([ 0.01, 25, 60, 100, 150, 200, 250, 300])
logKRxn = calcRxnlogK(T = T, dbaccessdic = db_dic, P = 'T', X = 0.634, Specie = 'Plagioclase')#densityextrap = True)
logK, Rxn = logKRxn.logK, logKRxn.Rxn
logKRxnph = calcRxnlogK(T = T, P = 'T', X = 0.634, Specie = 'Plagioclase',
dbaccessdic = db_dic, densityextrap = True, Al_Si = 'Arnórsson_Stefánsson')
# output in EQ36 format
outputfmt(fid, logK, Rxn, dataset = 'EQ36')
# output in GWB format
outputfmt(fid, logK, Rxn, dataset = 'GWB')
# output in PHREEQC format
outputfmt(fid, logKRxnph.logK, logKRxnph.Rxn, *TK_range, dataset = 'PHREEQC')
# output in Pflotran format
outputfmt(fid, logK, Rxn, dataset = 'Pflotran')
# output in ToughReact format
outputfmt(fid, logK, Rxn, dataset = 'ToughReact')
fid.close()
Example: Calculate new reaction equilibrium constants and export to specific formats
#%% specify the direct access thermodynamic database
db_dic = db_reader(dbaccess = os.path.join(os.getcwd(), 'database', 'speq21.dat'),
sourcedb = os.path.join(os.getcwd(), 'database', 'thermo_latest.tdat'), sourceformat = 'GWB')
Temp = np.array([25, 80, 140, 190, 240, 300, 350, 400])
TK_range = Temp + 273.15
folder_to_save = 'output'
if os.path.exists(os.path.join(os.getcwd(), folder_to_save)) == False:
os.makedirs(os.path.join(os.getcwd(), folder_to_save))
fid = open('./output/logK_test.txt', 'w')
db_dic.sourcedic['HWO4-'] = ['HWO4-', 2, '1.000', 'H+', '1.000', 'WO4--']
logK = calcRxnlogK(T = Temp, P = 250, Specie = 'HWO4-', dbaccessdic = db_dic.dbaccessdic, densityextrap = True,
sourcedic = db_dic.sourcedic, specielist = db_dic.specielist, Specie_class = 'aqueous').logK
Rxn = {}
Rxn['type'] = 'tungsten'
Rxn['name'] = 'HWO4-'
Rxn['formula'] = ps.dbaccessdic['HWO4-'][0]
Rxn['MW'] = calc_elem_count_molewt(ps.dbaccessdic['HWO4-'][0])[1]
Rxn['min'] = ps.dbaccessdic['HWO4-']
Rxn['V'] = Rxn['min'][5]
Rxn['source'] = Rxn['min'][1]
Rxn['spec'] = ['H+', 'WO4--']
Rxn['coeff'] = [1.000, 1.000]
Rxn['nSpec'] = len(Rxn['coeff'])
d = calc_elem_count_molewt(ps.dbaccessdic['HWO4-'][0])[0]
Rxn['elements'] = [item for sublist in [['%s' % v,k] for k,v in d.items()] for item in sublist]
# output in EQ36 format
outputfmt(fid, logK, Rxn, dataset = 'EQ36')
# output in GWB format
outputfmt(fid, logK, Rxn, *TK_range, logK_form = 'polycoeffs', dataset = 'GWB')
# output in PHREEQC format
outputfmt(fid, logK, Rxn, *TK_range, dataset = 'PHREEQC')
# output in Pflotran format
outputfmt(fid, logK, Rxn, dataset = 'Pflotran')
# output in ToughReact format
outputfmt(fid, logK, Rxn, dataset = 'ToughReact')
fid.close()
Example: Calculate CO$_2$ activity and molality
T = np.array([ 0.010, 25 , 60, 100, 150, 175, 200, 250])
P = 250*np.ones(np.size(T))
TK = convert_temperature(T, Out_Unit = 'K')
#%% Calculate CO2 activity and molality at ionic strength of 0.5M
# with Duan_Sun
log10_co2_gamma, mco2 = Henry_duan_sun(TK, P, 0.5)
co2_activity = 10**log10_co2_gamma
# with Drummond
log10_co2_gamma = drummondgamma(TK, 0.5)
co2_activityD = 10**log10_co2_gamma
for i in range(len(T)):
print('Fluid Temperature [C]: ', T[i])
print('Fluid Pressure [bar]: ', P[i])
print('Molality of CO2 in aqueous phase: ', mco2.ravel()[i])
print('Activity of CO2 in aqueous phase (Duan_Sun): ', co2_activity.ravel()[i])
print('Activity of CO2 in aqueous phase (Drummond): ', co2_activityD.ravel()[i])
print('\n')
Fluid Temperature [C]: 0.01
Fluid Pressure [bar]: 250.0
Molality of CO2 in aqueous phase: 1.9318448998874207
Activity of CO2 in aqueous phase (Duan_Sun): 1.1291811186477334
Activity of CO2 in aqueous phase (Drummond): 1.1340282640704162
Fluid Temperature [C]: 25.0
Fluid Pressure [bar]: 250.0
Molality of CO2 in aqueous phase: 1.4446373503231484
Activity of CO2 in aqueous phase (Duan_Sun): 1.117606348338049
Activity of CO2 in aqueous phase (Drummond): 1.1228777554452154
Fluid Temperature [C]: 60.0
Fluid Pressure [bar]: 250.0
Molality of CO2 in aqueous phase: 1.1670235810670366
Activity of CO2 in aqueous phase (Duan_Sun): 1.1097025985737983
Activity of CO2 in aqueous phase (Drummond): 1.1184645553679928
Fluid Temperature [C]: 100.0
Fluid Pressure [bar]: 250.0
Molality of CO2 in aqueous phase: 1.0886091317050277
Activity of CO2 in aqueous phase (Duan_Sun): 1.109458812741491
Activity of CO2 in aqueous phase (Drummond): 1.1250331593722136
Fluid Temperature [C]: 150.0
Fluid Pressure [bar]: 250.0
Molality of CO2 in aqueous phase: 1.1736136778845458
Activity of CO2 in aqueous phase (Duan_Sun): 1.1191485747389718
Activity of CO2 in aqueous phase (Drummond): 1.1457708279040804
Fluid Temperature [C]: 175.0
Fluid Pressure [bar]: 250.0
Molality of CO2 in aqueous phase: 1.2783771607594785
Activity of CO2 in aqueous phase (Duan_Sun): 1.1276219378379773
Activity of CO2 in aqueous phase (Drummond): 1.1602093887224851
Fluid Temperature [C]: 200.0
Fluid Pressure [bar]: 250.0
Molality of CO2 in aqueous phase: 1.4208468773780791
Activity of CO2 in aqueous phase (Duan_Sun): 1.1385618727912967
Activity of CO2 in aqueous phase (Drummond): 1.1769259206651324
Fluid Temperature [C]: 250.0
Fluid Pressure [bar]: 250.0
Molality of CO2 in aqueous phase: 1.77261489074432
Activity of CO2 in aqueous phase (Duan_Sun): 1.1700221965116315
Activity of CO2 in aqueous phase (Drummond): 1.2163347543456382
Example: Calculate water activity
#%% Calculate Water activity, osmotic coefficient and NaCl mean activity coefficient at an ionic strength of 0.5M
aw, phi, mean_act = Helgeson_activity(T, P, 0.5, Dielec_method = 'JN91')
for i in range(len(T)):
print('Fluid Temperature [C]: ', T[i])
print('Fluid Pressure [bar]: ', P[i])
print('Water activity: ', aw.ravel()[i])
print('Water osmotic coefficient: ', phi.ravel()[i])
print('NaCl mean activity coefficient: ', mean_act.ravel()[i])
print('\n')
Fluid Temperature [C]: 0.01
Fluid Pressure [bar]: 250.0
Water activity: 0.9832827062135284
Water osmotic coefficient: 0.9357947724582291
NaCl mean activity coefficient: 0.7028482566450444
Fluid Temperature [C]: 25.0
Fluid Pressure [bar]: 250.0
Water activity: 0.9834680138300679
Water osmotic coefficient: 0.925334742394683
NaCl mean activity coefficient: 0.6831227245014528
Fluid Temperature [C]: 60.0
Fluid Pressure [bar]: 250.0
Water activity: 0.9834674900419101
Water osmotic coefficient: 0.925364305805263
NaCl mean activity coefficient: 0.6731415110785025
Fluid Temperature [C]: 100.0
Fluid Pressure [bar]: 250.0
Water activity: 0.9836234739653248
Water osmotic coefficient: 0.9165610286207091
NaCl mean activity coefficient: 0.6470401472320056
Fluid Temperature [C]: 150.0
Fluid Pressure [bar]: 250.0
Water activity: 0.9839617696261588
Water osmotic coefficient: 0.8974734053857932
NaCl mean activity coefficient: 0.601472728601265
Fluid Temperature [C]: 175.0
Fluid Pressure [bar]: 250.0
Water activity: 0.9841776700481395
Water osmotic coefficient: 0.8852951070804821
NaCl mean activity coefficient: 0.5748272882994092
Fluid Temperature [C]: 200.0
Fluid Pressure [bar]: 250.0
Water activity: 0.9844315601276349
Water osmotic coefficient: 0.870977343196605
NaCl mean activity coefficient: 0.5454072718625895
Fluid Temperature [C]: 250.0
Fluid Pressure [bar]: 250.0
Water activity: 0.9850902029390286
Water osmotic coefficient: 0.8338513425870471
NaCl mean activity coefficient: 0.4767054502731833
Example: Calculate salt water properties
from pygcc.pygcc_utils import Driesner_NaCl, concentration_converter
water_salt = Driesner_NaCl(T = 200, P = 500, xNaCl = concentration_converter(val = 0.5, In_Unit='x_mol_kgsol', Out_Unit='x_mol'))
print('Pressure of vapor + liquid + halite coexistence [bar]: ', water_salt.PVLH)
print('Composition of halite-saturated liquid (halite liquidus) [-]: ', water_salt.xL_NaCl)
print('Composition of halite-saturated vapor [-]: ', water_salt.xV_NaCl)
print('Composition of liquid at vapor + liquid coexistence [-]: ', water_salt.xVL_Liq)
print('Composition of vapor at vapor+liquid coexistence [-]: ', water_salt.xVL_Vap)
print('Liquid NaCl density [kg/m3]: ', water_salt.rho)
print('molar volume [mol/m3]: ', water_salt.vm)
print('Specific enthalpy of an H2O–NaCl solution [J/kg]: ', water_salt.H)
print('viscosity [Pa-s]: ', water_salt.mu)
print('Isobaric heat capacity [J/kg/K]: ', water_salt.Cp)
Pressure of vapor + liquid + halite coexistence [bar]: 11.151719135466353
Composition of halite-saturated liquid (halite liquidus) [-]: 0.12609895126323659
Composition of halite-saturated vapor [-]: nan
Composition of liquid at vapor + liquid coexistence [-]: [-0.30583232]
Composition of vapor at vapor+liquid coexistence [-]: [nan]
Liquid NaCl density [kg/m3]: 0.918850325064366
molar volume [mol/m3]: 20.010808061663727
Specific enthalpy of an H2O–NaCl solution [J/kg]: 856.75482704851
viscosity [Pa-s]: 0.00014914486904138454
Isobaric heat capacity [J/kg/K]: 4161.461344041091