Hybrid Computational Strategy for Predicting CO2 Solubilities in

Subscriber access provided by Kaohsiung Medical University

Article

Hybrid Computational Strategy for Predicting CO Solubilities in Reactive Ionic Liquids 2

Quintin R. Sheridan, Ryan G. Mullen, Tae Bum Lee, Edward J. Maginn, and William F. Schneider J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.8b02095 • Publication Date (Web): 11 May 2018 Downloaded from http://pubs.acs.org on May 16, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Hybrid Computational Strategy for Predicting CO2 Solubilities in Reactive Ionic Liquids Quintin R. Sheridan, Ryan Gotchy Mullen, Taebum Lee, Edward J. Maginn,∗ and William F. Schneider∗ Department of Chemical and Biomolecular Engineering, The University of Notre Dame Notre Dame, Indiana, 46556 USA E-mail: [email protected]; [email protected]

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 31

Abstract A combination of ab initio calculations and classical molecular dynamics simulations was used to calculate the free energy of reacting an aprotic heterocyclic anion ionic liquid with CO2 . The overall reaction was broken into a series of steps using a thermodynamic cycle to calculate the free energy of the gas phase reaction and the free energy contributions of solvation environment effects, which make comparable contributions to the total free energy of reaction. CO2 absorption isotherms that agree reasonably well with experimental data were calculated using a derived expression for the free energy of reaction as a function of temperature, pressure, and the extent of reaction.

Introduction Ionic liquids (ILs) have gained much attention in recent years as solvents that could potentially outperform aqueous amines in CO2 separations. 1–5 ILs offer a wide range of chemical tunability and often have negligible vapor pressure and high thermal stability, all of which make them excellent candidate solvents for CO2 separations. 5–7 Recently, a class of aprotic heterocyclic anion (AHA) ILs have been synthesized that react stoichiometrically with CO2 : 8–12

[C]+ [A]− + CO2 −−→ [C]+ [A−CO2 ]−

(1)

significantly enhancing the CO2 solubility over non-reacting, physical solvents. ILs can be prepared from many combinations of anions and cations, which makes it difficult to screen candidate ILs experimentally. It is thus desirable to have computational methods that are able to predict CO2 solubilities in CO2 -reactive ILs prior to synthesis and experimental evaluation. Methods for computing the physical solubility of CO2 in ILs are well developed. The physical interactions between a CO2 molecule and the components of an IL are generally 2


Page 3 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


well described by appropriately parameterized classical force fields. Solubilities can be computed indirectly from classical molecular dynamics (CMD) trajectories using free energy methods to extract excess chemical potentials and Henry’s Law constants. 13,14 Physical absorption isotherms can be calculated directly using Gibbs ensemble Monte Carlo simulations to determine equilibrium compositions of liquid and vapor phases over a range of thermodynamic states. 15–17 When CO2 chemically reacts with an IL, however, solubilities become much more difficult to predict, as the reaction equilibria involves chemical bonding that is best described by quantum chemical methods. The first computational attempts to determine relative CO2 solubilities in CO2 -reactive ILs were based on the reaction of a single isolated anion with CO2 . 8–10,18–21 The gas-phase enthalpy of reaction was used to rank the relative CO2 solubilities as a function of the IL anion. While correlations with observation are generally good, there are exceptions. 8 Further, anion-only models are by definition insensitive to cation choice, ignore the influence of ion solvation on CO2 binding energies, and provide no information about temperatureand composition-dependent free energies of absorption. Further, they are unable to capture any secondary chemistry involving the anions and cations together. 22,23 A number of approaches have been used to incorporate the effect of the IL solvent on CO2 solubilities in reactive ILs. Firaha et al. used the COnductor-like Screening MOdel for Realistic Solvents (COSMO-RS) solvation model to determine the geometries and free energies of implicitly solvated anions with and without CO2 and thus to compute a free energy of CO2 absorption. 24 The COSMO-RS predicted reaction free energies were in qualitative agreement with experimental CO2 solubilities and in many cases quantitative agreement with free energies obtained by fitting isotherm models to experimental data. Implicit solvation models, however, can have difficulty describing solvents with micro-heterogeneous charged regions that influence solubility. 25–27 Further, they provide no information about composition-dependent absorption. The deficiencies of an implicit solvation model can in principle be overcome by including

