A Quantum Mechanically Derived Force Field to Predict CO

‡Curtin Institute for Computation, The Institute for Geoscience Research (TIGeR),. Department of Chemistry, Curtin University, PO Box U1987, Perth, ...
0 downloads 0 Views 2MB Size
Subscriber access provided by LAURENTIAN UNIV

Article 2

A Quantum Mechanically Derived Force Field to Predict CO Adsorption on Calcite {10.4} in An Aqueous Environment Alessandro Silvestri, Akin Budi, Evren Ataman, Mats H. M. Olsson, Martin Peter Andersson, Susan Louise Svane Stipp, Julian D. Gale, and Paolo Raiteri

J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.7b06700 • Publication Date (Web): 06 Oct 2017 Downloaded from http://pubs.acs.org on October 11, 2017

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 free 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 accessible to all readers and 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.

The Journal of Physical Chemistry C 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 21

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

A Quantum Mechanically Derived Force Field to Predict CO2 Adsorption on Calcite {10.4} in an Aqueous Environment A. Silvestri,† A. Budi,† E. Ataman,† M. H. M. Olsson,† M. P. Andersson,† S. L. S. Stipp,† J. D. Gale,‡ and P. Raiteri∗,‡ †Nano-Science Center, Department of Chemistry, University of Copenhagen, Universitetsparken 5, København Ø DK-2100, Denmark ‡Curtin Institute for Computation, The Institute for Geoscience Research (TIGeR), Department of Chemistry, Curtin University, PO Box U1987, Perth, WA 6845, Australia E-mail: [email protected]

Abstract

surface, nor contact binding of CO2 with Ca2+ . Rather, we saw a weak affinity for the surface of the ordered water layers. At 5 MPa and 323 K, we observed the nucleation of a CO2 droplet located above two structured water layers over the solid. It could not penetrate the structured water but remained bound to the second water layer for the first 10 ns of the simulation before eventually detaching and diffusing away.

Density functional theory (DFT) with semiempirical dispersion corrections (DFT-D2) has been used to calculate the binding energy of a CO2 molecule on the calcite {10.4} surface for different positions and orientations. This generated potential energy landscape was then used to parameterize a classical force field. From this, we used metadynamics (MTD) to derive free energy profiles at 300 and 350 K for CO2 binding to calcite, CO2 binding with Ca2+ and pairing of two CO2 molecules, all for aqueous conditions. We subsequently performed classical molecular dynamics (MD) simulations of CO2 and water on the {10.4} surface at pressures and temperatures relevant for CO2 geological storage. Density profiles show characteristic structured water layering at the calcite surface and two distinct phases of water and CO2 . We have also calculated the densities of the CO2 rich and water rich phases and thereby determined the mutual solubilities. For all the pressures and temperatures in the studied range, CO2 was unable to penetrate the ordered water layers and adsorb directly on the solid surface. This is further confirmed by the free energy profiles showing that, in the presence of water, there is neither direct adsorption to the {10.4}

Introduction Calcium carbonate (CaCO3 ) is a commonly occurring mineral in the natural environment. Living organisms use CaCO3 as a building block to grow biomineralized structures, such as shells. Coral polyps have shaped Earth’s largest single living structure, the Great Barrier Reef, a submerged CaCO3 formation that supports a vast diversity of life. CaCO3 forms an important part of the ocean floor sediment that is brought back to the surface by continental movements in the form of spectacular white chalk cliffs. Because these rocks are stable over millennia, calcium carbonate represents a natural material that sequesters carbon. CaCO3 is formed by the reaction of carbon dioxide (CO2 ) with dissolved calcium in seawater, releasing protons. Formation and dissolution of CaCO3

ACS Paragon Plus Environment

1

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 21

sults showed that {10.4} is the most stable surface under both wet and dry conditions and that the relative stability of calcite cleavage surfaces remains unchanged in the presence of CO2 –water mixtures. In a subsequent publication, Van Cuong et al. 47 used molecular dynamics simulations to study the effect of temperature variation and the presence of the {10.4} surface on the transport, adsorption and stability of CO2 and water. For a particular combination of force field parameters, the simulations led to unphysical results, showing a preferential adsorption of CO2 over water on {10.4}. The authors therefore stressed the importance of reliable force field parameters that correctly describe the CO2 –calcite interaction. Standard atomistic force fields are usually parameterized against the bulk properties of either solids or liquids and cannot be applied directly to simulations of interfaces. To further complicate things, the lack of experimental data describing CO2 adsorption on calcite limits the possibility of validating parameters for cross terms that result from direct application of standard combination rules. In this context, quantum mechanical approaches become the method of choice that can provide a detailed picture of the potential energy landscape at the surface. Specifically, the interaction energy of the molecule can be calculated at different surface sites and a set of binding energies is then used to fit a classical force field. In this work, we used Buckingham and Coulomb potentials to model the interaction between CO2 and calcite. Binding energy data were obtained by performing calculations for a three dimensional grid for CO2 centre of mass positions near the surface. A classical force field was parametrized using a subset of binding energies and validated against the overall quantum mechanical binding energy landscape. The force field was finally used to perform molecular dynamics and metadynamics simulations to answer a fundamental question: Can CO2 penetrate the structured layers of water at the calcite interface, displace it and bind directly to the surface?

