Module pychnosz.models
Equation of state models and water property models for CHNOSZ.
Sub-modules
pychnosz.models.archer_wang-
Archer & Wang (1990) dielectric constant correlation for water …
pychnosz.models.berman-
Berman mineral equations of state implementation …
pychnosz.models.dew-
DEW (Deep Earth Water) model implementation …
pychnosz.models.hkf_helpers-
Helper functions for HKF calculations, integrated from HKF_cgl.py …
pychnosz.models.iapws95-
IAPWS-95 water model implementation …
pychnosz.models.supcrt92_fortran-
SUPCRT92 water model with Fortran backend …
Functions
def Berman(name: str,
T: float | List[float] = 298.15,
P: float | List[float] = 1,
check_G: bool = False,
calc_transition: bool = True,
calc_disorder: bool = True) ‑> pandas.core.frame.DataFrame-
Expand source code
def Berman(name: str, T: Union[float, List[float]] = 298.15, P: Union[float, List[float]] = 1, check_G: bool = False, calc_transition: bool = True, calc_disorder: bool = True) -> pd.DataFrame: """ Calculate thermodynamic properties of minerals using Berman equations. Parameters ---------- name : str Name of the mineral T : float or list, optional Temperature in Kelvin (default: 298.15) P : float or list, optional Pressure in bar (default: 1) check_G : bool, optional Check consistency of G in data file (default: False) calc_transition : bool, optional Calculate polymorphic transition contributions (default: True) calc_disorder : bool, optional Calculate disorder contributions (default: True) Returns ------- pd.DataFrame DataFrame with columns T, P, G, H, S, Cp, V """ # Reference temperature and pressure Pr = 1 Tr = 298.15 # Make T and P the same length if isinstance(T, (int, float)): T = [T] if isinstance(P, (int, float)): P = [P] # Convert to list if numpy array (to avoid element-wise multiplication bug) if isinstance(T, np.ndarray): T = T.tolist() if isinstance(P, np.ndarray): P = P.tolist() ncond = max(len(T), len(P)) T = np.array(T * (ncond // len(T) + 1), dtype=float)[:ncond] P = np.array(P * (ncond // len(P) + 1), dtype=float)[:ncond] # Get parameters in the Berman equations # Start with thermodynamic parameters provided with CHNOSZ thermo_sys = thermo() if thermo_sys.Berman is None: raise RuntimeError("Berman data not loaded. Please run pychnosz.reset() first.") dat = thermo_sys.Berman.copy() # TODO: Handle user-supplied data file (thermo()$opt$Berman) # For now, just use the default data # Remove duplicates (only the first, i.e. most recent entry is kept) dat = dat.drop_duplicates(subset=['name'], keep='first') # Remove the multipliers on volume parameters vcols = ['v1', 'v2', 'v3', 'v4'] # columns with v1, v2, v3, v4 multexp = [5, 5, 5, 8] for i, col in enumerate(vcols): if col in dat.columns: dat[col] = dat[col] / (10 ** multexp[i]) # Which row has data for this mineral? matching_rows = dat[dat['name'] == name] if len(matching_rows) == 0: raise ValueError(f"Data for {name} not available in Berman database") dat_mineral = matching_rows.iloc[0] # Extract parameters for easier access GfPrTr = dat_mineral['GfPrTr'] HfPrTr = dat_mineral['HfPrTr'] SPrTr = dat_mineral['SPrTr'] VPrTr = dat_mineral['VPrTr'] k0 = dat_mineral['k0'] k1 = dat_mineral['k1'] k2 = dat_mineral['k2'] k3 = dat_mineral['k3'] k4 = dat_mineral['k4'] if not pd.isna(dat_mineral['k4']) else 0 k5 = dat_mineral['k5'] if not pd.isna(dat_mineral['k5']) else 0 k6 = dat_mineral['k6'] if not pd.isna(dat_mineral['k6']) else 0 v1 = dat_mineral['v1'] if not pd.isna(dat_mineral['v1']) else 0 v2 = dat_mineral['v2'] if not pd.isna(dat_mineral['v2']) else 0 v3 = dat_mineral['v3'] if not pd.isna(dat_mineral['v3']) else 0 v4 = dat_mineral['v4'] if not pd.isna(dat_mineral['v4']) else 0 # Transition parameters Tlambda = dat_mineral['Tlambda'] if not pd.isna(dat_mineral['Tlambda']) else None Tref = dat_mineral['Tref'] if not pd.isna(dat_mineral['Tref']) else None dTdP = dat_mineral['dTdP'] if not pd.isna(dat_mineral['dTdP']) else None l1 = dat_mineral['l1'] if not pd.isna(dat_mineral['l1']) else None l2 = dat_mineral['l2'] if not pd.isna(dat_mineral['l2']) else None # Disorder parameters Tmin = dat_mineral['Tmin'] if not pd.isna(dat_mineral['Tmin']) else None Tmax = dat_mineral['Tmax'] if not pd.isna(dat_mineral['Tmax']) else None d0 = dat_mineral['d0'] if not pd.isna(dat_mineral['d0']) else None d1 = dat_mineral['d1'] if not pd.isna(dat_mineral['d1']) else None d2 = dat_mineral['d2'] if not pd.isna(dat_mineral['d2']) else None d3 = dat_mineral['d3'] if not pd.isna(dat_mineral['d3']) else None d4 = dat_mineral['d4'] if not pd.isna(dat_mineral['d4']) else None Vad = dat_mineral['Vad'] if not pd.isna(dat_mineral['Vad']) else None # Get the entropy of the elements using the chemical formula # Get formula from OBIGT and calculate using entropy() function like in R CHNOSZ SPrTr_elements = 0 if thermo_sys.obigt is not None: obigt_match = thermo_sys.obigt[thermo_sys.obigt['name'] == name] if len(obigt_match) > 0: formula = obigt_match.iloc[0]['formula'] # Import entropy function and calculate SPrTr_elements properly from ..utils.formula import entropy SPrTr_elements = entropy(formula) # Check that G in data file follows Benson-Helgeson convention if check_G and not pd.isna(GfPrTr): GfPrTr_calc = HfPrTr - Tr * (SPrTr - SPrTr_elements) Gdiff = GfPrTr_calc - GfPrTr if abs(Gdiff) >= 1000: warnings.warn(f"{name}: GfPrTr(calc) - GfPrTr(table) is too big! == {round(Gdiff)} J/mol") ### Thermodynamic properties ### # Calculate Cp and V (Berman, 1988 Eqs. 4 and 5) # k4, k5, k6 terms from winTWQ documentation (doi:10.4095/223425) Cp = k0 + k1 * T**(-0.5) + k2 * T**(-2) + k3 * T**(-3) + k4 * T**(-1) + k5 * T + k6 * T**2 P_Pr = P - Pr T_Tr = T - Tr V = VPrTr * (1 + v1 * T_Tr + v2 * T_Tr**2 + v3 * P_Pr + v4 * P_Pr**2) # Calculate Ha (symbolically integrated using sympy - expressions not simplified) intCp = (T*k0 - Tr*k0 + k2/Tr - k2/T + k3/(2*Tr**2) - k3/(2*T**2) + 2.0*k1*T**0.5 - 2.0*k1*Tr**0.5 + k4*np.log(T) - k4*np.log(Tr) + k5*T**2/2 - k5*Tr**2/2 - k6*Tr**3/3 + k6*T**3/3) intVminusTdVdT = (-VPrTr + P*(VPrTr + VPrTr*v4 - VPrTr*v3 - Tr*VPrTr*v1 + VPrTr*v2*Tr**2 - VPrTr*v2*T**2) + P**2*(VPrTr*v3/2 - VPrTr*v4) + VPrTr*v3/2 - VPrTr*v4/3 + Tr*VPrTr*v1 + VPrTr*v2*T**2 - VPrTr*v2*Tr**2 + VPrTr*v4*P**3/3) Ha = HfPrTr + intCp + intVminusTdVdT # Calculate S (also symbolically integrated) intCpoverT = (k0*np.log(T) - k0*np.log(Tr) - k3/(3*T**3) + k3/(3*Tr**3) + k2/(2*Tr**2) - k2/(2*T**2) + 2.0*k1*Tr**(-0.5) - 2.0*k1*T**(-0.5) + k4/Tr - k4/T + T*k5 - Tr*k5 + k6*T**2/2 - k6*Tr**2/2) intdVdT = -VPrTr*(v1 + v2*(-2*Tr + 2*T)) + P*VPrTr*(v1 + v2*(-2*Tr + 2*T)) S = SPrTr + intCpoverT - intdVdT # Calculate Ga --> Berman-Brown convention (DG = DH - T*S, no S(element)) Ga = Ha - T * S ### Polymorphic transition properties ### if (Tlambda is not None and Tref is not None and not pd.isna(Tlambda) and not pd.isna(Tref) and np.any(T > Tref) and calc_transition): # Starting transition contributions are 0 Cptr = np.zeros(ncond) Htr = np.zeros(ncond) Str = np.zeros(ncond) # Eq. 9: Tlambda at P Tlambda_P = Tlambda + dTdP * (P - 1) # Eq. 8a: Cp at P Td = Tlambda - Tlambda_P Tprime = T + Td # With the condition that Tref < Tprime < Tlambda(1bar) iTprime = (Tref < Tprime) & (Tprime < Tlambda) # Handle NA values iTprime = iTprime & ~np.isnan(Tprime) if np.any(iTprime): Tprime_valid = Tprime[iTprime] Cptr[iTprime] = Tprime_valid * (l1 + l2 * Tprime_valid)**2 # We got Cp, now calculate the integrations for H and S iTtr = T > Tref if np.any(iTtr): Ttr = T[iTtr].copy() Tlambda_P_tr = Tlambda_P[iTtr].copy() Td_tr = Td[iTtr] if hasattr(Td, '__len__') else np.full_like(Ttr, Td) # Handle NA values Tlambda_P_tr[np.isnan(Tlambda_P_tr)] = np.inf # The upper integration limit is Tlambda_P Ttr[Ttr >= Tlambda_P_tr] = Tlambda_P_tr[Ttr >= Tlambda_P_tr] # Derived variables tref = Tref - Td_tr x1 = l1**2 * Td_tr + 2 * l1 * l2 * Td_tr**2 + l2**2 * Td_tr**3 x2 = l1**2 + 4 * l1 * l2 * Td_tr + 3 * l2**2 * Td_tr**2 x3 = 2 * l1 * l2 + 3 * l2**2 * Td_tr x4 = l2**2 # Eqs. 10, 11, 12 Htr[iTtr] = (x1 * (Ttr - tref) + x2/2 * (Ttr**2 - tref**2) + x3/3 * (Ttr**3 - tref**3) + x4/4 * (Ttr**4 - tref**4)) Str[iTtr] = (x1 * (np.log(Ttr) - np.log(tref)) + x2 * (Ttr - tref) + x3/2 * (Ttr**2 - tref**2) + x4/3 * (Ttr**3 - tref**3)) Gtr = Htr - T * Str # Apply the transition contributions Ga = Ga + Gtr Ha = Ha + Htr S = S + Str Cp = Cp + Cptr ### Disorder thermodynamic properties ### if (Tmin is not None and Tmax is not None and not pd.isna(Tmin) and not pd.isna(Tmax) and np.any(T > Tmin) and calc_disorder): # Starting disorder contributions are 0 Cpds = np.zeros(ncond) Hds = np.zeros(ncond) Sds = np.zeros(ncond) Vds = np.zeros(ncond) # The lower integration limit is Tmin iTds = T > Tmin if np.any(iTds): Tds = T[iTds].copy() # The upper integration limit is Tmax Tds[Tds > Tmax] = Tmax # Ber88 Eqs. 15, 16, 17 Cpds[iTds] = d0 + d1*Tds**(-0.5) + d2*Tds**(-2) + d3*Tds + d4*Tds**2 Hds[iTds] = (d0*(Tds - Tmin) + d1*(Tds**0.5 - Tmin**0.5)/0.5 + d2*(Tds**(-1) - Tmin**(-1))/(-1) + d3*(Tds**2 - Tmin**2)/2 + d4*(Tds**3 - Tmin**3)/3) Sds[iTds] = (d0*(np.log(Tds) - np.log(Tmin)) + d1*(Tds**(-0.5) - Tmin**(-0.5))/(-0.5) + d2*(Tds**(-2) - Tmin**(-2))/(-2) + d3*(Tds - Tmin) + d4*(Tds**2 - Tmin**2)/2) # Eq. 18; we can't do this if Vad == 0 (dolomite and gehlenite) if Vad is not None and not pd.isna(Vad) and Vad != 0: Vds = Hds / Vad # Include the Vds term with Hds Hds = Hds + Vds * (P - Pr) # Disordering properties above Tmax (Eq. 20) ihigh = T > Tmax if np.any(ihigh): Hds[ihigh] = Hds[ihigh] - (T[ihigh] - Tmax) * Sds[ihigh] Gds = Hds - T * Sds # Apply the disorder contributions Ga = Ga + Gds Ha = Ha + Hds S = S + Sds V = V + Vds Cp = Cp + Cpds ### (for testing) Use G = H - TS to check that integrals for H and S are written correctly Ga_fromHminusTS = Ha - T * S if not np.allclose(Ga_fromHminusTS, Ga, atol=1e-6): raise RuntimeError(f"{name}: incorrect integrals detected using DG = DH - T*S") ### Thermodynamic and unit conventions used in SUPCRT ### # Use entropy of the elements in calculation of G --> Benson-Helgeson convention (DG = DH - T*DS) Gf = Ga + Tr * SPrTr_elements # The output will just have "G" and "H" G = Gf H = Ha # Convert J/bar to cm^3/mol V = V * 10 return pd.DataFrame({ 'T': T, 'P': P, 'G': G, 'H': H, 'S': S, 'Cp': Cp, 'V': V })Calculate thermodynamic properties of minerals using Berman equations.
Parameters
name:str- Name of the mineral
T:floatorlist, optional- Temperature in Kelvin (default: 298.15)
P:floatorlist, optional- Pressure in bar (default: 1)
check_G:bool, optional- Check consistency of G in data file (default: False)
calc_transition:bool, optional- Calculate polymorphic transition contributions (default: True)
calc_disorder:bool, optional- Calculate disorder contributions (default: True)
Returns
pd.DataFrame- DataFrame with columns T, P, G, H, S, Cp, V
def G2logK(G, Tc)-
Expand source code
def G2logK(G, Tc): # Gas constant R is in cal/mol K return G / (-math.log(10) * 1.9872 * (273.15+Tc)) def OBIGT2eos(OBIGT, fixGHS=True, tocal=True, messages=True)-
Expand source code
def OBIGT2eos(OBIGT, fixGHS=True, tocal=True, messages=True): """ Convert OBIGT dataframe to equation of state parameters. This function processes the OBIGT thermodynamic database to prepare it for equation of state calculations. It handles energy unit conversions and optionally fills in missing G, H, or S values. Parameters ---------- OBIGT : pd.DataFrame OBIGT thermodynamic database fixGHS : bool, default True Fill in one missing value among G, H, S using thermodynamic relations tocal : bool, default True Convert energy units from Joules to calories messages : bool, default True Print informational messages (currently not used, reserved for future) Returns ------- pd.DataFrame Modified OBIGT dataframe with converted parameters """ OBIGT_out = OBIGT.copy() # Get column indices for named columns (to handle varying column positions) G_idx = OBIGT.columns.get_loc('G') H_idx = OBIGT.columns.get_loc('H') S_idx = OBIGT.columns.get_loc('S') Cp_idx = OBIGT.columns.get_loc('Cp') V_idx = OBIGT.columns.get_loc('V') omega_lambda_idx = OBIGT.columns.get_loc('omega.lambda') for i in range(0, OBIGT.shape[0]): # we only convert omega for aqueous species, not lambda for cgl species if tocal and OBIGT.iloc[i, :]["E_units"] == "J" and OBIGT.iloc[i, :]["state"] == "aq": # Convert G, H, S, Cp OBIGT_out.iloc[i, G_idx:Cp_idx+1] = OBIGT.iloc[i, G_idx:Cp_idx+1]/4.184 # Convert V through omega (includes omega for aq species) OBIGT_out.iloc[i, V_idx:omega_lambda_idx+1] = OBIGT.iloc[i, V_idx:omega_lambda_idx+1]/4.184 OBIGT_out.iloc[i, OBIGT.columns.get_loc('E_units')] = "cal" elif tocal and OBIGT.iloc[i, :]["E_units"] == "J": # Convert G, H, S, Cp OBIGT_out.iloc[i, G_idx:Cp_idx+1] = OBIGT.iloc[i, G_idx:Cp_idx+1]/4.184 # Convert V through c2.f (exclude omega.lambda for non-aq species) OBIGT_out.iloc[i, V_idx:omega_lambda_idx] = OBIGT.iloc[i, V_idx:omega_lambda_idx]/4.184 OBIGT_out.iloc[i, OBIGT.columns.get_loc('E_units')] = "cal" # fill in one of missing G, H, S # for use esp. by subcrt because NA for one of G, H or S # will preclude calculations at high T if fixGHS: # which entries are missing just one GHS_values = [OBIGT.iloc[i, G_idx], OBIGT.iloc[i, H_idx], OBIGT.iloc[i, S_idx]] imiss = [pd.isna(v) for v in GHS_values] if sum(imiss) == 1: ii = imiss.index(True) if ii == 0: # G is missing H = OBIGT_out.iloc[i, H_idx] S = OBIGT_out.iloc[i, S_idx] Selem = entropy(OBIGT_out.iloc[i, OBIGT_out.columns.get_loc('formula')]) T = 298.15 G = H - T*(S - Selem) OBIGT_out.iloc[i, G_idx] = G elif ii == 1: # H is missing G = OBIGT_out.iloc[i, G_idx] S = OBIGT_out.iloc[i, S_idx] Selem = entropy(OBIGT_out.iloc[i, OBIGT_out.columns.get_loc('formula')]) T = 298.15 H = G + T*(S - Selem) OBIGT_out.iloc[i, H_idx] = H elif ii == 2: # S is missing G = OBIGT_out.iloc[i, G_idx] H = OBIGT_out.iloc[i, H_idx] Selem = entropy(OBIGT_out.iloc[i, OBIGT_out.columns.get_loc('formula')]) T = 298.15 S = Selem + (H - G)/T OBIGT_out.iloc[i, S_idx] = S return OBIGT_outConvert OBIGT dataframe to equation of state parameters.
This function processes the OBIGT thermodynamic database to prepare it for equation of state calculations. It handles energy unit conversions and optionally fills in missing G, H, or S values.
Parameters
OBIGT:pd.DataFrame- OBIGT thermodynamic database
fixGHS:bool, defaultTrue- Fill in one missing value among G, H, S using thermodynamic relations
tocal:bool, defaultTrue- Convert energy units from Joules to calories
messages:bool, defaultTrue- Print informational messages (currently not used, reserved for future)
Returns
pd.DataFrame- Modified OBIGT dataframe with converted parameters
def available_properties(model: str = 'SUPCRT92') ‑> List[str]-
Expand source code
def available_properties(model: str = 'SUPCRT92') -> List[str]: """ Get list of available properties for a water model. Parameters ---------- model : str Water model name ('SUPCRT92', 'IAPWS95', or 'DEW') Returns ------- List[str] List of available property names """ model = model.upper() if model in ['SUPCRT92', 'SUPCRT']: from .supcrt92 import supcrt92_water return supcrt92_water.available_properties() elif model in ['IAPWS95', 'IAPWS']: from .iapws95 import iapws95_water return iapws95_water.available_properties() elif model == 'DEW': from .dew import dew_water return dew_water.available_properties() else: raise ValueError(f"Unknown water model: {model}")Get list of available properties for a water model.
Parameters
model:str- Water model name ('SUPCRT92', 'IAPWS95', or 'DEW')
Returns
List[str]- List of available property names
def calc_G_TP(OBIGT, Tc, P, water_model)-
Expand source code
def calc_G_TP(OBIGT, Tc, P, water_model): aq_out, H2O_Pt = hkf(property=["G"], parameters=OBIGT, T=273.15+Tc, P=P, contrib=["n", "s", "o"], H2O_props=["rho"], water_model=water_model) cgl_out = cgl(property=["G"], parameters=OBIGT, T=273.15+Tc, P=P) aq_col = pd.DataFrame.from_dict(aq_out, orient="index") cgl_col = pd.DataFrame.from_dict(cgl_out, orient="index") G_TP_df = pd.concat([aq_col, cgl_col], axis=1) G_TP_df.columns = ['aq','cgl'] OBIGT["G_TP"] = G_TP_df['aq'].combine_first(G_TP_df['cgl']) rows_added = 0 # add a row for water if "H2O" not in list(OBIGT["name"]): # Set the water model (without printing messages) water(water_model, messages=False) # water() returns a scalar when called with single property and scalar T, P # The result is in J/mol, need to convert to cal/mol by dividing by 4.184 G_water = water("G", T=Tc+273.15, P=P, messages=False) # Handle both scalar and DataFrame returns if isinstance(G_water, pd.DataFrame): G_water_cal = G_water.iloc[0]["G"] / 4.184 else: G_water_cal = float(G_water) / 4.184 OBIGT = pd.concat([OBIGT, pd.DataFrame({"name": "H2O", "tag": "nan", "G_TP": G_water_cal}, index=[OBIGT.shape[0]])], ignore_index=True) rows_added += 1 # add a row for protons if "H+" not in list(OBIGT["name"]): OBIGT = pd.concat([OBIGT, pd.DataFrame({"name": "H+", "tag": "nan", "G_TP": 0}, index=[OBIGT.shape[0]])], ignore_index=True) rows_added += 1 return OBIGT, rows_added def calc_logK(OBIGT_df, Tc, P, TP_i, water_model)-
Expand source code
def calc_logK(OBIGT_df, Tc, P, TP_i, water_model): OBIGT_TP, rows_added = calc_G_TP(OBIGT_df, Tc, P, water_model) dissrxn2logK_out = [] for i in OBIGT_TP.index: dissrxn2logK_out.append(dissrxn2logK(OBIGT_TP, i, Tc)) assert len(dissrxn2logK_out) == OBIGT_TP.shape[0] OBIGT_TP['dissrxn_logK_'+str(TP_i)] = dissrxn2logK_out # remove any rows added by calc_G_TP OBIGT_TP.drop(OBIGT_TP.tail(rows_added).index, inplace = True) return OBIGT_TP def cgl(property=None, parameters=None, T=298.15, P=1)-
Expand source code
def cgl(property = None, parameters = None, T = 298.15, P = 1): # calculate properties of crystalline, liquid (except H2O) and gas species Tr = 298.15 Pr = 1 # Convert T and P to arrays for vectorized operations T = np.atleast_1d(T) P = np.atleast_1d(P) # make T and P equal length if P.size < T.size: P = np.full_like(T, P[0] if P.size == 1 else P) if T.size < P.size: T = np.full_like(P, T[0] if T.size == 1 else T) n_conditions = T.size # initialize output dict out_dict = dict() # loop over each species # Iterate over each row by position to handle duplicate indices properly for i in range(len(parameters)): # Get the index label for this row k = parameters.index[i] # Get the row data by position (iloc) to avoid duplicate index issues PAR = parameters.iloc[i] if PAR["state"] == "aq": # For aqueous species processed by CGL, return NaN # (they should be processed by HKF instead) out_dict[k] = {p:float('NaN') for p in property} else: # OBIGT database stores G, H, S in calories (E_units = "cal") # CGL calculations use calories (integrals intCpdT, intCpdlnT, intVdP are in cal) # Results are output in calories and converted to J in subcrt.py at line 959 # Parameter scaling - SUPCRT92 data is already in correct units # PAR["a2.b"] = copy.copy(PAR["a2.b"]*10**-3) # PAR["a3.c"] = copy.copy(PAR["a3.c"]*10**5) # PAR["c1.e"] = copy.copy(PAR["c1.e"]*10**-5) # Check if this is a Berman mineral (columns 9-21 are all NA in R indexing) # In Python/pandas, we check the relevant thermodynamic parameter columns # NOTE: A mineral is only Berman if it LACKS standard thermodynamic data (G,H,S) # If G,H,S are present, use regular CGL even if heat capacity coefficients are all zero berman_cols = ['a1.a', 'a2.b', 'a3.c', 'a4.d', 'c1.e', 'c2.f', 'omega.lambda', 'z.T'] has_standard_thermo = pd.notna(PAR.get('G', np.nan)) and pd.notna(PAR.get('H', np.nan)) and pd.notna(PAR.get('S', np.nan)) all_coeffs_zero_or_na = all(pd.isna(PAR.get(col, np.nan)) or PAR.get(col, 0) == 0 for col in berman_cols) is_berman_mineral = all_coeffs_zero_or_na and not has_standard_thermo if is_berman_mineral: # Use Berman equations (parameters not in thermo()$OBIGT) from .berman import Berman try: # Berman is already vectorized - pass T and P arrays directly properties_df = Berman(PAR["name"], T=T, P=P) # Extract the requested properties as arrays values = {} for prop in property: if prop in properties_df.columns: # Get all values as an array prop_values = properties_df[prop].values # IMPORTANT: Berman function returns values in J/mol (Joules) # but CGL returns values in cal/mol (calories) # Convert Berman results from J/mol to cal/mol for consistency # Energy properties that need conversion: G, H, S, Cp # Volume (V) and other properties don't need conversion energy_props = ['G', 'H', 'S', 'Cp'] if prop in energy_props: # Convert J/mol to cal/mol by dividing by 4.184 prop_values = prop_values / 4.184 values[prop] = prop_values else: values[prop] = np.full(n_conditions, float('NaN')) except Exception as e: # If Berman calculation fails, fall back to NaN arrays values = {prop: np.full(n_conditions, float('NaN')) for prop in property} else: # Use regular CGL equations # in CHNOSZ, we have # 1 cm^3 bar --> convert(1, "calories") == 0.02390057 cal # but REAC92D.F in SUPCRT92 uses cm3bar_to_cal = 0.023901488 # cal # start with NA values values = dict() # a test for availability of heat capacity coefficients (a, b, c, d, e, f) # based on the column assignments in thermo()$OBIGT # Check for heat capacity coefficients, handling both NaN and non-numeric values # Heat capacity coefficients are at positions 14-19 (a1.a through c2.f) # Position 13 is V (volume), not a heat capacity coefficient has_hc_coeffs = False try: hc_values = list(PAR.iloc[14:20]) has_hc_coeffs = any([pd.notna(p) and p != 0 for p in hc_values if pd.api.types.is_numeric_dtype(type(p))]) # DEBUG if False and PAR["name"] == "rhomboclase": print(f"DEBUG for rhomboclase:") print(f" hc_values (iloc[14:20]): {hc_values}") print(f" has_hc_coeffs: {has_hc_coeffs}") except Exception as e: has_hc_coeffs = False if has_hc_coeffs: # we have at least one of the heat capacity coefficients; # zero out any NA's in the rest (leave lambda and T of transition (columns 20-21) alone) for i in range(14, 20): if pd.isna(PAR.iloc[i]) or not pd.api.types.is_numeric_dtype(type(PAR.iloc[i])): PAR.iloc[i] = 0.0 # calculate the heat capacity and its integrals (vectorized) Cp = PAR["a1.a"] + PAR["a2.b"]*T + PAR["a3.c"]*T**-2 + PAR["a4.d"]*T**-0.5 + PAR["c1.e"]*T**2 + PAR["c2.f"]*T**PAR["omega.lambda"] intCpdT = PAR["a1.a"]*(T - Tr) + PAR["a2.b"]*(T**2 - Tr**2)/2 + PAR["a3.c"]*(1/T - 1/Tr)/-1 + PAR["a4.d"]*(T**0.5 - Tr**0.5)/0.5 + PAR["c1.e"]*(T**3-Tr**3)/3 intCpdlnT = PAR["a1.a"]*np.log(T / Tr) + PAR["a2.b"]*(T - Tr) + PAR["a3.c"]*(T**-2 - Tr**-2)/-2 + PAR["a4.d"]*(T**-0.5 - Tr**-0.5)/-0.5 + PAR["c1.e"]*(T**2 - Tr**2)/2 # do we also have the lambda parameter (Cp term with adjustable exponent on T)? if pd.notna(PAR["omega.lambda"]) and PAR["omega.lambda"] != 0: # equations for lambda adapted from Helgeson et al., 1998 (doi:10.1016/S0016-7037(97)00219-6) if PAR["omega.lambda"] == -1: intCpdT = intCpdT + PAR["c2.f"]*np.log(T/Tr) else: intCpdT = intCpdT - PAR["c2.f"]*( T**(PAR["omega.lambda"] + 1) - Tr**(PAR["omega.lambda"] + 1) ) / (PAR["omega.lambda"] + 1) intCpdlnT = intCpdlnT + PAR["c2.f"]*(T**PAR["omega.lambda"] - Tr**PAR["omega.lambda"]) / PAR["omega.lambda"] else: # use constant heat capacity if the coefficients are not available (vectorized) # If Cp is NA/NaN, use 0 (matching R CHNOSZ behavior) Cp_value = PAR["Cp"] if pd.notna(PAR["Cp"]) else 0.0 Cp = np.full(n_conditions, Cp_value) intCpdT = Cp_value*(T - Tr) intCpdlnT = Cp_value*np.log(T / Tr) # in case Cp is listed as NA, set the integrals to 0 at Tr at_Tr = (T == Tr) intCpdT = np.where(at_Tr, 0, intCpdT) intCpdlnT = np.where(at_Tr, 0, intCpdlnT) # volume and its integrals (vectorized) if PAR["name"] in ["quartz", "coesite"]: # volume calculations for quartz and coesite qtz = quartz_coesite(PAR, T, P) V = qtz["V"] intVdP = qtz["intVdP"] intdVdTdP = qtz["intdVdTdP"] else: # for other minerals, volume is constant (Helgeson et al., 1978) V = np.full(n_conditions, PAR["V"]) # if the volume is NA, set its integrals to zero if pd.isna(PAR["V"]): intVdP = np.zeros(n_conditions) intdVdTdP = np.zeros(n_conditions) else: intVdP = PAR["V"]*(P - Pr) * cm3bar_to_cal intdVdTdP = np.zeros(n_conditions) # get the values of each of the requested thermodynamic properties (vectorized) for i,prop in enumerate(property): if prop == "Cp": values["Cp"] = Cp if prop == "V": values["V"] = V if prop == "E": values["E"] = np.full(n_conditions, float('NaN')) if prop == "kT": values["kT"] = np.full(n_conditions, float('NaN')) if prop == "G": # calculate S * (T - Tr), but set it to 0 at Tr (in case S is NA) Sterm = PAR["S"]*(T - Tr) Sterm = np.where(T == Tr, 0, Sterm) # DEBUG if False and PAR["name"] == "iron" and PAR.get("state") == "cr4": print(f"DEBUG G calculation for {PAR['name']} {PAR.get('state', 'unknown')}:") print(f" PAR['G'] = {PAR['G']}") print(f" PAR['S'] = {PAR['S']}") print(f" model = {PAR.get('model', 'unknown')}") print(f" Sterm[0] = {Sterm[0] if hasattr(Sterm, '__len__') else Sterm}") print(f" intCpdT[0] = {intCpdT[0] if hasattr(intCpdT, '__len__') else intCpdT}") print(f" T[0]*intCpdlnT[0] = {(T[0]*intCpdlnT[0]) if hasattr(intCpdlnT, '__len__') else T*intCpdlnT}") print(f" intVdP[0] = {intVdP[0] if hasattr(intVdP, '__len__') else intVdP}") G_calc = PAR['G'] - Sterm + intCpdT - T*intCpdlnT + intVdP print(f" G[0] (before subcrt conversion) = {G_calc[0] if hasattr(G_calc, '__len__') else G_calc}") values["G"] = PAR["G"] - Sterm + intCpdT - T*intCpdlnT + intVdP if prop == "H": values["H"] = PAR["H"] + intCpdT + intVdP - T*intdVdTdP if prop == "S": values["S"] = PAR["S"] + intCpdlnT - intdVdTdP out_dict[k] = values # species have to be numbered instead of named because of name repeats (e.g., cr polymorphs) return out_dict def compare_models(property: str,
T: float | numpy.ndarray = 298.15,
P: float | numpy.ndarray = 1.0) ‑> pandas.core.frame.DataFrame-
Expand source code
def compare_models(property: str, T: Union[float, np.ndarray] = 298.15, P: Union[float, np.ndarray] = 1.0) -> pd.DataFrame: """ Compare water property calculations across different models. Parameters ---------- property : str Property to compare T : float or array Temperature in Kelvin P : float or array Pressure in bar Returns ------- pd.DataFrame Comparison of property values from different models """ T = np.atleast_1d(np.asarray(T, dtype=float)) P = np.atleast_1d(np.asarray(P, dtype=float)) results = {} models = ['SUPCRT92', 'IAPWS95', 'DEW'] for model in models: try: if model == 'SUPCRT92': result = water_SUPCRT92(property, T, P) elif model == 'IAPWS95': # Convert bar to kPa for IAPWS95 result = water_IAPWS95(property, T, P * 100.0) # Convert units back if needed if property == 'rho': result = result / 1000.0 # kg/m³ to g/cm³ elif model == 'DEW': result = water_DEW(property, T, P) results[model] = result except Exception as e: print(f"Error with {model}: {e}") results[model] = np.full_like(T, np.nan) # Create DataFrame if len(T) == 1 and len(P) == 1: # Single point data = {model: [results[model]] if np.isscalar(results[model]) else results[model] for model in models} df = pd.DataFrame(data, index=[f"T={T[0]:.1f}K, P={P[0]:.1f}bar"]) else: # Multiple points data = {model: results[model] for model in models} index = [f"T={t:.1f}K, P={p:.1f}bar" for t, p in zip(T, P)] df = pd.DataFrame(data, index=index) return dfCompare water property calculations across different models.
Parameters
property:str- Property to compare
T:floatorarray- Temperature in Kelvin
P:floatorarray- Pressure in bar
Returns
pd.DataFrame- Comparison of property values from different models
def convert_cm3bar(value)-
Expand source code
def convert_cm3bar(value): return value*4.184 * 10 def dissrxn2logK(OBIGT, i, Tc)-
Expand source code
def dissrxn2logK(OBIGT, i, Tc): this_dissrxn = OBIGT.iloc[i, OBIGT.columns.get_loc('dissrxn')] if this_dissrxn == "nan": this_dissrxn = OBIGT.iloc[i, OBIGT.columns.get_loc('regenerate_dissrxn')] # print(OBIGT["name"][i], this_dissrxn) try: this_dissrxn = this_dissrxn.strip() split_dissrxn = this_dissrxn.split(" ") except: return float('NaN') coeff = [float(n) for n in split_dissrxn[::2]] species = split_dissrxn[1::2] try: G = sum([float(c*OBIGT.loc[OBIGT["name"]==sp, "G_TP"].iloc[0]) for c,sp in zip(coeff, species)]) except: G_list = [] for ii, sp in enumerate(species): G_TP = OBIGT.loc[OBIGT["name"]==sp, "G_TP"] if len(G_TP) == 1: G_list.append(float(coeff[ii]*OBIGT.loc[OBIGT["name"]==sp, "G_TP"])) else: ### check valid polymorph T # get polymorph entries of OBIGT that match mineral poly_df = copy.copy(OBIGT.loc[OBIGT["name"]==sp,:]) # ensure polymorph df is sorted according to cr, cr2, cr3... etc. poly_df = poly_df.sort_values("state") z_Ts = list(poly_df.loc[poly_df["name"]==sp, "z.T"]) last_t = float('-inf') appended=False for iii,t in enumerate(z_Ts): if Tc+273.15 > last_t and Tc+273.15 < t: G_list.append(float(coeff[ii]*list(poly_df.loc[poly_df["name"]==sp, "G_TP"])[iii])) appended=True if not appended and z_Ts[-1] == t: G_list.append(float(coeff[ii]*list(poly_df.loc[poly_df["name"]==sp, "G_TP"])[iii])) last_t = t G = sum(G_list) return G2logK(G, Tc) def get_water_models() ‑> List[str]-
Expand source code
def get_water_models() -> List[str]: """ Get list of available water models. Returns ------- List[str] List of available water model names """ return ['SUPCRT92', 'IAPWS95', 'DEW']Get list of available water models.
Returns
List[str]- List of available water model names
def gfun(rhohat, Tc, P, alpha, daldT, beta)-
Expand source code
def gfun(rhohat, Tc, P, alpha, daldT, beta): ## g and f functions for describing effective electrostatic radii of ions ## split from hkf() 20120123 jmd ## based on equations in ## Shock EL, Oelkers EH, Johnson JW, Sverjensky DA, Helgeson HC, 1992 ## Calculation of the Thermodynamic Properties of Aqueous Species at High Pressures ## and Temperatures: Effective Electrostatic Radii, Dissociation Constants and ## Standard Partial Molal Properties to 1000 degrees C and 5 kbar ## J. Chem. Soc. Faraday Trans., 88(6), 803-826 doi:10.1039/FT9928800803 # rhohat - density of water in g/cm3 # Tc - temperature in degrees Celsius # P - pressure in bars # Vectorized version - handle both scalars and arrays rhohat = np.atleast_1d(rhohat) Tc = np.atleast_1d(Tc) P = np.atleast_1d(P) alpha = np.atleast_1d(alpha) daldT = np.atleast_1d(daldT) beta = np.atleast_1d(beta) # Broadcast to same shape shape = np.broadcast_shapes(rhohat.shape, Tc.shape, P.shape, alpha.shape, daldT.shape, beta.shape) rhohat = np.broadcast_to(rhohat, shape) Tc = np.broadcast_to(Tc, shape) P = np.broadcast_to(P, shape) alpha = np.broadcast_to(alpha, shape) daldT = np.broadcast_to(daldT, shape) beta = np.broadcast_to(beta, shape) # Initialize output arrays g = np.zeros(shape) dgdT = np.zeros(shape) d2gdT2 = np.zeros(shape) dgdP = np.zeros(shape) # only rhohat less than 1 will give results other than zero mask = rhohat < 1 if not np.any(mask): return {"g": g, "dgdT": dgdT, "d2gdT2": d2gdT2, "dgdP": dgdP} # eta in Eq. 1 eta = 1.66027E5 # Table 3 ag1 = -2.037662 ag2 = 5.747000E-3 ag3 = -6.557892E-6 bg1 = 6.107361 bg2 = -1.074377E-2 bg3 = 1.268348E-5 # Work only with masked values Tc_m = Tc[mask] P_m = P[mask] rhohat_m = rhohat[mask] alpha_m = alpha[mask] daldT_m = daldT[mask] beta_m = beta[mask] # Eq. 25 ag = ag1 + ag2 * Tc_m + ag3 * Tc_m ** 2 # Eq. 26 bg = bg1 + bg2 * Tc_m + bg3 * Tc_m ** 2 # Eq. 24 g_m = ag * (1 - rhohat_m) ** bg # Table 4 af1 = 0.3666666E2 af2 = -0.1504956E-9 af3 = 0.5017997E-13 # Eq. 33 f = ( ((Tc_m - 155) / 300) ** 4.8 + af1 * ((Tc_m - 155) / 300) ** 16 ) * \ ( af2 * (1000 - P_m) ** 3 + af3 * (1000 - P_m) ** 4 ) # limits of the f function (region II of Fig. 6) ifg = (Tc_m > 155) & (P_m < 1000) & (Tc_m < 355) # Eq. 32 - apply f correction where ifg is True # Check for complex values f_is_real = ~np.iscomplex(f) apply_f = ifg & f_is_real g_m = np.where(apply_f, g_m - f.real, g_m) # at P > 6000 bar (in DEW calculations), g is zero 20170926 g_m = np.where(P_m > 6000, 0, g_m) ## now we have g at P, T # put the results in their right place (where rhohat < 1) g[mask] = g_m ## the rest is to get its partial derivatives with pressure and temperature ## after Johnson et al., 1992 # alpha - coefficient of isobaric expansivity (K^-1) # daldT - temperature derivative of coefficient of isobaric expansivity (K^-2) # beta - coefficient of isothermal compressibility (bar^-1) # Eqn. 76 d2fdT2 = (0.0608/300*((Tc_m-155)/300)**2.8 + af1/375*((Tc_m-155)/300)**14) * (af2*(1000-P_m)**3 + af3*(1000-P_m)**4) # Eqn. 75 dfdT = (0.016*((Tc_m-155)/300)**3.8 + 16*af1/300*((Tc_m-155)/300)**15) * \ (af2*(1000-P_m)**3 + af3*(1000-P_m)**4) # Eqn. 74 dfdP = -(((Tc_m-155)/300)**4.8 + af1*((Tc_m-155)/300)**16) * \ (3*af2*(1000-P_m)**2 + 4*af3*(1000-P_m)**3) d2bdT2 = 2 * bg3 # Eqn. 73 d2adT2 = 2 * ag3 # Eqn. 72 dbdT = bg2 + 2*bg3*Tc_m # Eqn. 71 dadT = ag2 + 2*ag3*Tc_m # Eqn. 70 # Convert complex to NaN d2fdT2 = np.where(np.iscomplex(d2fdT2), np.nan, np.real(d2fdT2)) dfdT = np.where(np.iscomplex(dfdT), np.nan, np.real(dfdT)) dfdP = np.where(np.iscomplex(dfdP), np.nan, np.real(dfdP)) # Initialize derivative arrays for masked region dgdT_m = np.zeros_like(g_m) d2gdT2_m = np.zeros_like(g_m) dgdP_m = np.zeros_like(g_m) # Calculate derivatives where alpha and daldT are not NaN alpha_valid = ~np.isnan(alpha_m) & ~np.isnan(daldT_m) if np.any(alpha_valid): # Work with valid subset av_idx = alpha_valid bg_av = bg[av_idx] rhohat_av = rhohat_m[av_idx] alpha_av = alpha_m[av_idx] daldT_av = daldT_m[av_idx] g_av = g_m[av_idx] ag_av = ag[av_idx] Tc_av = Tc_m[av_idx] dbdT_av = dbdT[av_idx] dadT_av = dadT[av_idx] # Handle log of (1-rhohat) safely with np.errstate(divide='ignore', invalid='ignore'): log_term = np.log(1 - rhohat_av) log_term = np.where(np.isfinite(log_term), log_term, 0) # Eqn. 69 dgadT = bg_av*rhohat_av*alpha_av*(1-rhohat_av)**(bg_av-1) + log_term*g_av/ag_av*dbdT_av D = rhohat_av # transcribed from SUPCRT92/reac92.f dDdT = -D * alpha_av dDdTT = -D * (daldT_av - alpha_av**2) Db = (1-D)**bg_av dDbdT = -bg_av*(1-D)**(bg_av-1)*dDdT + log_term*Db*dbdT_av dDbdTT = -(bg_av*(1-D)**(bg_av-1)*dDdTT + (1-D)**(bg_av-1)*dDdT*dbdT_av + \ bg_av*dDdT*(-(bg_av-1)*(1-D)**(bg_av-2)*dDdT + log_term*(1-D)**(bg_av-1)*dbdT_av)) + \ log_term*(1-D)**bg_av*d2bdT2 - (1-D)**bg_av*dbdT_av*dDdT/(1-D) + log_term*dbdT_av*dDbdT d2gdT2_calc = ag_av*dDbdTT + 2*dDbdT*dadT_av + Db*d2adT2 # Apply f correction where ifg is True ifg_av = ifg[av_idx] d2fdT2_av = d2fdT2[av_idx] dfdT_av = dfdT[av_idx] d2gdT2_calc = np.where(ifg_av, d2gdT2_calc - d2fdT2_av, d2gdT2_calc) dgdT_calc = g_av/ag_av*dadT_av + ag_av*dgadT # Eqn. 67 dgdT_calc = np.where(ifg_av, dgdT_calc - dfdT_av, dgdT_calc) dgdT_m[av_idx] = dgdT_calc d2gdT2_m[av_idx] = d2gdT2_calc # Calculate dgdP where beta is not NaN beta_valid = ~np.isnan(beta_m) if np.any(beta_valid): bv_idx = beta_valid bg_bv = bg[bv_idx] rhohat_bv = rhohat_m[bv_idx] beta_bv = beta_m[bv_idx] g_bv = g_m[bv_idx] dgdP_calc = -bg_bv*rhohat_bv*beta_bv*g_bv*(1-rhohat_bv)**-1 # Eqn. 66 ifg_bv = ifg[bv_idx] dfdP_bv = dfdP[bv_idx] dgdP_calc = np.where(ifg_bv, dgdP_calc - dfdP_bv, dgdP_calc) dgdP_m[bv_idx] = dgdP_calc # Put results back into full arrays dgdT[mask] = dgdT_m d2gdT2[mask] = d2gdT2_m dgdP[mask] = dgdP_m return {"g": g, "dgdT": dgdT, "d2gdT2": d2gdT2, "dgdP": dgdP} def hkf(property=None,
parameters=None,
T=298.15,
P=1,
contrib=['n', 's', 'o'],
H2O_props=['rho'],
water_model='SUPCRT92')-
Expand source code
def hkf(property=None, parameters=None, T=298.15, P=1, contrib = ["n", "s", "o"], H2O_props=["rho"], water_model="SUPCRT92"): # calculate G, H, S, Cp, V, kT, and/or E using # the revised HKF equations of state # H2O_props - H2O properties needed for subcrt() output # constants Tr = 298.15 # K Pr = 1 # bar Theta = 228 # K Psi = 2600 # bar # Convert T and P to arrays for vectorized operations T = np.atleast_1d(T) P = np.atleast_1d(P) # DEBUG if False: print(f"\nDEBUG HKF input:") print(f" T (K): {T}") print(f" P (bar): {P}") # make T and P equal length if P.size < T.size: P = np.full_like(T, P[0] if P.size == 1 else P) if T.size < P.size: T = np.full_like(P, T[0] if T.size == 1 else T) n_conditions = T.size # GB conversion note: handle error messages later # # nonsolvation, solvation, and origination contribution # notcontrib <- ! contrib %in% c("n", "s", "o") # if(TRUE %in% notcontrib) stop(paste("contrib must be in c('n', 's', 'o); got", c2s(contrib[notcontrib]))) # get water properties # rho - for subcrt() output and g function # Born functions and epsilon - for HKF calculations H2O_props += ["QBorn", "XBorn", "YBorn", "epsilon"] if water_model == "SUPCRT92": # using H2O92D.f from SUPCRT92: alpha, daldT, beta - for partial derivatives of omega (g function) H2O_props += ["alpha", "daldT", "beta"] elif water_model == "IAPWS95": # using IAPWS-95: NBorn, UBorn - for compressibility, expansibility H2O_props += ["alpha", "daldT", "beta", "NBorn", "UBorn"] elif water_model == "DEW": # using DEW model: get beta to calculate dgdP H2O_props += ["alpha", "daldT", "beta"] # DEBUG: Print T and P being passed to water if False: print(f"DEBUG HKF calling water():") print(f" T type: {type(T)}, T: {T}") print(f" P type: {type(P)}, P: {P}") print(f" H2O_props: {H2O_props}") H2O_PrTr = water(H2O_props, T=Tr, P=Pr) H2O_PT = water(H2O_props, T=T, P=P) # DEBUG: Print what water returned if False: print(f"DEBUG HKF water() returned:") print(f" H2O_PT type: {type(H2O_PT)}") if isinstance(H2O_PT, dict): print(f" H2O_PT keys: {H2O_PT.keys()}") print(f" epsilon: {H2O_PT.get('epsilon', 'NOT FOUND')}") # Handle dict output from water function def get_water_prop(water_dict, prop): """Helper function to get water property from dict or DataFrame""" if isinstance(water_dict, dict): return water_dict[prop] else: return water_dict.loc["1", prop] # Get epsilon values and handle potential zeros epsilon_PT = get_water_prop(H2O_PT, "epsilon") epsilon_PrTr = get_water_prop(H2O_PrTr, "epsilon") # Check for zero or very small epsilon values and warn if np.any(epsilon_PT == 0) or np.any(np.abs(epsilon_PT) < 1e-10): warnings.warn(f"HKF: epsilon at P,T is zero or very small: {epsilon_PT}. H2O_PT keys: {H2O_PT.keys() if isinstance(H2O_PT, dict) else 'not dict'}") with np.errstate(divide='ignore', invalid='ignore'): ZBorn = -1 / epsilon_PT ZBorn_PrTr = -1 / epsilon_PrTr # a class to store the result out_dict = {} # dictionary to store output for k in parameters.index: if parameters["state"][k] != "aq": out_dict[k] = {p:float('NaN') for p in property} else: sp = parameters["name"][k] # loop over each species PAR = copy.copy(parameters.loc[k, :]) PAR["a1.a"] = copy.copy(PAR["a1.a"]*10**-1) PAR["a2.b"] = copy.copy(PAR["a2.b"]*10**2) PAR["a4.d"] = copy.copy(PAR["a4.d"]*10**4) PAR["c2.f"] = copy.copy(PAR["c2.f"]*10**4) PAR["omega.lambda"] = copy.copy(PAR["omega.lambda"]*10**5) # substitute Cp and V for missing EoS parameters # here we assume that the parameters are in the same position as in thermo()$OBIGT # we don't need this if we're just looking at solvation properties (Cp_s_var, V_s_var) # GB conversion note: this block checks various things about EOS parameters. # for now, just set hasEOS to True hasEOS = True # delete this once the following block is converted to python # if "n" in contrib: # # put the heat capacity in for c1 if both c1 and c2 are missing # if all(is.na(PAR[, 18:19])): # PAR[, 18] = PAR["Cp"] # # put the volume in for a1 if a1, a2, a3 and a4 are missing # if all(is.na(PAR[, 14:17])): # PAR[, 14] = convert(PAR["V"], "calories") # # test for availability of the EoS parameters # hasEOS = any(!is.na(PAR[, 14:21])) # # if at least one of the EoS parameters is available, zero out any NA's in the rest # if hasEOS: # PAR[, 14:21][, is.na(PAR[, 14:21])] = 0 # compute values of omega(P,T) from those of omega(Pr,Tr) # using g function etc. (Shock et al., 1992 and others) omega = PAR["omega.lambda"] # omega_PrTr # its derivatives are zero unless the g function kicks in dwdP = np.zeros(n_conditions) dwdT = np.zeros(n_conditions) d2wdT2 = np.zeros(n_conditions) Z = PAR["z.T"] omega_PT = np.full(n_conditions, PAR["omega.lambda"]) if Z != 0 and Z != "NA" and PAR["name"] != "H+": # compute derivatives of omega: g and f functions (Shock et al., 1992; Johnson et al., 1992) rhohat = get_water_prop(H2O_PT, "rho")/1000 # just converting kg/m3 to g/cm3 # temporarily filter out Python's warnings about dividing by zero, which is possible # with the equations in the gfunction # Possible complex output is acounted for in gfun(). with warnings.catch_warnings(): warnings.simplefilter('ignore') g = gfun(rhohat, T-273.15, P, get_water_prop(H2O_PT, "alpha"), get_water_prop(H2O_PT, "daldT"), get_water_prop(H2O_PT, "beta")) # after SUPCRT92/reac92.f eta = 1.66027E5 reref = (Z**2) / (omega/eta + Z/(3.082 + 0)) re = reref + abs(Z) * g["g"] omega_PT = eta * (Z**2/re - Z/(3.082 + g["g"])) Z3 = abs(Z**3)/re**2 - Z/(3.082 + g["g"])**2 Z4 = abs(Z**4)/re**3 - Z/(3.082 + g["g"])**3 dwdP = (-eta * Z3 * g["dgdP"]) dwdT = (-eta * Z3 * g["dgdT"]) d2wdT2 = (2 * eta * Z4 * g["dgdT"]**2 - eta * Z3 * g["d2gdT2"]) # loop over each property w = float('NaN') for i,PROP in enumerate(property) : # over nonsolvation, solvation, or origination contributions - vectorized hkf_p = np.zeros(n_conditions) for icontrib in contrib : # various contributions to the properties if icontrib == "n": # nonsolvation ghs equations if PROP == "H": p_c = PAR["c1.e"]*(T-Tr) - PAR["c2.f"]*(1/(T-Theta)-1/(Tr-Theta)) p_a = PAR["a1.a"]*(P-Pr) + PAR["a2.b"]*np.log((Psi+P)/(Psi+Pr)) + \ ((2*T-Theta)/(T-Theta)**2)*(PAR["a3.c"]*(P-Pr)+PAR["a4.d"]*np.log((Psi+P)/(Psi+Pr))) p = p_c + p_a elif PROP == "S": p_c = PAR["c1.e"]*np.log(T/Tr) - \ (PAR["c2.f"]/Theta)*( 1/(T-Theta)-1/(Tr-Theta) + \ np.log( (Tr*(T-Theta))/(T*(Tr-Theta)) )/Theta ) p_a = (T-Theta)**(-2)*(PAR["a3.c"]*(P-Pr)+PAR["a4.d"]*np.log((Psi+P)/(Psi+Pr))) p = p_c + p_a elif PROP == "G": p_c = -PAR["c1.e"]*(T*np.log(T/Tr)-T+Tr) - \ PAR["c2.f"]*( (1/(T-Theta)-1/(Tr-Theta))*((Theta-T)/Theta) - \ (T/Theta**2)*np.log((Tr*(T-Theta))/(T*(Tr-Theta))) ) p_a = PAR["a1.a"]*(P-Pr) + PAR["a2.b"]*np.log((Psi+P)/(Psi+Pr)) + \ (PAR["a3.c"]*(P-Pr) + PAR["a4.d"]*np.log((Psi+P)/(Psi+Pr)))/(T-Theta) p = p_c + p_a # at Tr,Pr, if the origination contribution is not NA, ensure the solvation contribution is 0, not NA if not np.isnan(PAR["G"]): p = np.where((T==Tr) & (P==Pr), 0, p) # nonsolvation cp v kt e equations elif PROP == "Cp": p = PAR["c1.e"] + PAR["c2.f"] * ( T - Theta ) ** (-2) elif PROP == "V": p = convert_cm3bar(PAR["a1.a"]) + \ convert_cm3bar(PAR["a2.b"]) / (Psi + P) + \ (convert_cm3bar(PAR["a3.c"]) + convert_cm3bar(PAR["a4.d"]) / (Psi + P)) / (T - Theta) # elif PROP == "kT": # p = (convert(PAR["a2.b"], "cm3bar") + \ # convert(PAR["a4.d"], "cm3bar") / (T - Theta)) * (Psi + P) ** (-2) # elif PROP == "E": # p = convert( - (PAR["a3.c"] + PAR["a4.d"] / convert((Psi + P), "calories")) * \ # (T - Theta) ** (-2), "cm3bar") else: print("BAD") if icontrib == "s": # solvation ghs equations if PROP == "G": p = -omega_PT*(ZBorn+1) + omega*(ZBorn_PrTr+1) + omega*get_water_prop(H2O_PrTr, "YBorn")*(T-Tr) # at Tr,Pr, if the origination contribution is not NA, ensure the solvation contribution is 0, not NA if(np.isnan(PAR["G"])): p = np.where((T==Tr) & (P==Pr), 0, p) if PROP == "H": p = -omega_PT*(ZBorn+1) + omega_PT*T*get_water_prop(H2O_PT, "YBorn") + T*(ZBorn+1)*dwdT + \ omega*(ZBorn_PrTr+1) - omega*Tr*get_water_prop(H2O_PrTr, "YBorn") if PROP == "S": p = omega_PT*get_water_prop(H2O_PT, "YBorn") + (ZBorn+1)*dwdT - omega*get_water_prop(H2O_PrTr, "YBorn") # solvation cp v kt e equations if PROP == "Cp": p = omega_PT*T*get_water_prop(H2O_PT, "XBorn") + 2*T*get_water_prop(H2O_PT, "YBorn")*dwdT + T*(ZBorn+1)*d2wdT2 if PROP == "V": term1 = -convert_cm3bar(omega_PT) * get_water_prop(H2O_PT, "QBorn") term2 = convert_cm3bar(dwdP) * (-ZBorn - 1) p = term1 + term2 # DEBUG if False: print(f"\nDEBUG solvation V terms:") print(f" omega_PT: {omega_PT}") print(f" QBorn: {get_water_prop(H2O_PT, 'QBorn')}") print(f" dwdP: {dwdP}") print(f" ZBorn: {ZBorn}") print(f" term1 (-ω*QBorn): {term1}") print(f" term2 (dwdP*(-Z-1)): {term2}") print(f" total p: {p}") # TODO: the partial derivatives of omega are not included here here for kt and e # (to do it, see p. 820 of SOJ+92 ... but kt requires d2wdP2 which we don"t have yet) if PROP == "kT": p = convert_cm3bar(omega) * get_water_prop(H2O_PT, "NBorn") if PROP == "E": p = -convert_cm3bar(omega) * get_water_prop(H2O_PT, "UBorn") if icontrib == "o": # origination ghs equations if PROP == "G": p = PAR["G"] - PAR["S"] * (T-Tr) # don"t inherit NA from PAR$S at Tr p = np.where(T == Tr, PAR["G"], p) elif PROP == "H": p = np.full(n_conditions, PAR["H"]) elif PROP == "S": p = np.full(n_conditions, PAR["S"]) # origination eos equations (Cp, V, kT, E): senseless else: p = np.zeros(n_conditions) # accumulate the contribution hkf_p = hkf_p + p # DEBUG if False and PROP == "V": print(f"\nDEBUG HKF V calculation (species {k}, contrib={icontrib}):") print(f" T: {T}") print(f" P: {P}") print(f" contribution p: {p}") print(f" accumulated hkf_p: {hkf_p}") # species have to be numbered (k) instead of named because of name repeats in db (e.g., cr polymorphs) if i > 0: out_dict[k][PROP] = hkf_p else: out_dict[k] = {PROP:hkf_p} # DEBUG if False and PROP == "V": print(f"\nDEBUG HKF final V for species {k}: {hkf_p}") return(out_dict, H2O_PT) def quartz_coesite(PAR, T, P)-
Expand source code
def quartz_coesite(PAR, T, P): # the corrections are 0 for anything other than quartz and coesite if not PAR["name"] in ["quartz", "coesite"]: n = T.size if isinstance(T, np.ndarray) else 1 return(dict(G=np.zeros(n), H=np.zeros(n), S=np.zeros(n), V=np.zeros(n))) # Vectorized version T = np.atleast_1d(T) P = np.atleast_1d(P) # Tr, Pr and TtPr (transition temperature at Pr) Pr = 1 # bar Tr = 298.15 # K TtPr = 848 # K # constants from SUP92D.f aa = 549.824 ba = 0.65995 ca = -0.4973e-4 VPtTta = 23.348 VPrTtb = 23.72 Stran = 0.342 # constants from REAC92D.f VPrTra = 22.688 # VPrTr(a-quartz) Vdiff = 2.047 # VPrTr(a-quartz) - VPrTr(coesite) k = 38.5 # dPdTtr(a/b-quartz) #k <- 38.45834 # calculated in CHNOSZ: dPdTtr(info("quartz")) # code adapted from REAC92D.f qphase = PAR["state"].replace("cr", "") if qphase == "2": Pstar = P.copy() Sstar = np.zeros_like(T) V = np.full_like(T, VPrTtb) else: Pstar = Pr + k * (T - TtPr) Sstar = np.full_like(T, Stran) V = VPrTra + ca*(P-Pr) + (VPtTta - VPrTra - ca*(P-Pr))*(T-Tr) / (TtPr + (P-Pr)/k - Tr) # Apply condition: if T < TtPr below_transition = T < TtPr Pstar = np.where(below_transition, Pr, Pstar) Sstar = np.where(below_transition, 0, Sstar) if PAR["name"] == "coesite": VPrTra = VPrTra - Vdiff VPrTtb = VPrTtb - Vdiff V = V - Vdiff cm3bar_to_cal = 0.023901488 # Vectorized log calculation with np.errstate(divide='ignore', invalid='ignore'): log_term = np.log((aa + P/k) / (aa + Pstar/k)) log_term = np.where(np.isfinite(log_term), log_term, 0) GVterm = cm3bar_to_cal * (VPrTra * (P - Pstar) + VPrTtb * (Pstar - Pr) - \ 0.5 * ca * (2 * Pr * (P - Pstar) - (P**2 - Pstar**2)) - \ ca * k * (T - Tr) * (P - Pstar) + \ k * (ba + aa * ca * k) * (T - Tr) * log_term) SVterm = cm3bar_to_cal * (-k * (ba + aa * ca * k) * log_term + ca * k * (P - Pstar)) - Sstar # note the minus sign on "SVterm" in order that intdVdTdP has the correct sign return dict(intVdP=GVterm, intdVdTdP=-SVterm, V=V) def water(property: str | List[str] | None = None,
T: float | numpy.ndarray | List[float] = 298.15,
P: float | List[float] | numpy.ndarray | str = 1.0,
Psat_floor: float | None = 1.0,
model: str | None = None,
messages: bool = True) ‑> str | float | numpy.ndarray | Dict[str, Any]-
Expand source code
def water(property: Optional[Union[str, List[str]]] = None, T: Union[float, np.ndarray, List[float]] = 298.15, P: Union[float, np.ndarray, List[float], str] = 1.0, Psat_floor: Union[float, None] = 1.0, model: Optional[str] = None, messages: bool = True) -> Union[str, float, np.ndarray, Dict[str, Any]]: """ Calculate thermodynamic and electrostatic properties of liquid H2O. This is the main water function that provides the same interface as the R CHNOSZ water() function, with support for multiple water models. Parameters ---------- property : str, list of str, or None Properties to calculate. If None, returns current water model. If water model name (SUPCRT92, IAPWS95, DEW), sets the water model. Available properties depend on the water model used. T : float or array-like Temperature in Kelvin P : float, array-like, or "Psat" Pressure in bar, or "Psat" for saturation pressure Psat_floor : float or None Minimum pressure floor for Psat calculations (SUPCRT92 only) model : str, optional Override the default water model for this calculation messages : bool, default True Whether to print informational messages Returns ------- str, float, array, or dict Current water model name, single property value, array of values, or dictionary with calculated properties Examples -------- >>> import pychnosz >>> pychnosz.reset() >>> >>> # Get current water model >>> model = pychnosz.water() >>> print(model) # 'SUPCRT92' >>> >>> # Set water model >>> old_model = pychnosz.water('IAPWS95') >>> >>> # Calculate single property >>> density = pychnosz.water('rho', T=298.15, P=1.0) >>> >>> # Calculate multiple properties >>> props = pychnosz.water(['rho', 'epsilon'], T=298.15, P=1.0) >>> >>> # Temperature array >>> temps = np.array([273.15, 298.15, 373.15]) >>> densities = pychnosz.water('rho', T=temps, P=1.0) >>> >>> # Saturation pressure >>> psat = pychnosz.water('Psat', T=373.15) """ # Get thermo system thermo_system = thermo() # Ensure thermo is initialized before accessing/setting options # This prevents reset() from clearing options later if not thermo_system.is_initialized(): thermo_system.reset(messages=False) # Case 1: Query current water model if property is None: return thermo_system.get_option('water', 'SUPCRT92') # Case 2: Set water model if isinstance(property, str) and property.upper() in ['SUPCRT92', 'SUPCRT', 'IAPWS95', 'IAPWS', 'DEW']: old_model = thermo_system.get_option('water', 'SUPCRT92') # Normalize model name if property.upper() in ['SUPCRT92', 'SUPCRT']: new_model = 'SUPCRT92' elif property.upper() in ['IAPWS95', 'IAPWS']: new_model = 'IAPWS95' elif property.upper() == 'DEW': new_model = 'DEW' thermo_system.set_option('water', new_model) if messages: print(f"water: setting water model to {new_model}") return # Return None instead of the old model # Case 3: Calculate properties # Determine which model to use if model is not None: water_model = model.upper() else: water_model = thermo_system.get_option('water', 'SUPCRT92').upper() # Normalize model names if water_model in ['SUPCRT92', 'SUPCRT']: water_model = 'SUPCRT92' elif water_model in ['IAPWS95', 'IAPWS']: water_model = 'IAPWS95' elif water_model == 'DEW': water_model = 'DEW' else: warnings.warn(f"Unknown water model '{water_model}', using SUPCRT92") water_model = 'SUPCRT92' # Convert inputs T = np.atleast_1d(np.asarray(T, dtype=float)) if isinstance(P, str): P_input = P else: P_input = np.atleast_1d(np.asarray(P, dtype=float)) # Make T and P same length if len(P_input) < len(T): P_input = np.resize(P_input, len(T)) elif len(T) < len(P_input): T = np.resize(T, len(P_input)) # Call appropriate water model try: if water_model == 'SUPCRT92': result = _call_supcrt92(property, T, P_input, Psat_floor) elif water_model == 'IAPWS95': result = _call_iapws95(property, T, P_input, Psat_floor) elif water_model == 'DEW': result = _call_dew(property, T, P_input) else: raise ValueError(f"Unsupported water model: {water_model}") except Exception as e: raise WaterModelError(f"Error calculating water properties with {water_model} model: {e}") # Apply Psat rounding to match R CHNOSZ behavior # Round Psat values to 4 decimal places (round up to ensure liquid phase) result = _apply_psat_rounding(result, property) return resultCalculate thermodynamic and electrostatic properties of liquid H2O.
This is the main water function that provides the same interface as the R CHNOSZ water() function, with support for multiple water models.
Parameters
property:str, listofstr,orNone- Properties to calculate. If None, returns current water model. If water model name (SUPCRT92, IAPWS95, DEW), sets the water model. Available properties depend on the water model used.
T:floatorarray-like- Temperature in Kelvin
P:float, array-like,or"Psat"- Pressure in bar, or "Psat" for saturation pressure
Psat_floor:floatorNone- Minimum pressure floor for Psat calculations (SUPCRT92 only)
model:str, optional- Override the default water model for this calculation
messages:bool, defaultTrue- Whether to print informational messages
Returns
str, float, array,ordict- Current water model name, single property value, array of values, or dictionary with calculated properties
Examples
>>> import pychnosz >>> pychnosz.reset() >>> >>> # Get current water model >>> model = pychnosz.water() >>> print(model) # 'SUPCRT92' >>> >>> # Set water model >>> old_model = pychnosz.water('IAPWS95') >>> >>> # Calculate single property >>> density = pychnosz.water('rho', T=298.15, P=1.0) >>> >>> # Calculate multiple properties >>> props = pychnosz.water(['rho', 'epsilon'], T=298.15, P=1.0) >>> >>> # Temperature array >>> temps = np.array([273.15, 298.15, 373.15]) >>> densities = pychnosz.water('rho', T=temps, P=1.0) >>> >>> # Saturation pressure >>> psat = pychnosz.water('Psat', T=373.15) def water_SUPCRT92(property: str | List[str],
T: float | numpy.ndarray = 298.15,
P: float | str | numpy.ndarray = 1.0,
**kwargs) ‑> float | numpy.ndarray | Dict[str, Any]-
Expand source code
def water_SUPCRT92(property: Union[str, List[str]], T: Union[float, np.ndarray] = 298.15, P: Union[float, np.ndarray, str] = 1.0, **kwargs) -> Union[float, np.ndarray, Dict[str, Any]]: """ Calculate water properties using SUPCRT92 model. This function provides the same interface as the original Python implementation but with Fortran backend support for high accuracy. Parameters ---------- property : str or list of str Property or list of properties to calculate T : float or array Temperature in Kelvin P : float, array, or 'Psat' Pressure in bar, or 'Psat' for saturation pressure **kwargs Additional options (e.g., Psat_floor) Returns ------- float, array, or dict Calculated properties Examples -------- >>> # Basic usage >>> rho = water_SUPCRT92('rho', 298.15, 1.0) >>> print(f"Density: {rho:.3f} g/cm³") >>> # Multiple properties >>> props = water_SUPCRT92(['rho', 'epsilon'], 298.15, 1.0) >>> print(f"ρ = {props['rho']:.3f}, ε = {props['epsilon']:.1f}") >>> # Saturation pressure >>> Psat = water_SUPCRT92('Psat', 373.15, 'Psat') >>> print(f"Psat at 100°C: {Psat:.2f} bar") """ return supcrt92_water.calculate(property, T, P, **kwargs)Calculate water properties using SUPCRT92 model.
This function provides the same interface as the original Python implementation but with Fortran backend support for high accuracy.
Parameters
property:strorlistofstr- Property or list of properties to calculate
T:floatorarray- Temperature in Kelvin
P:float, array,or'Psat'- Pressure in bar, or 'Psat' for saturation pressure
**kwargs- Additional options (e.g., Psat_floor)
Returns
float, array,ordict- Calculated properties
Examples
>>> # Basic usage >>> rho = water_SUPCRT92('rho', 298.15, 1.0) >>> print(f"Density: {rho:.3f} g/cm³")>>> # Multiple properties >>> props = water_SUPCRT92(['rho', 'epsilon'], 298.15, 1.0) >>> print(f"ρ = {props['rho']:.3f}, ε = {props['epsilon']:.1f}")>>> # Saturation pressure >>> Psat = water_SUPCRT92('Psat', 373.15, 'Psat') >>> print(f"Psat at 100°C: {Psat:.2f} bar")
Classes
class SUPCRT92Water-
Expand source code
class SUPCRT92Water: """ SUPCRT92 water model with Fortran backend. This class provides an interface to the original SUPCRT92 Fortran subroutines for calculating water properties. If the Fortran library is not available, it falls back to Python approximations. The Fortran implementation gives exact compatibility with R CHNOSZ and includes all 23+ thermodynamic properties. """ def __init__(self): """ Initialize SUPCRT92 water model. SUPCRT92 requires the compiled FORTRAN interface for accuracy and compatibility with R CHNOSZ. No pure Python fallback is provided. """ if not HAS_FORTRAN: raise ImportError( f"SUPCRT92 water model requires compiled FORTRAN interface. " f"Error: {_fortran_error}. " f"Please compile the FORTRAN subroutines using: " f"python compile_fortran.py. " f"See setup_fortran_instructions.md for details." ) # Initialize Fortran interface (required) try: self._fortran_interface = get_h2o92_interface() except Exception as e: raise RuntimeError(f"Failed to initialize FORTRAN interface: {e}") from e def calculate(self, properties: Union[str, List[str]], T: Union[float, np.ndarray] = 298.15, P: Union[float, np.ndarray, str] = 1.0, **kwargs) -> Union[float, np.ndarray, Dict[str, Any]]: """ Calculate water properties using SUPCRT92 model. Parameters ---------- properties : str or list of str Property or list of properties to calculate T : float or array Temperature in Kelvin P : float, array, or 'Psat' Pressure in bar, or 'Psat' for saturation pressure **kwargs Additional options (e.g., Psat_floor) Returns ------- float, array, or dict Calculated properties """ # Handle input types if isinstance(properties, str): properties = [properties] single_prop = True else: single_prop = False # Convert inputs to arrays T = np.atleast_1d(np.asarray(T, dtype=float)) # Always use FORTRAN backend (no fallback) return self._calculate_fortran(properties, T, P, single_prop, **kwargs) def _calculate_fortran(self, properties: List[str], T: np.ndarray, P: Union[np.ndarray, str], single_prop: bool, **kwargs) -> Union[float, np.ndarray, Dict[str, Any]]: """Calculate properties using Fortran backend.""" # Handle pressure input if isinstance(P, str) and P == 'Psat': P_vals = 'Psat' else: P = np.atleast_1d(np.asarray(P, dtype=float)) if len(P) < len(T): P = np.resize(P, len(T)) elif len(T) < len(P): T = np.resize(T, len(P)) P_vals = P # Use batched calculation for better performance try: results = self._fortran_interface.calculate_properties_batch( T, P_vals, properties ) except Exception as e: warnings.warn(f"Batch Fortran calculation failed: {e}") # Fallback to individual calculations if batch fails results = {} for prop in properties: results[prop] = np.full_like(T, np.nan) for i in range(len(T)): T_i = T[i] # Skip invalid points if np.isnan(T_i): continue if isinstance(P_vals, str): P_i = P_vals else: P_i = P_vals[i] if np.isnan(P_i): continue try: # Call Fortran interface props_i = self._fortran_interface.calculate_properties( T_i, P_i, properties ) # Store results for prop in properties: if prop in props_i: results[prop][i] = props_i[prop] except Exception as e: warnings.warn(f"Fortran calculation failed at T={T_i:.1f}K, P={P_i}: {e}") continue # Handle Psat_floor for saturation pressure if 'Psat' in results and 'Psat_floor' in kwargs: Psat_floor = kwargs['Psat_floor'] if Psat_floor is not None: results['Psat'] = np.maximum(results['Psat'], Psat_floor) # Return results if single_prop: result = results[properties[0]] return result[0] if len(result) == 1 else result else: # Convert single-element arrays to scalars if appropriate if len(T) == 1: for key in results: if len(results[key]) == 1: results[key] = results[key][0] return results def available_properties(self) -> List[str]: """ Get list of available properties. Returns ------- List[str] List of available property names """ if self._use_fortran: return self._fortran_interface.available_properties() else: return self._python_backend.available_properties() @property def backend(self) -> str: """Get the active backend ('fortran' or 'python').""" return 'fortran' if self._use_fortran else 'python' @property def has_fortran(self) -> bool: """Check if Fortran backend is available.""" return HAS_FORTRANSUPCRT92 water model with Fortran backend.
This class provides an interface to the original SUPCRT92 Fortran subroutines for calculating water properties. If the Fortran library is not available, it falls back to Python approximations.
The Fortran implementation gives exact compatibility with R CHNOSZ and includes all 23+ thermodynamic properties.
Initialize SUPCRT92 water model.
SUPCRT92 requires the compiled FORTRAN interface for accuracy and compatibility with R CHNOSZ. No pure Python fallback is provided.
Instance variables
prop backend : str-
Expand source code
@property def backend(self) -> str: """Get the active backend ('fortran' or 'python').""" return 'fortran' if self._use_fortran else 'python'Get the active backend ('fortran' or 'python').
prop has_fortran : bool-
Expand source code
@property def has_fortran(self) -> bool: """Check if Fortran backend is available.""" return HAS_FORTRANCheck if Fortran backend is available.
Methods
def available_properties(self) ‑> List[str]-
Expand source code
def available_properties(self) -> List[str]: """ Get list of available properties. Returns ------- List[str] List of available property names """ if self._use_fortran: return self._fortran_interface.available_properties() else: return self._python_backend.available_properties()Get list of available properties.
Returns
List[str]- List of available property names
def calculate(self,
properties: str | List[str],
T: float | numpy.ndarray = 298.15,
P: float | str | numpy.ndarray = 1.0,
**kwargs) ‑> float | numpy.ndarray | Dict[str, Any]-
Expand source code
def calculate(self, properties: Union[str, List[str]], T: Union[float, np.ndarray] = 298.15, P: Union[float, np.ndarray, str] = 1.0, **kwargs) -> Union[float, np.ndarray, Dict[str, Any]]: """ Calculate water properties using SUPCRT92 model. Parameters ---------- properties : str or list of str Property or list of properties to calculate T : float or array Temperature in Kelvin P : float, array, or 'Psat' Pressure in bar, or 'Psat' for saturation pressure **kwargs Additional options (e.g., Psat_floor) Returns ------- float, array, or dict Calculated properties """ # Handle input types if isinstance(properties, str): properties = [properties] single_prop = True else: single_prop = False # Convert inputs to arrays T = np.atleast_1d(np.asarray(T, dtype=float)) # Always use FORTRAN backend (no fallback) return self._calculate_fortran(properties, T, P, single_prop, **kwargs)Calculate water properties using SUPCRT92 model.
Parameters
properties:strorlistofstr- Property or list of properties to calculate
T:floatorarray- Temperature in Kelvin
P:float, array,or'Psat'- Pressure in bar, or 'Psat' for saturation pressure
**kwargs- Additional options (e.g., Psat_floor)
Returns
float, array,ordict- Calculated properties
class WaterModelError (*args, **kwargs)-
Expand source code
class WaterModelError(Exception): """Exception raised for water model calculation errors.""" passException raised for water model calculation errors.
Ancestors
- builtins.Exception
- builtins.BaseException