3



more than a single ion within a quantum mechanical framework. Such explicit solvation models can capture specific ion-ion interactions and their influence on CO2 binding, 28 but translating these results into a full, composition-dependent isotherm is difficult. The IL environment can have an influence on CO2 binding, and such influence must be averaged over ion conformations and across the composition domain. While such sampling is readily accomplished within a classical framework, full quantum mechanical evaluation of these effects, through first-principles molecular dynamics, is prohibitively expensive. In general, the equilibrium fraction of chemically bound CO2 in an IL is determined by the free energy of reaction, ∆G(T, PCO2 , ξ), where T is temperature, PCO2 is the partial pressure of CO2 , and ξ is an extent of reaction, or fraction of anions bonded to CO2 . Free energy is a state function; it can be evaluated through an appropriately constructed thermodynamic cycle. 29–31 In this work, we evaluate a hybrid computational approach to predict ∆G(T, PCO2 , ξ), employing a thermodynamic cycle that combines the advantages of classical and quantum dynamics methods. Large-scale sampling of the bulk IL across the composition space is performed using classical molecular dynamics, the free energy of reaction of an ion-pair with CO2 is computed using a first-principles-computed potential of mean force, and these two steps are connected through classical molecular dynamics evaluations of IL vaporization free energies. The method is thus fully predictive. We demonstrate the method for the absorption of CO2 into an IL comprised of triethyl-octyl-phosphonium cations ([P2228 ]+ ) paired with 2-cyanopyrrolide anions ([2CNpyr]− ), as this is one of the most widely studied CO2 -reactive ILs (Figure 1). 8,9,32 We compute full CO2 absorption isotherms at 300, 333, and 360 K. We find that the gas phase reaction and solvation free energies make comparable contributions to the total free energy of reaction, that the unreacted and reacted IL act like a nearly ideal mixture, and that the differences in vaporization free energies are well described by the differences in vaporization enthalpies.

4


Page 4 of 31

Page 5 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


[P2228]+

[2CNpyr]-

[2CNpyr:CO2]-

Figure 1: Schematic illustrations of the triethyl-octyl-phosphonium cation, the 2cyanopyrrolide anion, and the 2-cyanopyrrolide carbamate.

Method Thermodynamic Cycle Figure 2 shows the overall thermodynamic cycle describing the absorption of CO2 into [P2228 ][2CNpyr] due to reaction of CO2 with the [2CNpyr]− anion. The free energy of the entire cycle is given by

g IL:CO2 ∆G(T, PCO2 , ξ) = ξ ∆GIL (T, P ) + ∆Glmix (T, ξ) (2) vap (T, P ) + ∆Grxn (T, PCO2 ) − ∆Gvap

where ∆G(T, PCO2 , ξ) is the composition-dependent free energy of chemical absorption of CO2 (Step 1), ∆GIL vap (T, P ) is the free energy to vaporize an unreacted cation/anion pair from the bulk IL to an ideal gas (Step 2), ∆Ggrxn (T, PCO2 ) is the free energy of reacting (T, P ) is the free energy a cation/anion pair with CO2 in the vacuum (Step 3), ∆GIL:CO2 vap to vaporize a cation/anion pair from the CO2 -reacted bulk IL to an ideal gas (Step 4), and ∆Glmix (T, ξ) is the free energy of mixing unreacted and CO2 -reacted IL to the desired composition (Step 5). Because the IL is essentially incompressible, we neglect in these expressions the pressure-dependence of the free energy of the bulk phases. An isotherm is defined by the value of ξ that minimizes ∆G(T, PCO2 , ξ) at a given T and PCO2 . Steps involving bulk IL, including Steps 2, 4, and 5, are directly amenable to classical simulation methods, while Step 3, linking unreacted and reacted IL, demands a first-principles treatment. In the following subsections, the methods used to calculate the free energy con-

5



Page 6 of 31

tribution of each piece of the thermodynamic cycle along with the assumptions used in each calculation are described. Further computational details are provided in the subsequent Simulation Details section. IL Pair

CO2

IL:CO2

(3)

(4) Liquid IL

(2)

(2’)

Liquid IL:CO2

(5) Liquid IL

CO2

Liquid Mixture

ΔGrxn

(1) Vapor Liquid