is an important link in the carbon cycle, regulating the composition of the ocean and its response to the increasing partial pressure of CO2 in the atmosphere. Mineral carbonation, the reaction of CO2 with silicate minerals to form stable carbonates, is therefore a promising option to sequester CO2 and mitigate the induced ocean acidification and the greenhouse gas effect that leads to climate change. 1 Matter et. al 2 demonstrated how, in contrast to the common belief that the conversion of CO2 to a stable mineral requires hundreds of years, mineral carbonation in basaltic rocks can occur rapidly, in less than two years. This shows the lack of a fundamental understanding of the carbonation process and of the reactions and interactions of CO2 that result from the interplay of intermolecular and surface forces at the mineral–water interface. Molecular simulations can be a useful tool to fill this knowledge gap and have been used to study, for example, the interactions of CO2 with clay minerals. 3–16 At atmospheric pressure and temperature, the thermodynamically most stable polymorph of CaCO3 is calcite. The {10.4} family of faces represents the most stable surfaces. Because crystal dissolution and growth are processes that occur in aqueous solution, the water– calcite interface has been studied extensively for more than two decades, using both experiment 17–24 and molecular simulations. 25–44 However, only a few computational studies have investigated the interface with a CO2 –water mixture. Kerisit et al. 45 used density functional theory calculations to study the surface termination of the {10.4} calcite surface in contact with a gaseous phase containing water and CO2 . They produced the surface phase diagram showing the ranges of water and carbon dioxide chemical potentials for which a particular termination is stable and showed that a hydrated, calcium-poor surface is the most stable form close to atmospheric conditions at high relative humidity. Kvamme et al. 46 used classical molecular dynamics to predict the excess surface energies of the two most stable calcite surfaces - the {10.4} and {10.0} - under dry and wet conditions. The wetting phase comprised both pure water and CO2 –water mixtures. Their re-

ACS Paragon Plus Environment

2

Page 3 of 21

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

Methods

A [2×2] surface cell was used to reduce CO2 – CO2 self-interaction. Using this cell, the closest distance between the oxygen atoms of two CO2 periodic images was 7.77 Å and the distance between the CO2 molecule and the periodic image of the calcite slab was never less than ∼19 Å. The CO2 –calcite surface interaction was sampled on a grid: 4 evenly spaced positions in the x direction, 3 on the y direction and 8 in the z direction. At each point, the CO2 molecule was rotated about two Cartesian axes from 0 to 150 degrees in increments of 30◦ . This led to a total of 3,456 configurations. The binding energy was calculated according to the following equation;

DFT calculations We used the periodic plane wave DFT method, as implemented in Quantum ESPRESSO, 48 to calculate the binding energy of a CO2 molecule on the calcite {10.4} surface. Electronic exchange correlation effects were modelled within the generalized gradient approximation (GGA) using the revPBE functional 49 and the projector augmented wave (PAW) method 50 was used to generate the pseudopotentials 51 describing the core electrons and nuclei. We included London dispersion corrections following the DFTD2 approach, 52,53 with modified C6 parameters for ionic solids, 54 which performed well for adsorption properties of molecules on the {10.4} calcite surface in vacuo. 55–57 The optimized bulk unit cell of calcite was taken from Ataman et al. 56 where the same calculation set up was used. The convergence of unit cell dimensions was checked with respect to both the k-point grid and the cutoff energy and the lattice parameters were converged to within 0.03 Å. The optimized unit cell dimensions were a = 5.06 Å and c = 17.25 Å. These accord well with previous ab initio calculations, 58 5.03 Å and 17.17 Å, and experimental values derived from X-ray diffraction, 59 4.99 Å and 17.06 Å. In the present study, the calcite slab geometry was optimized using a plane wave energy cutoff of 37 Ry and only a single k-point (Gamma). We modelled the calcite {10.4} surface using a [2×2] supercell consisting of 4 molecular layers, to which 25 Å of vacuum was added, resulting in a cell with dimensions of 16.40 Å × 10.12 Å × 35.80 Å. The first two molecular layers of the substrate were allowed to relax and the bottom two molecular layers were fixed to their bulk positions. Although the electronic structure of the lower two layers would deviate from that of the bulk, this approach provides approximate geometric constraints to the surface that mimic those of the underlying bulk. With this setting, the adsorption energies were converged with respect to the number of layers, the cutoff energy and the k-point, to within 0.02 eV (Supporting Information).

