Solving Chemical Absorption Equilibria using Free Energy and Quantum Chemistry Calculations: Methodology, Limitations, and New Open-Source Software

We developed an open-source chemical reaction equilibrium solver in Python (CASpy, https://github.com/omoultosEthTuDelft/CASpy) to compute the concentration of species in any reactive liquid-phase absorption system. We derived an expression for a mole fraction-based equilibrium constant as a function of excess chemical potential, standard ideal gas chemical potential, temperature, and volume. As a case study, we computed the CO2 absorption isotherm and speciation in a 23 wt % N-methyldiethanolamine (MDEA)/water solution at 313.15 K, and compared the results with available data from the literature. The results show that the computed CO2 isotherms and speciations are in excellent agreement with experimental data, demonstrating the accuracy and the precision of our solver. The binary absorptions of CO2 and H2S in 50 wt % MDEA/water solutions at 323.15 K were computed and compared with available data from the literature. The computed CO2 isotherms showed good agreement with other modeling studies from the literature while the computed H2S isotherms did not agree well with experimental data. The experimental equilibrium constants used as an input were not adjusted for H2S/CO2/MDEA/water systems and need to be adjusted for this system. Using free energy calculations with two different force fields (GAFF and OPLS-AA) and quantum chemistry calculations, we computed the equilibrium constant (K) of the protonated MDEA dissociation reaction. Despite the good agreement of the OPLS-AA force field (ln[K] = −24.91) with the experiments (ln[K] = −23.04), the computed CO2 pressures were significantly underestimated. We systematically investigated the limitations of computing CO2 absorption isotherms using free energy and quantum chemistry calculations and showed that the computed values of μiex are very sensitive to the point charges used in the simulations, which limits the predictive power of this method.

The following items are presented in this Supporting Information: • Derivation of an expression for mole fraction-based reaction equilibrium constants (section S1), • Detailed explanation of the input file for chemical reaction equilibrium solver (section S2), • Details of computing µ 0 i with quantum chemistry calculations and the JANAF tables (sections S3.1 and S3.2), • Details of computing µ ex i using Brick-CFCMC (section S3.3), • Simulation details including all force field parameters (section S4), • Details of accounting for CO 2 evaporation in sequential absorption of CO 2 and next H 2 S (section S5), • Derivation of an expression for the Henry constant of CO 2 in aqueous MDEA solutions (section S6).

S1 Derivation of an Expression for Mole Fraction-based Reaction Equilibrium Constants
Our aim is to derive an expression for the reaction equilibrium constant K as a function of the standard ideal gas chemical potential of species i (µ 0 i ), the excess chemical potential of species i (µ ex i ), and the absolute temperature T . Using the equilibrium condition and the definitions of chemical potential for solutes and the solvent (Eqs. 3 and 4 of the main text), we derive an expression for desired equilibrium constant. The equilibrium condition of reaction j can be expressed as 1 : where N species is the number of species involved in reaction j, µ i is the chemical potential of species i, and ν i,j is the stoichiometric coefficient of species i in reaction j. In our methodology, we assign positive stoichiometric coefficients to the reaction products while signed stoichiometric coefficients are assigned to the reactants. The chemical potential of solutes are calculated using 2 : where ρ i is the number density of the solute i in the solvent, ρ 0 is the reference number density of 1 molecule Å −3 , µ 0 i is the standard state ideal gas chemical potential 2,3 of species i, µ ex i is the excess chemical potential of species i, R is the ideal gas constant, and T is the absolute temperature. The chemical potential of the solvent is computed with ideal gas reference state using 1 : where ρ pure is the number density of the pure solvent and X s is the mole fraction of the solvent s in the solution (X s = N s /N total , N total is the sum of number of molecules of all species in the solution including the solvent). As an example, let us assume a reaction A + S ⇌ B + C where S is the solvent. Using Eq. (S1), Eq. (S2), and Eq. (S3), we can write 1 : For an arbitrary chemical reaction, this can be rewritten as: where N solute is the number of solutes (excluding the solvent) in the reaction, N i is the number molecules of the species i, ν s,j are the stoichiometric coefficient of the solvent in reaction j, and ν total,solute,j is the sum of the stoichiometric coefficients of all solutes in reaction j ( N solute i=1 ν i,j ). We define the equilibrium constant of reaction j as: The left side of the Eq. (S5) can be used to compute the desired equilibrium constant K ′ des as a function of µ 0 i , µ ex i , T , and V .
Using the values of µ 0 i (computed by quantum chemical calculations 2,4 ), µ ex i (computed by Monte Carlo simulations 2,4,5 ), T , V , and ρ pure , we can compute the equilibrium constant K ′ des for any reaction.
The equilibrium constants can also be defined with the number molecules of each species or the mole fraction of each species: where X i is the mole fraction of species i. Note that the summation in Eq. (S8) includes the solvent mole fraction and originates from defining the chemical potential of the species using a pure-liquid reference state 6 : where γ i is the activity coefficient of species i and µ * i is the reference chemical potential of the pure component i in the liquid phase. γ i is incorporated into K j in our calculations and assumed constant. To convert K ′ j to K j , we can use: The desired equilibrium constant K j (K j,des ) as a function of µ 0 , µ ex , T , and V can be computed as: This means that K j,des = K j at equilibrium. K j,des can be constant or it can be solved iteratively for changing values of X s . In our solver, we solve K j,des iteratively for changing values of X s . That means we compute new values of X s and K j,des in every iteration of the solver and this is continued until the difference between the new X s and the old X s no longer changes.

S2 Input File for the Chemical Reaction Equilibrium Solver
In this section, we explain the input file for our chemical reaction equilibrium solver in detail. For this purpose, we used the CO 2 /MDEA/water system as a case study. The reactions (Reactions R1-R4 of the main text) and the mass balance equations (Eqs. 11-14 of the main text) involved in this system are shown in the Methods section of the main text.
Note that the input file should be in the same directory with main.py and functions.py for solver to perform properly. An example input file for our solver is:  The lines in the input file represent the following: • Temperature: The absolute temperature.
• Number of Species: Number of species in the liquid phase.
• C0 (initial guess): Initial guess of the composition in the liquid phase. This is a list of initial concentrations of species in mol dm −3 . In our example, we used a lean solvent (only MDEA and water in the solution) as our initial guess. Note that none of the concentrations in the initial guess should be zero (due to the boundaries we use in our solver), instead, one can input a very low concentration. If any concentration in the initial guess is inputted zero or lower than zero (≤ 0), then it is changed by Although the solver works if the initial guess does not satisfy charge neutrality, we recommend an initial guess that satisfies charge neutrality for quicker results. The solver does not print a warning if charge neutrality is not satisfied by the initial guess and will continue to run.
• Names of species: Names of the species in the liquid phase.
• Charges: Net charges of the molecules/ions in the liquid phase.

S8
• Name of the solvent: The name of the solvent.
• Pure Density of the Solvent: Density of the pure solvent in mol dm −3 . The name of the solvent and the pure density of the solvent are used to compute the desired equilibrium constants of reactions using Eq. (S11).
• muˆ0 species: A list of the values of µ 0 i in kJ mol −1 for the species involved in reactions. Only used if the equilibrium constants are computed using Eq. (S11).
• muˆex species: A list of the values of µ ex i in kJ mol −1 for the species involved in reactions. Only used if the equilibrium constants are computed using Eq. (S11).
• Impose Ptotal and gas composition? (T=True or F=False): Are the total gas pressure and the gas composition imposed in the calculation? If True, the solver assumes an infinite gas phase and the speciations are computed for all Ptotal (in kPa) listed in the next line. If the total gas pressure and gas composition are imposed, no mass balance equation is used for the species in the gas phase since there is mass transfer from the infinite gas phase to the liquid phase.
• Ptotal / [kPa] (only used if Ptotal and gas composition is imposed): The list of total gas pressures. Used only if the total gas pressure and gas composition are imposed to compute the partial pressures of the species in the gas phase.
• Gas phase species: The names of the species in the gas phase.
• Gas phase composition (only used if Ptotal and gas composition is imposed): The composition of the gas phase. The values in this list are normalized so the values sum up to 1.
• muˆex gases / [kJ/mol]: The values of µ ex i for the gas phase species in the solvent.
• Ctotal,gas / [mol/dm3] (Total concentration of the gases in liquid phase): A list of concentrations on the gas phase species in the liquid phase in mol dm −3 . Only S9 used if the total gas pressure and gas composition are not imposed.
• Number of Reactions: Number of reactions in the liquid phase.
• Stoichiometry: The stoichiometric coefficients of all species for each reaction.
• ln(K) for reactions: The mole fraction-based equilibrium constants for each reaction in the liquid phase. In case the desired equilibrium constant should be computed using Eq. (S11), the input should be "QMMC".  can be used to compute molecular partition functions of isolated molecules. Molecular partition functions can be used to compute heat capacities, internal energies, or chemical potentials of species 2,9 . In this section, we explain how to obtain the standard ideal gas chemical potential of species using the Gaussian09 software 10 . For more detail on molecular partition functions, the reader is referred to Refs. 2,9 . The standard ideal gas chemical potential can be computed using the molecular partition function using 2,9 : where µ 0 i is the standard ideal gas chemical potential of species i, R is the ideal gas constant, T is the absolute temperature, q 0,i is the molecular partition function (excluding the translational part) of the molecule with the ground state energy of the molecule is taken as zero 2 , ρ 0 is the reference number density of 1 molecule Å −3 , Λ i is the thermal de Broglie wavelength of molecule i, and D 0,i is the atomization energy of molecule i, which is the energy required to break all bonds in the molecule 11 (D 0,i > 0). The atomization energy of molecule i can be computed using 2,10 : where N atoms,i is the number of atoms in the molecule i, y j is the number of atoms of type j It should be noted that Gaussian09 10 does not print the electronic energies of individual atoms when the energy of a molecule is computed. Therefore, the electronic energies of individual atoms must be computed separately. It is also important that the zero point energy is included in the electronic energy computed by Gaussian09, so the zero point energy should not be subtracted from the electronic energies of the individual atoms (Eq. (S13)).
Also, larger molecules such as MDEA and MEA has many different conformers with different ground state energies. A conformer search must be performed for these type of molecules to obtain the free energies of conformers. Although the conformers with similar free energies can be accounted using a "lumping" procedure 13 , we only use the conformer with minimum free energy since the differences between the free energy of the conformer with minimum free energy and the free energies of the other conformers are large (>> 1 k B T ). The standard S11 state ideal chemical potential of a molecule can be computed using Gaussian09 (with the Freq keyword). However, Gaussian09 uses a different reference state (P 0 = 1 bar) than Brick-CFCMC (ρ 0 = 1 molecule Å −3 ). The standard state ideal gas chemical potential at the Gaussian09 reference state (µ G i ) can be computed using: Also in Eq. (S14), q 0,i is used with the ground state energy of the molecule as a reference.
The value of the term from µ G i (Eq. (S14)) to µ 0 i (Eq. (S12)) can be performed using: As an example, we will compute the value of µ 0 i for water at 313.15 K. We optimized and computed free energy of water at 313.15 K using the G4 method. Below the input file to compute the value of µ 0 i for water using Gaussian09 10 can be found: We show a small part of the Gaussian09 10 output below.
The Gaussian09 10 output shows that the natural logarithm of the term Using Eq. (S13), we computed the atomization energy D 0,water as 916.06 kJ mol −1 . Using Eq. (S14), and Eq. (S15), we compute the value of µ 0 i of water as -965.23 kJ mol −1 at 313.15 K.

S3.2 Computing µ 0 i using the JANAF tables
Thermodynamic data sets such as the JANAF tables can also be used to compute µ 0 i of a molecule. The JANAF tables provide thermodynamic functions and parameters such as the Gibbs free energy, the enthalpy of formation, and heat capacity 7,8 as a function of temperature. As explained in the previous section, we use the ground state energy as the reference state in our calculations, so the Gibbs free energy values should be shifted to the enthalpy at T = 0 K. JANAF tables use the same reference state as Gaussian09 (P 0 = 1 bar). Therefore, we use the same symbol for the standard state ideal gas chemical potential computed using the JANAF tables as the standard state ideal gas chemical potential computed using Note that µ G i should be converted to µ 0 i using Eq. (S15) to use the correct reference state in Brick-CFCMC 2,5 or CASpy. The standard ideal gas chemical potential of a molecule can be computed in terms of the entries in JANAF tables using 2 : where G 0 i is the standard ideal gas Gibbs free energy of molecule i, H 0 i is the standard enthalpy of molecule i, T is the absolute temperature, and T r is the reference temperature S16 (T r = 298.15 K for the JANAF tables). The terms −[G 0 reported in the JANAF tables can be used to compute the terms and in Eq. (S16). The atomization energy D 0,i can also be computed using the JANAF tables using the difference between the enthalpy of formation of the molecule and the enthalpy of formation of the individual atoms 2 : where N atoms,i is the number of atoms in molecule i, y j denotes the number of atoms of type j in molecule i, and ∆ f H 0 is the enthalpy of formation as tabulated in the JANAF tables.
As an example, we compute µ 0 i of water at 313.15 K.
are reported as 188.87 J K −1 mol −1 and −9.904 kJ mol −1 for water, respec- implemented in Brick-CFCMC. The first method is the "probability" route and it uses the probability distribution of the interaction scaling factor λ at λ = 1 and λ = 0 to compute where p(λ i = 1) and p(λ i = 0) are the probabilities of λ = 1 and λ = 0, respectively.
This method requires a flat distribution of (observed) λ and this is obtained using a biasing function 2,5,12,15-17 . A rule of thumb for the flat distribution of λ is that the difference between the maximum and minimum probabilities should be lower than 20%. The second method is thermodynamic integration 5 . In thermodynamic integration, we use the average value of the derivative of the potential energy with respect to λ, ∂U ∂λ , and compute µ ex i using 4,5 : It is very challenging to obtain a flat probability distribution of λ in a single simulation for large and/or polar molecules 5 . As mentioned in the methodology section of the main text, we need to compute µ ex i of ionic and/or large molecules such as the hydronium ion (H 3 O + ) and the protonated MDEA ion (MDEAH + ). Using thermodynamic integration, a S18 flat probability distribution of λ is not required since ∂U ∂λ term can be computed from independent MC simulations at different and fixed λ values. Therefore, the thermodynamic integration as implemented in Brick-CFCMC is used to compute µ ex i in this study. The term ∂U ∂λ can only be computed for one charge-neutral fractional group 2,5 . The fractional group can consist of multiple molecules or ions. In this study, we computed the µ ex i of the reactants and reaction products of MDEAH + dissociation reaction. As the fractional group should be charge-neutral, an HCO − 3 ion is added to the fractional groups.

S4 Simulation Details
For thermodynamic integration, we used 50 equidistant and fixed values of λ with two additional simulations at λ = 10 −6 and λ = 1 − 10 −6 to increase the accuracy of thermodynamic integration. In these simulations, a simulation box of 300 water molecules and a single fractional group (see Table S1) was used. Initial configurations for these simulations were created S19 by randomly inserting molecules to a simulation box with a length of 20.8 Å. The random insertion of molecules into the simulation box causes atomic overlaps and these overlaps were eliminated by 10 3 initialization cycles in which only single molecule translation and single molecule rotation trial moves were performed. The initialization cycles were followed by For the computation of µ ex i of all species, the TIP3P 22 force field was used for (rigid) water molecules. We used the TIP3P 22 force field for water because µ ex i of water computed using TIP3P force field agrees much better with µ ex i measured using empirical data than the in Gaussian09 10 at HF/6-31G* level of theory. The electrostatic potentials computed by Gaussian09 were then fitted with Antechamber package 28 using a two-step Restrained Electrostatic Potential Surface (RESP) fitting method. The RESP fitted point charges were used as the point charges of GAFF. The parameters for the OPLS-AA force field with 1.14*CM1A point charges were generated using the LibParGen web server 29 . We used the TraPPE 30 and the force field from Kristóf and Lizsi 31 for CO 2 and H 2 S molecules, respectively. All force field parameters used in this study are listed in tables S1-S17. The force field parameters for unlike LJ interactions were computed using Lorentz-Berthelot mixing rules 18 (Table S6).
For flexible molecules (MDEA and MDEAH + ), the bond lengths are fixed. For flexible molecules, the bending potential was computed using 2 : where K is the bending constant, θ is the bending angle, θ 0 is the bending angle at equilibrium. For the GAFF 24 , the torsion potential was computed using 2,33 : where p 0 ..p 5 are the torsion constants and ϕ is the torsion angle. For the OPLS-AA force field, the torsion potential was computed using 2,25,26 : where K 1 ..K 4 are the torsion constants and ϕ is the torsion angle.         Table S15: Bending potential parameters, the equilibrium angle (θ 0 ) and the bending constant (K) of the MDEA molecule for the GAFF 24 . The bending potential was computed using Eq. (S20). Fig. S1 shows the atom types designation of MDEA for the GAFF 24 .
Bending  Figure S2: Schematic representation showing atom type designations of the MDEAH + ion for the GAFF 24 .   Table S19: Bending potential parameters, the equilibrium angle (θ 0 ) and the bending constant (K) of the MDEAH + ion for the GAFF 24 . The bending potential was computed using Eq. (S20). Fig. S2 shows the atom types designation of MDEAH + for the GAFF 24 .
Bending  Table S20: Torsion potential parameters of the MDEAH + ion for the GAFF 24 . The torsion potential was computed using Eq. (S21). Fig. S2 shows the atom types designation of MDEAH + for the GAFF 24 . The values of p 4 and p 5 are zero for all torsions in the MDEAH + ion.
Torsion  Figure Table S23: Bending potential parameters, the equilibrium angle (θ 0 ) and the bending constant (K) of the MDEA molecule for the OPLS-AA force field 25,26 . The bending potential was computed using Eq. (S20). Fig. S3 shows the atom types designation of MDEA for the OPLS-AA force field 25,26 .
Bending  Table S24: Torsion potential parameters of the MDEA molecule for the OPLS-AA force field 25,26 . The torsion potential was computed using Eq. (S22). Fig. S3 shows the atom types designation of MDEA for the OPLS-AA force field 25,26 .
Torsion  Figure S4: Schematic representation showing atom type designations of the MDEAH + ion for the OPLS-AA force field 25,26 . Table S25: Non-bonded interaction parameters for MDEAH + . For Lennard-Jones (LJ) interactions, the OPLS-AA force field 25,26 was used. For electrostatic interactions, 1.14*CM1A point charges 35 were used. The point charges of the MDEAH + ion sum up to -1. Fig. S4 shows the atom types designation of MDEAH + for the OPLS-AA force field 25,26 .   Table S27: Bending potential parameters, the equilibrium angle (θ 0 ) and the bending constant (K) of the MDEAH + ion for the OPLS-AA force field 25,26 . The bending potential was computed using Eq. (S20). Fig. S4 shows the atom types designation of MDEAH + for the OPLS-AA force field 25,26 .
Bending  Table S28: Torsion potential parameters of the MDEAH + ion for the OPLS-AA force field 25,26 . The torsion potential was computed using Eq. (S22). Fig. S4 shows the atom types designation of MDEAH + for the OPLS-AA force field 25,26 .
Torsion   36 , we have one more variable additional to the concentrations of the species in liquid phase (N i ), which is the amount of CO 2 evaporated to the gas phase (N CO 2 ,gas ). As a result of this, the CO 2 balance in the system has an additional term N CO 2 ,gas . The CO 2 balance in the H 2 S/CO 2 /MDEA/water system changes from Eq. 9 of the main text to: We also have an additional equation derived using the chemical equilibrium between the free CO 2 absorbed in the liquid solution and the CO 2 evaporated to the gas phase and the equilibrium condition (Eq. 1 of the main text): The chemical equilibrium between the free CO 2 absorbed in the liquid solution and the CO 2 evaporated to the gas phase can be shown as: which results in: where N CO 2 ,liquid represent the number of molecules of free CO 2 in liquid phase, V liquid is the volume of the liquid phase, V gas is the volume of the gas phase, and µ ex CO 2 ,liquid is the excess chemical potential of CO 2 in liquid phase. Note that we used a V liquid Vgas value of 0.3 in these calculations to mimic the conditions in the experiments by Dicko et al. 36 . In total, we have 12 variables and 12 equations to solve in computing the speciation of the H 2 S/CO 2 /MDEA/water system while accounting for the effect of evaporating CO 2 during sequential H 2 S absorption. We solved the speciation in this system using a numerical least squares solver for nonlinear equations. Fig. S5   CO 2 ) can be expressed as: where P CO 2 is the partial pressure of CO 2 in the gas phase and X CO 2 ,total is the total mole fraction of the free CO 2 , HCO − 3 , and CO 2− 3 (X CO 2 ,total = X CO 2 + X HCO − 3 + X CO 2− 3 ). There is already an expression to compute P CO 2 (Eq. 7 of the main text), so we need to derive an expression for X CO 2 ,total . As P CO 2 approaches 0, we assume that the solution is only composed of water and MDEA (X MDEA + X H 2 O = 1) and the net charge of the OH − ion is the only negative charge in the solution that balances the positive charge from the MDEAH + ion (X MDEAH + = X OH − ) i.e. MDEA is a weak base. By using the equilibrium constant of reaction R4 (K R4 ) and replacing X OH − with X MDEAH + , we obtain an expression for the mole fraction of the H 3 O + ion: By replacing the term X H 3 O + with Eq. (S27) in the equilibrium constant of reaction R3 (K R3 ), an expression to compute the mole fraction of MDEAH + is obtained: The mole fraction of free CO 2 in the solution can be computed using Eq. 7 of the main text and total number of molecules (N MDEA + N H 2 O ) as: where V is the volume of the liquid phase, µ ex CO 2 is the excess chemical potential of CO 2 , k B is the Boltzmann constant, and T is the absolute temperature. By using the equilibrium constant of reaction R1 (K R1 ), we obtain an expression to compute the mole fraction of HCO − 3 which is: By using the equilibrium constant of reaction R2 (K R2 ), we can compute the mole fraction of CO 2− 3 as: By summing Eq. (S29), Eq. (S30), and Eq. (S31) up, an expression for X CO 2 ,total is obtained: Because we show the CO 2 pressure as a function of the CO 2 loading in our isotherms (see Fig. 5 of the main text), we prefer to compute the Henry coefficient of CO 2 as: where α CO 2 ,total is the total loading of CO 2 in the solution in the units of mol where α CO 2 , α HCO − 3 , and α CO 2− 3 are the loadings of free CO 2 , HCO − 3 , and CO 2− 3 , respectively). The mole fractions computed using Eq. (S29),

S42
Eq. (S30), and Eq. (S31) can be converted to loading α i using: where K H i is the Henry constant of species i and K 0 is the arbitrary reference Henry constant to make the argument of the logarithm dimensionless. However, we feel that computing the S43 heat of absorption is out of the scope of this work.