[P2228]+

[2CNpyr]-

CO2

[2CNpyr:CO2]-

Figure 2: Thermodynamic free energy cycle for the reaction of CO2 in IL (1). Free energies are determined based upon the following steps: (2) vaporization of IL at constant T and P; (2’) a portion of the unreacted IL remains liquid and is mixed with reacted IL in step 5 (3) reaction of IL with CO2 in the gas phase; (4) condensation of the CO2 -reacted IL vapor at constant T and P; (5) mixing of the unreacted IL with the CO2 -reacted IL.

Vaporization and Condensation: Steps 2 and 4 The free energy changes in Steps 2 and 4 correspond to the free energy of vaporization/condensation of the unreacted and CO2 -reacted IL, respectively. The difference in the free energies of these two steps enters linearly into the cycle with the extent of reaction and therefore can be com-

6


Page 7 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


bined into one quantity:

IL:CO2 ∆∆Gvap (T ) = ∆GIL (T, P ) vap (T, P ) − ∆Gvap

(3)

Because ion pairs from both the unreacted and reacted ILs are vaporized to an ideal gas, the pressure-dependence of the vaporization difference is negligible. ∆∆Hvap (T ) and ∆∆Svap are defined similarly.

Gas-Phase Reaction: Step 3 An ion pair is the smallest charge neutral unit that can be extracted from the IL, the smallest chemical unit that can capture specific interactions between the reactive IL center and a charge-balancing cation, and is sufficiently small that it can be readily treated by firstprinciples methods. For these reasons, Step 3 is constructed as the free energy of reacting an isolated cation/anion pair with CO2 :

[P2228 ][2CNpyr](g) + CO2 (g) −−→ [P2228 ][2CNpyr−CO2 ](g)

(4)

This free energy includes contributions from the new chemical bond and from changes in the interaction of cation and anion, all of which must be sampled over the various degrees of freedom of the cation/anion pair. To perform this sampling, we used a potential of mean force (PMF) constraint force method, 33,34 in which a sequence of constant-temperature ab initio molecular dynamics (AIMD) simulations are performed at various constrained separations between the CO2 molecule and the anion reaction site. The separation between the CO2 carbon and the [2CNpyr]− nitrogen is a natural choice of reaction coordinate, and was held fixed at each separation by applying the necessary constraint force. The average constraint force was sampled at constraint lengths that span from the bound state to free CO2 . The work of separation is related to the average constraint force by

7



Page 8 of 31

dW (r) = −hF (r)i dr

(5)

where W (r) is the PMF and hF (r)i is the average value of the constraint force at a given constraint length, r. The free energy of the ion pair reaction is " ∆G◦rxn (T ) = −RT ln 4π

Z

r∗

r2 e−β[W (r)−W (rD )] dr

0

P◦ kB T

−∆ν # (6)

where r∗ is the distance cutoff beyond which CO2 and the ion pair are taken as unbound and rD is the longest distance sampled in the PMF calculation, at which W(rD ) and its first derivative are both zero. Results are insensitive to the exact choice of r∗ , as the integral is dominated by the minimum in W (r) and converges for all values of r∗ beyond which W (r) is non-negative. The complete derivation of eq 6 is shown in the Supporting Information (SI). In more conventional applications of first-principles methods, the free energies of each species are computed within the analytical ideal gas/rigid rotor/quantum harmonic oscillator approximations. These approximations are expected to work well for stiff modes bit to perform less well for low-frequency and anharmonic modes. As a result, these approximations are challenged to describe cation/anion pairs. For comparison purposes, we calculated the free energy of reaction 4 for the even simpler anion-only reaction:

[2 CNpyr]− (g) + CO2 (g) −−→ [2CNpyr − CO2 ]− (g)

(7)

using the analytical approximations.

Liquid Mixing: Step 5 The free energy of mixing

l l ∆Glmix (T, ξ) = ∆Hmix (T, ξ) − T ∆Smix (T, ξ)

8


(8)

Page 9 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


captures the nonlinear composition dependence of the overall reaction free energy. Mixing is a bulk property, evaluated most directly from CMD simulations on mixtures of various compositions. While the mixing free energy is difficult to compute, the enthalpy of mixing and the volume of mixing can be readily extracted from averages computed for MD trajectories over a range of compositions. While the entropy of mixing is also difficult to compute, it can be computed analytically for ideal solutions as a function of the extent of reaction

