Many-Body Potentials for Aqueous Be2+ Derived from ab Initio

Article pubs.acs.org/JPCB

Many-Body Potentials for Aqueous Be2+ Derived from ab Initio Calculations Nicolas D. Winter* Physical Sciences Department, Dominican University, River Forest, Illinois 60305, United States ABSTRACT: An effective three-body potential for the aqueous Be2+ ion has been constructed from a large number of high-level ab initio cluster calculations. The new potential was validated in subsequent molecular dynamics simulations of both gas phase ion−water clusters and bulk liquid. The structures of the first and second solvation shells were studied using radial distribution functions and angular distribution functions. The vibrational spectrum of Be2+ and first shell waters was examined by computing power spectra from the molecular dynamics simulations. The observed bands showed reasonable agreement with experimental spectroscopic frequencies. The potential of mean force for water exchange between the first and second solvation shells was calculated and the energy barrier for exchange was found to have improved agreement with experiment relative to two-body force fields. Examination of the solvation structure near the transition state yielded results consistent with an associative mechanism. including second hydration shell effects via the reaction field approach21 or explicit waters.25 MD simulations by Probst6 showed that assuming pairwise intermolecular interactions results in a six-water first hydration shell whereas including three-body effects brought the hydration number down to four consistent with X-ray diffraction experiments. More recently Gnanakaran et al.30 developed a two-body Be2+−H2O potential utilizing calculated London dispersion energies which correctly predicts the 4-fold coordination. Be2+ hydration has also been simulated using the Car−Parrinello approach.7,31,32 Ab initio molecular dynamics for a Be2+ atom in a periodic box of 31 water molecules showed that the ion was tetrahydrated.7 Due to its small size and double charge, Be2+ is able to highly polarize the surrounding solvent molecules. For aqueous ions, it has been shown that many-body effects play a crucial role in the ion−water interactions33−38 and how best to include these effects in a simulation depends on the types of properties being investigated. Many-body effects due to the redistribution of electronic charge density upon hydration can be captured to some extent by pairwise potentials which include a term to account for polarization of each atom/ion. Three-body potentials can explicitly incorporate many-body effects at the cost of longer computation times, however with ever increasing computational power this barrier becomes less relevant. Classical simulation techniques, such as MD, depend critically on the quality of the force field used. For structural or thermodynamic properties that can be calculated on a short time scale, the QM/MM method8 or the hydrated ion model39 is appropriate. However, to study longer time scale phenomena, such as the exchange of water molecules around a small, highly

1. INTRODUCTION The study of ions in solution has a long history.1,2 The interaction of metal ions with water is important for many processes occurring in biology, chemistry, and geology. Experimental techniques such as X-ray and neutron diffraction as well as infrared, Raman, and NMR spectroscopy have provided insight into the nature of ions in solution.3−5 Theoretical methods used to study ions in solution range from classical methods, such as molecular dynamics (MD),6 to first-principles methods, such as Car−Parrinello,7 as well as hybrid methods, such as quantum mechanical and molecular mechanical (QM/MM) simulations.8,9 The rates of exchange of a water molecule in the first coordination shell of a hydrated metal ion with one in the second coordination shell vary by approximately 20 orders of magnitude from the most labile to the most inert. The exchange rates for the alkaline earth metal ions vary by nearly 6 orders of magnitude with [Be(H2O)4]2+ being the least labile and [Ba(H2O)8]2+ the most labile.10 The study of aqueous Be2+ in particular is of academic, industrial, and medical interest11,12 as it has the highest surface charge density of any divalent ion. The behavior of aqueous Be2+ has been the subject of numerous experimental13−17 and theoretical studies.6,15,18−25 The reactivity of Be2+ in solution is generally controlled by the lability of the coordinated solvent molecules.26 NMR studies have been used to elucidate the rate and mechanism of solvent exchange around Be2+.16,27 In particular, Merbach applied highpressure NMR techniques to study the exchange of various solvents around Be2+ and deduced an associative mechanism for water exchange around hydrated Be2+ by measuring the pressure dependence of the exchange rate and determining the activation volume.28,29 Previous electronic structure calculations on hydrated beryllium were performed on small clusters15,20 with some © 2016 American Chemical Society

Received: August 20, 2016 Revised: October 21, 2016 Published: November 7, 2016 12371

DOI: 10.1021/acs.jpcb.6b08419 J. Phys. Chem. B 2016, 120, 12371−12378

Article

