How to calculate the activity coefficient using CHNOSZ?

[Zqice]

Hello everyone,

I am currently trying to reproduce the results from the following paper:

Wu, H., Liu, L., Dick, J. M., and Liu, H. (2023)
Speciation and Thermodynamic Properties of Pb–HS Complexes in Sulfur-Rich Hydrothermal Fluids: Insights from Ab Initio Molecular Dynamics Simulations
ACS Earth and Space Chemistry, 7(12), 2393–2409
DOI: 10.1021/acsearthspacechem.3c00183

I am having difficulty reproducing the part of the paper where the activity-coefficient correction was calculated using CHNOSZ.

In short, the paper used distance-constrained ab initio molecular dynamics (AIMD) combined with thermodynamic integration (TI) to obtain Gibbs free energies for a series of ligand-dissociation reactions of Pb–HS complexes.

My understanding is that these Gibbs free energies are the uncorrected values (i.e., before activity correction). The authors then applied the B-dot extension of the Debye–Hückel theory to calculate activity coefficients, and used these to obtain the activity-corrected Gibbs free energies. According to the paper, this activity-coefficient calculation was done using CHNOSZ, and this is the part I am currently trying to reproduce.

More specifically, for a simulation box with the composition:

Pb(HS)3^- + HS^- + 2 Na^+ + 55 H2O

one of the reactions considered is:

Pb(HS)3^- + H2O = Pb(HS)2(H2O) + HS^-

If I label the reactants and products as A, B, C, and D, then the corrected Gibbs free energy can be calculated as:

ΔrG° = ΔrG + RT ln[(C_A γ_A · C_B γ_B) / (C_C γ_C · C_D γ_D)]

I selected one specific condition from the paper: 200 °C and 20 bar.

At this condition, the paper reports:

Uncorrected Gibbs free energy: −18.7 ± 7.5 kJ/mol

Activity-corrected Gibbs free energy: −17.1 ± 7.5 kJ/mol

To calculate the activity coefficients in CHNOSZ for this system, I used the following R code

library(CHNOSZ)

# 1) Select the non-ideal model (Bdot is also the default, but I specify it explicitly here)
nonideal("Bdot")

# 2) Set temperature, pressure, and ionic strength
T <- 200
P <- 20
IS <- 2.02

# 3) Calculate activity coefficients for the species
sp <- c("Pb(HS)3-", "Pb(HS)2", "HS-", "Na+")
sres <- subcrt(sp, T = T, P = P, IS = IS)

out <- lapply(seq_along(sp), function(i) {
  data.frame(
    species = sp[i],
    loggam = sres$out[[i]]$loggam,
    gamma  = 10^(sres$out[[i]]$loggam)
  )
})
do.call(rbind, out)

This gives the following results:

  species      loggam     gamma
1 Pb(HS)3- -0.26463477 0.5437074
2  Pb(HS)2 -0.01552359 0.9648869
3      HS- -0.32743155 0.4705096
4      Na+ -0.29338742 0.5088767

However, when I substitute these values into the equation above to calculate ΔrG°, I obtain −18.0 kJ/mol, which does not match the value reported in the paper (−17.1 ± 7.5 kJ/mol).

My Python script for the calculation is:

import math

def calc_deltaG_theta(deltaG, T, C_A, C_B, C_C, C_D, gamma_A, gamma_B, gamma_C, gamma_D):
    """
    Calculate:
    ΔrG° = ΔrG + RT ln[(C_A*γ_A * C_B*γ_B) / (C_C*γ_C * C_D*γ_D)]

    Parameters
    ----------
    deltaG : float
        Input ΔrG in kJ/mol
    T : float
        Temperature in K
    C_A, C_B, C_C, C_D : float
        Concentrations of A, B, C, D
    gamma_A, gamma_B, gamma_C, gamma_D : float
        Activity coefficients of A, B, C, D

    Returns
    -------
    deltaG_theta : float
        Calculated ΔrG° in kJ/mol
    """
    R = 8.314462618  # J/(mol·K)

    numerator = (C_A * gamma_A) * (C_B * gamma_B)
    denominator = (C_C * gamma_C) * (C_D * gamma_D)

    if numerator <= 0 or denominator <= 0:
        raise ValueError("The argument of the logarithm must be greater than 0.")

    correction_J = R * T * math.log(numerator / denominator)
    correction_kJ = correction_J / 1000.0

    deltaG_theta = deltaG + correction_kJ
    return deltaG_theta