l,ideal ∆Smix (ξ) = −R [ξlnξ + (1 − ξ)ln(1 − ξ)] .

(9)

Ideal solutions are mixtures for which the intensive properties are the sum of the product of the component properties multiplied by the component mole fractions. The validity of the ideal solution model is confirmed in the results section by showing that the enthalpy of mixing and the volume of mixing are negligible.

Simulation Details CMD Simulations The enthalpy of vaporization, enthalpy of mixing, and volume of mixing were calculated from separate CMD simulations of the bulk IL and vapor phases. All of the CMD simulations were run with GROMACS 4.5.5 35–37 using initial structures generated by Packmol. 38 Full details of the force field parameters for the ions can be found in the SI and are taken from our previous work. 39 For all simulations, a 1 fs time step was used and all covalent bonds to hydrogen were constrained using the LINCS algorithm. 40 The vapor phase was simulated as a single ion pair in a non-periodic box. Coulomb and Lennard-Jones (LJ) interactions were calculated using a cutoff of 20 nm. LJ parameters ij and σij were calculated using Lorentz-Berthelot combining rules. 41 The vapor systems were equilibrated for 10 ns using a Nosé-Hoover thermostat 42,43 followed by 25 ns production run

9



while sampling energies every ps. Ten simulations with different initial configurations were run in order to get better statistics and estimate uncertainties. The standard deviation of the total energy between the replicate runs was used as an uncertainty estimate for the vapor phase enthalpy. The bulk IL phase was simulated using 200 ion pairs with compositions ranging from ξ = 0 to 1 in 0.1 intervals. Uncertainties in bulk liquid properties are based on block averaging of five equal space intervals of the CMD trajectories. Both LJ and Coulomb interactions used a switch potential with a 1.2 nm cutoff and a switch applied at 1.15 nm. Long range electrostatic interactions were calculated using particle mesh Ewald summation. 44 Long-range corrections to the LJ potential were applied to energy and pressure. Following an initial steepest descent energy minimization, the systems were annealed to 700 K in the NVT ensemble for 500 ps. Then an equilibration was performed for 8 ns in the NPT ensemble with a Berendsen thermostat and barostat. 45 A final NPT equilibration using Nosé-Hoover thermostat 42,43 and a Parrinello-Rahman barostat 46 was performed for 2 ns. After equilibration, a production run of 4 ns was run sampling energies every ps using the same run parameters as the final equilibration.

Vaporization Free Energy Calculations The free energy of vaporization of a cation/anion pair was calculated using the Bennett Acceptance Ratio (BAR) method as implemented in GROMACS version 4.5.5. 47 The method uses a coupling parameter, λ, that scales the interaction between one “solute” molecule and the surrounding “solvent” molecules. The details of the BAR method and similar free energy methods can be found in Bennett’s original paper and several other works. 47–49 Due to software limitations, a cation/anion pair could not be uncoupled together from the bulk IL while maintaining their interactions with each other. Instead, a lone cation was uncoupled from the bulk IL using one series of simulations, and a separate set were run to uncouple an anion. The solvation free energy of each ion was calculated using a two step 10


Page 10 of 31

Page 11 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


approach. In the first step, the electrostatic contribution to the solvation free energy was calculated. The partial charges of the “solute” ion were scaled by a different value of λ on the interval 1 ≥ λ ≥ 0 in each of 15 simulations. In a second step, the van der Waals contribution to the solvation free energy was calculated. Charges on the “solute” ion were set to zero to guard against particle overlap that might effect the simulation stability and convergence. Similar to the electrostatic step, the LJ interactions were scaled by a different λ in each of 15 simulations. The solvation free energies computed in this way do not account for the free energy of bringing the uncoupled cation and anion together in vacuum and thus overestimate IL:CO2 (T, P ). These missing quantities are likely dominated by elec∆GIL vap (T, P ) and ∆Gvap

trostatic energies that, beyond a molecular radius, will depend only on the net charge of the ion and not the charge distribution or molecular configuration. The free energy benefit from combining the cation with an unreacted anion will largely cancel the free energy penalty from separating the cation from a reacted anion in the thermodynamic cycle. The net contribution to ∆∆Gvap is expected to be negligible.