The Journal of Physical Chemistry B charged ion like Be2+, a technique, such as classical MD, utilizing an analytical force field is more practical. Spangberg developed effective three-body potential functions for the interaction of water with Li+, Na+, Mg2+, and Al3+.40,41 The Li+−H2O potential function was demonstrated to provide a reasonable prediction of the zero pressure water exchange rate.42 Rustad and Stack43 used this force field to determine the potential of mean force (PMF) and rate constant for water exchange around Li+ at different pressures and from this data made an estimate of the activation volume. Recently Dang44 published a similar analysis of the ligand exchange around aqueous Be2+ and Mg2+ using classical MD and an empirically parametrized two-body polarizable force field. Currently the only available analytical three-body potential for Be2+-water is from the work of Probst.6 However, only ∼150 ab initio points were included in the fitting of the three-body portion of the potential and clusters with more than two waters were not included. In this article we report the development of an effective three-body analytical potential to describe the interactions between Be2+ and SPC/E water with the ultimate goal of using this potential to study kinetics and thermodynamics of water exchange using long classical MD simulations. The parameters in the potential are derived from a fit to high level ab initio calculations including Be2+(H2O)n clusters with up to six waters. The remainder of the article is organized as follows. In Section 2.1 we describe the ab initio calculations and the development of the analytical potential functions. In Section 2.2 we outline the MD simulations used to test the potentials. In Section 3.1 we outline results for the MD simulations of Be2+(H2O)n gas phase clusters and in Section 3.2 we discuss the results for MD simulations of Be2+ in bulk water. We summarize the article with concluding remarks in Section 4.

USR(ion − wat er) = UTIE − UCoulomb(ion − wat er) − Uwat er − wat er (3)

where all solvent−solvent interactions have been collected into one term. The short-range ion−solvent interaction energy can be expanded as follows, where i and j denote an ion or a water molecule interaction site:47 USR(ion − water) =

ij

∑ Ui i

ijk

(4)

ijkl

Truncation of this series after the second term results in a three-body potential. The resulting ion−solvent energies are then fitted to an analytical equation. If fewer terms are included in the analytical form than are present in the “‘true’” total interaction energy, the potentials derived are “effective potentials”. To be consistent with the other metal ion−water many-body potentials developed by Spangberg41 the empirical SPC/E48 model was used to describe the water−water interactions and the water geometry. The SPC/E water model is widely used and reproduces various bulk water properties reasonably well. The total interaction energy is computed for ion−water clusters of various sizes ranging from 1 to 6 waters. All ab initio calculations in this study were performed using Spartan ’14 for Linux.49 All calculations utilized the 6-311+G(d,p) basis set and the basis set superposition error (BSSE) correction. Since electron correlation plays a significant role50 for the water− water interactions we include MP2 correlation corrections determined using a frozen core approximation. All internal degrees of freedom of the water molecules were relaxed during the optimization of cluster geometries while they were kept rigid in all subsequent single point energy calculations needed for the construction of potentials. The difference in interaction energy between relaxed and rigidized clusters was found to be small, on the order of 0.1 kJ/mol for the water dimer.50 Spangberg40 performed the bulk of the calculations for the construction of the potentials at the HF/SBKJC++(d,p) level with a subsequent correction to the MP2(FC)/TZV++(d,p)/ CP level using five points on each curve computed using both methods. A single point energy calculation of a Be2+(H2O)6 cluster at the MP2(FC)/6-311+G(d,p) level in Spartan ‘14 took approximately 2 min on a 2.67 GHz Intel Xeon processor so this time-saving procedure was deemed unnecessary and all points were calculated at this level. The resulting Be−O distances for all the clusters are listed in Table 1 and are similar to previous ab initio results.23 A table summarizing Be−O distances calculated using different methods and basis sets can be found in the paper by Pye.51 Also listed in Table 1 are the total interaction energy per water and successive binding energies for each cluster size which are in good agreement with other ab initio literature values.19,52,53

2. METHODOLOGY 2.1. Derivation of Analytical Potentials. To create the effective three-body potentials we follow closely the methodology of Spangberg40 which is an extension of the “single molecule detachment-interaction energy” (SMD-IE) method developed by Periole et al.45 used to extract effective two-body potentials for aqueous ions. The procedure will be summarized below and for additional details the reader is referred to ref 40. The total interaction energy (TIE) of an ion−water cluster can be computed as UTIE = Ucluster −

∑ Uij + ∑ Uijk + ∑ Uijkl + ...

(1)

