Modeling Differential Enthalpy of Absorption of CO2 with Piperazine as a Function of Temperature

Temperature-dependent correlations for equilibrium constants (ln K) and heat of absorption (ΔHabs) of different reactions (i.e., deprotonation, double deprotonation, carbamate formation, protonated carbamate formation, dicarbamate formation) involved in the piperazine (PZ)/CO2/H2O system have been calculated using computational chemistry based ln K values input to the Gibbs–Helmholtz equation. This work also presents an extensive study of gaseous phase free energy and enthalpy for different reactions using composite (G3MP2B3, G3MP2, CBS-QB3, and G4MP2) and density functional theory [B3LYP/6-311++G(d,p)] methods. The explicit solvation shell (ESS) model and SM8T solvation free energy coupled with gaseous phase density functional theory calculations give temperature-dependent reaction equilibrium constants for different reactions. Calculated individual and overall reaction equilibrium constants and enthalpies of different reactions involved in CO2 absorption in piperazine solution are compared against experimental data, where available, in the temperature range 273.15–373 K. Postcombustion CO2 capture (PCC) is a temperature swing absorption–desorption process. The enthalpy of the solution directly correlates with the steam requirement of the amine regeneration step. Temperature-dependent correlations for ln K and ΔHabs calculated using computational chemistry tools can help evaluate potential PCC solvents’ thermodynamics and cost-efficiency. These correlations can also be employed in thermodynamic models (e.g., e-UNIQUAC, e-NRTL) to better understand postcombustion CO2 capture solvent chemistry.


Description
Pages Gaseous Phase Gibbs free energy of piperazine species, using G3MP2B3, G3MP2, G4MP2, CBS-QB3 and DFT level of theories, at 298 K. S3 Gaseous Phase enthalpy of piperazine species, using G3MP2B3, G3MP2, G4MP2, CBS-QB3 and DFT level of theories, at 298 K. S3 Gaseous Phase Gibbs free energy and enthaply of bicarbonate (HCO 3 -), water (H 2 O), CO 2 and H 3 O + using G3MP2B3, G3MP2, G4MP2, CBS-QB3 and DFT level of theories, at 298 K. S4 Underlying data for Table 2 for piperazine species: energy of solute and cluster in gas phase, thermal corrections to the energy and entropy of solute and cluster. S4 Quantitative details of optimized structures of the first solvation shell from molecular simulations for a set of 100 initial geometries of Piperazine. S5 -S7 Molecular Geometries of ESS clusters for various Piperazine species S8 -S12 Simulation Details S13 Explanation of QM/PB Continuum Solvent Model S14 -S15 Constant terms utilized in calculations S16 Full citation for References used in manuscript S16 -S17     19.120977000 11.935583000 8.491351000

Simulation Details
The molecular dynamics (MD) simulations were isothermal-isobaric simulations (nPT) with periodic boundary conditions at 298K and 1 bar. The particle-mesh Ewald procedure was used to handle long-range electrostatic interactions. Temperature was controlled by Langevin dynamics, while pressure was controlled by weak coupling to an external bath. The time-step in the simulations was set to 0.002 ps. Solute molecule geometries were fully constrained during simulations. Solute coordinates were restrained with a harmonic potential. The restraining weight was set at 5 kcal/mol Å2. The nonbonded cutoff was set to 8 Å for all simulations. Bonds length involving hydrogen atoms were maintained by use of a SHAKE algorithm. The systems were equilibrated in nPT simulations for 800 ps before geometries were extracted. Cluster geometries were extracted from the MD trajectory every 2 ps. These clusters consisted of the solute and the 5 closest solvent molecules.
Water solvent was represented with the TIP3P force field1. The number of solvent molecules in simulations for a given solute was taken to be 100. The solutes geometries were obtained from HF/6-31+G (d) calculations in vacuum. Solute charges were calculated with the CM2 model2.
The force field Lennard-Jones parameters were drawn from force fields reported in the scientific literature. We selected the force fields that appeared to have the most detailed parameterization for a given solute. For many ionic species no special force field have been developed, and is these cases we drew on atomic parameters for neutral species. For atoms where no force field stood out in terms of parameterization we drew on the GAFF force field3. GAFF would be our recommended default when no solute specific force field is available. Below are given the Lennar d Jones parameters utilized in the present work.

QM/PB Continuum Solvent Model
In the PB models, the solution is characterized as a charge distribution of the solute in a cavity mimicking the molecular shape which surrounded by a continuous dielectric representing the solvent. The solvent is polarized by the solute and the solvent polarization generates an electrostatic field called the reaction field. The electrostatic potential of the reaction field can be described by the Poisson-Boltzmann equation: The (r) is the charge distribution of the solute calculated via semi empirical QM methods. The (r) is the electrostatic potential that is to be obtained. The (r) is the dielectric function describing the dielectric discontinuity and it usually only has two values: out for the region outside of the cavity and in for the space within the (molecular) cavity. The (r) is a modified Debye-Huckel parameter that reflects the salt concentration and temperature. When the charge density on the solute is low and the ionic strength is low, the term can be approximated via a linear term yielding the so-called linear Poisson-Boltzmann equation: Even for the linearized PB equation, analytical solutions can only be obtained for systems defined by a simple dielectric boundary, e.g., a sphere. For a boundary as complex as that of a protein or DNA molecule, the finite difference approach (among others) substitutes for the analytical differential equation solution yielding a numerical solution of the PB equation.
The electrostatic potential of the reaction field generated by a fixed set of (r) will, in turn, polarize the solute charge distribution and generate a new '(r) (i.e. distort the gas phase wave function). Thus the gas phase solute Hamiltonian, H 0 is perturbed by a potential energy operator coming from the interaction between the solute and the reaction field: The equation 4 is an integral over all the virtual surface charges, and r' is the coordinates of the surface charges. The virtual surface charges (σ) are calculated from the converged electrostatic potential of the reaction field. And the full set of virtual surface charges should generate an electrostatic field which is identical to the reaction field. Thus, the perturbation (Vint) due to the reaction field could be accurately calculated by equation 4 The G rf is the electrostatic interaction of the solute charge distribution (including core and electrons) with the electrostatic potential generated by the reaction field. It is more efficient and accurate to calculate the G rf using the virtual surface charge, (see equation 6), instead of using the electrostatic potential directly. = * + (7) Figure 1 summarizes the origin of these three energy components of the solvation free energy The components of the solvation free energy in the QM/PB model. G wfd is the penalty originating from the polarization of the solute's wave function in vacuum relative to the wave function in solution. G np models non-polar effects, including the entropy to make the molecular cavity in the solvent. G rf is the electrostatic interaction between the polarized solute and the surrounding solvent.