Reaction Free Energy from Potential of Mean Force The PMF simulations were performed using Car-Parrinello molecular dynamics (CPMD). 50 For each PMF simulation, a single cation with an anion-CO2 complex was placed in a 15 Å cubic box with periodic boundary conditions. Calculations employed the PBE exchangecorrelation functional 51 with an empirical van der Waals correction applied using Grimme’s damped dispersion model. 52 The Kohn-Sham orbitals were expanded in a plane wave basis with a cutoff of 30 Ry. An electronic mass of 400 a.u. and a time step of 0.12 fs were used. To reduce the computational demands of sampling the high frequency C−H modes, hydrogen atoms were replaced with deuterium. Temperature was controlled using a NoséHoover thermostat. The constraint force was sampled along the reaction coordinate from a distance between the [2CNpyr]− ring nitrogen and the CO2 carbon from 1.1 to 2.8 Å in 0.1 Å 11



increments, from 2.8 to 4.2 Å in 0.2 Å increments, and from 4.2-5.8 Å in 0.4 Å increments. At each distance interval, a constraint was placed on the C−N distance and the system was equilibrated for 5 ps followed by a 5 ps production run where the constraint force was sampled every ten time steps. Two simulations where run at each constraint distance starting from different initial configurations in order to sample a diverse set of ion configurations. The resulting constraint forces at each constraint length were then averaged. The PMF profiles were calculated by trapezoid integration of the average constraint forces over the constraint distances using eq 5.

Reaction Free Energy from Partition Functions Partition functions were calculated for each molecule in reaction 7 with Gaussian09 53 using the PBE functional 51 for consistency with the PMF calculations and the aug-cc-pvdz basis set. 54 Computed geometries and reaction energies for the [2CNpyr]− anion with CO2 are consistent with previous results. 8

Results and Discussion Vaporization and Condensation: Steps 2 and 4 We first consider Steps 2 and 4 in the cycle, corresponding to vaporization of bulk unreacted and CO2 -reacted IL. The full set of vaporization and solvation results are provided in Tables S14 and S15. Table 1 summarizes results for ∆∆Gvap (T ) computed using the Bennett Acceptance Ratio (BAR) approach, in which a cation/anion pair “solute” was gradually deleted from the bulk IL during the course of an MD simulation. 47 Values are on the order of −10 kJ mol−1 , showing that the free energy of vaporization of the reacted IL is greater than the unreacted which serves to promote CO2 absorption. Absolute vaporization free energies decrease with temperature but differences between the reacted and unreacted ILs are nearly temperature-independent and on the order of the simulation uncertainties. 12


Page 12 of 31

Page 13 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1: Computed thermodynamic properties of steps 2, 4, and 3, in kJ mol−1 . ∆∆Svap is given in J mol−1 K−1

300

Temperature (K) 333

360

Step 4 − 2 ∆∆GBAR vap (T ) −8.60 ± 2.28 −10.64 ± 1.60 −6.85 ± 1.15 ∆∆Hvap (T ) −6.72 ± 0.82 −7.79 ± 0.48 −7.80 ± 0.34 ∆∆Svap (T ) 6.3 8.6 −2.6 Step 3 ∆G◦,PMF 0.9 ± 2.9 7.3 ± 0.6 9.7 ± 1.6 rxn -1.7 3.1 7.1 ∆G◦,anion rxn To understand the physical origins of ∆∆Gvap (T ), we also computed ∆∆Hvap from CMD simulations 55–58 using the averaged energies of the bulk phases and of ion pairs simulated in the vacuum (Table 1). The enthalpy of vaporization is larger for the CO2 -reacted IL than for the unreacted IL. This difference is due to increased van der Waals interactions with the addition of the CO2 moiety to the anion. The enthalpies of vaporization from 300 to 360 K ranged from 153 to 165 kJ mol−1 and 161 to 172 kJ mol−1 for unreacted and CO2 reacted IL, respectively. These values are comparable to both experimental measurements and calculations on similar ILs, which vary from 100 to 200 kJ mol−1 . 56,59 As shown in Table 1, enthalpy differences are of similar order of magnitude to ∆∆Gvap (T ), indicating that the majority of the difference in vaporization free energies can be attributed to an increase in cohesive energy of the larger CO2 -bound IL over the unreacted IL. The difference between the enthalpy of vaporization and BAR-computed Gibbs energy of vaporization corresponds to the difference in vaporization entropies, −T ∆∆Svap . Dividing out temperature yields the entropy differences shown in Table 1. Values are near zero and vary non-systematically with temperature. The calculated values for ∆∆Hvap and ∆∆GBAR are similar; there is only a statistivap cal difference in these values greater than one standard deviation at 333 K. The fact that ∆∆GBAR vap ≈ ∆∆Hvap suggests that ∆∆Hvap is a good first approximation to ∆∆Gvap . Calculating ∆∆Hvap only requires two liquid simulations and two vapor phase simulations whereas 13