where Ui is the total energy of ion or molecule i. Ui is computed using the full basis set of the cluster to obtain the counterpoisecorrected interaction energy.46 The total interaction energy within the ion−water cluster can divided into separate contributions:

Table 1. Geometries and Energies of Be2+(H2O)n Complexes Optimized at the MP2(FC)/6-311+G(d,p) Level

UTIE = UCoulomb(ion − water) + UCoulomb(water − water) (2)

n

symmetry

r(Be−O) (Å)

TIE/n (kJ/mol)

ΔEbinding (kJ/mol)

ΔEmany‑body (kJ/mol)

where the ion−solvent and solvent−solvent Coulomb interactions are separated from the non-Coulombic short-range (SR) interactions. Higher-order electrostatic interactions are included in the short-range solvent−solvent and ion−solvent interaction terms. The ion−water short-range interaction energy can be obtained from

1 2 3 4 5 6

C2v C1 C1 C1 C1 D2

1.507 1.529 1.580 1.650 1.648 (3.375) 1.644 (3.425)

−582.3 −532.6 −459.3 −393.9 −341.1 −304.7

−582.3 −486.9 −315.9 −202.2 −121.3 −112.0

−2.2 28.7 80.6 153.1 208.0

+ USR(ion − water) + USR(water − water)

12372

DOI: 10.1021/acs.jpcb.6b08419 J. Phys. Chem. B 2016, 120, 12371−12378

Article

The Journal of Physical Chemistry B The last column in Table 1 lists ΔEmany‑body which gives the difference in binding energy per water in the cluster and the binding energy of water in a Be2+−H2O complex with the same Be−O distance minus the water−water repulsion: ΔEmany − body

N

U3body =

2

2

E(F + cos(θij))3 e−G(riO+ r jO) (8)

i=1 j=i+1

where θij is the angle between the water oxygen in molecule i, the ion and the water oxygen in molecule j. Note that the form of the two-body potential function uses different powers of r for the ion-oxygen interaction vs the ion-hydrogen interaction. This was done to maximize the quality of the fit to the ab initio data. The function reported here gave the best fit of the different combinations of powers of r that were tested. The parameters AO(H),bO(H), CO(H), DO(H), E, F, and G were fitted to the ab initio energies (with the SPC/E water−water and ion− water Coulomb energies subtracted) using a weighted nonlinear least-squares scheme as implemented by the NonlinearModelFit function in Mathematica.54 Values for all parameters are listed in Table 2.

E(Be2 +(H 2O)n ) = − E(Be2 + − H 2O) − E((H 2O)n )/n n

(5)

This value quantifies the error introduced by assuming pairwise additivity of the interactions and the importance of including many-body terms.19 For clusters with three or more waters, this error becomes increasingly large and the positive sign indicates assuming pairwise interactions (without any polarization correction) overestimates the binding energy. After geometry optimization of the cluster the water molecule geometry was modified to correspond to the SPC/ E geometry by modifying the positions of the hydrogen atoms.40 The OH distances were changed to 1.000 Å and the H−O−H angle changed to 109.47° by moving the hydrogen atoms in the plane of the molecule. Neither the position of the oxygen atom nor the direction of the vector between the oxygen atom and the point midway between the hydrogen atoms were modified. The “rigidized” clusters were used in all single-point energy calculations. For each cluster several curves of the total interaction energy as a function of ion−water oxygen distance were computed. One water molecule was moved toward and away from the ion along the Be−O vector with all the other water molecules held fixed. The following orientations were explored for each cluster size:40 (1) The ion−oxygen vector coinciding with the water dipole vector. (2) The ion−oxygen vector coinciding with the negative of the water dipole vector. (3) The ion−oxygen vector coinciding with a water O−H vector with the hydrogen pointing toward the ion. (4) The water molecule pointing a “lone pair” toward the ion, meaning the ion-oxygen vector forms a 45° angle with the water bisector vector defined as the vector between the oxygen atom and the point halfway between the two hydrogens. Additional potential energy curves were computed, where for each curve the internal geometry of the cluster apart from the nonfixed water molecule was varied by perturbing slightly the locations and orientations of the remaining fixed water molecules to more extensively probe the potential energy surface. For the Be2+(H2O)n clusters, 116, 3600, 750, 257, and 600 single-point energies were computed for n = 1, 2, 4, 5, and 6 respectively, resulting in a total of 5323 points. The analytical form of the short-range ion−water potential can be written as the sum of a two-body and a three-body term: Uion − water = U2body + U3body

N