Ebin = Eslab+mol − Eslab − Emol ,

(1)

where Eslab+mol represents the energy of the slab + molecule configuration, Eslab represents the energy of the isolated calcite slab and Emol , the energy of the isolated CO2 molecule. The energy of the isolated CO2 molecule was calculated in a box of the same size as the one used for the slab calculations.

Force field fitting The interaction between atom i of the CO2 molecule and atom j of the calcite surface was described by a potential with a Buckingham pairwise nonbonded term and a Coulomb term and, in atomic units, it takes the form: 

rij U (rij ) = Aij exp − ρij





Cij qi qj , + rij6 rij

(2)

where A, ρ and C represent the constants of the Buckingham potential. Here, C corresponds to the C6 term referred to earlier. To derive our force field, we started from a reasonable set of parameters and calculated the binding energy for the full set of configurations. According to the Boltzmann distribution, the fraction of configurations with an energy above 100 kJ/mol is on the order of 10−18 at ambient temperature. These high energy configurations would therefore not be explored in MD simulations and were excluded from further consideration.

ACS Paragon Plus Environment

3

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

A total of 20 configurations consisting of those for which the prediction was poor and the one that gave the lowest predicted binding state were used as an initial training set. We used the same weights for all the configurations except for the minimum energy configuration, for which we chose a 10,000 times larger weighting factor. The large weighting factor on the minimum energy structure ensures that this is reproduced, however it is insufficient to constrain all the fitting parameters to a unique solution and the remaining 19 configurations were used to drive the fitting of the curvature of the binding energy curve. Atomic partial charges were taken from Harris and Yung’s CO2 EPM2 model 60 and from Raiteri et al. 61 for calcite and were not modified in the fitting procedure. The parameters Aij and ρij were adjusted using least squares fitting, as implemented in the GULP package. 62 The attractive C term is included only for the O–O interaction. One possible strategy would have been to take the C6 parameters straight from DFT-D2. However, because these are damped at short range, while in the typical force field forms this is not the case, we decided to explore varying the C coefficient. This term was initially kept constant throughout the Aij and ρij fitting process. However, after an initial fit of the repulsive Buckingham parameters, the starting value of C was adjusted manually until good agreement with the DFT energies was achieved.

Page 4 of 21

molecules as rigid entities. It was previously shown that the incorporation of flexible bonds does not change the computed phase boundaries appreciably. 66 Therefore, because we used flexible models for calcite and water, we incorporated the bond stretching and angle bending terms from Cygan et al. 4 The use of standard combination rules to get the cross term interaction parameters has been questioned by Vlcek et al. 67 They optimized the cross term interaction parameters of water and CO2 , described respectively by the rigid SPC/E and EPM2 models, to accurately reproduce experimental mutual solubilities in fluid– fluid phase equilibria at conditions relevant for CO2 geological sequestration. Myshakin et al. 6 later studied the intercalation of CO2 in hydrated Na–montmorillonite using the flexible SPC water model 68 and the Cygan version of the EPM2 model. 4 They used standard combination rules for the cross term interaction parameters and showed how the radial distribution functions resulting from their simulations closely resembled those of Vlcek et al. 67 They concluded that the deficiencies in the force field were not related to the use of combination rules for determining the cross term parameters but might be associated with the inability of a 3point charge model to correctly represent the electrostatic potential of CO2 . We therefore derived the water–CO2 cross terms by applying Lorentz–Berthelot mixing rules; ǫij =

Molecular dynamics simulations and

Classical molecular dynamics simulations were performed using the LAMMPS package. 63 Calcite was modelled with the thermodynamically consistent force field of Raiteri et al. 61 The SPC/Fw model 64 was used for water. This force field has been shown to be transferable to high temperatures 61 and was used to study diffusion processes at the calcite {10.4} interface under water supercritical conditions. 65 For CO2 , we used the elementary physical model (EPM2), 60 which is able to correctly describe the phase boundaries in the liquid–vapor equilibrium and the critical properties of pure CO2 . 66 Simulations using the rigid EPM2 model treat CO2