calculation of ∆∆GBAR vap requires 120 liquid simulations. These results suggest that one may choose to use the less computationally intensive approximate calculation of ∆∆Hvap for screening larger data sets and to use ∆∆GBAR vap to investigate only the most promising candidates.

Gas-Phase Reaction: Step 3 We next consider Step 3 in the cycle, the reaction of CO2 with IL. Figure 3a reports the mean and standard deviation (computed as half the difference between the two trials) of the averaged constraint forces from the two independent trials. Figure 3b reports the corresponding mean and standard deviation of the integrated PMFs from these two independent trials. The free energy profiles have minima at the equilibrium bond length near 1.55 Å, labeled a. A representative snapshot at this constraint distance is shown in Figure 4a. The CO2 is oriented towards the anion and the O−C−O angle is approximately 133◦ . As the CO2 is drawn from the anion it climbs a steep free energy hill that maximizes at approximately 2.1 Å, reflecting a free energy barrier for CO2 desorption from the ion pair. Figure 5 reports the mean O−C−O angle along one of the computed PMFs. The angle increases nearly linearly with PMF distance up to the point that CO2 is liberated from the cation/anion pair, beyond which it approaches a constant value representative of that of free CO2 , or about 174◦ . A representative snapshot from point c at 3.5 Å is shown in Figure 4c. PMF-computed free energies of reaction 4, computed as averages of the two trials, are reported in Table 1. Results increase from about 1 to about 10 kJ mol−1 with increasing temperature, consistent with an increasing entropic penalty for CO2 binding with increasing temperature. Uncertainties are reported as half the difference between the two independent trials and are less than 3 kJ mol−1 at all temperatures. This error analysis and the physically reasonable temperature dependence both support the reliability of the calculations. Clearly the overall robustness of the PMF approach depends on the quality of the sampling. The average ∆G◦,PMF (T ) are of comparable magnitude but opposite in sign to ∆∆Gvap (T ), so rxn 14


Page 14 of 31

100

0

a

c

b

−100 1

2

3

4

5

6

2CNpyr−CO2 Bond Length (Å)

(a) 20 b

10 −1

PMF (kJ mol )

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Constraint Force (kJ mol−1 Å−1)

Page 15 of 31

c

0 −10 −20 −30

a 1

2 3 4 5 2CNpyr−CO2 Bond Length (Å)

6

(b)

Figure 3: Average and standard deviation of PMF-computed constraint forces (a) and integrated potential of mean force (b) at 300 K from two independent trials. that these two contributions counter-balance one another. As shown in Table 1, ∆G◦rxn calculated for a single anion binding with CO2 and free energies computed from partition functions are systematically 3 to 4 kJ mol−1 lower than the PMF-computed values. To leading order, the free energy of reaction should have a linear temperature dependence assuming that the enthalpy of reaction and entropy of reaction are approximately constant over the temperature range studied. An approximately linear temperature dependence is seen in calculated reaction free energies from both models. The magnitude of the Gibbs free energy of reaction in the gas phase is comparable to the sol15



(a)

(b)

(c)

Figure 4: Representative snapshops from the 300 K PMF trajectories. Configurations a, b, c correspond to the points marked in Figures 3a and 3b.

16


Page 16 of 31

Page 17 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5: Average CO2 angle from 300 K PMF trajectories. Error bars show the standard deviation of the CO2 angle at each constraint length. vation environment effects. In previous work, solvation effects were estimated to contribute

Hybrid Computational Strategy for Predicting CO2 Solubilities in

Recommend Documents