∑ ∑

Table 2. Effective Three-Body Potential Parameters for Be2+ in SPC/E Water parameter

value

units

AO bO CO DO AH bH CH DH E F G

69063.3 3.49425 2279.29 −1790.88 36018.3 5.87503 327.601 −37.8169 0.870429 6.17802 0.187647

(kJ mol−1) (Å−1) (kJ mol−1 Å4) (kJ mol−1 Å6) (kJ mol−1) (Å−1) (kJ mol−1 Å2) (kJ mol−1 Å3) (kJ mol−1) (Å−2)

A comparison of the effective three-body potential fit to the ab initio points is shown in the top panel of Figure 1. The bottom panel of Figure 1 shows the total interaction energies for the same points calculated with the effective two-body

(6)

where the form of the two-body term is of the Buckingham type, N

U2body =

⎛

∑ ⎜⎜AOe−bOriO − i=1

⎝

CO D − 6O + ri4O riO

2

∑ AHe−bHriHk − k=1

CH D ⎞ − 3H ⎟⎟ 2 ri Hk ri Hk ⎠

Figure 1. Comparison of energies calculated using different analytical potential expressions vs ab initio energies. The top panel shows the quality of the fit using the effective three-body potential developed in this work. For comparison the bottom panel shows the potential energy values from the effective two-body potential of Periole.45

(7)

where rik is the distance between the ion and the corresponding atom in the water molecule. The three-body term has the same form as in ref 40: 12373

DOI: 10.1021/acs.jpcb.6b08419 J. Phys. Chem. B 2016, 120, 12371−12378

Article

The Journal of Physical Chemistry B potential for Be2+−H2O developed by Periole et al.45 with the water−water interactions described by the TIP3P water potential.55 The three-body potential provides a closer fit to the ab initio points than the two-body potential. The two-body potential was derived from ab initio calculations which neglected electron correlation and the BSSE correction and only included one cluster size, namely Be2+(H2O)4. Although not included in Figure 1, we also calculated the total interaction energies for all points using the three-body potential of Probst as well as the empirically parametrized two-body polarizable potential of Dang (with and without the polarizable term). Both potentials provided a poor fit to the ab initio points. This is not surprising since in the case of Probst, only clusters containing two waters amounting to 150 total points were used in fitting the three-body potential. We found that the inclusion of the polarizable term in Dang’s potential had a negligible effect on the quality of the fit to the ab initio points. This suggests that the addition of polarizability to a two-body potential may not be sufficient to describe the many body effects. 2.2. Molecular Dynamics Simulations. MD simulations were performed both for gas-phase clusters and liquid systems in order to test the new effective three-body potential. Cubic periodic boundary conditions were employed in the bulk liquid simulations with interactions smoothly switched to zero beyond spherical cutoffs at half the box length for non-Coulombic pair interactions and 6.0 Å for three-body interactions, respectively. Ewald lattice sums56 were used to calculate the Coulomb interactions. No cutoffs were used for gas-phase simulations of ion−water clusters. All MD simulations in this work were done with a 1.0 fs time step using the velocity version of the Verlet algorithm. The water molecules were kept rigid using the RATTLE57 version of the SHAKE58 method of constraint dynamics. All simulations were performed in the NVT ensemble at 298 K using the Langevin thermostat of Bussi,59,60 except NVE dynamics were used for calculation of power spectra. The system included one Be2+ ion and 550 water molecules in a cubic box of length 25.34 Å (ρH2O = 1.01 g/cm3) to be consistent with the MD simulations of Dang.44 At this density the average system pressure was < P > = −300 ± 600 bar. All simulations were duplicated for one Be2+ ion and 216 water molecules at the same density to test for size dependence. With the exception of the PMF in Figure 6, no significant differences were observed in the results for the smaller water box, so these results are not shown. The PMF was calculated as follows. The reaction coordinate is defined to be the distance between Be2+ and the oxygen atom in the water molecule. Each PMF was calculated via weighted histogram analysis method (WHAM)61 using 49 windows with centers at RC spaced 0.1 Å apart between 1.1 and 5.9 Å.43 Biasing was done with a harmonic spring, 1/2ksp(R − RC)2, with ksp = 140 kJ mol−1 Å−2. Histograms of the Be−O distance for each window were accumulated over 2 ns. W(r) was calculated on 500 bins between 1.0 and 6.0 Å. The WHAM equations were solved to a tolerance of 0.00001. Errors were estimated via Monte Carlo bootstrap analysis and were found to be