ǫi ǫj

(3)

σi + σj , (4) 2 where ǫ and σ represent the Lennard-Jones parameters. To check the validity of using combination rules to obtain the cross terms, we calculated the solvation free energy of CO2 in water using the free energy perturbation technique (FEP), 69,70 as implemented in LAMMPS. The interactions between the solute and solvent molecules were gradually switched off in 40 stages of 5 ns each. The first 20 stages for the electrostatic interactions were followed by 20 stages for the van der Waals terms. We σij =

ACS Paragon Plus Environment

4

Page 5 of 21

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

included the finite size 71 and standard state 72 corrections. The volume correction required for the NPT simulations was negligible for our system. The computed value at 300 K was 1.15 kJ/mol, which compares well with the experimental value of 1.00 kJ/mol. 73 Finally, the interactions between CO2 and calcite were described by our new set of force field parameters. A cutoff distance of 9.0 Å was used for the dispersion interactions in CO2 , water and between CO2 and water. As for the dispersion interactions of calcite with both CO2 and water, we used the tapering function of Mei et al.; 74 x= and

r − rt rc − rt

ing the NPT ensemble at 5, 10 and 25 MPa and for each pressure, at 298, 323 and 373 K. From these first simulations, we calculated the average cell volume and used it to perform NVT simulations for another 20 ns. Atomic configurations were stored every 1,000 steps (1 ps) for trajectory analysis.

Metadynamics Pairing free energies in water were calculated with the metadynamics technique, 78 using PLUMED 2.0, 79 which provides a plugin for free energy calculations that can be interfaced with LAMMPS. We calculated the CO2 – CO2 and CO2 –Ca2+ pairing free energies and the CO2 –calcite adsorption free energy in water using as collective variables, the distance from the carbon atom of CO2 and i) the carbon atom of the second CO2 molecule, ii) the calcium ion and iii) a calcium ion from the bottom layer of the calcite slab. The free energy profiles were constructed by running metadynamics simulations with Gaussians laid every 1 ps, with an initial height equal to 2.5 kJ/mol. The Gaussian exponents were 0.1 Å along the interatomic distance. To ensure convergence of the final free energy, we made 30 parallel simulations using the multiple walkers algorithm 80 for a total simulation time of 150 ns and we employed the well tempered technique 81 with a bias factor of 5 to progressively reduce the heights of the Gaussians and ensure a smooth convergence of the free energy profile. For the binding of CO2 on calcite, we used the same simulation box described above. For the CO2 – CO2 and CO2 –Ca2+ pairing, we used a cubic box with side length of 49.72 Å containing 4,136 water molecules. All metadynamics simulations were performed at 300 and 350 K.

(5)

  1 f (x) = (1 + 3x + 6x2 )(1 − x)3   0

if r < rt if rt ≤ r ≤ rc , if r > rc (6) where we used 6 Å for rt and 9 Å for rc . A cutoff of 10 Å was used for the real space contribution to the electrostatic interaction. The reciprocal space electrostatics were calculated using the PPPM algorithm 75 with an accuracy of 10−5 . The simulation box for the calcite-water interface had starting dimensions of 49.32 Å × 48.55 Å × 150.00 Å. The calcite slab was approximately 30 Å thick and contained 1,440 CaCO3 units distributed in 12 layers. We used 1,000 CO2 molecules and 4,000 water molecules for all simulations. For a 4:1 proportion of water to CO2 , two distinct phases were formed during each simulation. The temperature was controlled via a Nosé–Hoover chain of thermostats 76 of length 5 with a relaxation time of 0.1 ps. The cell volume was controlled with the Parrinello-Rahman algorithm 77 with a relaxation time of 1 ps, the cell angles were kept at 90 ◦ while the cell lengths were allowed to fluctuate independently. Periodic boundary conditions were applied in all directions. The equations of motion were integrated with a time step of 1 fs. All systems were energy minimized before starting the molecular dynamics simulations. We initially performed our simulations for 20 ns us-

Results and discussion DFT binding energy landscape First, we analyse the features of the binding energy landscape that are predicted by DFT. Figure 1 shows three x-y scans. The distance between the molecule and the surface is defined

ACS Paragon Plus Environment

5

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 6 of 21

low energy configurations that would be explored during MD simulations to within a deviation of ± 2.5 kJ/mol, which is comparable to the magnitude of thermal fluctuations per degree of freedom.