deltaG = -18.7   # kJ/mol
T = 473.15       # K

C_A = 1.009244
C_B = 1.0
C_C = 1.009244
C_D = 1.009244

gamma_A = 0.5437074
gamma_B = 1.0
gamma_C = 0.9648869
gamma_D = 0.4705096

result = calc_deltaG_theta(deltaG, T, C_A, C_B, C_C, C_D, gamma_A, gamma_B, gamma_C, gamma_D)
print(f"ΔrG° = {result:.6f} kJ/mol")

So my question is:

Am I using CHNOSZ incorrectly, or could there be another reason why my calculated ΔrG° does not match the value reported in the paper?

I would greatly appreciate any guidance or suggestions. Thank you very much!

[jedick]

Hi! Thanks for your question. It looks like you are calculating the activity coefficients correctly. However, the calculation of ΔG° should only use the activity coefficients of species in the reaction, not other species like Na⁺ (although all ionic species contribute to the ionic strength).

I tried this to calculate ΔG° for the reaction Pb(HS)₂ + HS^- = Pb(HS)₃^-. Because - (the minus sign) can't be used in variable names, I used _ instead (e.g. gamma_HS_). Here I used rounded values of 1 mol/kg for concentrations and 2 mol/kg for ionic strength (adapted from your code).

# Non-standard Gibbs energy of reaction (J/mol)
Delta_G <- -18700
# Ionic strength used for calculating activity coefficients (mol/kg)
IS <- 2
# Temperature and pressure
T <- 200
P <- 20
# Calculate activity coefficients of charged species
# NOTE: loggam is the base-10 logarithm of activity coefficient
gamma_HS_ <- 10^subcrt("HS-", T = T, P = P, IS = IS)$out[[1]]$loggam
gamma_PbHS3_ <- 10^subcrt("Pb(HS)3-", T = T, P = P, IS = IS)$out[[1]]$loggam
# Use one for activity coefficient of uncharged species
gamma_PbHS2 <- 1
# Use one molal concentration for all species
C_HS_ <- 1
C_PbHS2 <- 1
C_PbHS3_ <- 1
# Gas constant (J/K/mol) and temperature (K)
R <- 8.314
TK <- T + 273.15
# Calculate Delta G0
numerator <- (C_PbHS2 * gamma_PbHS2) * (C_HS_ * gamma_HS_)
denominator <- (C_PbHS3_ * gamma_PbHS3_)
# NOTE: log() is the natural logarithm
Delta_G0 <- Delta_G + R * TK * log(numerator/denominator)
print(round(Delta_G0))

This gives -19267 J/mol for ΔG° of the reaction. That's different from both the value you calculated (-18.0 kJ/mol) and the value in the paper (-17.1 kJ/mol). Increasing the ionic strength to 4 mol/kg (to account for 2 mol/kg Na⁺ ions in the simulation box) pushes the ΔG° more negative (-19395 J/mol).

I checked that the above code uses values of gamma_HS_ and gamma_PbHS3_ that are very close to the ones you calculated. So the issue probably isn't with the calculation of activity coefficients, but how they're used in the equation to calculate ΔG°.

At the moment I don't have more insight into the calculations used for the paper. Perhaps one of the other coauthors would be able to help.

[Zqice]

Dear Prof. Dick,
Thank you very much for your detailed reply and for checking this so carefully.
Your explanation is very helpful. I now understand that the main issue is probably not the calculation of the activity coefficients themselves, but rather how they are used in the equation to calculate ΔG°. I will look into this part more carefully.
Thank you again for your time and kind help.
Best regards,
Qi Zhang