Source code for pygcc.pygcc_utils

#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Created on Wed Mar 17 16:02:22 2021

@author: Adedapo Awolayo and Ben Tutolo, University of Calgary

Copyright (c) 2020 - 2021, Adedapo Awolayo and Ben Tutolo, University of Calgary

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.

"""
from .read_db import db_reader
from .water_eos import iapws95, ZhangDuan, water_dielec, readIAPWS95data, convert_temperature
from .species_eos import heatcap, supcrtaq
from .solid_solution import solidsolution_thermo
from .clay_thermocalc import calclogKclays
import warnings
import re
import os
import json
import numpy as np
import pandas as pd
import time
from scipy.optimize import fsolve, curve_fit
from scipy.linalg import lu_factor, lu_solve
from scipy.interpolate import splev, splrep, Rbf
import inspect
import math
from collections import OrderedDict
np.random.seed(4321)
warnings.filterwarnings("ignore", message="divide by zero encountered")
warnings.filterwarnings("ignore", message="invalid value encountered")

eps = 2.220446049250313e-16
J_to_cal = 4.184
IAPWS95_COEFFS = readIAPWS95data()

[docs]def calc_elem_count_molewt(formula, **kwargs): """ This function calculates the molecular mass and the elemental composition of a substance given by its chemical formula. This is modified from https://github.com/cgohlke/molmass/blob/master/molmass/molmass.py Parameters ---------- formula : string chemical formula Elementdic : dict dictionary, containing the atomic mass of element database Returns ---------- elements : dict dictionary of elemental composition and their respective number of atoms molewt : float Calculated Molecular Weights [g/mol] Usage: ---------- [elements, molewt] = calc_elem_count_molewt(formula) Examples of valid formulas are "H2O", "[2H]2O", "CH3COOH", "EtOH", "CuSO4:5H2O", "(COOH)2", "AgCuRu4(H)2[CO]12{PPh3}2", "CGCGAATTCGCG", and, "MDRGEQGLLK" . Examples -------- >>> elements, molewt = calc_elem_count_molewt("CuSO4:5H2O") >>> elements, molewt {'O': 9.0, 'H': 10.0, 'S': 1, 'Cu': 1}, 249.6773 """ kwargs = dict({"Elementdic": None }, **kwargs) Elementdic = kwargs['Elementdic'] if Elementdic is None: periodic_table = open(os.path.join(os.path.dirname(os.path.abspath(__file__)), 'PeriodicTableJSON.json'), encoding='utf8') data = json.load(periodic_table) Elementdic = {data['elements'][x]['symbol'] : pd.DataFrame([data['elements'][x]['name'], data['elements'][x]['atomic_mass']], index = ['name', 'mass']).T for x in range(len(data['elements']))} periodic_table.close() validchars = set('([{<123456789ABCDEFGHIKLMNOPRSTUVWXYZ') validchars |= set(']})>0abcdefghiklmnoprstuy') elements = {} ele = '' # parsed element num = 0 # number level = 0 # parenthesis level counts = [1] # parenthesis level multiplication i = len(formula) while i: i -= 1 char = formula[i] if char in '([{<': level -= 1 elif char in ')]}>': if num == 0: num = 1 level += 1 if level > len(counts) - 1: counts.append(0) counts[level] = num * counts[level - 1] num = 0 elif char.isdigit(): j = i while i and (formula[i - 1].isdigit() or formula[i - 1] == '.'): i -= 1 num = float(formula[i : j + 1]) elif char.islower(): ele = char elif char.isupper(): ele = char + ele if num == 0: num = 1 j = i number = num * counts[level] if ele in elements.keys(): elements[ele] = number + elements[ele] else: elements[ele] = number ele = '' num = 0 elif char == ':': if num == 0: num = 1 for k in elements.keys(): elements[k] = elements[k]*num molewt = 0 for symbol in elements: molewt += Elementdic[symbol].mass[0] * elements[symbol] return elements, molewt
[docs]def importconfile(filename, *Rows): """ This function imports numeric data from a text file as column vectors. [Var1, Var2] = importconfile(filename) Reads data from text file filename for the default selection. [Var1, Var2] = importconfile(filename, StartRow, EndRow) Reads data from rows StartRow through EndRow of text file filename. Example: [Var1, Var2] = importconfile('100bar.con', 6, 13); """ # %% Initialize variables. if (len(Rows) == 0): startRow = 5 endRow = np.inf else: startRow, endRow = Rows # %% Open the text file. fileID = open(filename,'r'); df=fileID.readlines() Var1 = [] Var2 = [] for block in range(startRow,len(df)): if not df[block].startswith('\n'): dataArray = df[block].split() Var1.append(float(dataArray[0])) Var2.append(float(dataArray[1])) Var1 = np.asarray(Var1) Var2 = np.asarray(Var2) # %% Close the text file. fileID.close() # %% return Var1, Var2
[docs]def var_name(var): callers_local_vars = inspect.currentframe().f_back.f_locals.items() print(str([k for k, v in callers_local_vars if v is var][0])+': '+str(var))
[docs]def feval(funcName, *args): return eval(funcName)(*args)
[docs]def roundup_tenth(x): return int(math.ceil(x / 10.0)) * 10
[docs]def roundup_hundredth(x): return int(math.ceil(x / 100)) * 100
[docs]def read_specific_lines(file, lines_to_read): lines = set(lines_to_read) last = max(lines) for n, line in enumerate(file): if n + 1 in lines: yield line if n + 1 > last: return if not line: continue
[docs]def derivative(f, a, method = 'central', h = 0.001): '''Compute the derivative of f, f'(a) with step size h. Parameters ---------- f : function Vectorized function of one variable a : number Compute derivative at x = a method : string Difference formula: 'forward', 'backward' or 'central' h : number Step size in difference formula Returns ------- float Difference formula: central: f(a + h) - f(a - h))/2h forward: f(a + h) - f(a))/h backward: f(a) - f(a-h))/h ''' if method == 'central': return (f(a + h) - f(a - h))/(2*h) elif method == 'forward': return (f(a + h) - f(a))/h elif method == 'backward': return (f(a) - f(a - h))/h else: raise ValueError("Method must be 'central', 'forward' or 'backward'.")
[docs]def info(name, dic): """ This function checks for naming convention of the species in the direct-access or source database Parameters ---------- name : string species name dic : dict dictionary of species from direct-access or source database Returns ---------- lst : list resulting search of all species with the input name Examples -------- >>> from pygcc.pygcc_utils import db_reader >>> ps = db_reader() # utilizes the default direct-access database, speq21 >>> info('ss_', ps.dbaccessdic) ['ss_Anorthite', 'ss_Albite_high', 'ss_K-feldspar', 'ss_Ferrosilite', 'ss_Enstatite', 'ss_Clinoenstatite', 'ss_Hedenbergite', 'ss_Diopside', 'ss_Forsterite', 'ss_Fayalite', 'ss_Annite', 'ss_Phlogopite', 'ss_Anorthite', 'ss_Albite_high', 'ss_K-feldspar', 'ss_Ferrosilite', 'ss_Enstatite', 'ss_Clinoenstatite', 'ss_Hedenbergite', 'ss_Diopside', 'ss_Forsterite', 'ss_Fayalite', 'ss_Annite', 'ss_Phlogopite'] """ lst = [i for i in list(dic.keys()) if i.startswith(name)==True] # starts with name lst = lst + [i for i in list(dic.keys()) if i.__contains__(name)] # contains name return lst
#%%------------------------------------------------------------------------------------
[docs]def drummondgamma(TK, I): """ This function models solubility of CO2 gas in brine using Drummond equation Parameters ---------- TK : float, vector Temperature [K] I : float, vector Ionic strength Returns ---------- log10_gamma : float, vector co2 aqueous activity coefficients in log10 Examples -------- >>> log10_gamma = drummondgamma( 500, 0.5) >>> log10_gamma 0.0781512920184902 """ # A=20.244 # B=-0.016323 C=-1.0312 # D=-3629.7 E=0.4445 F=0.0012806 G=255.9 H=-0.001606 #gamma=(((C + F*TK + (G/TK))*I - (E + H*TK)*(I / (I + 1)))/np.log(10)) log10_gamma = np.log10(np.exp((C + F*TK + (G/TK))*I - (E + H*TK)*(I / (I + 1)))) return log10_gamma
[docs]def Henry_duan_sun(TK, P, I): """ This function evaluates the solubility of CO2 gase in brine using Duan_Sun Formulation Parameters ---------- TK : float, vector Temperature [K] P : float, vector Pressure [bar] I : float Ionic strength Returns ---------- log10_co2_gamma : float, vector co2 aqueous activity coefficients in log10 mco2 : float, vector co2 aqueous molalities Usage ---------- log10_co2_gamma, mco2 = Henry_duan_sun( TK, P, I) Examples -------- >>> log10_co2_gamma, mco2 = Henry_duan_sun( 500, 250, 0.5) >>> log10_co2_gamma, mco2 0.06202532, 1.6062557 References ---------- (1) Duan, Z. and Sun, R., 2003. An improved model calculating CO2 solubility in pure water and aqueous NaCl solutions from 273 to 533 K and from 0 to 2000 bar. Chemical geology, 193(3-4), pp.257-271. """ if (np.ndim(TK) == 0): TK = np.array(TK).ravel() if (np.ndim(P) == 0): P = np.array(P).ravel() if (np.ndim(I) == 0): I = np.array(I).ravel() if np.size(TK) < np.size(P): TK = TK*np.ones_like(P) if np.size(P) < np.size(TK): P = P*np.ones_like(TK) # calculate CO2 fugacity and partial pressure # vco2fugacity = np.vectorize(co2fugacity) [fCO2, denCO2, SI, pCO2] = co2fugacity(TK, P) # Equation B1 of Duan & Sun model to calculate vapor pressure of water c1 = -38.640844 c2 = 5.894842 c3 = 59.876516 c4 = 26.654627 c5 = 10.637097 Pc = 220.85 ## bar Tc = 647.29 ## K t = (TK - Tc) / Tc PH2O = Pc * TK/Tc*(1 + c1*(-t)**1.9 + c2*t + c3*t**2 + c4*t**3 + c5*t**4) yCO2 = (P - PH2O) / P ## gas phase molar fraction of CO2 mNaCl = 2 * I/ ((1)**2 + (-1)**2) par_mu = np.array([28.9447706, -0.0354581768, -4770.67077, 0.0000102782768, 33.8126098, 0.0090403714, -0.00114934031, -0.307405726, -0.0907301486, 0.000932713393, 0]) par_lambda = np.array([-0.411370585, 0.000607632013, 97.5347708, 0, 0, 0, 0, -0.0237622469, 0.0170656236, 0, 0.0000141335834]) par_xi = np.array([0.000336389723, -0.000019829898, 0, 0, 0, 0, 0, 0.0021222083, -0.00524873303, 0, 0]) fTP = [1*np.ones([1, len(TK)]), TK, 1/TK, TK**2, 1/(630-TK), P, P*np.log(TK), P/TK, P/(630-TK), (P/(630-TK))**2, TK*np.log(P)] fTP = np.vstack(fTP) muCO2 = np.sum(par_mu.reshape(-1,1) * fTP, 0) # mu0/RT lambdaCO2Na = np.sum(par_lambda.reshape(-1,1) * fTP, 0) # lambda_CO2-Na Pitzer 2nd order int. param. xiCO2NaCl = np.sum(par_xi.reshape(-1,1) * fTP, 0) # zeta_CO2-Na-Cl Pitzer 3rd order int. param. # activity coef. aqueous co2 # Honoring limits for Duan and Sun - P-T-X range (0– 2000 bar, 0–260°C, 0–4.3 m NaCl) lngamco2 = np.zeros([len(I), len(TK)]); mco2 = np.zeros([len(I), len(TK)]) for j in range(len(I)): for i in range(len(TK)): if (TK[i] <= convert_temperature( 260, Out_Unit = 'K' ) ) and (I[j] <= 4.3) and (P[i] <= 2000): lngamco2[j, i] = 2*lambdaCO2Na[i]*mNaCl[j] + xiCO2NaCl[i]*mNaCl[j]**2 mco2[j, i] = P[i] * yCO2[i]/np.exp(muCO2[i] - np.log(fCO2[i]/P[i]) + \ 2*lambdaCO2Na[i]*mNaCl[j] + xiCO2NaCl[i]*mNaCl[j]**2) else: lngamco2[j, i] = 0 mco2[j, i] = 0 log10_co2_gamma = np.log10(np.exp(lngamco2)) # ((lngamco2)/np.log(10)) # return log10_co2_gamma, mco2
[docs]def co2fugacity(TK, P, poy = True): """ This function computes the fugacity and density of CO2 by Duan and Sun 2003 Also Calculate the Saturation Index SI and partial pressure of CO2(g) at any given T, P using the Duan equation of state. A Poynting correction factor is also applied. Parameters ---------- TK : float, vector Temperature [K] P : float, vector pressure [bar] poy : bool Option to apply a Poynting correction factor [True is the default] Returns ---------- fCO2 : float, vector co2 fugacity [bar] denCO2 : float, vector co2 density [g/cm3] SI : float, vector co2 saturation index pCO2 : float, vector co2 partial pressure [bar] Usage ---------- [fCO2, denCO2, SI, pCO2] = co2fugacity( TK, P) Examples -------- >>> fCO2, denCO2, SI, pCO2 = co2fugacity( np.array([400,420]), np.array([150, 200])) >>> fCO2 array([114.50387423, 150.87672517]) >>> denCO2 array([2.6741242 , 3.33134397]) >>> SI array([1.73291553, 1.84597591]) >>> pCO2 array([54.78127571, 71.0710152]) """ if np.ndim(TK) == 0 : TK = np.array(TK).ravel() else: TK = np.ravel(TK) if np.ndim(P) == 0 : P = np.array(P).ravel() else: P = np.ravel(P) Rgas = 0.0831446261815324 #bar*L/mol/K PcCO2 = 73.825 ## bar TcCO2 = 31.05 + 273.15 ## K VcCO2 = Rgas * TcCO2 / PcCO2 # L/mol Pr = P / PcCO2 # unitless Tr = TK / TcCO2 # unitless xmwc = 44.0098 # g/mol a1 = 0.0899288497 a2 = -0.494783127 a3 = 0.0477922245 a4 = 0.0103808883 a5 = -0.0282516861 a6 = 0.0949887563 a7 = 0.00052060088 a8 = -0.000293540971 a9 = -0.00177265112 a10 = -0.0000251101973 a11 = 0.0000893353441 a12 = 0.0000788998563 a13 = -0.0166727022 a14 = 1.398 a15 = 0.0296 # Equation A1 Z = np.zeros([len(TK), 1]).ravel() for k in range(len(TK)): zfun = lambda Z : -Z + (1 + (a1 + a2 / Tr[k]**2 + a3 / Tr[k]**3) / (Z * Tr[k]) * Pr[k] + (a4 + a5 / Tr[k]**2 + a6 / Tr[k]**3) / (Z * Tr[k])**2 * Pr[k]**2 + (a7 + a8 / Tr[k]**2 + a9 / Tr[k]**3) / (Z * Tr[k])**4 * Pr[k]**4 + (a10 + a11 / Tr[k]**2 + a12 / Tr[k] ** 3) / (Z * Tr[k])**5 * Pr[k]**5 + a13 / Tr[k]**3 / (Z * Tr[k])**2 * Pr[k]**2 * (a14 + a15 / (Z * Tr[k])**2 * Pr[k]**2)* np.exp(-a15 / (Z * Tr[k])**2 * Pr[k]**2)) # set initial guess from ideal gas law Videal = Rgas*TK[k]/P[k] Zguess = Videal/VcCO2*Pr[k]/Tr[k] Z[k] = fsolve(zfun, Zguess, xtol=1.0e-10) Vr = Z * Tr / Pr V = Vr * VcCO2 ## L / mol phi = np.exp(Z - 1 - np.log(Z) + (a1 + a2 / (Tr ** 2) + a3 / (Tr ** 3)) / Vr + \ (a4 + a5 / (Tr ** 2) + a6 / (Tr ** 3)) / (2 * Vr ** 2) + \ ((a7 + a8 / (Tr ** 2) + a9 / (Tr ** 3)) / (4 * Vr ** 4) + \ (a10 + a11 / (Tr ** 2) + a12 / (Tr ** 3)) / (5 * Vr ** 5) + \ a13 / (2 * Tr ** 3 * a15) * (a14 + 1 - (a14 + 1 + a15 / Vr ** 2) * \ np.exp(-a15 / Vr ** 2)))) # print(phi) #-----fugacity fCO2 = phi * P #bar #-----density denCO2 = xmwc*1e-3 / V # g/cm^3 Poy = 1 Patm = P / 1.01325 if (poy): R = 0.082057366080960 ## L atm /K/mol Vm = np.where(V < 0, 32.0e-3, V) ## L / mol Poy = np.exp(-(Patm - 1)*Vm/R/TK) pCO2 = phi*Patm*Poy # atm SI = np.log10(pCO2) pCO2 = pCO2 * 1.01325 # bar return fCO2, denCO2, SI, pCO2
[docs]def gamma_correlation(TC, P, method = None): """ This function calculates the CO2 activity correlation coefficients at given temperature T and pressure P Parameters ---------- TC : float, vector Temperature [°C] P : float, vector Pressure [bar] method : string specify the activity model [``'Duan_Sun'`` or ``'Drummond'``] Returns ---------- cco2 : float, vector co2 correlation coefficients Usage ---------- cco2 = gamma_correlation( TC, P) Examples -------- >>> log10_co2_gamma, mco2 = Henry_duan_sun( 500, 250, 0.5) >>> log10_co2_gamma, mco2 0.06202532, 1.6062557 References ---------- (1) Segal Edward Drummond, 1981, Boiling and Mixing of Hydrothermal Fluids: Chemical Effects on Mineral Precipitation, page 19 (2) Wolery, T. J., Lawrence Livermore National Laboratory, United States Dept. of Energy, 1992. EQ3/6: A software package for geochemical modeling of aqueous systems: package overview and installation guide (version 7.0) """ if np.ndim(TC) == 0: TC = np.array(TC).ravel() if np.ndim(P) == 0: P = np.array(P).ravel() TK = convert_temperature( TC, Out_Unit = 'K' ) # assign ionic strength from 0 to 3 N = 100 I = np.zeros([N, 1]) for i in range(N): I[i] = 3 * i /(N - 1) A = np.zeros([3, 3]) B = np.zeros([3, 1]) cco2 = np.zeros([4, len(TK)]) if method is not None: if method.lower() == 'duan_sun': ccoef = Henry_duan_sun(TK, P, I)[0] elif method.lower() == 'drummond': ccoef = drummondgamma(TK, I) else: ccoef = drummondgamma(TK, I) for i in range(len(TK)): for ii in range(3): for jj in range(3): A[ii, jj] = np.sum(I**((ii+1) + (jj+1))) B[ii] = np.sum(ccoef[:, i]*I.ravel()**(ii+1)) Coef = lu_solve(lu_factor(A), B) cco2[:, i] = np.concatenate([Coef.ravel(),np.array([0])]) return cco2
[docs]def Helgeson_activity(TC, P, I, Dielec_method = None, **rhoEDB): """ This function calculates the solute activity coefficient, solvent osmotic coefficient, and solvent activity at given temperature and pressure using equations 298, 190 and 106 in Helgeson, Kirkham and Flowers, 1981, A.J.S. p.1249-1516 Parameters: ---------- TC : Temperature [°C] P : pressure [bar] I : ionic strength Dielec_method : specify either 'FGL97' or 'JN91' or 'DEW' as the method to calculate dielectric constant (optional), if not specified default - 'JN91' rhoEDB : dictionary of water properties like density (rho), dielectric factor (E) and Debye–Hückel coefficients (optional) Returns: ---------- aw : solvent activity phi : solvent osmotic coefficient mean_act : solute activity coefficient Usage ------- [aw, phi, mean_act] = Helgeson_activity( TC, P, I) """ if np.ndim(TC) == 0: TC = np.array(TC).ravel() if np.ndim(P) == 0: P = np.array(P).ravel() if np.ndim(I) == 0 : I = np.array(I).ravel() if np.size(TC) < np.size(P): TC = TC*np.ones_like(P) if np.size(P) < np.size(TC): P = P*np.ones_like(TC) I = I.reshape(-1,1) Dielec_method = 'JN91' if Dielec_method is None else Dielec_method if rhoEDB.__len__() != 0: Ah = rhoEDB['Ah'].ravel() Bh = rhoEDB['Bh'].ravel() else: water = water_dielec(T = TC, P = P, Dielec_method = Dielec_method) Ah, Bh = water.Ah, water.Bh z = [-1.0, 1.0] mwH2O = 18.01528/1000 # kg/mol R = 1.9872041 # cal/mol/K mstar = 2*I # total solute in solution mtj = I # total concentration of ion j mchr= 2*I # total solute excluding neutral species loggamma = [0]*2; summt = [0]*2 rej = [1.810, 1.910] # Rej of Cl- and Na+ bijl = np.zeros(np.size(TC)); bi = np.zeros(np.size(TC)) # bil and bihat correlation for k in range(len(TC)): if TC[k] < 350: Psat = iapws95(T = TC[k], P = 'T').P else: Psat = [0] if P[k] == Psat: #%% P = Psat Region x = np.array([ 25., 50., 75., 100., 125., 150., 175., 200., 225., 250., 275., 300., 325.]) y = np.array([ 2.47, 2.15, 1.79, 1.39, 0.93, 0.41, -0.18, -0.85, -1.64, -2.57, -3.71, -5.21, -7.32]) fun = splrep(x, y) bhat = splev(TC[k], fun) y = np.array([-9.77, -5.59, -2.43, 0.23, 2.44, 4.51, 6.38, 8.11, 9.73, 11.29, 11.71, 14.15, 15.49]) fun = Rbf(x, y) bil = fun(TC[k]) elif P[k] != Psat and P[k] < 1000: isoline = P[k] x = np.array([ 25., 50., 75., 100., 125., 150., 175., 200., 225., 250., 275., 300., 325.]) y = np.array([ 2.47, 2.15, 1.79, 1.39, 0.93, 0.41, -0.18, -0.85, -1.64, -2.57, -3.71, -5.21, -7.32]) fun = splrep(x, y) xpred = np.linspace(0, 500, 100) ypred1 = splev(xpred, fun) ypred1 = np.where(ypred1 < -17.12, -17.12, ypred1) x = np.arange(25,525,25) y = np.array([ 2.58, 2.28, 1.95, 1.58, 1.18, 0.73, 0.25, 0.28, -0.86, -1.5 , -2.2 , -2.99, -3.88, -4.91, -6.11, -7.56, -9.31, -11.43, -14.01, -17.12]) fun = splrep(x, y) ypred2 = splev(xpred, fun) newy = (1000-isoline)/1000*ypred1 + (isoline/1000)*ypred2 fun = splrep(xpred, newy) bhat = splev(TC[k], fun) x = np.array([ 25., 50., 75., 100., 125., 150., 175., 200., 225., 250., 275., 300., 325.]) y = np.array([-9.77, -5.59, -2.43, 0.23, 2.44, 4.51, 6.38, 8.11, 9.73, 11.29, 11.71, 14.15, 15.49]) fun = Rbf(x, y) xpred = np.linspace(0, 500, 100) ypred1 = fun(xpred) ypred1 = np.where(ypred1 > 24.03, 24.03, ypred1) x = np.arange(25,525,25) y = np.array([ -9.19, -5.07, -2.63, -0.07, 1.71, 4.13, 5.91, 7.48, 9.27, 10.8 , 12.24, 13.68, 15.16, 16.48, 17.8 , 19.12, 20.33, 21.55, 22.86, 24.03]) fun = splrep(x, y) ypred2 = splev(xpred, fun) newy = (1000-isoline)/1000*ypred1 + (isoline/1000)*ypred2 fun = splrep(xpred, newy) bil = splev(TC[k], fun) elif P[k] >= 1000 and P[k] <= 2000: isoline = P[k] - 1000 x = np.arange(25,525,25) y = np.array([ 2.58, 2.28, 1.95, 1.58, 1.18, 0.73, 0.25, 0.28, -0.86, -1.5 , -2.2 , -2.99, -3.88, -4.91, -6.11, -7.56, -9.31, -11.43, -14.01, -17.12]) fun = splrep(x, y) xpred = np.linspace(0, 500, 100) ypred1 = splev(xpred, fun) y = np.array([ 2.66, 2.37, 2.06, 1.72, 1.35, 0.95, 0.51, 0.04, -0.46, -1. , -1.57, -2.19, -2.85, -3.58, -4.38, -5.26, -6.24, -7.34, -8.56, -9.88]) fun = splrep(x, y) ypred2 = splev(xpred, fun) newy = (1000-isoline)/1000*ypred1 + (isoline/1000)*ypred2 fun = splrep(xpred, newy) bhat = splev(TC[k], fun) y = np.array([ -9.19, -5.07, -2.63, -0.07, 1.71, 4.13, 5.91, 7.48, 9.27, 10.80, 12.24, 13.68, 15.16, 16.48, 17.8 , 19.12, 20.33, 21.55, 22.86, 24.03]) fun = splrep(x, y) xpred = np.linspace(0, 500, 100) ypred1 = splev(xpred, fun) y = np.array([ -9.12, -5.52, -2.65, -0.1 , 1.83, 4.12, 5.88, 7.66, 9.2 , 10.73, 12.19, 13.66, 15.1 , 16.42, 17.79, 19.11, 20.35, 21.61, 22.84, 23.95]) fun = splrep(x, y) ypred2 = splev(xpred, fun) newy = (1000-isoline)/1000*ypred1 + (isoline/1000)*ypred2 fun = splrep(xpred, newy) bil = splev(TC[k], fun) elif P[k] == 3000: #%% P = 3000 Region x = np.arange(25,525,25) y = np.array([ 2.72, 2.45, 2.15, 1.83, 1.48, 1.11, 0.72, 0.29, -0.15, -0.63, -1.13, -1.67, -2.24, -2.86, -3.52, -4.24, -5.03, -5.87, -6.78, -7.73]) fun = splrep(x, y) bhat = splev(TC[k], fun) y = np.array([ -9.51, -5.6 , -2.46, 0.13, 2.14, 4.28, 6.16, 8.18, 9.59, 11.15, 12.62, 13.99, 15.41, 16.76, 18.1 , 19.4 , 20.66, 21.85, 23.09, 24.32]) fun = splrep(x, y) bil = splev(TC[k], fun) elif P[k] == 4000: #%% P = 4000 Region x = np.arange(25,525,25) y = np.array([ 2.77, 2.51, 2.22, 1.91, 1.58, 1.23, 0.86, 0.47, 0.05, -0.38, -0.84, -1.33, -1.85, -2.4 , -2.99, -3.63, -4.32, -5.06, -5.84, -6.63]) fun = splrep(x, y) bhat = splev(TC, fun) y = np.array([-10.44, -5.68, -2.15, 0.71, 3.05, 5.13, 7.05, 9.05, 10.4 , 11.92, 13.39, 14.83, 16.23, 17.52, 18.81, 20.13, 21.38, 22.56, 23.74, 24.94]) fun = splrep(x, y) bil = splev(TC[k], fun) elif P[k] >= 5000: #%% P = 5000 Region and beyond x = np.arange(25,525,25) y = np.array([ 2.82, 2.56, 2.28, 1.99, 1.67, 1.33, 0.98, 0.6 , 0.21, -0.2 , -0.63, -1.09, -1.58, -2.1 , -2.65, -3.25, -3.89, -4.57, -5.28, -6. ]) fun = splrep(x, y) bhat = splev(TC[k], fun) y = np.array([-11.76, -5.87, -1.68, 1.39, 4.21, 6.26, 8.13, 10.3, 11.58, 13.15, 14.54, 16.01, 17.36, 18.65, 19.83, 21.23, 22.38, 23.64, 24.78, 25.95]) fun = splrep(x, y) bil = splev(TC[k], fun) # bihat NaCl (b NaCl = bhat/(2.303RT)) from Table 29 (bi here) in kg/mol * 1e+3 # b Na+Cl- from Table 30 (bil here) in kg/mol * 1e+2 bihat = bhat*1e-3 # [kg/mol] bi[k] = bihat/(np.log(10)*R*convert_temperature( TC[k], Out_Unit = 'K' )) #[kg/cal] bijl[k] = bil*1e-2 # activity and osmotic coefficients # for j in range(len(I)): for i in range(2): zabsi = np.abs(z[i]) if i < 1: rex = 1.81 else: rex = 1.91 azero = 2*(rej[i] + zabsi*rex)/(zabsi + 1) omega = 1.66027e5*z[i]**2/rej[i] lambdaa = 1 + Bh*azero*I**0.5 loggamma[i] = - Ah * z[i]**2 * I**0.5 / lambdaa - np.log10(1 + mwH2O * mstar) + \ (omega*bi + bijl - 0.19*(np.abs(z[i]) - 1))*I summt[i] = mtj * ( (Ah*z[i]**2/((azero*Bh)**3*I)) *\ (lambdaa - 1/lambdaa - 2*np.log(lambdaa)) + \ (-np.log10(1 + mwH2O * mstar)/(mwH2O*mstar)) - \ 0.5*(omega*bi*I + (bijl - 0.19 *(np.abs(z[i]) - 1)) * \ mchr*0.5) ) mean_act = 10**((loggamma[0]+loggamma[1])/2) phi = -np.log(10)*(summt[0]+summt[1])/mstar phi = np.where(np.isnan(phi), 0, phi) aw = np.exp(-phi*mstar*mwH2O) return aw, phi, mean_act
[docs]def aw_correlation(TC, P, Dielec_method = None, **rhoEDB): """ Calculates the water activity correlation coefficients at given temperature and pressure Parameters: ---------- TC : Temperature [°C] P : pressure [bar] Dielec_method : specify either 'FGL97' or 'JN91' or 'DEW' as the method to calculate dielectric constant (optional), if not specified default - 'JN91' rhoEDB : dictionary of water properties like density (rho), dielectric factor (E) and Debye–Hückel coefficients (optional) Returns: ---------- ch20 : water activity correlation coefficients Usage ------- [ch20] = aw_correlation( TC, P) """ Dielec_method = 'JN91' if Dielec_method is None else Dielec_method if rhoEDB.__len__() != 0: rho = rhoEDB['rho'].ravel() E = rhoEDB['E'].ravel() Ah = rhoEDB['Ah'].ravel() Bh = rhoEDB['Bh'].ravel() else: if Dielec_method.upper() == 'DEW': rho = ZhangDuan(T = TC, P = P).rho else: rho = iapws95(T = TC, P = P).rho water = water_dielec(T = TC, P = P, Dielec_method = Dielec_method) E, Ah, Bh = water.E, water.Ah, water.Bh rhoEDB = {'rho': rho, 'E': E, 'Ah': Ah, 'Bh': Bh} if np.ndim(TC) == 0 | np.ndim(P) == 0: TC = np.array(TC).ravel() P = np.array(P).ravel() # assign ionic strength from 0 to 3 N = 100 Is = np.linspace(0, 6, N).reshape(-1,1) mwH2O = 18.01528/1000 # kg/mol # %% Water Activity aw = Helgeson_activity(TC, P, Is, Dielec_method = Dielec_method, **rhoEDB)[0] ch20 = np.zeros([4, len(TC)]) x0 = [1.454, 0.02236, 9.380e-3, -5.362e-4] for i in range(len(TC)): if np.sum(np.isnan(aw[:, i])) == 0: Is_0 = Is [1:] # avoid zero values Ahi = Ah[i] ch20func = lambda Is_0, *x: (-2*Is_0*mwH2O * \ (1 - (np.log(10)*Ahi/(x[0]**3*Is_0)) * \ ((1 + x[0]*np.sqrt(Is_0)) - 2*np.log(1 + x[0]*np.sqrt(Is_0)) - \ ( 1/(1 + x[0]*np.sqrt(Is_0)) ) ) + \ (x[1]*Is_0/2) + (2/3*x[2]*Is_0**2) + (3/4*x[3]*Is_0**3) )) ch20[:, i], pcov = curve_fit(ch20func, Is_0.ravel(), np.log(aw[1:, i]).ravel(), p0=x0, maxfev = 1000000) else: ch20[:, i] = [500]*4 return ch20
#%%------------------------------------------------------------------------------------
[docs]class calcRxnlogK(): """ This class implemetation calculates logK values for any reaction with the option for extrapolation where rho < 350kg/m3 \n Parameters ---------- T : float, vector Temperature [°C] P : float, vector Pressure [bar] Specie : string specify the species for logK calculation, either the Product species of any reaction or solid solutions or clay like 'AnAb' or 'AbOr' or 'FoFa' or 'EnFe' or 'DiHedEnFe' or 'clay' \n Specie_class : string, optional specify the class of species like 'aqueous', 'minerals', 'liquids', or 'gases' elem : list list containing nine parameters with clay names and elements compositions with the following format ['Montmorillonite_Lc_MgK', 'Si', 'Al', 'FeIII', 'FeII', 'Mg', 'K', 'Na', 'Ca', 'Li'] \n dbaccessdic : dict direct-acess database dictionary rhoEGextrap : dict dictionary of water properties like density (rho), dielectric factor (E) and Gibbs Energy for density region 350-550kg/m3 group : string specify the structural layering of the phyllosilicate, for layers composed of ``1 tetrahedral + 1 octahedral sheet (1:1 layer)`` - specify ``'7A'``, ``2 tetrahedral + 1 octahedral sheet (2:1 layer)`` - specify ``'10A'``, or the latter with ``a brucitic sheet in the interlayer (2:1:1 layer)`` - specify ``'14A'`` (optional), if not specified, default is '10A' for smectites, micas, et cetera \n ClayMintype : string specify either 'Smectite' or 'Chlorite' or 'Mica' as the clay type, if not specified default - 'Smectites' X : float volume fractions of any (Anorthite, Albite, Forsterite, Enstatite) or mole fraction of Mg cpx_Ca : float number of moles of Ca in formula unit (=1 for Di, Hed), must be greater than zero sourcedic : dict source database reactions dictionary specielist : list of list, optional source database species grouped into categories [element, basis, redox, aqueous, minerals, gases, oxides] Dielec_method : string specify either 'FGL97' or 'JN91' or 'DEW' as the method to calculate dielectric constant, default is 'JN91' \n heatcap_method : string specify either 'SUPCRT' or 'Berman88' or 'HP11' or 'HF76' as the method to calculate thermodynamic properties of any mineral or gas, default is 'SUPCRT' \n ThermoInUnit : string specify either 'cal' or 'KJ' as the input units for species properties (optional), particularly used to covert KJ data to cal by supcrtaq function if not specified default - 'cal' Al_Si : string specify either 'pygcc' or 'Arnórsson_Stefánsson' as the input to express Al and Si species in solid solution (optional), 'Arnórsson_Stefánsson' expresses them as 'Al(OH)4-' and 'H4SiO4(aq)', respectively while pygcc uses 'Al3+' and 'SiO2(aq)', respectively if not specified default - 'pygcc' sourceformat : string specify the source database format, either 'GWB' or 'EQ36' densityextrap : float, vector specify the extrapolation option for density-logK, 'Yes'/True or 'No'/False Returns ------- log_K : float, vector logarithmic K value(s) \n dGrxn : float, vector Total reaction Gibbs energy [cal/mol] \n Usage ---------- The general usage of water_dielec is as follows: \n (1) For water dielectric properties at any Temperature and Pressure: \n calclogK = calcRxnlogK(T = T, P = P, Dielec_method = 'JN91', **kwargs), \n where T is temperature in celsius and P is pressure in bar (2) For water dielectric properties at any Temperature and density : \n calclogK = calcRxnlogK(T = T, rho = rho, Dielec_method = 'JN91', **kwargs), \n where T is temperature in celsius and rho is density in kg/m³ (3) For water dielectric properties at any Temperature and Pressure on steam saturation curve: \n calclogK = calcRxnlogK(T = T, P = 'T', Dielec_method = 'JN91', **kwargs), \n where T is temperature in celsius, followed with a quoted character 'T' to reflect steam saturation pressure \n calclogK = calcRxnlogK(P = P, T = 'P', Dielec_method = 'JN91', **kwargs), \n where P is pressure in bar, followed with a quoted character 'P' to reflect steam saturation temperature Examples -------- >>> ps = db_reader(sourcedb = './default_db/thermo.com.dat', sourceformat = 'gwb', dbaccess = './default_db/speq21.dat') >>> calclogK = calcRxnlogK(T = 100, P = 50, Specie = 'H2S(aq)', Specie_class = 'aqueous', dbaccessdic = ps.dbaccessdic, sourcedic = ps.sourcedic, specielist = ps.specielist) >>> calclogK.logK, calclogK.dGrxn -6.4713, 11049.3168 """ kwargs = {"T": None, "Specie": None, "Specie_class": None, "ThermoInUnit": 'cal', "P": None, "group": None, "X": None, "cpx_Ca": None, "elem": None, 'rhoEGextrap': None, "sourcedic": None, "specielist": None, 'dbaccessdic': None, "Dielec_method": None, "sourceformat": None, 'heatcap_method': None, "densityextrap": None, 'rhoEG': None, 'ClayMintype': 'Smectite', "Al_Si": 'pygcc'} def __init__(self, **kwargs): self.kwargs = calcRxnlogK.kwargs.copy() self.__calc__(**kwargs)
[docs] def __calc__(self, **kwargs): self.kwargs.update(kwargs) """initialization """ self.TC = self.kwargs["T"]; self.P = self.kwargs["P"] self.sourceformat = self.kwargs['sourceformat']; self.specielist = self.kwargs['specielist']; self.sourcedic = self.kwargs['sourcedic']; self.Dielec_method = self.kwargs['Dielec_method']; self.rhoEG = self.kwargs['rhoEG']; self.ThermoInUnit = self.kwargs['ThermoInUnit'] self.rhoEGextrap = self.kwargs['rhoEGextrap']; self.Specie = self.kwargs['Specie'] self.Specie_class = self.kwargs['Specie_class']; self.elem = self.kwargs['elem'] self.cpx_Ca = self.kwargs['cpx_Ca']; self.group = self.kwargs['group'] self.X = self.kwargs['X']; self.heatcap_method = self.kwargs['heatcap_method'] self.ClayMintype = self.kwargs['ClayMintype']; self.Al_Si = self.kwargs["Al_Si"] self.densityextrap = 'No' if (self.kwargs['densityextrap'] is None or self.kwargs['densityextrap'] is False) else 'Yes' if self.kwargs['densityextrap'] is True else self.kwargs['densityextrap'] self.Dielec_method = 'JN91' if self.Dielec_method is None else self.Dielec_method self.heatcap_method = 'SUPCRT' if self.heatcap_method is None else self.heatcap_method self.sourceformat = 'GWB' if self.sourceformat is None else self.sourceformat if self.kwargs['dbaccessdic'] == None: self.dbaccess_dir = './default_db/speq21.dat' self.dbaccess_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), self.dbaccess_dir) self.dbaccessdic = db_reader(dbaccess = self.dbaccess_dir).dbaccessdic else: self.dbaccessdic = self.kwargs['dbaccessdic'] if (type(self.P) == str) or (type(self.TC) == str): if self.P == 'T': self.P = iapws95(T = self.TC).P self.P[np.isnan(self.P) | (self.P < 1)] = 1.0133 elif self.TC == 'P': self.TC = iapws95(P = self.P).TC if np.ndim(self.TC) == 0 : self.TC = np.array(self.TC).ravel() elif np.size(self.TC) == 2: self.TC = np.array([roundup_tenth(j) if j != 0 else 0.01 for j in np.linspace(self.TC[0], self.TC[-1], 8)]) if np.size(self.P) <= 2: self.P = np.ravel(self.P) self.P = self.P[0]*np.ones(np.size(self.TC)) if self.rhoEG is None: if self.Dielec_method.upper() == 'DEW': water = ZhangDuan(T = self.TC, P = self.P) else: water = iapws95(T = self.TC, P = self.P) rho, dGH2O = water.rho, water.G E = water_dielec(T = self.TC, P = self.P, Dielec_method = self.Dielec_method).E self.rhoEG = {'rho': rho, 'E': E, 'dGH2O': dGH2O} if self.densityextrap.lower() == 'yes': subBornptrs = self.rhoEG['rho'] < 350 self.nonsubBornptrs = self.rhoEG['rho'] >= 350 rhoEG_nonsubBornptrs = {'rho': self.rhoEG['rho'][self.nonsubBornptrs], 'E': self.rhoEG['E'][self.nonsubBornptrs], 'dGH2O': self.rhoEG['dGH2O'][self.nonsubBornptrs]} if self.rhoEGextrap is None: # Calculate the rho E G for density extrapolation method here so we have it below self.rhoEGextrap = {} if any(subBornptrs): for i, j in enumerate(zip(self.TC[subBornptrs], self.P[subBornptrs])): rhoextrap = np.linspace(350, 550, 3) Pextrap = iapws95(T = j[0], rho = rhoextrap).P if self.Dielec_method.upper() != 'DEW' else ZhangDuan(T = j[0], rho = rhoextrap).P Textrap = j[0]*np.ones(np.size(Pextrap)) dGH2O = iapws95(T = Textrap, P = Pextrap).G if self.Dielec_method.upper() != 'DEW' else ZhangDuan(T = Textrap, P = Pextrap).G E = water_dielec(T = Textrap, P = Pextrap, Dielec_method = self.Dielec_method).E rhoextrap = np.around(rhoextrap, 3) self.rhoEGextrap['%d_%d' % (j[0], j[1])]= {'rho': rhoextrap,'E': E, 'dGH2O': dGH2O, 'Textrap': Textrap, 'Pextrap': Pextrap} if self.densityextrap.lower() == 'no': if self.Specie.lower().startswith(('plagio', 'oliv', 'pyroxe', 'alk-')) or self.Specie.lower() in ['cpx', 'clay']: self.logK, self.Rxn = self.AllRxnslogK( self.TC, self.P, self.rhoEG) else: self.logK, self.dGrxn, dGP, dGRs = self.AllRxnslogK( self.TC, self.P, self.rhoEG) self.nonsubBornptrs = [False]*len(self.TC) # Required to shut off density extrapolation prompt elif self.densityextrap.lower() == 'yes': self.logK = np.nan*np.zeros(len(self.TC)) self.dGrxn = np.nan*np.zeros(len(self.TC)) if any(subBornptrs): self.logK[subBornptrs] = self.densitylogKextrap( self.TC[subBornptrs], self.P[subBornptrs], self.rhoEGextrap) if any(self.nonsubBornptrs): if self.Specie.lower().startswith(('plagio', 'oliv', 'pyroxe', 'alk-')) or self.Specie.lower() in ['cpx', 'clay']: self.logK[self.nonsubBornptrs], self.Rxn = self.AllRxnslogK( self.TC[self.nonsubBornptrs], self.P[self.nonsubBornptrs], rhoEG_nonsubBornptrs) else: self.logK[self.nonsubBornptrs], self.dGrxn[self.nonsubBornptrs], dGP, dGRs = self.AllRxnslogK( self.TC[self.nonsubBornptrs], self.P[self.nonsubBornptrs], rhoEG_nonsubBornptrs)
[docs] def AllRxnslogK( self, TC, P, rhoEG): """ This function calculates logK values of all reactions including solid-solution and clay minerals \n Parameters ---------- TC : temperature [°C] \n P : pressure [bar] \n rhoEG : dictionary of water properties like density (rho), dielectric factor (E) and Gibbs Energy (optional) \n Returns ------- logK : logarithmic K value(s) \n dGrxn : Total reaction Gibbs energy [cal/mol] \n dGP : Product specie Gibbs energy [cal/mol] \n dGRs : Reactant species Gibbs energy [cal/mol] \n Rxn : dict The calculated dictionary of reaction thermodynamic properties has the following properties: * type: solid-solution mineral type, [-] * name: solid-solution name, [K] * formula: solid-solution mineral formula, [-] * MW: Molecular weight, [g/mol] * min: solid-solution mineral properties, ['formula', 'source date', dG[cal/ml], dH[cal/mol], S[cal/mol-K], V[cm3/mol], a[cal/mol-K], b[10^3 cal/mol/K^2], c[10^-5 cal/mol/K]] * spec: list of species, [-] * coeff: list of corresponding coefficients of species above, [-] * nSpec: Total number of species, [-] * V: Molar volume, [cm3/mol] * source: Source of thermo data, [kJ/kg·K] * elements: list of elements and their total numbers, [-] Usage ------- The general usage of calcRxnlogK without the optional argument is as follows: \n (1) Not on steam saturation curve: \n [logK, dGrxn, dGP, dGRs] = AllRxnslogK(TC, P, Prod, dbaccessdic, sourcedic, specielist), \n where TC is temperature in celsius and P is pressure in bar; (2) On steam saturation curve: \n [logK, dGrxn, dGP, dGRs] = AllRxnslogK(TC, 'T', Prod, dbaccessdic, sourcedic, specielist), \n where TC is temperature in celsius, followed with a quoted char 'T' \n [logK, dGrxn, dGP, dGRs] = AllRxnslogK(P, 'P', Prod, dbaccessdic, sourcedic, specielist), \n where P is pressure in bar, followed with a quoted char 'P'. (3) Meanwhile, usage with any specific dielectric constant method ('FGL97') for condition not on steam saturation curve is as follows. Default method is 'JN91' \n [logK, dGrxn, dGP, dGRs] = AllRxnslogK( TC, P, Prod, dbaccessdic, sourcedic, specielist, Dielec_method = 'FGL97') """ if self.Specie.lower().startswith(('plagio', 'oliv', 'pyroxe', 'alk-')): ss = solidsolution_thermo(X = self.X, T = TC, P = P, Dielec_method = self.Dielec_method, dbaccessdic = self.dbaccessdic, solidsolution_type = self.Specie, ThermoInUnit = self.ThermoInUnit, rhoEG = rhoEG, Al_Si = self.Al_Si) return ss.logK, ss.Rxn elif self.Specie.lower() == 'cpx': ss = solidsolution_thermo(cpx_Ca = self.cpx_Ca, X = self.X, T = TC, P = P, dbaccessdic = self.dbaccessdic, Dielec_method = self.Dielec_method, solidsolution_type = 'cpx', ThermoInUnit = self.ThermoInUnit, rhoEG = rhoEG, Al_Si = self.Al_Si) return ss.logK, ss.Rxn elif self.Specie.lower() == 'clay': logK, Rxn = calclogKclays(TC, P, *self.elem, dbaccessdic = self.dbaccessdic, group = self.group, Dielec_method = self.Dielec_method, ThermoInUnit = self.ThermoInUnit, ClayMintype = self.ClayMintype, **rhoEG) return logK, Rxn else: logK, dGrxn, dGP, dGRs = self.RxnlogK( TC, P, self.Specie, rhoEG) return logK, dGrxn, dGP, dGRs
[docs] def RxnlogK( self, TC, P, Prod, rhoEG): """ This function calculates logK values of any reaction \n Parameters ---------- TC : temperature [°C] \n P : pressure [bar] \n Prod : Product species of the reaction \n dbaccessdic : direct-acess database dictionary \n sourcedic : source database reactions dictionary \n specielist : source database species grouped into [element, basis, redox, aqueous, minerals, gases, oxides] \n Dielec_method : specify either 'FGL97' or 'JN91' or 'DEW' as the method to calculate dielectric constant, default is 'JN91' \n sourceformat: source database format, either 'GWB' or 'EQ36', default is 'GWB' heatcap_method : specify either 'SUPCRT' or 'Berman88' or 'HP11' or 'HF76' as the method to calculate thermodynamic properties of any mineral, default is 'SUPCRT' \n rhoEG : dictionary of water properties like density (rho), dielectric factor (E) and Gibbs Energy (optional) \n Returns ------- logK : logarithmic K value(s) \n dGrxn : Total reaction Gibbs energy [cal/mol] \n dGP : Product specie Gibbs energy [cal/mol] \n dGRs : Reactant species Gibbs energy [cal/mol] \n Usage ------- The general usage of calcRxnlogK without the optional argument is as follows: \n (1) Not on steam saturation curve: \n [logK, dGrxn, dGP, dGRs] = RxnlogK(TC, P, Prod, dbaccessdic, sourcedic, specielist), \n where TC is temperature in celsius and P is pressure in bar; (2) On steam saturation curve: \n [logK, dGrxn, dGP, dGRs] = RxnlogK(TC, 'T', Prod, dbaccessdic, sourcedic, specielist), \n where TC is temperature in celsius, followed with a quoted char 'T' \n [logK, dGrxn, dGP, dGRs] = RxnlogK(P, 'P', Prod, dbaccessdic, sourcedic, specielist), \n where P is pressure in bar, followed with a quoted char 'P'. (3) Meanwhile, usage with any specific dielectric constant method ('FGL97') for condition not on steam saturation curve is as follows. Default method is 'JN91' \n [logK, dGrxn, dGP, dGRs] = RxnlogK( TC, P, Prod, dbaccessdic, sourcedic, specielist, Dielec_method = 'FGL97') """ R = 1.9872041 # cal/mol/K dGH2O = rhoEG['dGH2O'].ravel() TK = convert_temperature( TC, Out_Unit = 'K' ) if self.sourceformat.upper() == 'EQ36': rxnspecies = [j for k,j in enumerate(self.sourcedic[Prod]) if k not in [2,3]] elif self.sourceformat.upper() == 'GWB': rxnspecies = self.sourcedic[Prod] method = 'SUPCRT' if (self.heatcap_method != 'HP11' and (Prod.endswith(('(g)', ',g')) or self.Specie_class == 'gases')) else self.heatcap_method # print(self.heatcap_method, Prod, method) if Prod == 'e-' or Prod == 'eh': dGP = 0 elif Prod == 'H2O': dGP = dGH2O elif Prod in ['Hydroxyapatite', 'Fluorapatite', 'Ankerite', 'Acmite', 'Molybdenite', 'Molybdite'] or Prod.startswith('ss_'): if self.heatcap_method != 'HP11': #and self.dbaccessdic[Prod][1].split()[0] in ['H&P2011', 'R&H95'] dGP = heatcap( T = TC, P = P, Species_ppt = self.dbaccessdic[Prod], method = 'HF76').dG else: dGP = heatcap( T = TC, P = P, Species_ppt = self.dbaccessdic[Prod], Species = Prod, method = 'HP11').dG elif self.specielist is None: if (Prod.endswith(('(aq)','+','-')) or Prod[-1].isdigit()): dGP = supcrtaq(TC, P, self.dbaccessdic[Prod.replace('(CH3COO)', '(Ac)').replace('CH3COO', '(Ac)')], Dielec_method = self.Dielec_method, ThermoInUnit = self.ThermoInUnit, **rhoEG) else: dGP = heatcap( T = TC, P = P, Species_ppt = self.dbaccessdic[Prod], Species = Prod, method = method).dG elif (Prod in self.specielist[4] + self.specielist[5] + self.specielist[6]) or (self.Specie_class in ['minerals', 'liquids', 'gases']): # dGP = heatcap( T = TC, P = P, Species_ppt = self.dbaccessdic[Prod], Species = Prod, method = method).dG else: dGP = supcrtaq(TC, P, self.dbaccessdic[Prod.replace('(CH3COO)', '(Ac)').replace('CH3COO', '(Ac)')], Dielec_method = self.Dielec_method, ThermoInUnit = self.ThermoInUnit, **rhoEG) total_reactants = int(len(rxnspecies[2:])/2) dGRs = 0 for i in range(total_reactants): R_coeff = float(rxnspecies[2 + 2*i]) R_specie = rxnspecies[4 + 2*i - 1] method = 'SUPCRT' if (self.heatcap_method != 'HP11' and R_specie.endswith(('(g)', ',g')) ) else self.heatcap_method # print(self.heatcap_method, R_specie, method) if R_specie == 'e-' or R_specie == 'eh': dGR = 0 elif R_specie == 'H2O': dGR = dGH2O elif R_specie in ['Hydroxyapatite', 'Fluorapatite', 'Ankerite', 'Acmite', 'Molybdenite', 'Molybdite'] or R_specie.startswith('ss_'): if self.heatcap_method != 'HP11': dGR = heatcap( T = TC, P = P, Species_ppt = self.dbaccessdic[R_specie], method = 'HF76').delG else: dGP = heatcap( T = TC, P = P, Species_ppt = self.dbaccessdic[R_specie], Species = R_specie, method = 'HP11').dG elif self.specielist is None: if (R_specie.endswith(('(aq)','+','-')) or R_specie[-1].isdigit()): dGR = supcrtaq(TC, P, self.dbaccessdic[R_specie.replace('(CH3COO)', '(Ac)').replace('CH3COO', '(Ac)')], Dielec_method = self.Dielec_method, ThermoInUnit = self.ThermoInUnit, **rhoEG) else: dGR = heatcap( T = TC, P = P, Species_ppt = self.dbaccessdic[R_specie], Species = R_specie, method = method).dG elif (R_specie in self.specielist[4] + self.specielist[5] + self.specielist[6]) or R_specie.endswith(('(g)', ',g')): dGR = heatcap( T = TC, P = P, Species_ppt = self.dbaccessdic[R_specie], Species = R_specie, method = method).dG else: dGR = supcrtaq(TC, P, self.dbaccessdic[R_specie.replace('(CH3COO)', '(Ac)').replace('CH3COO', '(Ac)')], Dielec_method = self.Dielec_method, ThermoInUnit = self.ThermoInUnit, **rhoEG) dGRs = dGRs + R_coeff*dGR dGrxn = - dGP + dGRs logK = (-dGrxn/R/(TK)/np.log(10)) return logK, dGrxn, dGP, dGRs
[docs] def densitylogKextrap(self, TC, P, rhoEGextrap ): """ This function calculates logK values extrapolation for conditions where rho < 350kg/m3 \n Parameters ---------- TC : temperature [°C] \n P : pressure [bar] \n rhoEGextrap : dictionary of water properties like density (rho), dielectric factor (E) and Gibbs Energy for density region 350-550kg/m3 \n Returns ------- logK : extrapolated logarithmic K value(s) \n Usage ------- [logK] = densitylogKextrap(TC, P, 'H2S(aq)', dbaccessdic = dbaccessdic, sourcedic = sourcedic, specielist = specielist), \n """ length = len(TC) logK = np.nan*np.ones(np.size(TC)) for i in range(length): #need a for loop because calculation is T-specific rhotarget = iapws95(T = TC[i], P = P[i]).rho rhotarget = rhotarget/1000 # kg/m3 => g/cm^3 Pextrap = rhoEGextrap['%d_%d' % (TC[i], P[i])]['Pextrap'] Textrap = rhoEGextrap['%d_%d' % (TC[i], P[i])]['Textrap'] rhoextrap = rhoEGextrap['%d_%d' % (TC[i], P[i])]['rho'] rhoextrap = rhoextrap/1000 # kg/m3 => g/cm^3 rhoextrap = np.around(rhoextrap, 3) rhoEGextrap_only = rhoEGextrap['%d_%d' % (TC[i], P[i])] logrho = np.log10(rhoextrap) logKextrap = self.AllRxnslogK( Textrap, Pextrap, rhoEGextrap_only)[0] p = np.polyfit(logrho, logKextrap, 1) logK[i] = np.polyval(p, np.log10(rhotarget)) return logK
[docs]def outputfmt(fid, logK, Rxn, *T, dataset = None, logK_form = None): """ This function writes logK and Rxn data to any file using GWB, EQ36, Pflotran and ToughReact format Parameters ---------- fid : string file ID logK : float, vector logarithmic K value(s) Rxn : dict dictionary of reaction thermodynamic properties T : float, vector Temperature value(s), optional, required when 'polycoeffs' is specified for logK_form dataset : string specify the dataset format, either 'GWB', 'EQ36', 'Pflotran' or 'ToughReact' logK_form : string specify the format of logK either as a set of eight values one for each of the dataset’s principal temperatures, or blocks of polynomial coefficients, [values, polycoeffs], default is 'a set of eight values' (optional) Returns ------- Output data to the file with filename described in fid with any format mentioned above. Examples -------- >>> fid = open('./logK_details.txt', 'w') >>> logKRxn = calcRxnlogK(T = 100, P = 'T', X = 0.634, Specie = 'Plagioclase', densityextrap = True) >>> # output in EQ36 format >>> outputfmt(fid, logKRxn.logK, logKRxn.Rxn, dataset = 'EQ36') >>> fid.close() """ if len(T) == 0: TK = 25*np.ones(np.size(logK)) else: TK = np.asarray(T).ravel() logK_form = 'values' if logK_form is None else logK_form.lower() # %% Open the text file. if dataset.lower() == 'gwb': fid.writelines("%s " % Rxn['name']) fid.writelines( "%s= " % list(Rxn.keys())[0]) if Rxn['type'].find('plag') == 0: fid.writelines( "plagioclase\n") else: fid.writelines( "%s\n" % Rxn['type']) fid.writelines( " %s= " % list(Rxn.keys())[2]) fid.writelines( "%s\n" % Rxn['formula']) fid.writelines( " mole vol.= %1.3f cc" % Rxn['V']) fid.writelines( " mole wt.= %1.4f g\n" % Rxn['MW']) fid.writelines( " %s species in reaction\n" % Rxn['nSpec']) for i in range(len(Rxn['spec'])): i = i + 1 fid.writelines( "%9.4f " % Rxn['coeff'][i-1]) fid.writelines( "%-9s " % Rxn['spec'][i-1]) if (i % 3 == 0) | (i % 6 == 0) | (i == len(Rxn['spec'])): fid.writelines( "\n") if logK_form.lower() == 'polycoeffs': Tr = 298.15 logKfunc = lambda TK, *x: x[0] + x[1]*(TK - Tr) + x[2]*(TK**2 - Tr**2) + x[3]*((1/TK) - (1/Tr)) + \ x[4]*((1/TK**2) - (1/Tr**2)) + x[5]*np.log(TK/Tr) x0 = [-31.9605, 20.6576, 3.73497e-2, -9.01862, 6.0111, 2.5] logKcorr = curve_fit(logKfunc, TK[logK!=500].ravel(), logK[logK!=500].ravel(), p0 = x0, maxfev = 1000000)[0] fid.writelines(' a= %15.9f ' % logKcorr[0] + 'b= %15.9f ' % logKcorr[1] + \ 'c= %15.6e\n' % logKcorr[2]) fid.writelines(' d= %15.6f ' % logKcorr[3] + 'e= %15.5f ' % logKcorr[4] + \ 'f= %15.8f \n' % logKcorr[5]) fid.writelines(' TminK= %-15.2f ' % np.min(TK) + 'TmaxK= %-7.2f\n' % np.max(TK)) else: for i in range(len(logK)): i = i + 1 if (i == 1) | (i == 5) | (i == 9) | (i == 13) | (i == 17): fid.writelines(" %9.4f" % logK[i-1]) else: fid.writelines(" %9.4f" % logK[i-1]) if (i % 4 == 0) | (i == 8) | (i == len(logK)): fid.writelines( "\n") fid.writelines( "* gflag = 1 [reported delG0f used]\n" ) if Rxn['type'].find('plag') == 0: fid.writelines( "* extrapolation algorithm: supcrt92 [92joh/oel]\n" ) else: fid.writelines( "* extrapolation algorithm: supcrt92/water95\n" ) if Rxn['type'].find('serp') == 0: fid.writelines( "* reference-state data source = Blanc et al 2015\n" ) elif 'source' in Rxn: fid.writelines( "* reference-state data source = %s\n" % Rxn['source']) else: fid.writelines( "* reference-state data source = supcrt92?\n" ) fid.write( "* delG0f = %8.3f kcal/mol\n" % (Rxn['min'][2]/1000) ) if Rxn['type'] == 'Smectites': fid.writelines( "* delH0f = %8.3f kcal/mol\n" % (Rxn['min'][3]/1000) ) else: fid.writelines( "* delH0f = NaN kcal/mol\n") fid.writelines( "* S0PrTr = %8.3f cal/mol\n" % Rxn['min'][4]) fid.writelines( "\n") elif dataset.lower() == 'eq36': fid.writelines('%-25s %s \n' % (Rxn['name'], Rxn['formula'])) fid.writelines(' sp.type = solid\n') fid.writelines('* EQ3/6 = com, alt, sup\n') fid.writelines(' revised = 01-Jan-2020\n') fid.writelines('* mol.wt. =%8.3f g/mol\n' % Rxn['MW']) fid.writelines(' V0PrTr = %8.3f cm**3/mol [source: %s ]\n' % (Rxn['V'], Rxn['source'])) fid.writelines('****\n') fid.writelines( " %s element(s):\n" % int(len(Rxn['elements'])/2)) for i in range(len(Rxn['elements'])): i = i + 1 if (i == 1) | (i == 7) | (i == 13): fid.writelines( " %9.4f " % float(Rxn['elements'][i - 1])) elif i % 2 != 0: fid.writelines( "%9.4f " % float(Rxn['elements'][i - 1])) else: fid.writelines( "%-9s " % (Rxn['elements'][i - 1])) if (i % 6 == 0) | (i == len(Rxn['elements'])): fid.writelines( "\n") fid.writelines('****\n') fid.writelines( " %s species in reaction:\n" % (Rxn['nSpec'] + 1)) fid.writelines( " %9.4f " % (-1)) fid.writelines( " %-21s " % Rxn['name']) for i in range(len(Rxn['spec'])): i = i + 1 fid.writelines( " %9.4f " % Rxn['coeff'][i-1]) fid.writelines( " %-21s " % Rxn['spec'][i-1]) if (i % 2 != 0) | (i == len(Rxn['spec'])): fid.writelines( "\n") fid.writelines('*\n') fid.writelines('**** logK grid [T, P @ Miscellaneous parameters]\n') for i in range(len(logK)): i = i + 1 if (i == 1) | (i == 5) | (i == 9) | (i == 13) | (i == 17): fid.writelines( " %9.4f" % logK[i-1]) else: fid.writelines( " %9.4f" % logK[i-1]) if (i % 4 == 0) | (i == len(logK)): fid.writelines( "\n") fid.writelines( "* gflag = 1 [reported delG0f used]\n" ) if Rxn['type'].find('plag') == 0: fid.writelines( "* extrapolation algorithm: supcrt92 [92joh/oel]\n" ) else: fid.writelines( "* extrapolation algorithm: supcrt92/water95\n" ) if Rxn['type'].find('serp') == 0: fid.writelines( "* ref-state data [source: Blanc et al 2015 ]\n" ) elif 'source' in Rxn: fid.writelines( "* ref-state data [source: %s ]\n" % Rxn['source']) else: fid.writelines( "* ref-state data [source: supcrt92? ]\n" ) fid.writelines( "* delG0f = %8.3f kcal/mol\n" % (Rxn['min'][2]/1000) ) if Rxn['type'] == 'Smectites': fid.writelines( "* delH0f = %8.3f kcal/mol\n" % (Rxn['min'][3]/1000) ) else: fid.writelines( "* delH0f = NaN kcal/mol\n") fid.writelines( "* S0PrTr = %8.3f cal/(mol*K)\n" % Rxn['min'][4]) fid.writelines( "* Cp coefficients [source: %s ]\n" % Rxn['source']) fid.writelines( "* T**0 = %11.8e \n" % (Rxn['min'][6]) ) fid.writelines( "* T**1 = %11.8e \n" % (Rxn['min'][7])) if Rxn['min'][8] < 1: fid.writelines( "* T**-2 = %12.8e \n" % (Rxn['min'][8])) else: fid.writelines( "* T**-2 = %11.8e \n" % (Rxn['min'][8])) if len(Rxn['min']) > 10: fid.writelines( "* T**-0.5 = %11.8e \n" % (Rxn['min'][9])) fid.writelines( "* T**2 = %11.8e \n" % (Rxn['min'][10])) fid.writelines( "+" + "-"*68 + "\n") elif dataset.lower() == 'pflotran': list_logk = ' '.join(str("%9.4f" % e) for e in list(logK)) Rxns_lst = ' '.join([ "%8.4f" % Rxn['coeff'][i]+' ' + "'%s'" % Rxn['spec'][i] for i in range(len(Rxn['spec']))]) info = "'%s'" % Rxn['name'] + ' ' + "%7.3f" % Rxn['V'] + ' ' + str(Rxn['nSpec']) + ' ' + \ Rxns_lst + ' ' + list_logk + ' ' + "%8.4f" % Rxn['MW'] fid.writelines('%s\n' % info) elif dataset.lower() == 'toughreact': logKfunc = lambda TK, *x: x[0]*np.log(TK) + x[1] + x[2]*TK + x[3]*TK**(-1) + x[4]*TK**(-2) x0 = [-3.19605e1, 2.06576e2, 3.73497e-2, -9.01862e3, 6.0111e5] logKcorr = curve_fit(logKfunc, TK[logK!=500].ravel(), logK[logK!=500].ravel(), p0 = x0, maxfev = 1000000)[0] list_logk = ' '.join(str("%9.4f" % e) for e in list(logK)) list_logKcorr = ' '.join(str("%.5e" % e) for e in list(logKcorr)) Rxns_spec = [j.replace('++++', '+4') if j.endswith('++++',0) else j.replace('+++', '+3') if j.endswith('+++',0) else j.replace('++', '+2') if j.endswith('++',0) else j.replace('----', '-4') if j.endswith('----',0) else j.replace('---', '-3') if j.endswith('---',0) else j.replace('--', '-2') if j.endswith('--',0) else j for j in Rxn['spec']] Rxn_lst = ' '.join([ "%8.4f" % Rxn['coeff'][i]+' ' + "'%s'" % Rxns_spec[i] for i in range(len(Rxn['spec']))]) info = "%-32s" % Rxn['name'] + "%8.3f" % Rxn['MW'] + " %7.3f" % Rxn['V'] +\ ' ' + str(Rxn['nSpec']) + ' ' + Rxn_lst info = "'%s'" % info[:len(info.split()[0])] + info[len(info.split()[0]) + 2:] fid.writelines('%s\n' % info) info = '%-35s' % Rxn['name'] + ' ' + list_logk info = "'%s'" % info[:len(info.split()[0])] + info[len(info.split()[0]) + 2:] fid.writelines('%s\n' % info) info = '%-35s' % Rxn['name'] + ' ' + list_logKcorr.replace('e','E') info = "'%s'" % info[:len(info.split()[0])] + info[len(info.split()[0]) + 2:] fid.writelines('%s\n' % info) return
[docs]class write_database(): """ Class to write the new database for either GWB, EQ3/6, ToughReact, Pflotran into a new folder called "output" \n Parameters ---------- T : float, vector Temperature [°C] \n P : float, vector Pressure [bar] \n cpx_Ca : string number of moles of Ca in solid solution of clinopyroxene (optional) if it is ommitted solid solution of clinopyroxene will not be included, 0 < nCa >=1 \n solid_solution : string, bool specify the inclusion of solid-solution [Yes/True or No/False], default is 'No' \n clay_thermo : float, vector specify the inclusion of clay thermodynamic properties [Yes/True or No/False], default is 'No' \n logK_form : float, vector specify the format of logK either as a set of eight values one for each of the dataset’s principal temperatures, or blocks of polynomial coefficients, [values, polycoeffs] default is 'a set of eight values' \n densityextrap : float, vector specify the utilization of density extrapolation [Yes/True or No/False], default is 'Yes' \n dbaccess : string direct-access database filename and location (optional) \n dbBerman_dir : string filename and location of the Berman mineral database (optional) \n dbHP_dir : string filename and location of the supcrtbl mineral and gas database, optional dbaccessformat : string, optional specify the direct-access/sequential-access database format, either 'speq' or 'supcrtbl', default is 'speq' sourcedb : string source database filename and location (optional) \n sourceformat : string source database format, either 'GWB' or 'EQ36', default is 'GWB' sourcedb_codecs : string specify the name of the encoding used to decode or encode the sourcedb file, optional objdb : string new database filename and location (optional) \n co2actmodel : string co2 activity model equation [Duan_Sun or Drummond] (optional), if not specified, default is 'Drummond' \n Dielec_method : string specify either 'FGL97' or 'JN91' or 'DEW' as the method to calculate dielectric constant, default is 'JN91' (optional) \n heatcap_method : string specify either 'SUPCRT' or 'Berman88' or 'HP11' or 'HF76' the method to calculate thermodynamic properties of any mineral, default is 'SUPCRT' \n ThermoInUnit : string specify either 'cal' or 'KJ' as the input units for species properties (optional), particularly used to covert KJ data to cal by supcrtaq function if not specified default - 'cal' dataset : string specify the dataset format, either 'GWB', 'EQ36', 'Pflotran' or 'ToughReact', default is old GWB database ['GWB'] (optional) \n print_msg : string, bool print debug message [True or False], default is False \n Returns ------- Output the new database to an ASCII file with filename described in 'objdb' if specified. Usage ------- With any Temperature and Pressure: (1) General format with default dielectric constant and CO2 activity model and exclusions of solid solutions for GWB \n write_database(T = T, P = P, cpx_Ca = nCa, dataset = 'GWB', sourceformat = 'GWB') \n (2) Inclusion of solid solutions and exclusion of solid solution of clinopyroxene and clay thermo \n write_database(T = T, P = P, solid_solution = 'Yes', clay_thermo = 'Yes', dataset = 'GWB', sourceformat = 'GWB') \n (3) Inclusion of all solid solutions and clay thermo with \emph{'Duan_Sun'} CO2 activity model and 'FGL97' dielectric constant calculation \n write_database(T = T, P = P, cpx_Ca = nCa, solid_solution = 'Yes', clay_thermo = 'Yes', co2actmodel = 'Duan_Sun', Dielec_method = 'FGL97', dataset = 'GWB', sourceformat = 'GWB') \n With any Temperature or Pressure on the steam saturation curve: (4) General format with default dielectric constant and CO2 activity model \n write_database(T = T, P = 'T', cpx_Ca = nCa, dataset = 'GWB', sourceformat = 'GWB'), \n where T is temperature in celsius and P is assigned a quoted character 'T' to reflect steam saturation pressure \n write_database(P = P, T = 'P', cpx_Ca = nCa, dataset = 'GWB', sourceformat = 'GWB'), \n where P is pressure in bar and T is assigned a quoted character 'P' to reflect steam saturation temperature Examples -------- >>> write_database(T=np.array([ 0.010, 25., 60., 100., 150., 200., 250., 300.]), P=200, cpx_Ca = 0.5, sourceformat = 'GWB', solid_solution = True, clay_thermo = True, dataset = 'GWB') \n >>> write_database(T=np.array([ 0.010, 25., 60., 100., 150., 200., 250., 300.]), P=200, cpx_Ca = 0.5, sourceformat = 'GWB', solid_solution = True, clay_thermo = True, dataset = 'Pflotran', sourcedb = './pygcc/default_db/thermo.2021.dat') \n """ kwargs = {"T": None, "P": None, "cpx_Ca": None, "solid_solution": None, "clay_thermo": None, "logK_form": None, 'dbBerman_dir': None, 'dbHP_dir': None, "dbaccess": None, "sourcedb": None, "objdb": None, "ThermoInUnit": 'cal', "co2actmodel": None, "Dielec_method": None, "heatcap_method": None, "dataset": None, "sourceformat": None, 'densityextrap': None, "sourcedb_codecs": None, "dbaccessformat": 'speq', "print_msg": False } def __init__(self, **kwargs): self.kwargs = write_database.kwargs.copy() self.__calc__(**kwargs)
[docs] def __calc__(self, **kwargs): self.kwargs.update(kwargs) self.T = self.kwargs["T"] self.P = self.kwargs["P"] if self.kwargs["dbaccess"] is None: self.dbaccess = './default_db/speq21.dat' self.dbaccess = os.path.join(os.path.dirname(os.path.abspath(__file__)), self.dbaccess) else: self.dbaccess = self.kwargs["dbaccess"] self.dbaccessformat = self.kwargs["dbaccessformat"] self.dbBerman_dir = self.kwargs['dbBerman_dir'] self.dbHP_dir = self.kwargs['dbHP_dir'] self.dataset = self.kwargs["dataset"] if self.kwargs["sourceformat"] == None: if self.dataset.upper() == 'GWB': self.sourceformat = 'GWB' elif self.dataset.upper() == 'EQ36': self.sourceformat = 'EQ36' else: self.sourceformat = self.kwargs["sourceformat"] if self.sourceformat.lower() == 'gwb': # options for all default database included with PyGCC if self.kwargs["sourcedb"] == 'thermo.com': self.sourcedb = './default_db/thermo.com.dat' self.sourcedb = os.path.join(os.path.dirname(os.path.abspath(__file__)), self.sourcedb) elif self.kwargs["sourcedb"] == 'thermo.2021': self.sourcedb = './default_db/thermo.2021.dat' self.sourcedb = os.path.join(os.path.dirname(os.path.abspath(__file__)), self.sourcedb) elif self.kwargs["sourcedb"] == 'thermo_latest': self.sourcedb = './default_db/thermo_latest.tdat' self.sourcedb = os.path.join(os.path.dirname(os.path.abspath(__file__)), self.sourcedb) elif self.kwargs["sourcedb"] == 'thermo_cemdata_mar': self.sourcedb = './default_db/thermo_cemdata_mar.tdat' self.sourcedb = os.path.join(os.path.dirname(os.path.abspath(__file__)), self.sourcedb) elif self.kwargs["sourcedb"] is None: self.sourcedb = './default_db/thermo.com.tdat' self.sourcedb = os.path.join(os.path.dirname(os.path.abspath(__file__)), self.sourcedb) else: self.sourcedb = self.kwargs["sourcedb"] elif self.sourceformat.lower() == 'eq36': if self.kwargs["sourcedb"] == 'data0' or self.kwargs["sourcedb"] is None: self.sourcedb = './default_db/data0.dat' self.sourcedb = os.path.join(os.path.dirname(os.path.abspath(__file__)), self.sourcedb) else: self.sourcedb = self.kwargs["sourcedb"] self.objdb = self.kwargs["objdb"] self.co2actmodel = self.kwargs["co2actmodel"] self.ThermoInUnit = self.kwargs['ThermoInUnit'] self.solid_solution = 'No' if (self.kwargs["solid_solution"] is None or self.kwargs['solid_solution'] is False) else 'Yes' if self.kwargs['solid_solution'] is True else self.kwargs["solid_solution"] self.clay_thermo = 'No' if (self.kwargs["clay_thermo"] is None or self.kwargs['clay_thermo'] is False) else 'Yes' if self.kwargs['clay_thermo'] is True else self.kwargs["clay_thermo"] self.densityextrap = 'Yes' if (self.kwargs["densityextrap"] is None or self.kwargs['densityextrap'] is True) else 'No' if self.kwargs['densityextrap'] is False else self.kwargs["densityextrap"] self.Dielec_method = 'JN91' if self.kwargs['Dielec_method'] is None else self.kwargs['Dielec_method'] self.heatcap_method = 'HP11' if self.dbHP_dir is not None else 'Berman88' if self.dbBerman_dir is not None else 'SUPCRT' if self.kwargs['heatcap_method'] is None else self.kwargs['heatcap_method'] self.logK_form = 'values' if self.kwargs["logK_form"] is None else self.kwargs['logK_form'] self.cpx_Ca = 0 if self.kwargs["cpx_Ca"] is None else self.kwargs["cpx_Ca"] self.sourcedb_codecs = self.kwargs["sourcedb_codecs"] self.dbr = db_reader(dbaccess = self.dbaccess, dbBerman_dir = self.dbBerman_dir, dbHP_dir = self.dbHP_dir, dbaccessformat = self.dbaccessformat, sourcedb = self.sourcedb, sourceformat = self.sourceformat, sourcedb_codecs = self.sourcedb_codecs) # condition to add Dimer from Sverjensky et al. 2014 if self.dbr.header_ref[self.dbr.dbaccessdic['SiO2(aq)'][1].strip(' ref:').split(' ')[0]].startswith('Sverjensky, D. A., Harrison, B., & Azzolini, D., 2014'): if 'Si2O4(aq)' not in self.dbr.dbaccessdic.keys() or 'Si2O4(aq)' not in self.dbr.sourcedic.keys(): warnings.warn('Warning: you are 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.') if (type(self.P) == str)or (type(self.T) == str): if self.P == 'T': if np.ndim(self.T) == 0: self.T = np.ravel(self.T) elif np.size(self.T) == 2: if self.T[-1] > 400: # for critical region of water (350 - 400C) self.T = np.array([0.01 if (x == 0)&(self.T[0] == 0) else self.T[0] if (x == 0)&(self.T[0] != 0) else roundup_hundredth(self.T[0] + x*(self.T[-1] - self.T[0])/(8 - 1)) if 350 <= (self.T[0] + x*(self.T[-1] - self.T[0])/(8 - 1)) <= 400 else round(self.T[-1]) if x == 7 else roundup_tenth(self.T[0] + x*(self.T[-1] - self.T[0])/(8 - 1)) for x in range(8)]) elif self.T[-1] > 350: self.T = np.array([0.01 if (x == 0)&(self.T[0] == 0) else self.T[0] if (x == 0)&(self.T[0] != 0) else roundup_hundredth(self.T[0] + x*(self.T[-1] - self.T[0])/(8 - 1)) if 350 <= (self.T[0] + x*(self.T[-1] - self.T[0])/(8 - 1)) <= 400 else roundup_tenth(self.T[0] + x*(self.T[-1] - self.T[0])/(8 - 1)) for x in range(8)]) else: self.T = np.array([roundup_tenth(j) if j != 0 else 0.01 for j in np.linspace(self.T[0], self.T[-1], 8)]) self.P = iapws95(T = self.T).P self.P[np.isnan(self.P) | (self.P < 1)] = 1.0133 elif self.T == 'P': self.T = iapws95(P = self.P).TC if np.ndim(self.T) == 0: self.T = np.ravel(self.T) elif np.size(self.T) == 2: if self.T[-1] > 400: # for critical region of water (350 - 400C) self.T = np.array([0.01 if (x == 0)&(self.T[0] == 0) else self.T[0] if (x == 0)&(self.T[0] != 0) else roundup_hundredth(self.T[0] + x*(self.T[-1] - self.T[0])/(8 - 1)) if 350 <= (self.T[0] + x*(self.T[-1] - self.T[0])/(8 - 1)) <= 400 else round(self.T[-1]) if x == 7 else roundup_tenth(self.T[0] + x*(self.T[-1] - self.T[0])/(8 - 1)) for x in range(8)]) elif self.T[-1] > 350: self.T = np.array([0.01 if (x == 0)&(self.T[0] == 0) else self.T[0] if (x == 0)&(self.T[0] != 0) else roundup_hundredth(self.T[0] + x*(self.T[-1] - self.T[0])/(8 - 1)) if 350 <= (self.T[0] + x*(self.T[-1] - self.T[0])/(8 - 1)) <= 400 else roundup_tenth(self.T[0] + x*(self.T[-1] - self.T[0])/(8 - 1)) for x in range(8)]) else: self.T = np.array([roundup_tenth(j) if j != 0 else 0.01 for j in np.linspace(self.T[0], self.T[-1], 8)]) if np.size(self.P) <= 2: if np.ndim(self.P) == 0: self.P = np.ravel(self.P) self.P = self.P[0]*np.ones(np.size(self.T)) if self.dataset.upper() == 'GWB': self.write_GWBdb(self.T, self.P) elif self.dataset.upper() == 'EQ36': self.write_EQ36db(self.T, self.P ) elif self.dataset.lower() == 'pflotran': self.write_pflotrandb(self.T, self.P ) elif self.dataset.lower() == 'toughreact': self.write_ToughReactdb(self.T, self.P ) delimiters = "/", "\\" patterns = '|'.join('(?<={})'.format(re.escape(delim)) for delim in delimiters) self.msg = 'Database for %s generated successfully using %s dielectric constant, \ %s %s source database' % (self.dataset.upper(), self.Dielec_method, re.split(patterns, self.sourcedb)[-1], self.sourceformat) if self.cpx_Ca != 0: self.msg += ', full solid solution included' elif self.solid_solution == 'Yes': self.msg += ', solid solution included with cpx excluded' if self.clay_thermo == 'Yes': self.msg += ', clay thermodynamics included' if self.kwargs["print_msg"] == True: print(self.msg)
[docs] def write_GWBdb(self, T, P ): """ This function writes the new GWB database into a new folder called "output" \n Parameters ---------- T : temperature [°C] \n P : pressure [bar] \n Returns ------- Outputs the new database to an ASCII file with filename described in 'objdb'. \n Usage ------- Example: (1) General format with default dielectric constant and CO2 activity model and exclusions of solid solutions \n write_GWBdb(T, P) \n (2) Inclusion of solid solutions and clay thermo and exclusion of solid solution of clinopyroxene \n write_GWBdb(T, P) \n (3) Inclusion of all solid solutions and clay thermo with \emph{'Duan_Sun'} CO2 activity model and 'JN91' dielectric constant calculation \n write_GWBdb(T, P) \n """ nCa_cpx = self.cpx_Ca; logK_form = self.logK_form solid_solution = self.solid_solution; clay_thermo = self.clay_thermo sourcedb = self.sourcedb objdb = self.objdb; Dielec_method = self.Dielec_method co2actmodel = self.co2actmodel; sourceformat = self.sourceformat heatcap_method = self.heatcap_method; densityextrap = self.densityextrap dbaccessdic, dbname, sourcedic, specielist = self.dbr.dbaccessdic, self.dbr.dbaccess, self.dbr.sourcedic, self.dbr.specielist MWdic, act_param, chargedic = self.dbr.MWdic, self.dbr.act_param, self.dbr.chargedic sourcedb_codecs = self.dbr.sourcedb_codecs if self.sourcedb_codecs is None else self.sourcedb_codecs activity_model = act_param['activity_model'] if sourceformat.upper() == 'GWB': Mineraltype, fugacity_info = self.dbr.Mineraltype, self.dbr.fugacity_info # fugacity_model = fugacity_info['fugacity_model'] elif sourceformat.upper() == 'EQ36': block_info, Elemlist = self.dbr.block_info, self.dbr.Elemlist if sourceformat.upper() == 'EQ36': dataset = 'tdat' dataset_format = 'apr20' else: dataset = sourcedb.split('.')[-1] dataset_format = act_param['dataset_format'] logK_form = 'values' if (logK_form is None) | (dataset_format in ['oct94', 'jul17']) else logK_form.lower() if Dielec_method.upper() == 'DEW': water = ZhangDuan(T = T, P = P) rho, dGH2O, dHH2O, SH2O = water.rho, water.G, np.nan*np.ones(len(T)), np.nan*np.ones(len(T)) else: water = iapws95(T = T, P = P) rho, dGH2O, dHH2O, SH2O = water.rho, water.G, water.H, water.S TK = convert_temperature( T, Out_Unit = 'K' ) if dataset_format in ['jan19', 'apr20', 'mar21'] and logK_form.lower() == 'polycoeffs': Tr = 298.15 logKfunc = lambda TK, *x: x[0] + x[1]*(TK - Tr) + x[2]*(TK**2 - Tr**2) + x[3]*((1/TK) - (1/Tr)) + \ x[4]*((1/TK**2) - (1/Tr**2)) + x[5]*np.log(TK/Tr) x0 = [-31.9605, 20.6576, 3.73497e-2, -9.01862, 6.0111, 2.5] if os.path.exists(os.path.join(os.getcwd(), 'output/GWB')) == False: os.makedirs(os.path.join(os.getcwd(), 'output/GWB')) periodic_table = open(os.path.join(os.path.dirname(os.path.abspath(__file__)), 'PeriodicTableJSON.json'), encoding='utf8') data = json.load(periodic_table) Element = {data['elements'][x]['symbol'] : pd.DataFrame([data['elements'][x]['name'], data['elements'][x]['atomic_mass']], index = ['name', 'mass']).T for x in range(len(data['elements']))} periodic_table.close() fid = open(sourcedb, 'r', encoding = sourcedb_codecs) missing_species = [] elemspeclist = [ symbol for x in specielist[0] for symbol, item in Element.items() if item.name[0][:5] == x[:5] ] if sourceformat.upper() == 'GWB' else specielist[0] form_del = [1] if sourceformat.upper() == 'GWB' else [1, 3, 4] all_species_source = [[i]+k for i, k in sourcedic.items() if i not in (['eh', 'e-', 'H2O']) or (dataset_format == 'mar21' and i not in specielist[6])] # all_species_source = [x.split()[0] if (dataset_format == 'mar21' and k in specielist[6]) # else [i]+k # for i, k in sourcedic.items() for x in k if i not in (['eh', 'e-', 'H2O']) or (dataset_format == 'mar21' and k in specielist[6] and x.strip('\n') and x.split()[0] not in ['a0', '*'] ) ] all_species_source = [[k for j, k in enumerate(all_species_source[i]) if (j not in form_del and k not in elemspeclist and str(k).strip('0123456789.- ') != '') ] if (i <= len(specielist[0])) else [k for j, k in enumerate(all_species_source[i]) if (j not in form_del and str(k).strip('0123456789.- ') != '') ] for i in range(len(all_species_source)) ] for num in range(len(all_species_source)): # if num < len(all_species_source): if dataset_format == 'mar21': lst = [v for v in all_species_source[num] if v not in (specielist[6] + ['eh', 'e-', 'H2O']) ] else: lst = [v for v in all_species_source[num] if v not in (['eh', 'e-', 'H2O']) ] bool_miss = [x.replace('(CH3COO)', '(Ac)').replace('CH3COO', '(Ac)') not in dbaccessdic.keys() for x in lst ] if any(bool_miss): sublist = [i for (i, v) in zip(lst, bool_miss) if v ] if lst[0] not in sublist: missing_species.append([lst[0]] + sublist) else: missing_species.append(sublist) missingfile = open(os.path.join(os.path.abspath("."), 'output', 'GWB', 'spxNotFound.txt'), 'w+') #missing_species = [i for i in missing_species if len(i)<1] for line in missing_species: if len(line) > 0: missingfile.writelines(line[0]) missingfile.writelines('\n') for i in range(len(line)): missingfile.writelines(' %s' % line[i]) missingfile.writelines('\n') missingfile.close() missing_species = [item for sublist in missing_species for item in sublist] missing_species = [i for n, i in enumerate(missing_species) if i not in missing_species[:n]] elem_avail = list(np.unique([k for i, j in enumerate([[i] + k for i, k in sourcedic.items() if i not in missing_species]) for l, k in enumerate(j) if i < len(specielist[1]) and str(k).strip('0123456789.- ') != '' and not str(k).endswith(("+", "-", '(aq)', '(g)')) and k not in ['O2', 'H2O'] and len(k) < 3])) if objdb == None: objdb = 'thermo.%sbars' % int(P[0]) logKnan_alert = False # timestr = '.' + time.strftime("%d%b%y_%H%M") fout = open(os.path.join(os.path.abspath("."),'output', 'GWB', objdb + '.' + dataset), 'w+') # + timestr dbname2 = 'and supcrtbl.dat' if self.dbHP_dir is not None else 'and berman.dat' if self.dbBerman_dir is not None else '' if sourceformat.upper() != 'EQ36': s = fid.readline() fout.writelines(s) s = fid.readline() fout.writelines(s[:16] + dataset_format + '\n') if dataset_format in ['jul17', 'jan19', 'apr20', 'mar21']: s = fid.readline() fout.writelines(s) s = fid.readline() fout.writelines(s) s = fid.readline() if s.strip(' \n*') != '': if s[:27].strip(' \n*:') == 'THERMODYNAMIC DATABASE': fout.writelines(s[:27] + dbname + '\n') else: fout.writelines('* THERMODYNAMIC DATABASE: ' + dbname + '\n') else: fout.writelines('* THERMODYNAMIC DATABASE: ' + dbname + '\n') s = fid.readline() if s[:15].strip(' \n*:') == 'generated by': fout.writelines(s[:15] + ': pyGeochemCalc, ' + time.ctime() + '\n') else: fout.writelines('* generated by' + ': pyGeochemCalc, ' + time.ctime() + '\n') s = fid.readline() fout.writelines(s) s = fid.readline() fout.writelines(s) else: fout.writelines('dataset of thermodynamic data for gwb programs \n' + \ 'dataset format: ' + dataset_format + '\n' + \ 'activity model: %s \n' % activity_model + \ 'fugacity model: tsonopoulos \n' + \ '* THERMODYNAMIC DATABASE: ' + dbname + ' ' + dbname2 + '\n' +\ '* generated by: pyGeochemCalc, ' + time.ctime() + '\n' +\ '* Output package: gwb \n' + \ '* Data set: com \n') if dataset_format in ['oct94', 'jul17']: fout.writelines('* Note: coefficients for calculating the activity coefficients \n' + \ '* for CO2 are based on Ref:\n' + \ '* S.E.Drummond,1981. Boiling and Mixing of Hydrothermal\n' + \ '* Fluids: Chemical Effects on Mineral Precipitation.\n') elif dataset_format in ['jan19', 'apr20', 'mar21'] and logK_form.lower() == 'polycoeffs': fout.writelines('* \n' + \ '* This thermo data file uses the polynomial expression of the logK values: \n' + \ '* \n' + \ '* log10 K(TK) = a + b*(TK-Tr) + c*(TK^2-Tr^2) + d*(1/TK-1/Tr) + e*(1/TK^2-1/Tr^2) + f*ln(TK/Tr) \n' + \ '* TK is the Temperature in Kelvin, Tr is the relative temperature (T = 298.15 K). \n') #skip lines till temperature rows for i in range(1000) : s = fid.readline() if s.strip('\n').strip('* ') in ['temperatures', 'temperatures (degC)', 'Temperature grid (degC)']: break fout.writelines( "\n") if s.strip('\n').strip('* ') == 'temperatures': if sourceformat.upper() != 'EQ36': fout.writelines(s[:-1] + '(degC)\n') else: fout.writelines('* ' + s[:-1] + '(degC)\n') elif s.strip('\n').strip('* ') == 'Temperature grid (degC)': fout.writelines('* temperatures (degC)\n') else: fout.writelines(s) for i in range(len(T)): i = i + 1 if (i == 1) | (i == 5): fout.writelines( " %9.4f" % T[i-1]) else: fout.writelines( " %9.4f" % T[i-1]) if (i % 4 == 0) | (i == len(T)): fout.writelines( "\n") #skip lines till pressure rows for i in range(1000) : s = fid.readline() if s.strip('\n').strip('* ') in ['pressures', 'pressures (bar)', 'Pressure grid (bars)']: break if s.strip('\n').strip('* ') == 'pressures': if sourceformat.upper() != 'EQ36': fout.writelines(s[:-1] + '(bar)\n') else: fout.writelines('* ' + s[:-1] + '(bar)\n') elif s.strip('\n').strip('* ') == 'Pressure grid (bars)': fout.writelines('* pressures (bar)\n') else: fout.writelines(s) for i in range(len(P)): i = i + 1 if (i == 1) | (i == 5): fout.writelines( " %9.4f" % P[i-1]) else: fout.writelines( " %9.4f" % P[i-1]) if (i % 4 == 0) | (i == len(P)): fout.writelines( "\n") #%% Calculation for debye huckel and bdot and water properties waterdielc = water_dielec(T = T, P = P, Dielec_method = Dielec_method) E, Adh, Bdh, bdot = waterdielc.E, waterdielc.Ah, waterdielc.Bh, waterdielc.bdot rhoEG = {'rho': rho, 'E': E, 'dGH2O': dGH2O} rhoEDB = {'rho': rho, 'E': E, 'Ah': Adh, 'Bh': Bdh} # Calculate the rho E G for density extrapolation method here so we have it below rhoEGextrap = {} if any(rhoEG['rho'] < 350): subBornptrs = rhoEG['rho'] < 350 for i, j in enumerate(zip(T[subBornptrs], P[subBornptrs])): rhoextrap = np.linspace(350, 550, 3) Pextrap = iapws95(T = j[0], rho = rhoextrap).P if Dielec_method.upper() != 'DEW' else ZhangDuan(T = j[0], rho = rhoextrap).P Textrap = j[0]*np.ones(np.size(Pextrap)) dGH2O = iapws95(T = Textrap, P = Pextrap).G if Dielec_method.upper() != 'DEW' else ZhangDuan(T = Textrap, P = Pextrap).G E = water_dielec(T = Textrap, P = Pextrap, Dielec_method = Dielec_method).E rhoextrap = np.around(rhoextrap, 3) rhoEGextrap['%d_%d' % (j[0], j[1])]= {'rho': rhoextrap,'E': E, 'dGH2O': dGH2O, 'Textrap': Textrap, 'Pextrap': Pextrap} #skip lines till adh rows # 'cco2' in s.strip('\n') 'log k for eh reaction' or 'Eh reaction: logKr' in s.strip('\n') for i in range(100) : s = fid.readline() if any(re.findall(r'|'.join(('(adh)', 'Debye-Huckel A_gamma')), s.strip('\n'), re.IGNORECASE)): break if sourceformat.upper() != 'EQ36': fout.writelines(s) else: fout.writelines('* debye huckel a (adh)\n') if (dataset_format in ['jan19', 'apr20', 'mar21']) & (logK_form.lower() == 'polycoeffs'): Adhcorr = curve_fit(logKfunc, TK[Adh!=500].ravel(), Adh[Adh!=500].ravel(), p0 = x0, maxfev = 1000000)[0] fout.writelines(' a= %15.9f ' % Adhcorr[0] + 'b= %15.9f ' % Adhcorr[1] + \ 'c= %15.6e\n' % Adhcorr[2]) fout.writelines(' d= %15.6f ' % Adhcorr[3] + 'e= %15.5f ' % Adhcorr[4] + \ 'f= %15.8f \n' % Adhcorr[5]) else: for i in range(len(Adh)): i = i + 1 if (i == 1) | (i == 5): fout.writelines( " %9.4f" % Adh[i-1]) else: fout.writelines( " %9.4f" % Adh[i-1]) if (i % 4 == 0) | (i == len(Adh)): fout.writelines( "\n") #skip lines till bdh rows for i in range(100): s = fid.readline() if any(re.findall(r'|'.join(('(bdh)', 'Debye-Huckel B_gamma')), s.strip('\n'), re.IGNORECASE)): break if sourceformat.upper() != 'EQ36': fout.writelines(s) else: fout.writelines('* debye huckel b (bdh)\n') if (dataset_format in ['jan19', 'apr20', 'mar21']) & (logK_form.lower() == 'polycoeffs'): Bdhcorr = curve_fit(logKfunc, TK[Bdh!=500].ravel(), Bdh[Bdh!=500].ravel(), p0 = x0, maxfev = 1000000)[0] fout.writelines(' a= %15.9f ' % Bdhcorr[0] + 'b= %15.9f ' % Bdhcorr[1] + \ 'c= %15.6e\n' % Bdhcorr[2]) fout.writelines(' d= %15.6f ' % Bdhcorr[3] + 'e= %15.5f ' % Bdhcorr[4] + \ 'f= %15.8f \n' % Bdhcorr[5]) else: for i in range(len(Bdh)): i = i + 1 if (i == 1) | (i == 5): fout.writelines( " %9.4f" % Bdh[i-1]) else: fout.writelines( " %9.4f" % Bdh[i-1]) if (i % 4 == 0) | (i == len(Bdh)): fout.writelines( "\n") #skip lines till bdot rows if activity_model in ['debye-huckel', 'b-dot']: for i in range(100) : s = fid.readline() if any(re.findall(r'|'.join(('bdot', 'B-dot')), s.strip('\n'), re.IGNORECASE)): break if sourceformat.upper() != 'EQ36': fout.writelines(s) else: if activity_model in ['debye-huckel', 'b-dot']: fout.writelines('* bdot\n') if activity_model in ['debye-huckel', 'b-dot']: for i in range(len(bdot)): i = i + 1 if (i == 1) | (i == 5): fout.writelines( " %9.4f" % bdot[i-1]) else: fout.writelines( " %9.4f" % bdot[i-1]) if (i % 4 == 0) | (i == len(bdot)): fout.writelines( "\n") #%% Calculation for co2 fitting coefs cco2 = gamma_correlation(T, P, co2actmodel) cco2 = cco2.T #.ravel() for j in range(4): fout.writelines('* c co2 %s\n' % (j+1)) for i in range(len(cco2)): i = i + 1 if (i == 1) | (i == 5): fout.writelines( " %9.4f" % cco2[i-1, j]) else: fout.writelines( " %9.4f" % cco2[i-1, j]) if (i % 4 == 0) | (i == len(cco2)): fout.writelines( "\n") #skip lines till c h2o rows if sourceformat.upper() != 'EQ36': for i in range(100) : s = fid.readline() # print(s) if s[:7] == '* c h2o': break #%% Calculation for h2o fitting coefs ch2o = aw_correlation(T, P, Dielec_method = Dielec_method, **rhoEDB) ch2o = ch2o.T #.ravel() for j in range(4): fout.writelines('* c h2o %s\n' % (j+1)) for i in range(len(cco2)): i = i + 1 if (i == 1) | (i == 5): fout.writelines( " %9.4f" % ch2o[i-1, j]) else: fout.writelines( " %9.4f" % ch2o[i-1, j]) if (i % 4 == 0) | (i == len(ch2o)): fout.writelines( "\n") if dataset_format == 'oct94': # #Copy and replace c h2o values for i in range(100) : s = fid.readline() if s[:7] == '* log k': break # fout.writelines(s) #%% Calculations for "log k for eh" rows fout.writelines(s) logK = calcRxnlogK( T = T, P = P, Specie = 'eh', dbaccessdic = dbaccessdic, sourcedic = sourcedic, specielist = specielist, Dielec_method = Dielec_method, rhoEG = rhoEG, sourceformat = sourceformat, densityextrap = densityextrap, rhoEGextrap = rhoEGextrap).logK logK = np.where(np.isnan(logK), 500, logK) # set abitrary 500 to nan values for i in range(len(logK)): i = i + 1 if (i == 1) | (i == 5): fout.writelines( " %9.4f" % logK[i-1]) else: fout.writelines( " %9.4f" % logK[i-1]) if (i % 4 == 0) | (i == len(logK)): fout.writelines( "\n") #%% Calculations for "log k for o2" #skip lines till "log k for o2" rows for i in range(50) : s = fid.readline() # print(s) if s[:14] == '* log k for o2': break fout.writelines(s) logK = calcRxnlogK( T = T, P = P, Specie = 'O2(g)', dbaccessdic = dbaccessdic, sourcedic = sourcedic, specielist = specielist, Dielec_method = Dielec_method, rhoEG = rhoEG, sourceformat = sourceformat, densityextrap = densityextrap, Specie_class = 'gases', heatcap_method = heatcap_method, rhoEGextrap = rhoEGextrap).logK logK = np.where(np.isnan(logK), 500, logK) # set abitrary 500 to nan values for i in range(len(logK)): i = i + 1 if (i == 1) | (i == 5): fout.writelines( " %9.4f" % logK[i-1]) else: fout.writelines( " %9.4f" % logK[i-1]) if (i % 4 == 0) | (i == len(logK)): fout.writelines( "\n") #%% Calculations for "log k for h2" #skip lines till "log k for h2" rows for i in range(50) : s = fid.readline() # print(s) if s[:14] == '* log k for h2': break fout.writelines(s) logK = calcRxnlogK( T = T, P = P, Specie = 'H2(g)', dbaccessdic = dbaccessdic, sourcedic = sourcedic, Specie_class = 'gases', specielist = specielist, Dielec_method = Dielec_method, sourceformat = sourceformat, densityextrap = densityextrap, rhoEG = rhoEG, heatcap_method = heatcap_method, rhoEGextrap = rhoEGextrap).logK logK = np.where(np.isnan(logK), 500, logK) # set abitrary 500 to nan values for i in range(len(logK)): i = i + 1 if (i == 1) | (i == 5): fout.writelines( " %9.4f" % logK[i-1]) else: fout.writelines( " %9.4f" % logK[i-1]) if (i % 4 == 0) | (i == len(logK)): fout.writelines( "\n") #%% Calculations for "log k for n2" #skip lines till "log k for n2" rows for i in range(50) : s = fid.readline() # print(s) if s[:14] == '* log k for n2': break fout.writelines(s) logK = calcRxnlogK( T = T, P = P, Specie = 'N2(g)', dbaccessdic = dbaccessdic, sourcedic = sourcedic, specielist = specielist, Dielec_method = Dielec_method, rhoEG = rhoEG, sourceformat = sourceformat, densityextrap = densityextrap, Specie_class = 'gases', heatcap_method = heatcap_method, rhoEGextrap = rhoEGextrap).logK logK = np.where(np.isnan(logK), 500, logK) # set abitrary 500 to nan values for i in range(len(logK)): i = i + 1 if (i == 1) | (i == 5): fout.writelines( " %9.4f" % logK[i-1]) else: fout.writelines( " %9.4f" % logK[i-1]) if (i % 4 == 0) | (i == len(logK)): fout.writelines( "\n") fout.writelines( "\n") else: fout.writelines('\n') if sourceformat.upper() == 'EQ36': fout.writelines( "\n") f_ionsize = open(os.path.join(os.path.dirname(os.path.abspath(__file__)), 'ion_size.txt'), 'r') Rd = f_ionsize.readlines() Rd = Rd[1:] ion_sizedic = {Rd[x].split()[0] : Rd[x].split()[1] for x in range(len(Rd))} f_ionsize.close() # Needed for rebalancing aqueous and mineral reactions in terms of O2(aq) instead of O2(g) dic = {'O2(g)' : ['', 1, '1.0000', 'O2(aq)']} logK_rebal = calcRxnlogK( T = T, P = P, Specie = 'O2(g)', dbaccessdic = dbaccessdic, sourcedic = dic, specielist = specielist, Dielec_method = Dielec_method, rhoEG = rhoEG, sourceformat = sourceformat, densityextrap = densityextrap, Specie_class = 'gases', heatcap_method = heatcap_method, rhoEGextrap = rhoEGextrap).logK logK_rebal = np.where(np.isnan(logK_rebal), 500, logK_rebal) # set abitrary 500 to nan values #%% Elements #skip lines till "elements" rows for i in range(5000) : s = fid.readline() if s.rstrip('\n').lstrip('0123456789.- ') == 'elements': break fout.writelines( " %s elements" % len(elem_avail)) fout.writelines( "\n\n") #copy and paste lines till "basis" rows for i in range(1000) : s = fid.readline() s_mod = s.split()[0] if s.strip('\n') != '' else s if s.rstrip('\n').lstrip('0123456789.- ') == 'basis species': break if sourceformat.upper() != 'EQ36': if s_mod[:5] in [Element[j].name[0][:5] for j in elem_avail]: fout.writelines(s) else: if s_mod in elem_avail and not s.startswith('+---'): fout.writelines('%-15s (%-2s) mole wt.= %8.4f g\n' % (Element[s_mod].name[0], s_mod, float(s.split()[1]))) else: continue fout.writelines( "\n-end-\n\n") #only used for counter counter = 0 for j in specielist[1]: if (j not in missing_species): counter +=1 fout.writelines( " %s basis species\n\n" % counter) #%% Basis reactions #skip lines till "redox couples" rows for i in range(2000): s = fid.readline() if s.rstrip('\n').lstrip('0123456789.- ') in ['redox couples', 'auxiliary basis species']: break for j in specielist[1]: if (j not in missing_species): k = j.replace('(CH3COO)', '(Ac)').replace('CH3COO', '(Ac)') if sourceformat.upper() != 'EQ36': fout.writelines('%s\n' % j) else: fout.writelines('%s\n' % j.replace('O2(g)', 'O2(aq)')) if sourceformat.upper() != 'EQ36': fout.writelines('%s\n' % chargedic[j]) else: if j == 'O2(g)': ionsize = [float(re.sub('[^0123456789\.]', '', x)) for x in block_info['O2(g)_b'] if x.strip('* ').startswith('DHazero')] else: ionsize = [float(re.sub('[^0123456789\.]', '', x)) for x in block_info[j] if x.strip('* ').startswith('DHazero')] if any(ionsize) and any([MWdic[j]]): fout.writelines(' charge= %d ion size= %.1f A mole wt.= %8.4f g\n' % (float(chargedic[j].split()[-1]), ionsize[0], MWdic[j])) else: formula = j if sourcedic[j][0] == '' else sourcedic[j][0] formula = formula.rstrip('(aq)(g)') ionsize = 3 if j.endswith('(aq)') else float(ion_sizedic[j]) if j in ion_sizedic.keys() else 500 fout.writelines(' charge= %d ion size= %.1f A mole wt.= %8.4f g\n' % (float(chargedic[j].split()[-1]), ionsize, calc_elem_count_molewt(formula, Elementdic = Element)[-1]) ) if sourceformat.upper() != 'EQ36': fout.writelines( " %s elements in species\n" % sourcedic[j][1]) Rxn = sourcedic[j][2:] else: fout.writelines( " %s elements in species\n" % int(len(Elemlist[j])/2) ) Rxn = Elemlist[j] for i in range(len(Rxn)): i = i + 1 if (i == 1) | (i == 7): fout.writelines( " %9.4f " % float(Rxn[i - 1])) elif i % 2 != 0: fout.writelines( "%9.4f " % float(Rxn[i - 1])) else: if sourceformat.upper() != 'EQ36': fout.writelines( "%-9s " % (Rxn[i - 1])) else: fout.writelines( "%-9s " % (Rxn[i - 1]).replace('O2(g)', 'O2(aq)')) if (i % 6 == 0) | (i == len(Rxn)): fout.writelines( "\n") fout.writelines( "* gflag = 1 [reported delG0f used]\n" ) fout.writelines( "* extrapolation algorithm: [92joh/oel]\n" ) if j != 'H2O': ref = dbaccessdic[k][1].split(' ')[0] ref = ref.split(':')[1] if ('ref' in ref) or ('REF' in ref) else ref dG, dH, S = dbaccessdic[k][2], dbaccessdic[k][3], dbaccessdic[k][4] else: dG, dH, S, ref = dGH2O[0], dHH2O[0], SH2O[0], 'iapws95/' + Dielec_method if Dielec_method in ['FGL97', 'JN91'] else 'DEW' fout.writelines( "* reference-state data source = %s\n" % ref) fout.writelines( "* delG0f = %8.3f kcal/mol\n" % (dG/1000) ) fout.writelines( "* delH0f = %8.3f kcal/mol\n" % (dH/1000)) fout.writelines( "* S0PrTr = %8.3f cal/(mol*K)\n" % S) fout.writelines( "\n") else: continue fout.writelines( "-end-\n\n") #only used for counter counter = 0 for j in specielist[2]: rxnlst = [b for a, b in enumerate(sourcedic[j]) if a not in [0, 1]] # remove formula and specie number rxnlst = [v for x, v in enumerate(rxnlst) if x % 2 != 0 ] # remove all coefficients if (j not in missing_species) and (len([k for k in rxnlst if k not in missing_species]) == len(rxnlst)): counter +=1 fout.writelines( " %s redox couples\n\n" % counter) #%% Redox reactions #skip lines till "aqueous species" rows for i in range(2000): s = fid.readline() # print(s) if s.rstrip('\n').lstrip('0123456789.- ') == 'aqueous species': break if sourceformat.upper() != 'EQ36': speclst = specielist[2] else: speclst = [ x for x in specielist[2] if x != 'O2(aq)'] for j in speclst: rxnlst = [b for a, b in enumerate(sourcedic[j]) if a not in [0, 1]] # remove formula and specie number rxnlst = [v for x, v in enumerate(rxnlst) if x % 2 != 0 ] # remove all coefficients if (j not in missing_species) and (len([i for i in rxnlst if i not in missing_species]) == len(rxnlst)): k = j.replace('(CH3COO)', '(Ac)').replace('CH3COO', '(Ac)') if dataset_format == 'oct94': fout.writelines('%s\n' % j) if sourcedic[j][0] != '': fout.writelines('* formula= %s\n' % sourcedic[j][0]) else: fout.writelines('* formula= %s\n' % dbaccessdic[k][0]) elif dataset_format in ['jul17', 'jan19', 'apr20', 'mar21']: if sourcedic[j][0] != '': fout.writelines('%-30s %s %s\n' % (j, 'formula=', sourcedic[j][0])) else: fout.writelines('%-30s %s %s\n' % (j, 'formula=', dbaccessdic[k][0])) if sourceformat.upper() != 'EQ36': fout.writelines('%s\n' % chargedic[j]) else: ionsize = [float(re.sub('[^0123456789\.]', '', x)) for x in block_info[j] if x.strip('* ').startswith('DHazero')] if any(ionsize) and any([MWdic[j]]): fout.writelines(' charge= %d ion size= %.1f A mole wt.= %8.4f g\n' % (float(chargedic[j].split()[-1]), ionsize[0], MWdic[j])) else: formula = j if sourcedic[j][0] == '' else sourcedic[j][0] formula = formula.rstrip('(aq)') ionsize = 3 if j.endswith('(aq)') else float(ion_sizedic[j]) if j in ion_sizedic.keys() else 500 fout.writelines(' charge= %d ion size= %.1f A mole wt.= %8.4f g\n' % (float(chargedic[j].split()[-1]), ionsize, calc_elem_count_molewt(formula, Elementdic = Element)[-1]) ) if sourceformat.upper() != 'EQ36': fout.writelines( " %s species in reaction\n" % sourcedic[j][1]) else: fout.writelines( " %s species in reaction\n" % (sourcedic[j][1] - 1) ) Rxn = sourcedic[j][2:] if sourceformat.upper() != 'EQ36' else sourcedic[j][4:] for i in range(len(Rxn)): i = i + 1 if (i == 1) | (i == 7): fout.writelines( " %9.4f " % float(Rxn[i - 1])) elif i % 2 != 0: fout.writelines( "%9.4f " % float(Rxn[i - 1])) else: fout.writelines( "%-9s " % (Rxn[i - 1])) if (i % 6 == 0) | (i == len(Rxn)): fout.writelines( "\n") logK = calcRxnlogK( T = T, P = P, Specie = j, dbaccessdic = dbaccessdic, sourcedic = sourcedic, specielist = specielist, Dielec_method = Dielec_method, rhoEG = rhoEG, sourceformat = sourceformat, densityextrap = densityextrap, Specie_class = 'aqueous', rhoEGextrap = rhoEGextrap, ThermoInUnit = self.ThermoInUnit) if densityextrap.lower() == 'yes': #any(np.isnan(logK)): if all(logK.nonsubBornptrs) == True: # if all densities are >= 350 logKnan_alert = False # turn off the prompts for using Density extrapolation else: logKnan_alert = True else: logKnan_alert = False # turn off the prompts for using Density extrapolation logK = logK.logK logK = np.where(np.isnan(logK), 500, logK) # set abitrary 500 to nan values if (dataset_format in ['jan19', 'apr20', 'mar21']) & (logK_form.lower() == 'polycoeffs'): logKcorr = curve_fit(logKfunc, TK[logK!=500].ravel(), logK[logK!=500].ravel(), p0 = x0, maxfev = 1000000)[0] fout.writelines(' a= %15.9f ' % logKcorr[0] + 'b= %15.9f ' % logKcorr[1] + \ 'c= %15.6e\n' % logKcorr[2]) fout.writelines(' d= %15.6f ' % logKcorr[3] + 'e= %15.5f ' % logKcorr[4] + \ 'f= %15.8f \n' % logKcorr[5]) fout.writelines(' TminK= %-15.2f ' % np.min(TK) + 'TmaxK= %-7.2f\n' % np.max(TK)) else: for i in range(len(logK)): i = i + 1 if (i == 1) | (i == 5): fout.writelines( " %9.4f" % logK[i-1]) else: fout.writelines( " %9.4f" % logK[i-1]) if (i % 4 == 0) | (i == len(logK)): fout.writelines( "\n") fout.writelines( "* gflag = 1 [reported delG0f used]\n" ) fout.writelines( "* extrapolation algorithm: supcrt92 [92joh/oel]\n" ) if j != 'H2O': ref = dbaccessdic[k][1].split(' ')[0] ref = ref.split(':')[1] if ('ref' in ref) or ('REF' in ref) else ref dG, dH, S = dbaccessdic[k][2], dbaccessdic[k][3], dbaccessdic[k][4] else: dG, dH, S, ref = dGH2O[0], dHH2O[0], SH2O[0], 'iapws95/' + Dielec_method if Dielec_method in ['FGL97', 'JN91'] else 'DEW' fout.writelines( "* reference-state data source = %s\n" % ref) fout.writelines( "* delG0f = %8.3f kcal/mol\n" % (dG/1000) ) fout.writelines( "* delH0f = %8.3f kcal/mol\n" % (dH/1000)) fout.writelines( "* S0PrTr = %8.3f cal/(mol*K)\n" % S) fout.writelines( "\n") else: continue fout.writelines( "-end-\n\n") if logKnan_alert == True: warnings.warn('Some temperature and pressure points are out of aqueous species HKF eqns regions of applicability, hence, density extrapolation has been applied') #only used for counter counter = 0 for j in specielist[3]: rxnlst = [b for a, b in enumerate(sourcedic[j]) if a not in [0, 1]] # remove formula and specie number rxnlst = [v for x, v in enumerate(rxnlst) if x % 2 != 0 ] # remove all coefficients if (j not in missing_species) and (len([k for k in rxnlst if k not in missing_species])==len(rxnlst)): counter +=1 fout.writelines( " %s aqueous species\n\n" % counter) #%% Aqueous reactions # skip lines till "minerals" rows for i in range(20000): s = fid.readline() if s.rstrip('\n').lstrip('0123456789.- ') in ['minerals', 'solids']: break for j in specielist[3]: rxnlst = [b for a, b in enumerate(sourcedic[j]) if a not in [0, 1]] # remove formula and specie number rxnlst = [v for x, v in enumerate(rxnlst) if x % 2 != 0 ] # remove all coefficients if (j not in missing_species) and (len([i for i in rxnlst if i not in missing_species]) == len(rxnlst)): k = j.replace('(CH3COO)', '(Ac)').replace('CH3COO', '(Ac)') # print(k) if dataset_format == 'oct94': fout.writelines('%s\n' % j) if sourcedic[j][0] != '': fout.writelines('* formula= %s\n' % sourcedic[j][0]) else: fout.writelines('* formula= %s\n' % dbaccessdic[k][0]) elif dataset_format in ['jul17', 'jan19', 'apr20', 'mar21']: if sourcedic[j][0] != '': fout.writelines('%-30s %s %s\n' % (j, 'formula=', sourcedic[j][0])) else: fout.writelines('%-30s %s %s\n' % (j, 'formula=', dbaccessdic[k][0])) if sourceformat.upper() != 'EQ36': fout.writelines('%s\n' % chargedic[j]) else: ionsize = [float(re.sub('[^0123456789\.]', '', x)) for x in block_info[j] if x.strip('* ').startswith('DHazero')] if any(ionsize) and any([MWdic[j]]): fout.writelines(' charge= %d ion size= %.1f A mole wt.= %8.4f g\n' % (float(chargedic[j].split()[-1]), ionsize[0], MWdic[j])) else: formula = j if sourcedic[j][0] == '' else sourcedic[j][0] formula = formula.rstrip('(aq)') ionsize = 3 if j.endswith('(aq)') else float(ion_sizedic[j]) if j in ion_sizedic.keys() else 500 fout.writelines(' charge= %d ion size= %.1f A mole wt.= %8.4f g\n' % (float(chargedic[j].split()[-1]), ionsize, calc_elem_count_molewt(formula, Elementdic = Element)[-1] ) ) if sourceformat.upper() != 'EQ36': fout.writelines( " %s species in reaction\n" % sourcedic[j][1]) else: fout.writelines( " %s species in reaction\n" % (sourcedic[j][1] - 1) ) Rxn = sourcedic[j][2:] if sourceformat.upper() != 'EQ36' else sourcedic[j][4:] for i in range(len(Rxn)): i = i + 1 if (i == 1) | (i == 7) | (i == 13): fout.writelines( " %9.4f " % float(Rxn[i - 1])) elif i % 2 != 0: fout.writelines( "%9.4f " % float(Rxn[i - 1])) else: if sourceformat.upper() != 'EQ36': fout.writelines( "%-9s " % (Rxn[i - 1])) else: fout.writelines( "%-9s " % (Rxn[i - 1]).replace('O2(g)', 'O2(aq)')) if (i % 6 == 0) | (i % 12 == 0) | (i == len(Rxn)): fout.writelines( "\n") logK = calcRxnlogK( T = T, P = P, Specie = j, dbaccessdic = dbaccessdic, sourcedic = sourcedic, specielist = specielist, Dielec_method = Dielec_method, rhoEG = rhoEG, sourceformat = sourceformat, densityextrap = densityextrap, Specie_class = 'aqueous', rhoEGextrap = rhoEGextrap, ThermoInUnit = self.ThermoInUnit).logK logK = np.where(np.isnan(logK), 500, logK) # set abitrary 500 to nan values if sourceformat.upper() == 'EQ36' and 'O2(g)' in Rxn: coeff_O2 = float(Rxn[Rxn.index('O2(g)') - 1]) logK = np.where(logK != 500, logK + coeff_O2*logK_rebal, logK) if (dataset_format in ['jan19', 'apr20', 'mar21']) & (logK_form.lower() == 'polycoeffs'): logKcorr = curve_fit(logKfunc, TK[logK!=500].ravel(), logK[logK!=500].ravel(), p0 = x0, maxfev = 1000000)[0] fout.writelines(' a= %15.9f ' % logKcorr[0] + 'b= %15.9f ' % logKcorr[1] + \ 'c= %15.6e\n' % logKcorr[2]) fout.writelines(' d= %15.6f ' % logKcorr[3] + 'e= %15.5f ' % logKcorr[4] + \ 'f= %15.8f \n' % logKcorr[5]) fout.writelines(' TminK= %-15.2f ' % np.min(TK) + 'TmaxK= %-7.2f\n' % np.max(TK)) else: for i in range(len(logK)): i = i + 1 if (i == 1) | (i == 5): fout.writelines( " %9.4f" % logK[i-1]) else: fout.writelines( " %9.4f" % logK[i-1]) if (i % 4 == 0) | (i == len(logK)): fout.writelines( "\n") fout.writelines( "* gflag = 1 [reported delG0f used]\n" ) fout.writelines( "* extrapolation algorithm: supcrt92 [92joh/oel]\n" ) if j != 'H2O': ref = dbaccessdic[k][1].split(' ')[0] ref = ref.split(':')[1] if ('ref' in ref) or ('REF' in ref) else ref dG, dH, S = dbaccessdic[k][2], dbaccessdic[k][3], dbaccessdic[k][4] else: dG, dH, S, ref = dGH2O[0], dHH2O[0], SH2O[0], 'iapws95/' + Dielec_method if Dielec_method in ['FGL97', 'JN91'] else 'DEW' fout.writelines( "* reference-state data source = %s\n" % ref) fout.writelines( "* delG0f = %8.3f kcal/mol\n" % (dG/1000) ) fout.writelines( "* delH0f = %8.3f kcal/mol\n" % (dH/1000)) fout.writelines( "* S0PrTr = %8.3f cal/(mol*K)\n" % S) fout.writelines( "\n") else: continue fout.writelines( "-end-\n\n") #%% free electron for tdat dataset format if dataset_format in ['jul17', 'jan19', 'apr20', 'mar21']: if sourceformat.upper() != 'EQ36': speclst = specielist[4] else: speclst = ['e-'] fout.writelines( " %s free electron\n\n" % len(speclst)) for j in speclst: fout.writelines('%s\n' % j) if sourceformat.upper() != 'EQ36': fout.writelines('%s\n' % chargedic[j]) fout.writelines( " %s species in reaction\n" % sourcedic[j][1]) else: fout.writelines(' charge= %d ion size= %.1f A mole wt.= %8.4f g\n' % (-1, 0, 0) ) fout.writelines( " %s species in reaction\n" % (sourcedic[j][1] - 1) ) Rxn = sourcedic[j][2:] if sourceformat.upper() != 'EQ36' else sourcedic[j][4:] for i in range(len(Rxn)): i = i + 1 if (i == 1) | (i == 7) | (i == 13): fout.writelines( " %9.4f " % float(Rxn[i - 1])) elif i % 2 != 0: fout.writelines( "%9.4f " % float(Rxn[i - 1])) else: fout.writelines( "%-9s " % (Rxn[i - 1])) if (i % 6 == 0) | (i % 12 == 0) | (i == len(Rxn)): fout.writelines( "\n") speclist = specielist #None if dataset_format == 'mar21' else specielist logK = calcRxnlogK( T = T, P = P, Specie = j, dbaccessdic = dbaccessdic, sourcedic = sourcedic, specielist = speclist, Dielec_method = Dielec_method, rhoEG = rhoEG, sourceformat = sourceformat, densityextrap = densityextrap, heatcap_method = heatcap_method, rhoEGextrap = rhoEGextrap).logK logK = np.where(np.isnan(logK), 500, logK) # set abitrary 500 to nan values if (dataset_format in ['jan19', 'apr20', 'mar21']) & (logK_form.lower() == 'polycoeffs'): logKcorr = curve_fit(logKfunc, TK[logK!=500].ravel(), logK[logK!=500].ravel(), p0 = x0, maxfev = 1000000)[0] fout.writelines(' a= %15.9f ' % logKcorr[0] + 'b= %15.9f ' % logKcorr[1] + \ 'c= %15.6e\n' % logKcorr[2]) fout.writelines(' d= %15.6f ' % logKcorr[3] + 'e= %15.5f ' % logKcorr[4] + \ 'f= %15.8f \n' % logKcorr[5]) fout.writelines(' TminK= %-15.2f ' % np.min(TK) + 'TmaxK= %-7.2f\n' % np.max(TK)) else: for i in range(len(logK)): i = i + 1 if (i == 1) | (i == 5): fout.writelines( " %9.4f" % logK[i-1]) else: fout.writelines( " %9.4f" % logK[i-1]) if (i % 4 == 0) | (i == len(logK)): fout.writelines( "\n") fout.writelines( "-end-\n\n") #only used for counter if clay_thermo is not None: if clay_thermo.lower() == 'yes': fclay = open(os.path.join(os.path.dirname(os.path.abspath(__file__)), 'clay_elements.dat'), 'r') Rd = fclay.readlines() Rd = [j.replace('-','_').strip('\n') for j in Rd] counter = len(Rd) else: counter = 0 else: # default no clay thermo calculation clay_thermo = 'no' counter = 0 solidsolution_no = 11 if solid_solution is not None: if solid_solution.lower() == 'yes': if nCa_cpx is None: nCa = 0 counter = counter + 3*solidsolution_no # no cpx else: nCa = nCa_cpx if nCa > 0: counter = counter + 4*solidsolution_no else: counter = counter + 3*solidsolution_no # no cpx else: counter = counter + 0 else: # default no solid solution counter = counter + 0 solid_solution = 'no' if sourceformat.upper() != 'EQ36': speclst = specielist[5] else: speclst = specielist[4] + specielist[5] for j in speclst: rxnlst = [b for a, b in enumerate(sourcedic[j]) if a not in [0, 1]] # remove formula and specie number rxnlst = [v for x, v in enumerate(rxnlst) if x % 2 != 0 ] # remove all coefficients if (j not in missing_species) and (len([k for k in rxnlst if k not in missing_species]) == len(rxnlst)): if (solid_solution.lower() == 'no') & (clay_thermo.lower() == 'no'): counter += 1 elif (solid_solution.lower() == 'yes') & (clay_thermo.lower() == 'no'): if j not in ['Anorthite', 'Albite', 'Forsterite', 'Fayalite', 'Enstatite', 'Ferrosilite']: counter += 1 else: continue elif (solid_solution.lower() == 'no') & (clay_thermo.lower() == 'yes'): if j not in [ Rd[h].split(',')[0] for h in range(len(Rd))]: counter += 1 else: continue else: if j not in ['Anorthite', 'Albite', 'Forsterite', 'Fayalite', 'Enstatite', 'Ferrosilite'] + [ Rd[h].split(',')[0] for h in range(len(Rd))]: counter += 1 else: continue if (solid_solution.lower() == 'yes') and (nCa == 1): fout.writelines( " %s minerals\n\n" % (counter - 2)) else: fout.writelines( " %s minerals\n\n" % (counter)) # print('-> Processing mineral output') #%% Mineral reactions #skip lines till "gases" rows for i in range(20000): s = fid.readline() if (dataset_format == 'mar21') and (s.rstrip('\n').lstrip('0123456789.- ') == 'solid solutions'): break elif dataset_format != 'mar21' and s.rstrip('\n').lstrip('0123456789.- ') == 'gases': break if solid_solution.lower() == 'yes': mineralcount = 0 fnlist = ['plagio', 'olivine', 'pyroxene', 'cpx'] if nCa > 0 else ['plagio', 'olivine', 'pyroxene'] for fn in fnlist: for nX in np.round(np.linspace(1, 0, solidsolution_no), 1): if fn != 'cpx': ss = calcRxnlogK(X = nX, T = T, P = P, Dielec_method = Dielec_method, rhoEG = rhoEG, dbaccessdic = dbaccessdic, Specie = fn, densityextrap = densityextrap, ThermoInUnit = self.ThermoInUnit, rhoEGextrap = rhoEGextrap) else: ss = calcRxnlogK(cpx_Ca = nCa, X = nX, T = T, P = P, Dielec_method = Dielec_method, rhoEG = rhoEG, dbaccessdic = dbaccessdic, Specie = fn, densityextrap = densityextrap, ThermoInUnit = self.ThermoInUnit, rhoEGextrap = rhoEGextrap) logK, Rxn = ss.logK, ss.Rxn logK = np.where(np.isnan(logK), 500, logK) # set abitrary 500 to nan values outputfmt(fout, logK, Rxn, TK, dataset = 'GWB', logK_form = logK_form) mineralcount += 1 # clay minerals if clay_thermo.lower() == 'yes': for i in range(len(Rd)): ss = calcRxnlogK(T = T, P = P, Specie = 'Clay', elem = Rd[i].split(','), dbaccessdic = dbaccessdic, ThermoInUnit = self.ThermoInUnit, rhoEG = rhoEG, rhoEGextrap = rhoEGextrap, densityextrap = densityextrap) logK, Rxn = ss.logK, ss.Rxn logK = np.where(np.isnan(logK), 500, logK) # set abitrary 500 to nan values outputfmt(fout, logK, Rxn, TK, dataset = 'GWB', logK_form = logK_form) for j in speclst: rxnlst = [b for a, b in enumerate(sourcedic[j]) if a not in [0, 1]] # remove formula and specie number rxnlst = [v for x, v in enumerate(rxnlst) if x % 2 != 0 ] # remove all coefficients if (j not in missing_species) and (len([k for k in rxnlst if k not in missing_species])==len(rxnlst)): if solid_solution.lower() == 'yes' and j in ['Anorthite', 'Albite', 'Forsterite', 'Fayalite', 'Enstatite', 'Ferrosilite']: continue elif solid_solution.lower() == 'yes' and (nCa == 1) and j in ['Diopside', 'Hedenbergite']: continue elif clay_thermo.lower() == 'yes' and j in [ Rd[h].split(',')[0] for h in range(len(Rd))]: continue else: k = j.replace('(CH3COO)', '(Ac)').replace('CH3COO', '(Ac)') fout.writelines("%s " % j) if sourceformat.upper() != 'EQ36': fout.writelines( "type= %s\n" % Mineraltype[j]) else: fout.writelines( "\n" ) if sourcedic[j][0] != '': fout.writelines(' formula= %s\n' % sourcedic[j][0]) else: fout.writelines(' formula= %s\n' % dbaccessdic[k][0]) fout.writelines(" mole vol.= %1.3f cc" % dbaccessdic[k][5]) if MWdic[j] != []: fout.writelines(" mole wt.= %1.4f g\n" % MWdic[j]) else: formula = k if sourcedic[j][0] == '' else sourcedic[j][0] formula = formula.rstrip('(aq)(am)')#.rstrip('+2').rstrip('+3').rstrip('+4') fout.writelines(" mole wt.= %1.4f g\n" % calc_elem_count_molewt(formula, Elementdic = Element)[-1] ) if sourceformat.upper() != 'EQ36': fout.writelines( " %s species in reaction\n" % sourcedic[j][1]) else: fout.writelines( " %s species in reaction\n" % (sourcedic[j][1] - 1) ) Rxn = sourcedic[j][2:] if sourceformat.upper() != 'EQ36' else sourcedic[j][4:] for i in range(len(Rxn)): i = i + 1 if (i == 1) | (i == 7) | (i == 13): fout.writelines("%9.4f " % float(Rxn[i - 1])) elif i % 2 != 0: fout.writelines("%9.4f " % float(Rxn[i - 1])) else: if sourceformat.upper() != 'EQ36': fout.writelines( "%-9s " % (Rxn[i - 1])) else: fout.writelines( "%-9s " % (Rxn[i - 1]).replace('O2(g)', 'O2(aq)')) if (i % 6 == 0) | (i % 12 == 0) | (i == len(Rxn)): fout.writelines( "\n") logK = calcRxnlogK( T = T, P = P, Specie = j, dbaccessdic = dbaccessdic, sourcedic = sourcedic, specielist = specielist, Dielec_method = Dielec_method, rhoEG = rhoEG, sourceformat = sourceformat, densityextrap = densityextrap, Specie_class = 'minerals', heatcap_method = heatcap_method, rhoEGextrap = rhoEGextrap).logK logK = np.where(np.isnan(logK), 500, logK) # set abitrary 500 to nan values if sourceformat.upper() == 'EQ36' and 'O2(g)' in Rxn: coeff_O2 = float(Rxn[Rxn.index('O2(g)') - 1]) logK = np.where(logK != 500, logK + coeff_O2*logK_rebal, logK) if (dataset_format in ['jan19', 'apr20', 'mar21']) & (logK_form.lower() == 'polycoeffs'): logKcorr = curve_fit(logKfunc, TK[logK!=500].ravel(), logK[logK!=500].ravel(), p0 = x0, maxfev = 1000000)[0] fout.writelines(' a= %15.9f ' % logKcorr[0] + 'b= %15.9f ' % logKcorr[1] + \ 'c= %15.6e\n' % logKcorr[2]) fout.writelines(' d= %15.6f ' % logKcorr[3] + 'e= %15.5f ' % logKcorr[4] + \ 'f= %15.8f \n' % logKcorr[5]) fout.writelines(' TminK= %-15.2f ' % np.min(TK) + 'TmaxK= %-7.2f\n' % np.max(TK)) else: for i in range(len(logK)): i = i + 1 if (i == 1) | (i == 5) | (i == 9): fout.writelines(" %9.4f" % logK[i-1]) else: fout.writelines(" %9.4f" % logK[i-1]) if (i % 4 == 0) | (i == len(logK)): fout.writelines( "\n") fout.writelines( "* gflag = 1 [reported delG0f used]\n" ) if dbaccessdic[k][0] == 'nan': fout.writelines( "* extrapolation algorithm: supcrt92/water95\n" ) else: fout.writelines( "* extrapolation algorithm: supcrt92 [92joh/oel]\n" ) ref = dbaccessdic[k][1].split(' ')[0] ref = ref.split(':')[1] if ('ref' in ref) or ('REF' in ref) else ref fout.writelines( "* reference-state data source = %s\n" % ref ) fout.writelines( "* delG0f = %8.3f kcal/mol\n" % (dbaccessdic[k][2]/1000/J_to_cal if heatcap_method.lower() == 'berman88' else dbaccessdic[k][2]/1000) ) fout.writelines( "* delH0f = %8.3f kcal/mol\n" % (dbaccessdic[k][3]/1000/J_to_cal if heatcap_method.lower() == 'berman88' else dbaccessdic[k][3]/1000)) fout.writelines( "* S0PrTr = %8.3f cal/(mol*K)\n" % (dbaccessdic[k][4]/J_to_cal if heatcap_method.lower() == 'berman88' else dbaccessdic[k][4])) fout.writelines( "\n") else: continue fout.writelines( "-end-\n\n") #%% Solid solutions for GWB internal calculation #only used for counter counter = 0; ss_list = [] if dataset_format == 'mar21': for j in specielist[6]: # print(j) lst = [x for x in sourcedic[j] if x.strip('\n')] minlst = [x.split()[0] for x in lst if x.split()[0] not in ['a0', '*']] # print(minlst) if all(x not in missing_species for x in minlst): counter +=1 ss_list.append(j) fout.writelines( " %s solid solutions\n\n" % counter ) for i in range(20000): s = fid.readline() if s.rstrip('\n').lstrip('0123456789.- ') == 'gases': break for j in ss_list: for k in range(len(sourcedic[j])): fout.writelines( sourcedic[j][k] ) if sourcedic[j][-1] != '\n': fout.writelines( "\n") fout.writelines( "-end-\n\n") #only used for counter counter = 0 spec = specielist[6] if dataset_format != 'mar21' else specielist[7] for j in spec: rxnlst = [b for a, b in enumerate(sourcedic[j]) if a not in [0, 1]] # remove formula and specie number rxnlst = [v for x, v in enumerate(rxnlst) if x % 2 != 0 ] # remove all coefficients if (j not in missing_species) and (len([k for k in rxnlst if k not in missing_species])==len(rxnlst)): counter +=1 fout.writelines( " %s gases\n\n" % counter ) #%% Gas reactions #skip lines till "oxides" rows for i in range(20000): s = fid.readline() # print(s) if s.rstrip('\n').lstrip('0123456789.- ') in ['oxides', 'solid solutions']: break for j in spec: rxnlst = [b for a, b in enumerate(sourcedic[j]) if a not in [0, 1]] # remove formula and specie number rxnlst = [v for x, v in enumerate(rxnlst) if x % 2 != 0 ] # remove all coefficients if (j not in missing_species) and (len([k for k in rxnlst if k not in missing_species])==len(rxnlst)): k = j.replace('(CH3COO)', '(Ac)').replace('CH3COO', '(Ac)') fout.writelines("%s\n" % j) if MWdic[j] != []: fout.writelines( " mole wt.= %1.4f g\n" % MWdic[j]) else: formula = j.rstrip('(g)') fout.writelines(" mole wt.= %1.4f g\n" % calc_elem_count_molewt(formula, Elementdic = Element)[-1] ) if sourceformat.upper() == 'EQ36': fugacity_info = {'fugacity_chi': {'CH4(g)': ' chi= -537.779 1.54946 -.000927827 1.20861 -.00370814 3.33804e-6\n', 'CO2(g)': ' chi= -1430.87 3.598 -.00227376 3.47644 -.0104247 8.46271e-6\n', 'H2(g)': ' chi= -12.5908 .259789 -7.2473e-5 .00471947 -2.69962e-5 2.15622e-8\n', 'H2O(g)': ' chi= -6191.41 14.8528 -.00914267 -66.3326 .18277 -.00013274\n'}, 'fugacity_Pcrit': {'Ar(g)': ' Pcrit= 48.7 bar Tcrit= 150.8 K omega= .001\n', 'CH4(g)': ' Pcrit= 46.0 bar Tcrit= 190.4 K omega= .011\n', 'CO2(g)': ' Pcrit= 73.8 bar Tcrit= 304.1 K omega= .239\n', 'H2(g)': ' Pcrit= 13.0 bar Tcrit= 33.2 K omega= -.218\n', 'H2O(g)': ' Pcrit= 221.2 bar Tcrit= 647.3 K omega= .344 a=-.0109 b= 0.0\n', 'H2S(g)': ' Pcrit= 89.4 bar Tcrit= 373.2 K omega= .097\n', 'He(g)': ' Pcrit= 2.27 bar Tcrit= 5.19 K omega= -.365\n', 'N2(g)': ' Pcrit= 33.9 bar Tcrit= 126.2 K omega= .039\n', 'NH3(g)': ' Pcrit= 113.5 bar Tcrit= 405.5 K omega= .250\n', 'O2(g)': ' Pcrit= 50.4 bar Tcrit= 154.6 K omega= .025\n', 'SO2(g)': ' Pcrit= 78.8 bar Tcrit= 430.8 K omega= .256\n'}} if dataset_format in ['jul17', 'jan19', 'apr20', 'mar21']: if j in fugacity_info['fugacity_chi'].keys(): fout.writelines("%s" % fugacity_info['fugacity_chi'][j]) if j in fugacity_info['fugacity_Pcrit'].keys(): fout.writelines("%s" % fugacity_info['fugacity_Pcrit'][j]) if sourceformat.upper() != 'EQ36': fout.writelines( " %s species in reaction\n" % sourcedic[j][1]) else: fout.writelines( " %s species in reaction\n" % (sourcedic[j][1] - 1) ) Rxn = sourcedic[j][2:] if sourceformat.upper() != 'EQ36' else sourcedic[j][4:] for i in range(len(Rxn)): i = i + 1 if (i == 1) | (i == 7) | (i == 13): fout.writelines( " %9.4f " % float(Rxn[i - 1])) elif i % 2 != 0: fout.writelines( "%9.4f " % float(Rxn[i - 1])) else: if sourceformat.upper() != 'EQ36': fout.writelines( "%-9s " % (Rxn[i - 1])) else: fout.writelines( "%-9s " % (Rxn[i - 1]).replace('O2(g)', 'O2(aq)')) if (i % 6 == 0) | (i % 12 == 0) | (i == len(Rxn)): fout.writelines( "\n") logK = calcRxnlogK( T = T, P = P, Specie = j, dbaccessdic = dbaccessdic, sourcedic = sourcedic, specielist = specielist, Dielec_method = Dielec_method, rhoEG = rhoEG, sourceformat = sourceformat, densityextrap = densityextrap, Specie_class = 'gases', heatcap_method = heatcap_method, rhoEGextrap = rhoEGextrap).logK logK = np.where(np.isnan(logK), 500, logK) # set abitrary 500 to nan values if sourceformat.upper() == 'EQ36' and 'O2(g)' in Rxn: coeff_O2 = float(Rxn[Rxn.index('O2(g)') - 1]) logK = np.where(logK != 500, logK + coeff_O2*logK_rebal, logK) if (dataset_format in ['jan19', 'apr20', 'mar21']) & (logK_form.lower() == 'polycoeffs'): TK = convert_temperature( T, Out_Unit = 'K' ) logKcorr = curve_fit(logKfunc, TK[logK!=500].ravel(), logK[logK!=500].ravel(), p0 = x0, maxfev = 1000000)[0] fout.writelines(' a= %15.9f ' % logKcorr[0] + 'b= %15.9f ' % logKcorr[1] + \ 'c= %15.6e\n' % logKcorr[2]) fout.writelines(' d= %15.6f ' % logKcorr[3] + 'e= %15.5f ' % logKcorr[4] + \ 'f= %15.8f \n' % logKcorr[5]) fout.writelines(' TminK= %-15.2f ' % np.min(TK) + 'TmaxK= %-7.2f\n' % np.max(TK)) else: for i in range(len(logK)): i = i + 1 fout.writelines( " %9.4f" % logK[i-1]) if (i % 4 == 0) | (i == len(logK)): fout.writelines( "\n") ref = dbaccessdic[k][1].split(' ')[0] ref = ref.split(':')[1] if ('ref' in ref) or ('REF' in ref) else ref fout.writelines( "* gflag = 1 [reported delG0f used]\n" ) fout.writelines( "* extrapolation algorithm: supcrt92 [92joh/oel]\n" ) fout.writelines( "* reference-state data source = %s\n" % ref ) fout.writelines( "* delG0f = %8.3f kcal/mol\n" % (dbaccessdic[k][2]/1000) ) fout.writelines( "* delH0f = %8.3f kcal/mol\n" % (dbaccessdic[k][3]/1000)) fout.writelines( "* S0PrTr = %8.3f cal/(mol*K)\n" % dbaccessdic[k][4]) fout.writelines( "\n") else: continue fout.writelines( "-end-\n\n") if sourceformat.upper() != 'EQ36': #only used for counter counter = 0 spec = specielist[7] if dataset_format != 'mar21' else specielist[8] for j in spec: rxnlst = [b for a, b in enumerate(sourcedic[j]) if a not in [0, 1]] # remove formula and specie number rxnlst = [v for x, v in enumerate(rxnlst) if x % 2 != 0 ] # remove all coefficients if (len([k for k in rxnlst if k not in missing_species])==len(rxnlst)): counter +=1 fout.writelines( " %s oxides\n\n" % counter ) #%% Oxides reactions #skip lines till "references" rows for i in range(20000): s = fid.readline() # print(s) if s.strip(' \n*').lstrip('0123456789.- ').startswith(("references", 'virial coefficients', 'Virial coefficients', 'SIT epsilon coefficients', 'Pitzer parameters')): break for j in spec: rxnlst = [b for a, b in enumerate(sourcedic[j]) if a not in [0, 1]] # remove formula and specie number rxnlst = [v for x, v in enumerate(rxnlst) if x % 2 != 0 ] # remove all coefficients if (len([k for k in rxnlst if k not in missing_species]) == len(rxnlst)): fout.writelines("%s\n" % j) fout.writelines( " mole wt.= %1.4f g\n" % MWdic[j]) fout.writelines( " %s species in reaction\n" % sourcedic[j][1]) Rxn = sourcedic[j][2:] for i in range(len(Rxn)): i = i + 1 if (i == 1) | (i == 7) | (i == 13): fout.writelines( " %9.4f " % float(Rxn[i - 1])) elif i % 2 != 0: fout.writelines( "%9.4f " % float(Rxn[i - 1])) else: fout.writelines( "%-9s " % (Rxn[i - 1])) if (i % 6 == 0) | (i % 12 == 0) | (i == len(Rxn)): fout.writelines( "\n") fout.writelines( "\n") else: continue else: fout.writelines( " 0 oxides\n\n" ) #%% Pitzer parameters if activity_model != 'h-m-w': fout.writelines( "-end-\n\n") elif activity_model == 'h-m-w': fout.writelines( "-end-\n*\n") if sourceformat.upper() == 'GWB': fout.writelines(s) for i in range(15000): s = fid.readline() fout.writelines(s) if not s.strip(' \n').startswith("*"): break for i in range(15000): s = fid.readline() if any(s.lstrip().rstrip('\n').startswith(x) for x in act_param['act_list']): rowlst = []; rowlst.append(s) for i in range(150): s = fid.readline() if s.lstrip().rstrip('\n').startswith(""): break rowlst.append(s) if all([x not in missing_species for x in rowlst[0].rstrip('\n').split()]): fout.writelines(rowlst) fout.writelines( "\n") else: fout.writelines(s) elif sourceformat.upper() == 'EQ36': fout.writelines("* Pitzer parameters are represented by the 25C-centric four-term temperature \n" + \ "* function given by: \n" + "* \n" + \ "* x(T) = a1 + a2*(1/T - 1/298.15) + a3*ln(T/298.15) + a4*(T - 298.15) \n" + "* \n" + \ "* where T is temperature in Kelvin and a1 through a4 denote the temperature \n" + \ "* function fitting coefficients for the temperature-dependent Pitzer \n" + \ "* parameters. The conversion of non-standard or expanded forms of Pitzer \n" + \ "* interaction parameters recently adopted by several workers for highly \n" + \ "* soluble salts to the standard form currently embedded in EQ3/6 Version 8.0 \n" + \ "* was conducted using the approach described in 02Rard/Wij. This conversion \n" + \ "* imposes usage limits on these parameters within a valid range of temperature \n" + \ "* and ionic strength. \n" + "* \n" + \ "* In GWB (Version 9.0.1 or lower) Pitzer parameters are represented by the \n" + \ "* 25C-centric five-term temperature function given by: \n" + "* \n" + \ "* val = val25 + c1*(Tk-Tr) + c2*(1/Tk-1/Tr) + c3*ln(Tk/Tr) + c4(Tk^2-Tr^2) \n" + "* \n" + \ "* So the last temperature term (c4(Tk^2-Tr^2)) will be set to zero \n" + \ "* since no such term is available in the data0.ypf.R0. \n" + "* \n" + \ "* In GWB (Version 9.0.2 or higher) Pitzer parameters are represented by the \n" + \ "* 25C-centric six-term temperature function given by: \n" + \ "* val = val25 + c1*(Tk-Tr) + c2*(1/Tk-1/Tr) + c3*ln(Tk/Tr) + c4(Tk^2-Tr^2) + c5(1/Tk^2-1/Tr^2) \n" + \ "* So the last two temperature terms (c4(Tk^2-Tr^2) and c5(1/Tk^2-1/Tr^2)) \n" + \ "* will be set to zero (or left blank) since no such term is available. \n\n") for k in act_param['alpha_beta'].keys(): if all([x not in missing_species for x in k.rstrip('\n').split()]): ks = k.rstrip('\n').split() fout.writelines('%-8s\n' % ks[0]) if len(ks) == 1 else fout.writelines('%-8s %-8s\n' % (ks[0], ks[1])) if len(ks) == 2 else fout.writelines('%-8s %-8s %-8s\n' % (ks[0], ks[1], ks[2])) lst = ['beta0', 'beta1', 'beta2', 'cphi', 'alpha1', 'alpha2'] for order in lst: fout.writelines(' %-6s = %s \n' % (order, act_param[order][k])) fout.writelines('\n') fout.writelines('-end- end of beta set, begin with theta set of 2nd virial coefficients \n\n') for k in act_param['theta'].keys(): if all([x not in missing_species for x in k.rstrip('\n').split()]): ks = k.rstrip('\n').split() fout.writelines('%-8s\n' % ks[0]) if len(ks) == 1 else fout.writelines('%-8s %-8s\n' % (ks[0], ks[1])) if len(ks) == 2 else fout.writelines('%-8s %-8s %-8s\n' % (ks[0], ks[1], ks[2])) fout.writelines(' %-6s = %s \n' % ('theta', act_param['theta'][k])) fout.writelines('\n') fout.writelines('-end- end of theta set, begin with lambda set \n\n') for k in act_param['lambda'].keys(): if all([x not in missing_species for x in k.rstrip('\n').split()]): ks = k.rstrip('\n').split() fout.writelines('%-8s\n' % ks[0]) if len(ks) == 1 else fout.writelines('%-8s %-8s\n' % (ks[0], ks[1])) if len(ks) == 2 else fout.writelines('%-8s %-8s %-8s\n' % (ks[0], ks[1], ks[2])) fout.writelines(' %-6s = %s \n' % ('lambda', act_param['lambda'][k])) fout.writelines('\n') fout.writelines('-end- end of lambda set, begin with psi set \n\n') for k in act_param['psi'].keys(): if all([x not in missing_species for x in k.rstrip('\n').split()]): ks = k.rstrip('\n').split() fout.writelines('%-8s\n' % ks[0]) if len(ks) == 1 else fout.writelines('%-8s %-8s\n' % (ks[0], ks[1])) if len(ks) == 2 else fout.writelines('%-8s %-8s %-8s\n' % (ks[0], ks[1], ks[2])) fout.writelines(' %-6s = %s \n' % ('psi', act_param['psi'][k])) fout.writelines('\n') fout.writelines('-end- end of psi set \n') fout.writelines('* references\n') fout.writelines('** Please copy references to here from the corresponding\n') fout.writelines('** sequential-access version of the direct-access SUPCRT database.\n* stop.\n\n') #%% close all files fid.close() fout.close() if clay_thermo.lower() == 'yes': fclay.close() return print('Success, your new GWB database is ready for download')
[docs] def write_EQ36db(self, T, P ): """ This function writes the new EQ3/6 database into a new folder called "output" \n Parameters ---------- T : temperature [°C] \n P : pressure [bar] \n Returns ------- Outputs the new database to an ASCII file with filename described in 'objdb'. \n Usage ------- Example: (1) General format with default dielectric constant and CO2 activity model and exclusions of solid solutions \n write_EQ36db(T, P ) \n (2) Inclusion of solid solutions and clay thermo and exclusion of solid solution of clinopyroxene \n write_EQ36db(T, P ) \n (3) Inclusion of all solid solutions and clay thermo with \emph{'Duan_Sun'} CO2 activity model and 'FGL97' dielectric constant calculation \n write_EQ36db(T, P ) \n """ nCa_cpx = self.cpx_Ca; solid_solution = self.solid_solution; clay_thermo = self.clay_thermo sourcedb = self.sourcedb objdb = self.objdb; Dielec_method = self.Dielec_method heatcap_method = self.heatcap_method; sourceformat = self.sourceformat densityextrap = self.densityextrap dbaccessdic, dbname, sourcedic, specielist = self.dbr.dbaccessdic, self.dbr.dbaccess, self.dbr.sourcedic, self.dbr.specielist if sourceformat.upper() == 'GWB': MWdic, act_param, chargedic = self.dbr.MWdic, self.dbr.act_param, self.dbr.chargedic dataset = '1kbu' periodic_table = open(os.path.join(os.path.dirname(os.path.abspath(__file__)), 'PeriodicTableJSON.json'), encoding='utf8') data = json.load(periodic_table) Element = {data['elements'][x]['symbol'] : pd.DataFrame([data['elements'][x]['name'], data['elements'][x]['atomic_mass']], index = ['name', 'mass']).T for x in range(len(data['elements']))} periodic_table.close() specielist[0] = [ symbol for x in specielist[0] for symbol, item in Element.items() if item.name[0][:4] == x[:4] ] elif sourceformat.upper() == 'EQ36': block_info, Elemlist, act_param = self.dbr.block_info, self.dbr.Elemlist, self.dbr.act_param dataset = sourcedb.split('.')[-1] if os.path.exists(os.path.join(os.getcwd(), 'output/EQ36')) == False: os.makedirs(os.path.join(os.getcwd(), 'output/EQ36')) logKnan_alert = False missing_species = [] all_species_source = [[i]+k for i, k in sourcedic.items() if i not in (['eh', 'e-', 'H2O']) ] all_species_source = [[k for j, k in enumerate(all_species_source[i]) if (j not in [1, 3] and k not in all_species_source[i][:j] and k not in specielist[0] and str(k).strip('0123456789.- ') != '') ] if (i <= len(specielist[0])) else [k for j, k in enumerate(all_species_source[i]) if (j not in [1, 3] and k not in all_species_source[i][:j] and str(k).strip('0123456789.- ') != '') ] for i in range(len(all_species_source)) ] for num in range(len(all_species_source)): # if num < len(all_species_source): lst = [v for v in all_species_source[num] if v not in (specielist[7] + ['eh', 'e-', 'H2O']) ] bool_miss = [x.replace('(CH3COO)', '(Ac)').replace('CH3COO', '(Ac)') not in dbaccessdic.keys() for x in lst ] if any(bool_miss): sublist = [i for (i, v) in zip(lst, bool_miss) if v ] if lst[0] not in sublist: missing_species.append([lst[0]] + sublist) else: missing_species.append(sublist) missingfile = open(os.path.join(os.path.abspath("."),'output', 'EQ36', 'spxNotFound.txt'), 'w') #missing_species = [i for i in missing_species if len(i)<1] for line in missing_species: if len(line) > 0: missingfile.writelines(line[0]) missingfile.writelines('\n') for i in range(len(line)): missingfile.writelines(' %s' % line[i]) missingfile.writelines('\n') missingfile.close() missing_species = [item for sublist in missing_species for item in sublist] elem_avail = list(set([k for i, j in enumerate([[i] + k for i, k in sourcedic.items() if i not in missing_species]) for l, k in enumerate(j) if i < len(specielist[1]) and str(k).strip('0123456789.- ') != '' and not str(k).endswith(("+", "-", '(aq)', '(g)')) and k not in ['O2', 'H2O'] and len(k) < 3])) all_species_source = [k for j in all_species_source for i, k in enumerate(j) ] all_species_source = list(set(all_species_source)) all_species_avail = [j for j in all_species_source if j not in missing_species] if objdb == None: objdb = 'data0.%s' % (int(P[0])) fout = open(os.path.join(os.path.abspath("."),'output', 'EQ36', objdb + '.%s' % dataset), 'w+') # + timestr fid = open(sourcedb, 'r') dbname2 = 'and supcrtbl.dat' if self.dbHP_dir is not None else 'and berman.dat' if self.dbBerman_dir is not None else '' if sourceformat.upper() != 'GWB': s = fid.readline() fout.writelines(s) for i in range(500) : s = fid.readline() if s.startswith('Generated', 0) | s.startswith('Data', 0) | (s.rstrip('\n') == ""): break fout.writelines('CII: ' + ' pyGeochemCalc.2021' + '\n') fout.writelines('Generated by: ' + ' pyGeochemCalc, ' + time.ctime() + '\n') fout.writelines('Output package: eq3\n' + 'Data set: ' + dbname + ' ' + dbname2 + '\n') else: fout.writelines('data0.com.RX ! dataset converted from gwb database \n' + \ 'CII: ' + ' pyGeochemCalc.2021' + '\n' + \ 'Generated by: ' + ' pyGeochemCalc, ' + time.ctime() + '\n' + \ 'Output package: eq3\n' + 'Data set: ' + dbname + ' ' + dbname2 + '\n') if Dielec_method.upper() == 'DEW': water = ZhangDuan(T = T, P = P) rho, dGH2O, dHH2O, SH2O = water.rho, water.G, np.nan*np.ones(len(T)), np.nan*np.ones(len(T)) else: water = iapws95(T = T, P = P) rho, dGH2O, dHH2O, SH2O = water.rho, water.G, water.H, water.S #copy and paste lines till temperature rows for i in range(2500): s = fid.readline() if sourceformat.upper() != 'GWB': fout.writelines(s) if s.startswith('+', 0): break else: if s.strip('\n').strip('* ') in ['temperatures', 'temperatures (degC)', 'Temperature grid (degC)']: break if sourceformat.upper() != 'GWB': s = fid.readline() fout.writelines(s) s = fid.readline() fout.writelines(s) s = fid.readline() fout.writelines(s) else: fout.write( "+" + "-"*68 + "\n") fout.write("Miscellaneous parameters\n") fout.write( "+" + "-"*68 + "\n") fout.write("Temperature limits (degC)\n") fout.writelines(' %9.4f %9.4f\n' % (T[0], T[-1])) if sourceformat.upper() != 'GWB': for i in range(50) : s = fid.readline() if s.strip('\n').strip('* ') in ['temperatures', 'temperatures (degC)', 'Temperature grid (degC)']: break fout.writelines(s) if sourceformat.upper() != 'GWB' else fout.writelines(s.strip('* ').split('(')[0] + '\n') for i in range(len(T)): i = i + 1 if (i == 1) | (i == 5): fout.writelines( " %9.4f" % T[i-1]) else: fout.writelines( " %9.4f" % T[i-1]) if (i % 4 == 0) | (i == len(T)): fout.writelines( "\n") for i in range(50) : s = fid.readline() if s.strip('\n').strip('* ') in ['pressures', 'Pressure grid (bars)', 'pressures (bar)']: break fout.writelines(s) if sourceformat.upper() != 'GWB' else fout.writelines(s.strip('* ').split('(')[0] + '\n') for i in range(len(P)): i = i + 1 if (i == 1) | (i == 5): fout.writelines( " %9.4f" % P[i-1]) else: fout.writelines( " %9.4f" % P[i-1]) if (i % 4 == 0) | (i == len(P)): fout.writelines( "\n") #%% Calculation for debye huckel and bdot and water properties waterdielc = water_dielec(T = T, P = P, Dielec_method = Dielec_method) E, Adh, Bdh, bdot = waterdielc.E, waterdielc.Ah, waterdielc.Bh, waterdielc.bdot rhoEG = {'rho': rho, 'E': E, 'dGH2O': dGH2O} # Calculate the rho E G for density extrapolation method here so we have it below rhoEGextrap = {} if any(rhoEG['rho'] < 350): subBornptrs = rhoEG['rho'] < 350 for i, j in enumerate(zip(T[subBornptrs], P[subBornptrs])): rhoextrap = np.linspace(350, 550, 3) Pextrap = iapws95(T = j[0], rho = rhoextrap).P if Dielec_method.upper() != 'DEW' else ZhangDuan(T = j[0], rho = rhoextrap).P Textrap = j[0]*np.ones(np.size(Pextrap)) dGH2O = iapws95(T = Textrap, P = Pextrap).G if Dielec_method.upper() != 'DEW' else ZhangDuan(T = Textrap, P = Pextrap).G E = water_dielec(T = Textrap, P = Pextrap, Dielec_method = Dielec_method).E rhoextrap = np.around(rhoextrap, 3) rhoEGextrap['%d_%d' % (j[0], j[1])]= {'rho': rhoextrap,'E': E, 'dGH2O': dGH2O, 'Textrap': Textrap, 'Pextrap': Pextrap} #skip lines till adh rows if act_param['activity_model'] == 'debye-huckel': for i in range(50) : s = fid.readline() if '(adh)' in s.strip('\n') or 'Debye-Huckel A_gamma' in s.strip('\n'): break fout.writelines(s) if sourceformat.upper() != 'GWB' else fout.writelines(s.strip('* ')) for i in range(len(Adh)): i = i + 1 if (i == 1) | (i == 5): fout.writelines( " %9.4f" % Adh[i-1]) else: fout.writelines( " %9.4f" % Adh[i-1]) if (i % 4 == 0) | (i == len(Adh)): fout.writelines( "\n") #skip lines till bdh rows for i in range(50) : s = fid.readline() if '(bdh)' in s.strip('\n') or 'Debye-Huckel B_gamma' in s.strip('\n'): break fout.writelines(s) if sourceformat.upper() != 'GWB' else fout.writelines(s.strip('* ')) for i in range(len(Bdh)): i = i + 1 if (i == 1) | (i == 5): fout.writelines( " %9.4f" % Bdh[i-1]) else: fout.writelines( " %9.4f" % Bdh[i-1]) if (i % 4 == 0) | (i == len(Bdh)): fout.writelines( "\n") #skip lines till bdot rows for i in range(50) : s = fid.readline() if 'bdot' in s.strip('\n') or 'B-dot' in s.strip('\n'): break fout.writelines(s) if sourceformat.upper() != 'GWB' else fout.writelines(s.strip('* ')) for i in range(len(bdot)): i = i + 1 if (i == 1) | (i == 5): fout.writelines( " %9.4f" % bdot[i-1]) else: fout.writelines( " %9.4f" % bdot[i-1]) if (i % 4 == 0) | (i == len(bdot)): fout.writelines( "\n") #skip lines till cco2 rows for i in range(50) : s = fid.readline() if 'cco2' in s.strip('\n') or 'c co2' in s.strip('\n'): break if sourceformat.upper() != 'GWB': fout.writelines(s) s = fid.readline() fout.writelines(s) s = fid.readline() fout.writelines(s) s = fid.readline() fout.writelines(s) else: fout.writelines('cco2 (coefficients for the Drummond (1981) polynomial) \n' + \ ' -1.0312 0.0012806 \n' + \ ' 255.9 0.4445 \n' + \ ' -0.001606 \n') elif act_param['activity_model'] == 'h-m-w': for i in range(50) : s = fid.readline() if any(re.findall(r'|'.join(('aphi', 'debye huckel')), s.strip('\n'), re.IGNORECASE)): break fout.writelines(s) Aphi = Adh*np.log(10)/3 for i in range(len(Aphi)): i = i + 1 if (i == 1) | (i == 5): fout.writelines( " %9.4f" % Aphi[i-1]) else: fout.writelines( " %9.4f" % Aphi[i-1]) if (i % 4 == 0) | (i == len(Aphi)): fout.writelines( "\n") for i in range(50) : s = fid.readline() if any(re.findall(r'|'.join(('log k for eh reaction', 'Eh reaction: logKr')), s.strip('\n').strip('*'), re.IGNORECASE)): break fout.writelines(s) if sourceformat.upper() != 'GWB' else fout.writelines(s.strip('* ')) #%% Calculations for "log k for eh" rows logK = calcRxnlogK( T = T, P = P, Specie = 'eh', dbaccessdic = dbaccessdic, sourcedic = sourcedic, specielist = specielist, Dielec_method = Dielec_method, rhoEG = rhoEG, sourceformat = sourceformat, densityextrap = densityextrap, rhoEGextrap = rhoEGextrap).logK logK = np.where(np.isnan(logK), 500, logK) # set abitrary 500 to nan values for i in range(len(logK)): i = i + 1 if (i == 1) | (i == 5): fout.writelines( " %9.4f" % logK[i-1]) else: fout.writelines( " %9.4f" % logK[i-1]) if (i % 4 == 0) | (i == len(logK)): fout.writelines( "\n") if sourceformat.upper() != 'GWB': for i in range(50) : s = fid.readline() if s.startswith('+', 0): break if act_param['activity_model'] == 'debye-huckel': if sourceformat.upper() != 'GWB': fout.writelines(s) s = fid.readline() fout.writelines(s) s = fid.readline() fout.writelines(s) s = fid.readline() fout.writelines(s) # copy and paste bdot parameters for i in range(3000) : s = fid.readline() s = s.replace(' acid','_acid').replace(' high','_high').replace(' low','_low') if 'acid' in s else s s_mod = s.split()[0] if s.startswith('+--', 0): break if s_mod in all_species_avail: fout.writelines(s) else: missing_species.append(s_mod) fout.writelines(s) s = fid.readline() fout.writelines(s) s = fid.readline() fout.writelines(s) else: fout.write( "+" + "-"*68 + "\n") fout.write("bdot parameters \n") fout.write( "+" + "-"*68 + "\n") fout.write("* species name azer0 neutral ion type \n") f_ionsize = open(os.path.join(os.path.dirname(os.path.abspath(__file__)), 'ion_size.txt'), 'r') Rd = f_ionsize.readlines() Rd = Rd[1:] Rd = list(OrderedDict.fromkeys(Rd)) Rd_dup = [x for n, x in enumerate(Rd) if x.split()[0] in [l.split()[0] for l in Rd[:n]]] Rd = [x for x in Rd if x not in Rd_dup] Rd = [x if x.split()[0] != 'O2(aq)' else x.replace('O2(aq)', 'O2(g) ') for x in Rd if x.split()[0] in chargedic.keys() ] f_ionsize.close() # for i in range(len(Rd)): fout.writelines(Rd[i]) fout.write( "+" + "-"*68 + "\n") elif act_param['activity_model'] == 'h-m-w': fout.writelines(s) s = fid.readline() fout.writelines(s) s = fid.readline() fout.writelines(s) for i in range(2000): s = fid.readline() if (s.rstrip('\n') == "elements"): break if any(s.lstrip().rstrip('\n').startswith(x) for x in act_param['act_list']): rowlst = []; rowlst.append(s) for i in range(50): s = fid.readline() if s.lstrip().rstrip('\n').startswith("+" + "-"*30): break rowlst.append(s) if all([x in all_species_avail for x in rowlst[0].rstrip('\n').split()]): fout.writelines(rowlst) fout.writelines( "+" + "-"*64 + "\n") else: #break fout.writelines(s) fout.writelines(s) s = fid.readline() fout.writelines(s) # copy and paste elements if sourceformat.upper() != 'GWB': counter = 0 for i in range(2000) : s = fid.readline() s_mod = s.split()[0] if s_mod in elem_avail: fout.writelines(s) counter = counter + 1 if s.startswith('+', 0): break fout.writelines(s) s = fid.readline() fout.writelines(s) s = fid.readline() fout.writelines(s) else: #skip lines till "elements" rows for i in range(5000) : s = fid.readline() if s.rstrip('\n').lstrip('0123456789.- ') == 'elements': break fout.writelines( "%s \n" % s.rstrip('\n').lstrip('0123456789.- ')) fout.writelines( "+" + "-"*68 + "\n") #copy and paste lines till "basis" rows for i in range(1000) : s = fid.readline() s_mod = s.split() if s.strip('\n') != '' else s.strip('\n') if s.rstrip('\n').lstrip('0123456789.- ') == 'basis species': break if len(s_mod) > 1: if s_mod[1].strip('()') in elem_avail: fout.writelines('%-2s %14s \n' % (s_mod[1].strip('()'), s_mod[-1])) fout.writelines( "+" + "-"*68 + "\n") fout.writelines( "%s \n" % s.rstrip('\n').lstrip('0123456789.- ')) fout.writelines( "+" + "-"*68 + "\n") #%% Basis reactions counter = 0 for j in specielist[1]: k = j.replace('(CH3COO)', '(Ac)').replace('CH3COO', '(Ac)') if j not in missing_species: fout.writelines('%s\n' % j.replace('O2(aq)', 'O2(g)')) if j == 'O2(aq)' else fout.writelines('%s\n' % j) if sourceformat.upper() != 'GWB': if j == 'O2(g)': fout.writelines(block_info['O2(g)_b']) else: fout.writelines(block_info[j]) else: filler = re.sub('[^-0123456789\.]', ' ', chargedic[j]).split() fout.writelines(' sp.type = basis \n') if j != 'O2(aq)' else fout.writelines(' sp.type = gas refstate \n') fout.writelines('* EQ3/6 = com, alt, sup, pit\n' +\ ' revised = - \n' +\ '* mol.wt. = %s g/mol\n' % filler[-1] +\ '* DHazero = %s \n' % filler[1] +\ ' charge = %s \n' % filler[0]) fout.writelines('****\n') if sourceformat.upper() != 'GWB': fout.writelines( " %s element(s):\n" % int(len(Elemlist[j])/2)) Rxn = Elemlist[j] else: Rxn = [[v,k] for k,v in dict( sorted(calc_elem_count_molewt(j.rstrip('(aq)(g)') )[0].items(), key=lambda x: x[0].lower()) ).items()] Rxn = [item for sublist in Rxn for item in sublist] fout.writelines( " %s element(s):\n" % int(len(Rxn)/2)) for i in range(len(Rxn)): i = i + 1 if (i == 1) | (i == 7): fout.writelines( " %9.4f " % float(Rxn[i - 1])) elif i % 2 != 0: fout.writelines( "%9.4f " % float(Rxn[i - 1])) else: fout.writelines( "%-9s " % (Rxn[i - 1])) if (i % 6 == 0) | (i == len(Rxn)): fout.writelines( "\n") fout.writelines('****\n') fout.writelines( "* gflag = 1 [reported delG0f used]\n" ) fout.writelines( "* extrapolation algorithm: supcrt92 [92joh/oel]\n" ) if j != 'H2O': ref = dbaccessdic[k][1].split(' ')[0] ref = ref.split(':')[1] if ('ref' in ref) or ('REF' in ref) else ref dG, dH, S = dbaccessdic[k][2], dbaccessdic[k][3], dbaccessdic[k][4] else: dG, dH, S, ref = dGH2O[0], dHH2O[0], SH2O[0], 'iapws95/' + Dielec_method fout.writelines( "* ref-state data [source: %s ]\n" % ref) fout.writelines( "* delG0f = %8.3f kcal/mol\n" % (dG/1000) ) fout.writelines( "* delH0f = %8.3f kcal/mol\n" % (dH/1000)) fout.writelines( "* S0PrTr = %8.3f cal/(mol*K)\n" % S) fout.writelines( "+" + "-"*68 + "\n") counter = counter + 1 else: continue fout.writelines( "auxiliary basis species\n") fout.writelines( "+" + "-"*68 + "\n") #%% Auxiliary Basis reactions specielist[2] = specielist[2] + ['O2(aq)'] if sourceformat.upper() == 'GWB' else specielist[2] for j in specielist[2]: k = j.replace('(CH3COO)', '(Ac)').replace('CH3COO', '(Ac)') rxnlst = [b for a, b in enumerate(sourcedic[j]) if a not in [0, 1, 2, 3]] rxnlst = [v for x, v in enumerate(rxnlst) if x % 2 != 0 ] # remove all coefficients if (j not in missing_species) and (len([i for i in rxnlst if i not in missing_species]) == len(rxnlst)): if sourcedic[j][0] != '': fout.writelines('%-25s %s \n' % (j, sourcedic[j][0])) else: fout.writelines('%-25s %s \n' % (j, dbaccessdic[k][0])) if sourceformat.upper() != 'GWB': fout.writelines(block_info[j]) else: filler = re.sub('[^-0123456789\.]', ' ', chargedic[j]).split() fout.writelines(' sp.type = aux\n' +\ '* EQ3/6 = com, alt, sup \n' +\ ' revised = - \n' +\ '* mol.wt. = %s g/mol\n' % filler[-1] +\ '* DHazero = %s \n' % filler[1] +\ ' charge = %s \n' % filler[0]) fout.writelines('****\n') if sourceformat.upper() != 'GWB': fout.writelines( " %s element(s):\n" % int(len(Elemlist[j])/2)) Elem = Elemlist[j] else: # print(j) formula = j.rstrip('(aq)(g)') if sourcedic[j][0] == '' else j if not j.endswith('(aq)') else dbaccessdic[k][0].rstrip('(aq)(g)') #.rstrip('(+0123456789)') Elem = [[v,k] for k,v in dict( sorted(calc_elem_count_molewt(formula)[0].items(), key=lambda x: x[0].lower()) ).items()] Elem = [item for sublist in Elem for item in sublist] fout.writelines( " %s element(s):\n" % int(len(Elem)/2)) for i in range(len(Elem)): i = i + 1 if (i == 1) | (i == 7): fout.writelines( " %9.4f " % float(Elem[i - 1])) elif i % 2 != 0: fout.writelines( "%9.4f " % float(Elem[i - 1])) else: fout.writelines( "%-9s " % (Elem[i - 1])) if (i % 6 == 0) | (i == len(Elem)): fout.writelines( "\n") fout.writelines('****\n') if sourceformat.upper() != 'GWB': fout.writelines( " %s species in aqueous dissociation reaction:\n" % sourcedic[j][1]) Rxn = sourcedic[j][2:] else: fout.writelines( " %s species in aqueous dissociation reaction:\n" % (sourcedic[j][1] + 1)) sourcedic[j] = [k if k != 'O2(aq)' else k.replace('O2(aq)', 'O2(g)') for k in sourcedic[j]] Rxn = ['-1.0000', 'O2(aq)', '1.0000', 'O2(g)'] if sourceformat.upper() == 'GWB' and j == 'O2(aq)' else ['-1.0000', '%s' % j] + sourcedic[j][2:] for i in range(len(Rxn)): i = i + 1 if (i == 1) | (i == 5) | (i == 9): fout.writelines( " %9.4f " % float(Rxn[i - 1])) elif i % 2 != 0: fout.writelines( "%9.4f " % float(Rxn[i - 1])) else: fout.writelines( " %-21s " % (Rxn[i - 1])) if (i % 4 == 0) | (i == len(Rxn)): fout.writelines( "\n") fout.writelines('*\n') fout.writelines('**** logK grid [T, P @ Miscellaneous parameters]\n') sourcedic[j] = ['', 1, '1.0000', 'O2(g)'] if sourceformat.upper() == 'GWB' and j == 'O2(aq)' else sourcedic[j] logK = calcRxnlogK( T = T, P = P, Specie = j, dbaccessdic = dbaccessdic, sourcedic = sourcedic, specielist = specielist, Dielec_method = Dielec_method, rhoEG = rhoEG, sourceformat = sourceformat, densityextrap = densityextrap, Specie_class = 'aqueous', rhoEGextrap = rhoEGextrap, ThermoInUnit = self.ThermoInUnit) if densityextrap.lower() == 'yes': if all(logK.nonsubBornptrs) == True: # if all densities are >= 350 logKnan_alert = False # turn off the prompts for using Density extrapolation else: logKnan_alert = True else: logKnan_alert = False # turn off the prompts for using Density extrapolation logK = logK.logK logK = np.where(np.isnan(logK), 500, logK) # set abitrary 500 to nan values for i in range(len(logK)): i = i + 1 if (i == 1) | (i == 5): fout.writelines( " %9.4f" % logK[i-1]) else: fout.writelines( " %9.4f" % logK[i-1]) if (i % 4 == 0) | (i == len(logK)): fout.writelines( "\n") fout.writelines( "* gflag = 1 [reported delG0f used]\n" ) fout.writelines( "* extrapolation algorithm: supcrt92 [92joh/oel]\n" ) if j != 'H2O': ref = dbaccessdic[k][1].split(' ')[0] ref = ref.split(':')[1] if ('ref' in ref) or ('REF' in ref) else ref dG, dH, S = dbaccessdic[k][2], dbaccessdic[k][3], dbaccessdic[k][4] else: dG, dH, S, ref = dGH2O[0], dHH2O[0], SH2O[0], 'iapws95/' + Dielec_method fout.writelines( "* ref-state data [source: %s ]\n" % ref) fout.writelines( "* delG0f = %8.3f kcal/mol\n" % (dG/1000) ) fout.writelines( "* delH0f = %8.3f kcal/mol\n" % (dH/1000)) fout.writelines( "* S0PrTr = %8.3f cal/(mol*K)\n" % S) fout.writelines( "+" + "-"*68 + "\n") else: continue fout.writelines( "aqueous species\n") fout.writelines( "+" + "-"*68 + "\n") if logKnan_alert == True: warnings.warn('Some temperature and pressure points are out of aqueous species HKF eqns regions of applicability, hence, density extrapolation has been applied') #%% Aqueous reactions for j in specielist[3]: k = j.replace('(CH3COO)', '(Ac)').replace('CH3COO', '(Ac)') rxnlst = [b for a, b in enumerate(sourcedic[j]) if a not in [0, 1, 2, 3]] rxnlst = [v for x, v in enumerate(rxnlst) if x % 2 != 0 ] if (j not in missing_species) and (len([i for i in rxnlst if i not in missing_species]) == len(rxnlst)): if sourcedic[j][0] != '': fout.writelines('%-25s %s \n' % (j, sourcedic[j][0])) else: fout.writelines('%-25s %s \n' % (j, dbaccessdic[k][0])) if sourceformat.upper() != 'GWB': fout.writelines(block_info[j]) else: filler = re.sub('[^-0123456789\.]', ' ', chargedic[j]).split() fout.writelines(' sp.type = aqueous\n' +\ '* EQ3/6 = com, alt, sup \n' +\ ' revised = - \n' +\ '* mol.wt. = %s g/mol\n' % filler[-1] +\ '* DHazero = %s \n' % filler[1] +\ ' charge = %s \n' % filler[0]) fout.writelines('****\n') if sourceformat.upper() != 'GWB': fout.writelines( " %s element(s):\n" % int(len(Elemlist[j])/2)) Elem = Elemlist[j] else: # print(j) filler = ('(aq)', '(But)', '(Prop)', '(Pent)', '(For)', '(Gly)', '(Glyc)', '(Lac)', 'Acetate') formula = j.rstrip('(aq)(g)') if sourcedic[j][0] == '' else j if not any([l in j for l in filler]) else re.sub('(aq)', '', dbaccessdic[k][0]) #.rstrip('(+0123456789)') Elem = [[v,k] for k,v in dict( sorted(calc_elem_count_molewt(formula)[0].items(), key=lambda x: x[0].lower()) ).items()] Elem = [item for sublist in Elem for item in sublist] fout.writelines( " %s element(s):\n" % int(len(Elem)/2)) for i in range(len(Elem)): i = i + 1 if (i == 1) | (i == 7): fout.writelines( " %9.4f " % float(Elem[i - 1])) elif i % 2 != 0: fout.writelines( "%9.4f " % float(Elem[i - 1])) else: fout.writelines( "%-9s " % (Elem[i - 1])) if (i % 6 == 0) | (i == len(Elem)): fout.writelines( "\n") fout.writelines('****\n') if sourceformat.upper() != 'GWB': fout.writelines( " %s species in aqueous dissociation reaction:\n" % sourcedic[j][1]) Rxn = sourcedic[j][2:] else: fout.writelines( " %s species in aqueous dissociation reaction:\n" % (sourcedic[j][1] + 1)) sourcedic[j] = [k if k != 'O2(aq)' else k.replace('O2(aq)', 'O2(g)') for k in sourcedic[j]] Rxn = ['-1.0000', '%s' % j] + sourcedic[j][2:] for i in range(len(Rxn)): i = i + 1 if (i == 1) | (i == 5) | (i == 9): fout.writelines( " %9.4f " % float(Rxn[i - 1])) elif i % 2 != 0: fout.writelines( "%9.4f " % float(Rxn[i - 1])) else: fout.writelines( " %-21s " % (Rxn[i - 1])) if (i % 4 == 0) | (i == len(Rxn)): fout.writelines( "\n") fout.writelines('*\n') fout.writelines('**** logK grid [T, P @ Miscellaneous parameters]\n') logK = calcRxnlogK( T = T, P = P, Specie = j, dbaccessdic = dbaccessdic, sourcedic = sourcedic, specielist = specielist, Dielec_method = Dielec_method, rhoEG = rhoEG, sourceformat = sourceformat, densityextrap = densityextrap, Specie_class = 'aqueous', rhoEGextrap = rhoEGextrap, ThermoInUnit = self.ThermoInUnit).logK logK = np.where(np.isnan(logK), 500, logK) # set abitrary 500 to nan values for i in range(len(logK)): i = i + 1 if (i == 1) | (i == 5): fout.writelines( " %9.4f" % logK[i-1]) else: fout.writelines( " %9.4f" % logK[i-1]) if (i % 4 == 0) | (i == len(logK)): fout.writelines( "\n") fout.writelines( "* gflag = 1 [reported delG0f used]\n" ) fout.writelines( "* extrapolation algorithm: supcrt92 [92joh/oel]\n" ) if j != 'H2O': ref = dbaccessdic[k][1].split(' ')[0] ref = ref.split(':')[1] if ('ref' in ref) or ('REF' in ref) else ref dG, dH, S = dbaccessdic[k][2], dbaccessdic[k][3], dbaccessdic[k][4] else: dG, dH, S, ref = dGH2O[0], dHH2O[0], SH2O[0], 'iapws95/' + Dielec_method fout.writelines( "* ref-state data [source: %s ]\n" % ref) fout.writelines( "* delG0f = %8.3f kcal/mol\n" % (dG/1000) ) fout.writelines( "* delH0f = %8.3f kcal/mol\n" % (dH/1000)) fout.writelines( "* S0PrTr = %8.3f cal/(mol*K)\n" % S) fout.writelines( "+" + "-"*68 + "\n") else: continue fout.writelines( "solids\n") fout.writelines( "+" + "-"*68 + "\n") #%% Mineral reactions solid_solution = 'no' if solid_solution is None else solid_solution clay_thermo = 'no' if clay_thermo is None else clay_thermo if solid_solution.lower() == 'yes': if nCa_cpx is None: nCa = 0 else: nCa = nCa_cpx solidsolution_no = 11 fnlist = ['plagio', 'olivine', 'pyroxene', 'cpx'] if nCa > 0 else ['plagio', 'olivine', 'pyroxene'] for fn in fnlist: for nX in np.round(np.linspace(1, 0, solidsolution_no), 1): if fn != 'cpx': ss = calcRxnlogK(X = nX, T = T, P = P, Dielec_method = Dielec_method, rhoEG = rhoEG, dbaccessdic = dbaccessdic, Specie = fn, densityextrap = densityextrap, ThermoInUnit = self.ThermoInUnit, rhoEGextrap = rhoEGextrap) else: ss = calcRxnlogK(cpx_Ca = nCa, X = nX, T = T, P = P, Dielec_method = Dielec_method, rhoEG = rhoEG, dbaccessdic = dbaccessdic, Specie = fn, densityextrap = densityextrap, ThermoInUnit = self.ThermoInUnit, rhoEGextrap = rhoEGextrap) logK, Rxn = ss.logK, ss.Rxn logK = np.where(np.isnan(logK), 500, logK) # set abitrary 500 to nan values outputfmt(fout, logK, Rxn, dataset = 'EQ36') # clay minerals if clay_thermo.lower() == 'yes': fclay = open(os.path.join(os.path.dirname(os.path.abspath(__file__)), 'clay_elements.dat'), 'r') Rd = fclay.readlines() Rd = [j.replace('-','_').strip('\n') for j in Rd] for i in range(len(Rd)): ss = calcRxnlogK(T = T, P = P, Specie = 'Clay', elem = Rd[i].split(','), dbaccessdic = dbaccessdic, ThermoInUnit = self.ThermoInUnit, rhoEG = rhoEG, rhoEGextrap = rhoEGextrap, densityextrap = densityextrap) logK, Rxn = ss.logK, ss.Rxn logK = np.where(np.isnan(logK), 500, logK) # set abitrary 500 to nan values outputfmt(fout, logK, Rxn, dataset = 'EQ36') # other minerals in the source database minlst = specielist[4] if sourceformat.upper() != 'GWB' else specielist[5] for j in minlst: k = j.replace('(CH3COO)', '(Ac)').replace('CH3COO', '(Ac)') rxnlst = [b for a, b in enumerate(sourcedic[j]) if a not in [0, 1, 2, 3]] rxnlst = [v for x, v in enumerate(rxnlst) if x % 2 != 0 ] if (j not in missing_species) and (len([k for k in rxnlst if k not in missing_species]) == len(rxnlst)): if solid_solution.lower() == 'yes' and j in ['Anorthite', 'Albite', 'Forsterite', 'Fayalite', 'Enstatite', 'Ferrosilite']: continue elif solid_solution.lower() == 'yes' and (nCa == 1) and j in ['Diopside', 'Hedenbergite']: continue elif clay_thermo.lower() == 'yes' and j in [ Rd[h].split(',')[0] for h in range(len(Rd))]: continue else: if sourcedic[j][0] != '': fout.writelines('%-25s %s \n' % (j, sourcedic[j][0])) else: fout.writelines('%-25s %s \n' % (j, dbaccessdic[k][0])) if sourceformat.upper() != 'GWB': fout.writelines(block_info[j]) else: fout.writelines(' sp.type = solid \n' +\ '* EQ3/6 = com, alt, sup \n' +\ ' revised = - \n' +\ '* mol.wt. = %s g/mol\n' % MWdic[j] +\ ' V0PrTr = %s cm**3/mol \n' % dbaccessdic[j][5] ) fout.writelines('****\n') if sourceformat.upper() != 'GWB': fout.writelines( " %s element(s):\n" % int(len(Elemlist[j])/2)) Elem = Elemlist[j] else: formula = re.sub('(s)', '', j) if sourcedic[j][0] == '' else re.sub('(s)', '', dbaccessdic[k][0]) Elem = [[v,k] for k,v in dict( sorted(calc_elem_count_molewt(formula)[0].items(), key=lambda x: x[0].lower()) ).items()] Elem = [item for sublist in Elem for item in sublist] fout.writelines( " %s element(s):\n" % int(len(Elem)/2)) for i in range(len(Elem)): i = i + 1 if (i == 1) | (i == 7) | (i == 13): fout.writelines( " %9.4f " % float(Elem[i - 1])) elif i % 2 != 0: fout.writelines( "%9.4f " % float(Elem[i - 1])) else: fout.writelines( "%-9s " % (Elem[i - 1])) if (i % 6 == 0) | (i == len(Elem)): fout.writelines( "\n") fout.writelines('****\n') if sourceformat.upper() != 'GWB': fout.writelines( " %s species in reaction:\n" % sourcedic[j][1]) Rxn = sourcedic[j][2:] else: fout.writelines( " %s species in reaction:\n" % (sourcedic[j][1] + 1)) sourcedic[j] = [k if k != 'O2(aq)' else k.replace('O2(aq)', 'O2(g)') for k in sourcedic[j]] Rxn = ['-1.0000', '%s' % j] + sourcedic[j][2:] for i in range(len(Rxn)): i = i + 1 if (i == 1) | (i == 5) | (i == 9) | (i == 13) | (i == 17) | (i == 21): fout.writelines( " %9.4f " % float(Rxn[i - 1])) elif i % 2 != 0: fout.writelines( "%9.4f " % float(Rxn[i - 1])) else: fout.writelines( " %-21s " % (Rxn[i - 1])) if (i % 4 == 0) | (i == len(Rxn)): fout.writelines( "\n") fout.writelines('*\n') fout.writelines('**** logK grid [T, P @ Miscellaneous parameters]\n') logK = calcRxnlogK( T = T, P = P, Specie = j, dbaccessdic = dbaccessdic, sourcedic = sourcedic, specielist = specielist, Dielec_method = Dielec_method, rhoEG = rhoEG, sourceformat = sourceformat, densityextrap = densityextrap, Specie_class = 'minerals', heatcap_method = heatcap_method, rhoEGextrap = rhoEGextrap).logK logK = np.where(np.isnan(logK), 500, logK) # set abitrary 500 to nan values for i in range(len(logK)): i = i + 1 if (i == 1) | (i == 5): fout.writelines( " %9.4f" % logK[i-1]) else: fout.writelines( " %9.4f" % logK[i-1]) if (i % 4 == 0) | (i == len(logK)): fout.writelines( "\n") fout.writelines( "* gflag = 1 [reported delG0f used]\n" ) fout.writelines( "* extrapolation algorithm: supcrt92 [92joh/oel]\n" ) if j != 'H2O': dG, dH, S = dbaccessdic[k][2], dbaccessdic[k][3], dbaccessdic[k][4] else: dG, dH, S, ref = dGH2O[0], dHH2O[0], SH2O[0], 'iapws95/' + Dielec_method fout.writelines( "* ref-state data [source: %s ]\n" % ref) fout.writelines( "* delG0f = %8.3f kcal/mol\n" % (dG/1000) ) fout.writelines( "* delH0f = %8.3f kcal/mol\n" % (dH/1000)) fout.writelines( "* S0PrTr = %8.3f cal/(mol*K)\n" % S) fout.writelines( "* Cp coefficients [source: %s ]\n" % ref) fout.writelines( "* T**0 = %11.8e \n" % (dbaccessdic[k][6]) ) fout.writelines( "* T**1 = %11.8e \n" % (dbaccessdic[k][7]*10**-3)) if dbaccessdic[j][8] < 1: fout.writelines( "* T**-2 = %12.8e \n" % (dbaccessdic[k][8]*10**5)) else: fout.writelines( "* T**-2 = %11.8e \n" % (dbaccessdic[k][8]*10**5)) fout.writelines( "* Tlimit = %7.2fC \n" % (dbaccessdic[k][9])) fout.writelines( "+" + "-"*68 + "\n") else: continue fout.writelines( "liquids\n") fout.writelines( "+" + "-"*68 + "\n") #%% Liquids reactions if sourceformat.upper() != 'GWB': for j in specielist