by the distance between the CO2 carbon atom and the z height relative to the first layer of Ca atoms. For the results shown in Figures 1a and 1b, the carbon atom of the CO2 molecule is 3.98 Å from the calcium atoms plane. When the molecule is perpendicular with respect to the surface (Figure 1a), the binding energy landscape has both attractive and repulsive sites. The attractive sites result from the interaction between the oxygen atom of CO2 and the calcium ion. The two repulsive sites result from the interaction between the oxygen atom of CO2 and the out-of-the-plane oxygen atoms of the carbonate units. When the CO2 molecule is parallel to the surface and aligned with the long axis of the calcite surface unit cell (Figure 1b), the two repulsive sites disappear at this height and the binding energy is weakly attractive. Figure 1c shows the adsorption site with the highest binding energy predicted in this study, -22.19 kJ/mol. This represents the only attractive site for an otherwise entirely repulsive scan for this particular configuration, where CO2 is close to the surface, at a height of ∼3 Å and oriented parallel to the short axis of the surface unit cell. This is about the distance at which we would expect the next layer of calcite molecules to be centred.

Table 1: Force field parameters for the CO2 -calcite interaction. CO2

calcite

A [eV]

ρ [Å]

O C O O

Ca O C O

55882.390340 431.436033 4140.137395 1578.843140

0.194676 0.300000 0.275682 0.279776

C [eV Å6 ] 0.0 0.0 0.0 20.284

Next, we considered three physically relevant configurations. Two of them are the adsorption configurations that resulted from the MD simulations of Van Cuong et al., 47 using two different force field combinations. These configurations are not part of the training set on which the force field parameters were fitted and can be used to validate the results. We therefore performed z dependent scans for: 1. the CO2 molecule perpendicular to the surface on top of the calcium atom; 2. the CO2 molecule parallel to the surface with the carbon atom on top of the above plane oxygen of a carbonate group;

Force field validation The resulting set of fitted parameters is shown in Table 1. To assess the accuracy of the force field, we estimated the root mean square deviation (rmsd) of the whole set of configurations with energy lower than 100 kJ/mol: v u N u1 X t δE 2 , (7) rmsd = N k=1 k

3. the converged adsorption geometry from Ataman et al. 56 Figure 3 shows the three configurations from top (left) and side (middle) views, along with the z scans predicted by DFT and the force field. Figure 3a shows that for the CO2 molecule perpendicular to the surface on top of the calcium atom, the force field predicts a binding energy that is too attractive. For the two other configurations (3b and c), the force field is in very good agreement with the DFT energies.

where N represents the total number of configurations and δEk stands for the difference between the DFT and force field energy for a particular configuration, k. The rmsd for these configurations is 7.91 kJ/mol. The corresponding scatter plot is shown in Figure 2. The force field overestimates the energy of repulsive configurations but predicts most of the attractive,

Free energy profiles Figure 4 shows the metadynamics reconstructed pairing free energy profiles for two

ACS Paragon Plus Environment

6

Page 7 of 21 z = 4 Å e = 0! q = 0!

z = 4 Å e = 0! q = 90!

!"#"$"%"""e"#"&'(!"""q"#"(!

"') "(

"*+,-"./012345

")(

ï') C O Ca

a

b

c

Figure 1: CO2 binding energy as a function of its positions on the calcite {10.4} surface unit cell (top row) for different orientations of the molecule (a, b and c). The sampled sites are separated by 2.05 Å in the x direction and 1.69 Å in the y direction. The center of each tile in the top row indicates a binding site and is colored according to the energy of the corresponding configuration; repulsive sites are shown in red and attractive sites in blue, with intensity of colour corresponding with level of energy. The dotted lines (top row) shows the periodic boundaries of the surface unit cell. The bottom row shows the orientation of the CO2 molecule used for the calculation of the binding energy in the corresponding panel on the first row; carbon atoms are represented in blue, oxygen atoms in red and calcium ions in orange.

CO2 molecules and Figure 5, the same relationship for a CO2 molecule and a calcium ion. Free energy is plotted as a function of the interatomic distance at two different temperatures with the corresponding analytical solution for the configurational entropy of the system, which constitutes the main contribution to the free energy, at large distances. To calculate the standard free energy of binding, we followed the procedure outlined by Raiteri et al., 61 which expands on the approach of Bjerrum and Fuoss, 82 where the dissociation constant of an ion pair is related to the integral of the potential of mean force φ(r):

100

75

50 EFF (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

The Journal of Physical Chemistry

25

0

-25

-50 -50

-25

0

25 50 EDFT (kJ/mol)

75

100

Kdiss = 4πC0

Figure 2: Scatter plot showing force field versus DFT energies for 2,990 configurations with DFT energy