Toward a Physically Motivated Force Field: Hydrogen Bond

Subscriber access provided by UNSW Library

Article

Toward a Physically-Motivated Force Field: Hydrogen Bond Directionality from a Symmetry-Adapted Perturbation Theory Perspective Maxim Tafipolsky, and Kay Ansorg J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.5b01057 • Publication Date (Web): 28 Jan 2016 Downloaded from http://pubs.acs.org on February 2, 2016

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.

Journal of Chemical Theory and Computation 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 28

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

Journal of Chemical Theory and Computation

Toward a Physically-Motivated Force Field: Hydrogen Bond Directionality from a Symmetry-Adapted Perturbation Theory Perspective Maxim Tafipolsky*, Kay Ansorg January 27, 2016 Institut für Physikalische und Theoretische Chemie, Universität Würzburg, Campus Hubland Nord, Emil-Fischer-Strasse 42, D-97074 Würzburg, Germany

Corresponding Author * Dr. Maxim Tafipolsky, Institut für Physikalische und Theoretische Chemie, Universität Würzburg, Campus Hubland Nord, Emil-Fischer-Strasse 42, D-97074 Würzburg, Germany Telefone: +49 9313186983, E-mail: [email protected]

Abstract It is argued here that the functional forms adopted in almost all popular force fields are too restrictive to allow for accurate and physics-based parameterization. Some important modifications are suggested based on symmetry-adapted intermolecular perturbation theory which directly separates the intermolecular interaction energy into four physically interpretable components: electrostatics, exchange-repulsion, dispersion, and induction. The exact electrostatic energy is approximated as a sum of the short range contribution (due to charge density penetration effects) included explicitly and the long range part (via distributed atomic multipoles), whereas the induction energy is evaluated by means of the distributed induced damped point dipole model. The dispersion energy is fitted to a simple analytical function and the exchange-repulsion contribution is approximated by the overlap of the valence-only electron charge densities of monomers. The water dimer is used to illustrate the approach and to discuss its potential and possible improvements. Analysis of the four main contributions to the binding energy allows for a deeper understanding of the hydrogen bond directionality. It is found that a notorious geometrical preference in the water dimer results mainly from large polarization contributions including induction and dispersion.

Keywords: Force field, AMOEBA, parameterization, genetic algorithm, charge penetration, electron density overlap, FF-SAPT

ACS Paragon Plus Environment 1


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

1

Page 2 of 28

Introduction Molecular simulations based on force fields are numerous and will remain very important in the near future. The

accuracy of the results, however, depends critically on the quality of the force field employed. This sounds trivial but nevertheless real progress in developing more accurate force fields is still slow. This conservatism is based mainly on the presumption (partly valid at the time when the force field development was started, in the 1980s) that the experimental data (mainly thermodynamic), used exclusively as reference, are scarce and therefore the number of parameters should be kept as small as possible. Additionally, within the paradigm of separated treatment of the intra- and intermolecular force fields, the functional forms are kept as simple as possible to speed up the calculations. Whereas the parameterization of the intramolecular force fields is well-understood and can nowadays be done routinely and sufficiently accurately,1 this is unfortunately not the case for the intermolecular force fields. For the parameterization of the intermolecular force fields, accurate ab initio results based on a supermolecular approach (giving the total interaction energy between monomers in a dimer) have been used to supplement (or replace) the available experimental data.2–38 But still, this latter strategy relies heavily on error compensation between various components of the force field. The situation has been changed dramatically with the advent of Symmetry-Adapted intermolecular Perturbation Theory (SAPT)39,40 which gives the interaction energy directly as a sum of four physically interpretable components: electrostatics, exchangerepulsion, dispersion, and induction. This powerful technique is still underappreciated in the community41,42 despite the fact that the parameterization of the intermolecular force fields from first-principles calculations exclusively based on intermolecular perturbation theory as reference was first successfully attempted in the 80s,43–46 and now started to emerge slowly.16,22,47–50 Recent developments and efficient implementations of SAPT51,52 enable highly accurate studies of intermolecular interactions at a level comparable to a state-of-the-art methods such as coupled-cluster theory, CCSD(T), but with a much reduced computational efforts. Moreover, it allows for a deeper understanding of different contributions to the intermolecular interactions thus giving a more detailed picture in comparison with the supermolecular approach where only the total interaction energy is available. Of course, only the total energy is a quantum-mechanical observable and all partitions into different components are ad hoc. Nevertheless, this partitioning is very useful and particularly suited if the development of an accurate force field is of concern. Here, we advocate some modifications to the popular intermolecular potentials needed on the way to a more robust and physically-motivated force field and illustrate the approach taking water as an example. To describe the short-range electrostatic contribution, we explicitly included the charge penetration term as is done before.21 In contrast to our previous work,21,53 we now fit the exchange-repulsion and dispersion components separately. As a starting point for our model potential developments, we use the AMOEBA54 polarizable force field that is capable of including non-pairwise-additive many body effects via induced atomic dipoles. We


Page 3 of 28

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 1: Three geometrical arrangements of the water dimer: linear (left), planar cyclic (middle) and bifurcated (right). implement a simple model for the van der Waals potential (with damped dispersion) proposed by Tang and Toennies55 some time ago. For the parameterization we use as reference published data on water dimers.56,57 We then compare our newly developed polarizable potential with the original AMOEBA model and also with highly accurate calculations. We seek for simple functional forms albeit physically motivated which can be readily implemented in standard force fields but, at the same time, keeping the number of fit parameters as small as possible. The flexibility of the water molecule is ignored in this work. We refer to the recent work of Jankowski et al.58 and of Babin et al.59 for a discussion on the importance of including the internal degrees of freedom in the intermolecular potential of water. The main objective of this work is not just to develop highly accurate force fields for water (see, for example, refs.48,49,59–64 ), but, in addition, to draw attention to some (partly known) deficiencies in popular simple potentials which definitely prevent further progress on the way to a more reliable and transferable force fields for atomistic simulations.65,66

2

Methods All 2510 water dimer configurations and their corresponding interaction energies were used as our training data set.57

These are based on SAPT(DFT) calculations, where the intramolecular correlation effects are evaluated at the DensityFunctional Theory (DFT) level. The vibrationally averaged monomer geometry of water56 (r(O−H) = 0.9716 Å, ∠H−O−H = 104.69º) was used and kept rigid in all dimer calculations. All data points were used with equal weight (no weighting scheme was applied). The quality of our new force field was evaluated by comparing its predictions for a test data set. We performed additional DFT-SAPT calculations for the test set which consists of four cuts through the water-water dimer potential energy surface including the most relevant dimer configurations.67–71 Three configurations are depicted in Figure 1 and will be called linear, planar cyclic and bifurcated water dimers, respectively.

Figure 1 near here The cuts were obtained by rigidly displacing the second water molecule along the (HO-)H...O vector for the linear and planar cyclic configurations. For the bifurcated structure, the distance between the two oxygen atoms was varied. In all cases the mutual orientation of the two water molecules was kept unchanged. A forth cut was added to the test set to investigate the ability of our force field to describe the angular dependence of the intermolecular energy in the linear dimer. The angle between the O...O line and the H-O-H bisector of the acceptor water molecule was varied keeping the



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 4 of 28

Figure 2: Definition of the angular variable for the linear water dimer. O...O distance between the two rigid monomers fixed at 2.91 Å (shown in Figure 2).

Figure 2 near here For our SAPT calculations, the dimer-centered basis sets, i.e. all atomic basis functions of both monomers were employed in both monomer and dimer calculations. To speed up the SAPT calculations, we employed the densityfitting variant72 of DFT-SAPT as implemented in the MOLPRO73,74 program package using the augmented correlation consistent aug-cc-pVQZ basis sets of Dunning75 (abbreviated as avqz, respectively) along with an appropriate combination of cc-pV(Q+1)Z JK and aug-cc-pVQZ MP2-fitting basis sets of Weigend et al.76,77 To produce accurate results, the asymptotic correction to the exchange-correlation potential was applied. We used the gradient-regulated asymptotic correction approach of Grüning et al.,78 to correct the PBE079,80 hybrid exchange and PW91 correlation functionals by employing the calculated ionization potentials of 0.4660 Hartree (exp: 0.463881 for water along with the HOMO energy of −0.3339 Hartree obtained using PBE019 functional. The corresponding shift parameters were approximated by the difference between the HOMO energy and the (negative) ionization potential of each monomer. Accurate evaluation of the dispersion (and its exchange counterpart) contribution is the most time-consuming and difficult part of the SAPT calculation.82,83 This particular term usually requires the use of extended basis sets.84 To check the accuracy of our SAPT results, the total intermolecular energies were obtained with the CCSD(T) method (only valence electrons were correlated) corrected for the basis set superposition error (BSSE) through the counterpoise correction,85 and utilizing the aug-cc-pVQZ basis set75 using the Gaussian program package.86 The deviations of the SAPT total interaction energies from their CCSD(T) counterparts were not exceeded 0.2 kcal/mol, see the Supporting Information. In a recent study by Korona87 it has been shown (see Table 4 there for the water dimer) that the PBE0 density functional is able to reproduce very closely the results from SAPT with intramonomer electron correlation described by coupled cluster theory limited to single and double excitations. The total interaction energy calculated within the SAPT framework, truncated at second order, is usually written as the sum of individual components containing the first (electrostatic) and the second (induction and dispersion) order interaction terms accompanied by their respective exchange counterparts (the superscripts 1 or 2 in parentheses refer to the order of the perturbation correction):


Page 5 of 28

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


(1)

(1)

(2)

(2)

(2)

(2)

SAPT Etot = Eelst + Eexch + Eind + Eexch−ind + Edisp + Eexch−disp + δ(HF)

(1)

As no higher than second-order terms are currently implemented in the MOLPRO program we used, third- and higherorder induction and exchange-induction contributions were estimated using the supermolecular approach at the HartreeFock level88 (counterpoise corrected for basis set superposition errors) and added to the sum of first- and second-order DFT-SAPT energy contributions:

(1)

(1)

(2)

(2)

SAPT δ(HF) = Etot (HF) − Eelst (HF) − Eexch (HF) − Eind (HF) − Eexch−ind (HF)

(2)

To get the dispersion energy, we combine two terms (dispersion and exchange-dispersion) together:

(2)

(2)

SAPT Edisp = Edisp + Eexch−disp

(3)

whereas the sum of induction, exchange-induction and higher-order (δ(HF)) contributions is used to represent the induction energy:

(2)

(2)

SAPT Eind = Eind + Eexch−ind + δ(HF)

(4)

Based on the analysis of the third-order induction energy in SAPT, Patkowski et al.89 have shown that the hybrid approach including δ(HF) for the polar molecules (like water) gives reliable results. We also note that, within the methodology of SAPT, a further separation of mutual perturbation of interacting monomers into internal polarization and external charge transfer is not done (see also recent developments by Misquitta90 ). In all standard force fields the so-called van der Waals term consists of two parts: repulsive and attractive. The “attractive” part is aimed at describing the dispersion contribution valid only at long range, whereas the “repulsive” part includes the exchange and Pauli repulsion terms and, in addition, absorbs other short-range contributions. The most important one is due to charge penetration that has electrostatic origin and is attractive. The remaining long-range part of the electrostatic energy is treated usually with point charges, bond dipoles or atomic multipoles. Since the intermolecular force fields are usually parameterized to reproduce the total interaction energy, some error compensation is unavoidable.



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 28

We follow a different approach because our objective is to develop an accurate force field with error compensation effects as small as possible using SAPT as the reference data. Hence, in our force field the long range part of the electrostatic energy is represented with atom-centered multipoles (up to and including quadrupoles) calculated with the Distributed Multipole Analysis (DMA) program91,92 (for more details, see54 ). We used the atom-centered multipoles calculated with the GDMA program directly (i.e. without any further refinement). Those were obtained from the wavefunctions calculated at the MP2(full)/aug-cc-pVTZ level of theory using the Gaussian program package.86 We used grid-based quadrature for partitioning the contributions to the charge density from diffuse basis functions (via the keyword “SWITCH 4.0”).92 Furthermore, as in our previous work, the electrostatic contribution due to charge penetration is included explicitly93 in order to estimate the short-range quantum effects that are not accounted for by the classical multipolar expansion valid only at long range. We should clearly stress here that our approach represents in some sense a paradigm shift: instead of refining some effective point charges to reproduce the long-range electrostatic interactions (this is done usually via fitting to mimic the electrostatic potential outside the atomic van der Waals region), we seek to reproduce the true electrostatic contribution both at short and long range. We evaluate this short-range contribution to the exact electrostatic energy by a sum of classical Coulomb interaction between spherical atomic charge densities by allowing their radial dependence to expand or contract (for more details, see21 ). To this end, the expansion-contraction parameters (one parameter for each atom type in our model), which reflect the tendency of atoms to expand or contract in a molecular environment, are used to scale the exponential parameter of the Slater-type orbitals which describe the atomic charge density (see also refs.94–96 ). These expansion-contraction parameters are the fitting parameters of our model. They reflect to some extent a deviation (albeit isotropic) of the valence charge density distributions around the nuclei in a molecule as compared to those in the unperturbed atoms. This short-range pairwise additive term is added to the long-range (multipolar) DMA (1)

contribution, and the expansion-contraction parameters were tuned to reproduce the SAPT electrostatic energy (Eelst , see Eq. 1). Their refined values were used to evaluate the integrals (see Eq. 4 in ref.

21

). To speed up the calculation

of the explicit charge penetration term, precalculated look-up tables can be utilized.97 For the ease of reproduction of our results, we have fitted the numerical values (integrals) using a five-parameter exponential function which is a sligthly modified function proposed by Ng et al.98 (for atom types a and b separated by a distance R):

Eab (R) = − exp(a0 R + a1 + a2 R−1 + a3 R−2 + a4 R−3 ),

(5)

where the coefficients a0 , ..., a4 for the O...O, O...H and H...H atomic pairs are given in Table 1 with the energy given in kcal/mol.

Table 1 near here


Page 7 of 28

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: Coefficients for the analytical fit of the charge penetration term (Eq. 5). Asymptotic standard errors are given in parentheses.99 a0 , 1/Å a1 a2 , Å a3 , Å2 a4 , Å3 O...O −3.88(1) 12.56(6) −4.1(1) 3.39(9) −1.59(3) O...H −4.32(2) 6.66(6) 0.27(6) 0.06(3) −0.108(5) H...H −13.85(5) 11.80(7) −2.13(4) 0.48(1) −0.059(1) The dispersion interactions are notoriously difficult to describe accurately with ab initio methods since higher level of electron correlation treatment and extended basis sets are needed. Therefore, more approximate but reliable and practicable models have been developed over the past years. In all standard force fields, however, this contribution is combined with the exchange-repulsion (including penetration!) part of the so-called van der Waals term. We advocate SAPT here the idea of parameterizing this particular contribution separately based on accurate SAPT reference data (Edisp ,

Eq. 3) on model systems small enough to be calculated reliably. This idea is not new, but to the best of our knowledge, has not been implemented routinely yet. Nevertheless, a number of successful attempts exist in the literature.22,27,45,82,100–103 In our approach, possible double counting effects (concerning dispersion contribution at short range) are not present by SAPT definition in contrast to some popular DFT plus dispersion approaches.104 To approximate the dispersion energy, Edisp ,

the London formula damped at short range is used:

FF−SAPT Edisp (R) = fdamp (βab R)

Cab . R6

(6)

To minimize the effects due to asymptotic divergence of the London term (∼ R−6 ) at short range, we decided to implement the damping function proposed long time ago by Tang and Toennies55 and used before (see, for example, refs.61,102,105–107 ):

fdamp (βab R) = 1 − e−βab R

6 X (βab R)n . n! n=0

(7)

The two parameters per atom type (βaa and Caa ) were fitted against the sum of the second-order dispersion and exchangeSAPT dispersion reference values, Edisp (Eq. 3), as has been done before.101 For the unlike atom pairs the following “combining

rules” are employed:

Cab =

p Caa Cbb ,

βab =

2βaa βbb . βaa + βbb

(8)

We have implemented both energy and analytic energy gradients for the damped dispersion term (Eq. 6) in a locally modified TINKER program package.108



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 8 of 28

Since the widely used potentials to describe the exchange-repulsion part (such as Lennard-Jones or Buckingham) lack anisotropy, which turned out to be very important, we adopted the overlap model used previously by Kim et al.109 and later refined by Wheatley and Price110 (see also refs.48,111 and references therein):

(9)

FF−SAPT Eexch−rep = K Sa,

where ˆ S=

ρA (r)ρB (r)d3 r,

with ρA(B) the unperturbed ground-state electron charge density of monomers A (B) with the proportionality constant, K, (1)

being an adjustable parameter. Below we will show, by comparison with the short-range first-order exchange, Eexch , that the use of the valence-only electron density (replacing the core electrons with effective core potentials taken from ref.112 and using the D95 basis set113 for the H atoms along with the PBE0 functional) can be a viable alternative to speed up the calculation of integrals in Eq. 9. To evaluate the induction contribution, we used the induced point dipole method as implemented in the AMOEBA force field. In this approach, molecular polarization is achieved via an interactive induction model with distributed atomic polarizabilities based on Thole’s damped interaction method in order to avoid a polarization catastrophe at very short range.54 Iterating the dipole−induced dipole interaction to self-consistency captures the many-body effects. The molecular polarizability is expressed as a sum over atomic isotropic dipole polarizabilities. All parameter refinements utilized the same genetic algorithm PIKAIA114 as done before21 with the root-mean-square deviation (rmsd) between the SAPT reference data and the corresponding part of the potential used as a fitness function.

3

Results and Discussion Below, we analyze in turn the four distinct and physically-motivated components originated from SAPT for our newly

developed force field for water. The electrostatic energy. The expansion-contraction parameters for the oxygen (0.86) and hydrogen (3.7) atoms of the water molecule were optimized in this work to reproduce the first-order electrostatic term. The refined values for the oxygen and hydrogen atoms in water reflect nicely the difference in electronegativity of the two atoms: whereas the refined kappa value for the oxygen atom is smaller than 1.0 (0.86) pointing to an expanded valence shell, the optimized kappa value for the hydrogen atom (3.7) deviates significantly from 1.0 indicating large contraction of its valence electron density. The long-range part of the electrostatic energy is described by Coulomb interactions between atomic multipoles up to


Page 9 of 28

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


and including quadrupole-quadrupole interactions. The calculated electrostatic energy as a function of the intermolecular separation is shown in Figure 3. It is clearly seen that the Coulomb energy due to the distributed multipoles appreciably underestimates the electrostatic component at shorter separations which is a manifestation of dominating role played by the charge penetration effects. The difference between the refined multipoles used in the AMOEBA force field (taken from water14.prm115 and labeled by “multipoles (AMOEBA)”) and those calculated from GDMA (labeled by “multipoles (GDMA)”) is also evident. If we include the short-range contribution, estimated by the method described in our previous work,21 the agreement with the reference data is significantly improved. The same conclusion has been reached in a recent study on water dimers by Kumar et al.,102 where charge penetration effects were approximated by a modified effective point charge model (see also refs.?,116 ). Similar improvements have been achieved in a recent screened point-charge approach developed by Wang and Truhlar.117

Figure 3 near here

It is worth noting that the discrepancy between the reference electrostatic energy and that calculated from distributed multipoles (up to quadrupole-quadrupole term) seen at the small intermolecular separations is entirely due to the missing charge penetration effects as shown in our previous work,21 and not to missing higher multipoles as sometimes claimed in the literature118,119 (see also the discussion in refs.33,91,120 ). Moreover, at short range, a description of the electrostatic interaction via multipole expansions becomes less and less valid anyway. For a careful discussion on the importance of the charge penetration effects we refer here to the book by Stone.121 To study hydrogen bond directionality, a more closer look at the angular dependence of the linear dimer is instructive (see Figure 3, bottom right). As can be seen from this Figure, the multipoles used in the AMOEBA force field (taken from water14.prm115 ) are not able to reproduce the reference electrostatic energy (from SAPT) even qualitatively. This is probably due to much smaller values of the dipole and quadrupole centered on the oxygen atom. At the same time the multipoles calculated by the GDMA procedure and used unchanged reproduce the angular dependence of the reference data fairly good. By adding the charge penetration term to the GDMA long-range multipole part we recover the reference data almost quantitatively. We note that point charges of −0.72 e (O) and 0.36 e (H) obtained by the fit (at MP2(full)/augcc-pVTZ level of theory; Merz-Kollman sampling scheme with ca. 6000 points per atom) to the electrostatic potential (ESP) also exhibit a wrong angular dependence due to unfavorable charge-charge interactions especially when the two hydrogens of the acceptor molecule are coming closer to the hydrogen atom of the donor molecule. The induction energy. Isotropic polarizabilities for oxygen (0.92 Å3 ) and hydrogen (0.539 Å3 ) in water were adopted from the AMOEBA force field (taken from water14.prm115 ). It turned out that to reproduce the induction energy component for water dimers (see Eq. 4) it was deemed necessary to use a value for the damping factor twice as large (0.78) as the recommended one (0.39), which is in contrast to a recent study of Cisneros.122 As can be seen from



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

0

−2 −4 −6 −8 −10 −12 −14 −16 −18

−5 −10 −15

SAPT FF−SAPT multipoles (AMOEBA) multipoles (GDMA)

−20 −25 1.6

0 −2 −4 −6 −8 −10 −12 −14 −16

1.8

2 2.2 2.4 RO...H (Å)

2.6

2.8

Page 10 of 28

SAPT FF−SAPT multipoles (AMOEBA) multipoles (GDMA) 1.8

2

2.2 2.4 RO...H (Å)

2.6

2.8

40

60

SAPT FF−SAPT multipoles (GDMA) multipoles (AMOEBA) charges (ESP fit)

−3 −4 −5 −6

SAPT FF−SAPT multipoles (AMOEBA) multipoles (GDMA) 2.4

2.6

2.8 3 3.2 RO...O (Å)

−7 −8 3.4

3.6

−9

−40

−20

0 20 α (deg)

Figure 3: Long-range multipole electrostatics with and without the optimized short-range contributions against the DFTSAPT reference first-order electrostatic energy (in kcal/mol) for the three water dimers (top left: linear, top right: planar cyclic, bottom left: bifurcated). Distances are in Angstrom and angle in degrees. The vertical broken line marks the approximate equilibrium distances. Bottom right: angular dependence for the linear water dimer with the O...O distance fixed at 2.91 Å. The vertical broken line marks the approximate equilibrium angle. See Figure 2 for the definition of the angle α.


Page 11 of 28

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


0

0

−2

−1

−4

−2

−6 −3

−8 SAPT FF−SAPT (damping=0.78) AMOEBA (damping=0.39)

−10 −12

1.6

1.8

2

2.2 2.4 RO...H (Å)

SAPT FF−SAPT (damping=0.78) AMOEBA (damping=0.39)

−4

2.6

2.8

−5

1.8

2

2.2 2.4 RO...H (Å)

2.6

2.8

0 −0.5

−1.6

−1 −1.5

−2

−2 SAPT FF−SAPT (damping=0.78) AMOEBA (damping=0.39)

−2.5 −3

2.4

2.6

2.8 3 3.2 RO...O (Å)

−2.4 3.4

3.6

−2.8

SAPT FF−SAPT (damping=0.78) AMOEBA (damping=0.39)

−40

−20

0 20 α (deg)

40

60

Figure 4: Induction energy (in kcal/mol) for the three water dimers (top left: linear, top right: planar cyclic, bottom left: bifurcated and bottom right: angular dependence for the linear water dimer with the O...O distance fixed at 2.91 Å, see Figures 1 and 2) from the induced dipole model against the reference values from SAPT. Distances are in Angstrom and angle in degrees. The vertical broken lines mark the approximate equilibrium distances/angle. Please note the difference in the scale of the y-axis. Figure 4, by changing this damping constant only, a drastic improvement is achieved for the linear water dimer, less so for the bifurcated configuration whereas almost no improvement is seen for the planar cyclic arrangement. The angular dependence for the linear water dimer is also shown in Figure 4.

Figure 4 near here Again, the agreement with the reference data is better for the model with the enlarged damping constant (0.78). We should note that the better match with the reference data is mainly due to enlarged oxygen dipoles and quadrupoles used unchanged from the GDMA calculations, whereas the newly refined parameters (atomic multipoles) for water? show almost no angular dependence of the induction energy. Of course, this effective parameter (damping constant) remedies to some extent the deficiency of the induced dipole model used to describe the induction energy component particularly seen at shorter range. It seems that this deficiency is partly due to the lack of anisotropy in describing atomic polarizabilities (mainly for oxygen),123 and, in addition, models based on induced higher order atomic multipoles would probably be necessary as well (see, for example, refs.27,48,60,61,90,124–126 ). We note that using this larger damping factor results in a



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 12 of 28

somewhat larger isotropic dipole polarizability ((αxx + αyy + αzz )/3) of the water molecule (1.62 Å3 ) as compared to the value of 1.52 Å3 obtained using the recommended damping factor (0.39) and also to the value of 1.44 Å3 calculated at the CCSD(T)/aug-cc-pVQZ level of theory (see also ref.125 ). We should also point out here that our reference data were calculated based on the so-called S 2 -approximation (or (2)

single-exchange of electrons) for the evaluation of the second-order exchange-induction term, (Eexch−ind ). Schäffer and Jansen127 have recently shown that at shorter intermolecular separations the S 2 -approximation can underestimate this energy component, and, as a consequence, the total induction energy can be overestimated, so that it becomes less attractive as shown in Figure 4. Recent developments in the DFT-SAPT methodology should clarify this point further.127 The dispersion energy. A comparison between the SAPT reference values for dispersion contributions and our model is given in Figure 5. From this Figure it is seen that the refined parameters reproduce closely the reference data for the three arrangements of the water dimers (Figure 1) and also for the angular dependence for the linear water dimer (Figure 2) as well.

Figure 5 near here Our refinements indicated that including a damping function (see Eq. 7) significantly improves the fit. The rootmean-square deviation drops significantly from 0.99 to 0.14 kcal/mol for water. At the same time we found almost no real improvements in the quality of the fit when combining rules (for unlike atom pairs) were not used.32 So, we decided to keep them as specified in Eq. 8 (because less parameters are needed in our model). We also did not include higher order terms in the asymptotic expansion (long-range correlation) since we believe that the errors introduced by other terms of our force field are definitely larger than those introduced by neglecting these higher terms. We also tested a model where dispersion interactions are included via damped R−6 terms between O atoms only as have been done in recent studies by DeFusco et al.100 and by Nicolini et al.103 Our refined parameters CO−O (∼1600 kcal/mol Å6 ) and βO−O (∼3.0 Å−1 ) were found to be larger than those used by Defusco et al.100 (CO−O = 1300 kcal/mol Å6 and βO−O = 2.23 Å−1 ) and by Nicolini et al.103 (CO−O ≈ 287 kcal/mol Å6 and βO−O = 0.762 Å−1 ). We should note that our results indicate that ignoring interactions involving hydrogen atoms appreciably worsened the quality of the fit. The root-mean-square deviation increased by an order of magnitude. Therefore it was deemed necessary to include all-atom dispersion terms. The parameters used in our model potential to evaluate the dispersion energy are given in Table 2.

Table 2 near here

Even though further validation studies are definitely needed to judge the transferability of the parameters, the results already shown are quite encouraging.


Page 13 of 28

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


0 −1 −2 −3 −4 −5 −6 −7

0 −1 −2 −3 −4 −5 −6 −7

SAPT FF−SAPT 1.6

1.8

2

2.2 2.4 RO...H (Å)

2.6

2.8

SAPT FF−SAPT 1.8

2

2.2 2.4 RO...H (Å)

2.6

2.8

40

60

−2.2 −1 −2

−2.3

−3 −2.4

−4 SAPT FF−SAPT

−5 2.4

2.6

2.8

3 RO...O (Å)

3.2

−2.5 3.4

SAPT FF−SAPT

3.6 −40

−20

0 20 α (deg)

Figure 5: Dispersion energy (in kcal/mol) for the three water dimers (top left: linear, top right: planar cyclic, bottom left: bifurcated and bottom right: angular dependence for the linear water dimer with the O...O distance fixed at 2.91 Å, see Figures 1 and 2) against the reference values from SAPT. Distances are in Angstrom and angle in degrees. The vertical broken lines mark the approximate equilibrium distances/angle.

Table 2: Parameters for the dispersion energy term (see Eq. 6). atom class C, kcal/mol Å6 β, 1/Å O (water) 673 3.439 H (water) 13 4.113



E(1)exch − K S (kcal/mol)

30 K S (kcal/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

25 20 15 10 5 0 0

5 10 15 20 25 30 E(1)exch (kcal/mol)

Page 14 of 28

10 5 0 −5 −10 0

5 10 15 20 25 30 E(1)exch (kcal/mol)

Figure 6: Scatter plots comparing scaled (K) intermolecular electron density overlap integrals (S) with their SAPT exchange-repulsion contributions (left) and their differences (right) for water dimers (in kcal/mol). As pointed out recently by Schäffer and Jansen,128 at shorter intermolecular separations the single-exchange or S 2 (2)

approximation (used here) tends to overestimate the second-order exchange-dispersion component (Eexch−disp ), and, as a consequence, the total dispersion energy term (see Eq. 3) can be underestimated. This is also true for the linear water dimer shown in Fig. 3 of ref.128 This is an important issue here since we parametrize all contributions separately and seek to avoid error compensation with other terms as much as possible. The advantage of the individual parameterization is that it allows to be improved as soon as a better theoretical description of each term is available. The exchange-repulsion energy. Apart from almost all previous studies, we use the first-order exchange-repulsion (1)

energy contribution calculated within SAPT (Eexch ) directly as our reference. We adopted the overlap model which was used in some previous studies48,110,129,130 to approximate the anisotropy of the exchange-repulsion contribution. It is the only term in our force field which is not a pairwise sum over atom-atom contributions between the two monomers but, instead, a monomer-based one. The first-order exchange-repulsion energy contributions for all 2510 water dimer configurations were used in a leastsquares fit to obtain both the scaling parameter (K) and the exponent (a), as defined in Eq. 9. Interestingly, our refined value for a (0.95) is the same as that found before by Elking et al.,130 where the exchange-repulsion reference values are calculated through CSOV decomposition using a much smaller number of dimer geometries. Since this value is very close to 1.0, we use this latter one to approximate the exchange-repulsion contribution in our force field. A fairly good correlation between the values obtained with the intermolecular electron density overlap model (Eq. 9 where a=1.0) and the SAPT reference data (whose values lie between 0 and 30 kcal/mol) can be judged from the scatter plots shown in Figure 6.

Figure 6 near here Figure 7 compares the first-order exchange-repulsion from SAPT with the scaled (K) intermolecular density overlap


Page 15 of 28

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


30

10

25

8

30

(1)

E exch (SAPT) overlap (avtz) overlap (sddall)

25

6

20

(1)

5

20

4

15

15

3

10

10

5

5

1.9

0

1.6

1.95

1.8

2

2

2.2 RO...H (Å)

30 25

5 4 3 2

20 15 10

2.4

2.6

2.8

2.25

1.8

2

9

E(1)exch (SAPT) overlap (avtz) overlap (sddall)

2.3

2.35

2.2 2.4 RO...H (Å)

2.6

2.8

E(1)exch (SAPT) overlap (avtz) overlap (sddall)

8.5 8 7.5

2.9 2.95

3

3.05 3.1 7

5 0

0

E exch (SAPT) overlap (avtz) overlap (sddall)

6.5 2.4

2.6

2.8

3 RO...O (Å)

3.2

3.4

3.6

6

−40

−20

0

20

α (deg)

40

60

Figure 7: The scaled intermolecular density overlap integral (in e2 /˚ A3 ) for the three water dimers (top left: linear, top right: planar cyclic, bottom left: bifurcated and bottom right: angular dependence for the linear water dimer with the O...O distance fixed at 2.91 Å, see Figures 1 and 2) against the reference values from SAPT (in kcal/mol, y-axis). Distances are in Angstrom and angle in degrees. The vertical broken lines mark the approximate equilibrium distances/angle. integrals (see Eq. 9) using both the all-electron (at PBE0/aug-cc-pVTZ level of theory, labeled by “overlap (avtz)”) and valence-only (at PBE0/SDD(O)+D95(H) level of theory, labeled by “overlap (ECP/sddall)”) electron densities. As before, the dependence of the exchange-repulsion energy upon the interatomic distance for the three arrangements of the water dimers is used to judge the quality of our force field. The intermolecular all-electron density overlap integrals calculated with the larger basis set (aug-cc-pVQZ) are found to be almost indistinguishable from those with the smaller one (augcc-pVTZ). Despite the fact that the scaling constants (K) are slightly different for all-electron (∼650) and valence-only (∼670) electron density overlap, both are able to reproduce the SAPT data quite reliably.

Figure 7 near here Two effects are responsible for the peculiar angular dependence of the exchange-repulsion contribution. It grows steadily away from α=0 deg. due to (1) an increase in the overlap between the hydrogen electron density and the lone pair(s) on oxygen and (2) a shortening of the H...H separations between the donor and the acceptor water monomers. It is very encouraging that the angular dependence for the linear water dimer is also qualitatively in line with our reference data (see



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 16 of 28

Figure 7). Our results also corroborate previous findings50,110,131 on the importance of the valence contributions for the determination of the intermolecular exchange-repulsion interaction (see also a discussion in refs.132,133 ). We should note, however, that an attempt to reproduce this angular dependence with the simple exponential function using the parameters fitted to match the corresponding first-order exchange-repulsion energy contributions was completely unsuccessful, see the Supporting Information. This points once again60 to the importance of including the anisotropy to treat this particular term in the force field accurately. Our results indicate that the overlap of the valence-only electron densities will do a good job with the core electrons of the oxygen atom replaced by the effective core potential. With the transferability of the force field parameters in mind, a logical step forward would be (1) to project the valence electron density of a monomer into its atomic contributions and (2) to approximate the first-order exchange-repulsion term between the monomers as an atomic pairwise additive density overlap model along with the appropriate scaling (see, for example, ref.130 ). These developments can be significantly facilitated by a recent study of Cisneros et al.,134 where the calculation of analytical forces for the overlap model has been reported. Total interaction energy. We emphasize, that in order to get the total interaction energy right with an error compensation as small as possible, it is necessary to reproduce accurately all its components, which is, of course, a great challenge. It should be stressed that we seek for a force field the parameters of which are found from a separate fit to the physically meaningful energy components given by SAPT. Our strategy is similar to that used previously by Torheyden and Jansen.48 This is to be contrasted with the approaches where the total energy is used directly to fit the parameters. In a recent study by Duke et al.135 the reduced variational space (RVS)136 method implemented at the Hartree-Fock level is used to produce the reference data to which the parameters of their model potential were fitted. This means that dispersion is not included explicitly and should be taken into account ad hoc. Despite some similarities between our approach and that used by Duke et al.135 (the Thole model for the induction contribution but using a different damping parameter and the overlap model for the exchange-repulsion term but employing the valence-only electron densities by us), the explicit treatment of the dispersion term in our potential is the important difference. It is crucial, for a potential to be useful and transferable, to develop an accurate and physically motivated representation of the dispersion contribution. It has been shown some time ago,137 that even for the water dimer the dispersion contribution is as important as the induction term. Moreover, even though the first-order electrostatic energy is the most attractive contributor to the total energy, it is calculated by utilizing the unperturbed molecular electron densities, thus violating the Pauli principle (with respect to the exchange of electrons between any two molecules). This term is, however, quenched to a large extent by the other first-order term, exchange-repulsion. For neutral species near equilibrium intermolecular geometries, the sum of these two first-order terms is usually small or can be even slightly repulsive (see below) and the second and higher-order terms (induction and dispersion) start to play a dominant role. Furthermore, our parameterization strategy is similar to that developed recently by McDaniel and Schmidt,22 where the DFT-SAPT energy contributions were fitted individually. Besides the different models used to describe the induction contribution (Drude oscillator was used in ref.22 ), our model


Page 17 of 28

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


3 2 1 0 −1 −2 −3 −4 −5 −6

8

SAPT AMOEBA FF−SAPT


6

−4.75

4

−5

2 0

−5.25 1.9

1.95

2

−2 −4 −6

1.6

1.8

2

2.2 RO...H (Å)

8

2.4

2.6

2.8

4

−2.5

2

−3

2

2.2 2.4 RO...H (Å)

2.6

0 20 α (deg)

40

2.8

−4


6

1.8

−4.2 −4.4 −4.6

0

−4.8

2.9 2.95 3 3.05 3.1

−5 SAPT

−2 −4

−5.2 AMOEBA 2.4

2.6

2.8

3

3.2

3.4

3.6

−5.4

RO...O (Å)

FF−SAPT

−40

−20

60

Figure 8: Total interaction energy (in kcal/mol) for the three water dimers (top left: linear, top right: planar cyclic, bottom left: bifurcated and bottom right: angular dependence for the linear water dimer with the O...O distance fixed at 2.91 Å, see Figures 1 and 2). Distances are in Angstrom and angle in degrees. potential, in addition, includes anisotropy for both the long-range electrostatics and exchange-repulsion terms. This last feature was deemed necessary for the accurate treatment of the water dimer (see, for example, ref.60 ). In contrast to the damped point charge model used by McDaniel and Schmidt,,22 our model potential describes electrostatics by distributed atomic multipoles up to quadrupoles. At the same time, to estimate the dispersion contribution, damped higher-order terms have been included in ref.22 In Figure 8 we compare the total interaction energies for the test set obtained from our force field with those calculated by DFT-SAPT. It shows that the agreement is improved in comparison to the AMOEBA force field. The angular dependence of the total energy for the linear water dimer, shown in Figure 8, is reproduced, however, less satisfactory. This results mainly from deficiencies in the used induction model (see Figure 4).

Figure 8 near here Analysis of the components of the binding energy using SAPT suggests that the global minimum configuration corresponding to the linear dimer is a consequence of a very subtle interplay between the two first-order contributions (1)

(1)

(electrostatic, Eelst , and exchange-repulsion, Eexch ). While electrostatics favors the geometrical arrangements away from



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 18 of 28

α=0 deg., where the O-H bond of one monomer points to the lone pair of the oxygen atom of another one (see Fig. 3), exchange-repulsion acts in the opposite direction (see Fig. 7) making such geometrical arrangements significantly less favorable. These two contributions, however, cancel each other almost perfectly (compare to ref.138 ). This is clearly seen in Figure 9, where the energy decomposition is provided with the angular dependence of each contribution plotted separately. A closer look at the separate contributions reveals that a particular configuration, corresponding to an almost linear OH...O arrangement (known as a fingerprint of the hydrogen bonding in water139 ), is a result of the attractive second-order (2)

(2)

(1)

(1)

contributions (Eind and Edisp ), being equally important, modulated by the sum of the first-order terms (Eelst + Eexch ). Moreover, the weak angular dependence of the sum of the first-order terms and of induction and dispersion contributions can nicely explain the ease of forming and breaking of hydrogen bonds in water. We should, therefore, stress that an accurate description of the polarization contributions, both induction and dispersion, is mandatory in order to explain geometrical preferences in all noncovalently interacting systems (see also ref.140 ).

Figure 9 near here Finally, since our work is a kind of proof-of-concept, we focused here on the water dimer (thus concentrating on 2-body effects only). The model parameters were refined using a large number of data points chosen randomly on the potential energy surface of water dimer using rigid monomers. The accuracy of the model can be judged by the success of prediction for a number of 1D cut-offs shown above. We note that our model predicts fairly accurately the distance dependence of the intermolecular energy in the linear and bifurcated configurations, while it is less satisfactory for the cyclic dimer. It is also very challenging to reproduce the angular dependence of each energy contribution in the linear dimer properly (rarely addressed in the literature!). We pointed to a number of deficiencies of our model and to some possible ways for its improvement as well. By a careful comparison of the outcome from our model potential with the respective energy contributions from DFT-SAPT, two main targets for further improvements can be identified, namely, exchange-repulsion and induction. Taken the approximations used in the DFT-SAPT implementation noted above into account, our model potential is still not without some error compensation, but even at the present state-of-the-art, SAPT seems to be a quite reliable method (among many other energy decomposition schemes) helping not only to uncover the physics behind the intermolecular interactions but also to develop accurate force fields for molecular simulations.

4

Conclusions In summary, we proposed here a simple redesign of the popular intermolecular force fields to describe non-bonded

interactions by explicitly including (1) a pairwise additive energy term to approximate the short-range contribution to the electrostatic energy due to charge (density) penetration effects and (2) a damped London term to estimate the dispersion


Page 19 of 28

2

2

Eelst Eexch Edisp Eind Etot

1.5 1 0.5

0

−1

0

−2

−0.5

−3

−1

−4 −40

−20

0 20 α (deg) ESAPT − EFF−SAPT (kcal/mol)

−1.5

Eelst + Eexch Edisp Eind Etot

1 ESAPT (kcal/mol)

∆ ESAPT (kcal/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


0.5 0.4 0.3 0.2 0.1 0 −0.1 −0.2 −0.3 −0.4 −0.5

40

−5

60

−40

−20

40

60

0 20 α (deg)

40

60


−40

−20

0 20 α (deg)

Figure 9: Angular dependence of the separate contributions (in kcal/mol) to the total interaction energy from SAPT (top) and their deviation from the FF-SAPT force field (bottom) for the linear water dimer with the O...O distance fixed at 2.91 Å. The vertical broken line marks the approximate equilibrium angle. Top left: The respective contributions are given relative to their values at 0 degrees. Top right: the absolute values are shown. See Figure 2 for the definition of the angle α.



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 20 of 28

interactions. The intermolecular valence-only electron density overlap model can be used to approximate the exchangerepulsion contribution. Our tests on the water dimers indicated that a more physically sound model potentials with only a small number of additional parameters can readily be constructed. The energy decomposition analysis based on the symmetry-adapted intermolecular perturbation theory allows the parametrization to be performed quite straightforwardly. Moreover, the deficiencies of the individual components of the model potential can be investigated on a term-by-term basis thus leading to a more accurate calculation of the intermolecular energies and forces. It is shown that utilizing a physically intuitive decomposition of the interaction energy into four main components, namely, electrostatics, exchange-repulsion, dispersion and induction, is very appealing for the parametrization of the intermolecular potential. Furthermore, we believe that, based on the physics of the interaction only, one can avoid introducing ad hoc terms such as “hydrogen” (or even “halogen”) bonding energy terms (see, for example, ref.141 ) and produce a set of fairly transferable parameters for next-generation force fields.

Acknowledgements The authors are grateful to the Deutsche Forschungsgemeinschaft (GRK1221, SFB 630) and the Volkswagen Stiftung for financial support. We would like to thank Prof. Krzysztof Szalewicz (Department of Physics and Astronomy, University of Delaware) for making the data sets on SAPT decomposition analysis available to me and Dr. Dennis Elking (OpenEye Scientific Software, Santa Fe) for providing the code to calculate the overlap integrals. We are also deeply indebted to Prof. Bernd Engels (University of Würzburg) for his continuous support throughout the course of this study and for many insightful discussions and fruitful suggestions.

Supporting Information Available: Plots comparing the total interaction energies from DFT-SAPT versus CCSD(T). This material is available free of charge via the Internet at http://pubs.acs.org. Note: The authors declare no competing financial interest.

References 1. Tafipolsky, M.; Schmid, R. J. Phys. Chem. B 2009, 113, 1341–1352. 2. Matsuoka, O.; Clementi, E.; Yoshimine, M. J. Chem. Phys. 1976, 64, 1351–1361. 3. Tsuzuki, S.; Uchimaru, T.; Tanabe, K.; Kuwajima, S. J. Phys. Chem. 1994, 98, 1830–1833. 4. Mooij, W. T. M.; van Duijneveldt, F. B.; van Duijneveldt-van de Rijdt, J. G. C. M.; van Eijck, B. P. J. Phys. Chem. A 1999, 103, 9872–9882.


Page 21 of 28

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


5. Špirko, V.; Engvist, O.; Soldán, P.; Selzle, H. L.; Schlag, E. W.; Hobza, P. J. Chem. Phys. 1999, 111, 572–582. 6. Mitchell, J. B. O.; Price, S. L. J. Phys. Chem. A 2000, 104, 10958–10971. 7. Hloucha, M.; Sum, A. K.; Sandler, S. I. J. Chem. Phys. 2000, 113, 5401–5406. 8. Engkvist, O.; Åstrand, P. O.; Karlström, G. Chem. Rev. 2000, 100, 4087–4108. 9. Donchev, A. G.; Ozrin, V. D.; Subbotin, M. V.; Tarasov, O. V.; Tarasov, V. I. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 7829–7834. 10. Podeszwa, R.; Bukowski, R.; Szalewicz, K. J. Phys. Chem. A 2006, 110, 10345–10354. 11. Li, X.; Volkov, A. V.; Szalewicz, K.; Coppens, P. Acta Cryst. D 2006, 62, 639–647. 12. Hellmann, R.; Bich, E.; Vogel, E. J. Chem. Phys. 2008, 128, 214303. 13. Szalewicz, K.; Leforestier, C.; van der Avoird, A. Chem. Phys. Lett. 2009, 482, 1–14. 14. Archambault, F.; Chipot, C.; Soteras, I.; Luque, F. J.; Schulten, K.; Dehez, F. J. Chem. Theory Comput. 2009, 5, 3022–3031. 15. Wang, B.; Truhlar, D. G. J. Chem. Theory Comput. 2010, 6, 3330–3342. 16. Totton, T. S.; Misquitta, A. J.; Kraft, M. J. Chem. Theory Comput. 2010, 6, 683–695. 17. Wu, J. C.; Piquemal, J.-P.; Chaudret, R.; Reinhardt, P.; Ren, P. J. Chem. Theory Comput. 2010, 6, 2059–2070. 18. Boese, A. D.; Forbert, H.; Masia, M.; Tekin, A.; Marx, D.; Jansen, G. Phys. Chem. Chem. Phys. 2011, 13, 14550– 14564. 19. Leforestier, C.; Tekin, A.; Jansen, G.; Herman, M. J. Chem. Phys. 2011, 135, 234306. 20. Taylor, D. E.; Rob, F.; Rice, B. M.; Podeszwa, R.; Szalewicz, K. Phys. Chem. Chem. Phys. 2011, 13, 16629–16636. 21. Tafipolsky, M.; Engels, B. J. Chem. Theory Comput. 2011, 7, 1791–1803. 22. McDaniel, J. G.; Schmidt, J. R. J. Phys. Chem. A 2013, 117, 2053–2066. 23. Yokogawa, D. J. Chem. Phys. 2012, 137, 204101. 24. Wang, L.-P.; Chen, J.; Van Voorhis, T. J. Chem. Theory Comput. 2013, 9, 452–460. 25. Misquitta, A. J.; Stone, A. J. J. Chem. Theory Comput. 2008, 4, 7–18. 26. Misquitta, A. J.; Stone, A. J.; Price, S. L. J. Chem. Theory Comput. 2008, 4, 19–32. 27. Holt, A.; Boström, J.; Karlström, G.; Lindh, R. J. Comput. Chem. 2010, 31, 1583–1591.



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 22 of 28

28. Wheatley, R. J.; Tulegenov, A. S.; Bichoutskaia, E. Int. Rev. Phys. Chem. 2004, 23, 151–185. 29. Hayes, J. M.; Greer, J. C.; Morton-Blake, D. A. J. Comput. Chem. 2004, 25, 1953–1966. 30. Kaminski, G. A.; Stern, H. A.; Berne, B. J.; Friesner, R. A. J. Phys. Chem. A 2004, 108, 621–627. 31. Arnautova, Y. A.; Jagielska, A.; Pillardy, J.; Scheraga, H. A. J. Phys. Chem. B 2003, 107, 7143–7154. 32. Bordner, A. J.; Cavasotto, C. N.; Abagyan, R. A. J. Phys. Chem. B 2003, 107, 9601–9609. 33. Gavezzotti, A. J. Phys. Chem. B 2003, 107, 2344–2353. 34. Manukyan, A.; Tekin, A. Phys. Chem. Chem. Phys. 2015, 17, 14685–14701. 35. Stone, A. J.; Misquitta, A. J. Int. Rev. Phys. Chem. 2007, 26, 193–222. 36. Fanourgakis, G. S.; Xantheas, S. S. J. Chem. Phys. 2008, 128, 074506. 37. Lao, K. U.; Herbert, J. M. J. Phys. Chem. A 2015, 119, 235–252. 38. Prampolini, G.; Livotto, P. R.; Cacelli, I. J. Chem. Theory Comput. 2015, 11, 5182–5196. 39. Jeziorski, B.; Moszynski, R.; Szalewicz, K. Chem. Rev. 1994, 94, 1887–1930. 40. Stone, A. J. The Theory of Intermolecular Forces; Oxford University Press: Oxford, 2013. 41. Phipps, M. J. S.; Fox, T.; Tautermann, C. S.; Skylaris, C.-K. Chem. Soc. Rev. 2015, 44, 3177–3211. 42. Schmidt, J. R.; Yu, K.; McDaniel, J. G. Acc. Chem. Res. 2015, 48, 548–556. 43. Singh, U. C.; Kollman, P. A. J. Chem. Phys. 1985, 83, 4033–4040. 44. Sokalski, W. A.; Lowrey, A. H.; Roszak, S.; Lewchenko, V.; Blaisdell, J.; Hariharan, P. C.; Kaufman, J. J. J. Comput. Chem. 1986, 7, 693–700. 45. Millot, C.; Stone, A. J. Mol. Phys. 1992, 77, 439–462. 46. Åstrand, P. O.; Wallqvist, A.; Karlström, G. J. Chem. Phys. 1994, 100, 1262–1273. 47. Mas, E. M.; Szalewicz, K.; Bukowski, R.; Jeziorski, B. J. Chem. Phys. 1997, 107, 4207–4218. 48. Torheyden, M.; Jansen, G. Mol. Phys. 2006, 104, 2101–2138. 49. Donchev, A. G.; Galkin, N. G.; Illarionov, A. A.; Khoruzhii, O. V.; Olevanov, M. A.; Ozrin, V. D.; Subbotin, M. V.; Tarasov, V. I. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 8613–8617. 50. Piquemal, J.-P.; Chevreau, H.; Gresh, N. J. Chem. Theory Comput. 2007, 3, 824–837.


Page 23 of 28

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


51. Szalewicz, K. WIREs Comput. Mol. Sci. 2012, 2, 254–272. 52. Jansen, G. WIREs Comput. Mol. Sci. 2014, 4, 127–144. 53. Ansorg, K.; Tafipolsky, M.; Engels, B. J. Phys. Chem. B 2013, 117, 10093–10102. 54. Ponder, J. W.; Wu, C.-J.; Ren, P.-Y.; Pande, V. S.; Chodera, J. D.; Schnieders, M. J.; Haque, I.; Mobley, D. L.; Lambrecht, D. S.; DiStasio, R. A., Jr.; et al., J. Phys. Chem. B 2010, 114, 2549–2564. 55. Tang, K. T.; Toennies, J. P. J. Chem. Phys. 1984, 80, 3726–3741. 56. Mas, E. M.; Bukowski, R.; Szalewicz, K.; Groenenboom, G. C.; Wormer, P. E. S.; van der Avoird, A. J. Chem. Phys. 2000, 113, 6687–6701. 57. Bukowski, R.; Szalewicz, K.; Groenenboom, G. C.; van der Avoird, A. J. Chem. Phys. 2008, 128, 094313. 58. Jankowski, P.; Murdachaew, G.; Bukowski, R.; Akin-Ojo, O.; Leforestier, C.; Szalewicz, K. J. Phys. Chem. A 2015, 119, 2940–2964. 59. Babin, V.; Leforestier, C.; Paesani, F. J. Chem. Theory Comput. 2013, 9, 5395–5403. 60. Millot, C.; Soetens, J. C.; Costa, M. T. C. M.; Hodges, M. P.; Stone, A. J. J. Phys. Chem. A 1998, 102, 754–770. 61. Wikfeldt, K. T.; Batista, E. R.; Vila, F. D.; Jónsson, H. Phys. Chem. Chem. Phys. 2013, 15, 16542–16556. 62. Morawietz, T.; Behler, J. J. Phys. Chem. A 2013, 117, 7356–7366. 63. Góra, U.; Cencek, W.; Podeszwa, R.; van der Avoird, A.; Szalewicz, K. J. Chem. Phys. 2014, 140, 194101. 64. Tröster, P.; Lorenzen, K.; Tavan, P. J. Phys. Chem. B 2014, 118, 1589–1602. 65. van Duijneveldt-van de Rijdt, J. G. C. M.; Mooij, W. T. M.; van Duijneveldt, F. B. Phys. Chem. Chem. Phys. 2003, 5, 1169–1180. 66. Ouyang, J. F.; Bettens, R. P. A. Chimia 2015, 69, 104–111. 67. Van Thiel, M.; Becker, E. D.; Pimentel, G. C. J. Chem. Phys. 1957, 27, 486–490. 68. Smith, B. J.; Swanton, D. J.; Pople, J. A.; Schaefer, H. F.; Radom, L. J. Chem. Phys. 1990, 92, 1240–1247. 69. Tschumper, G. S.; Leininger, M. L.; Hoffman, B. C.; Valeev, E. F.; Schaefer, H. F.; Quack, M. J. Chem. Phys. 2002, 116, 690–701. 70. Reinhardt, P.; Piquemal, J. P. Int. J. Quantum Chem. 2009, 109, 3259–3267. 71. Mukhopadhyay, A.; Cole, W. T. S.; Saykally, R. J. Chem. Phys. Lett. 2015, 633, 13–26.



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 24 of 28

72. Polly, R.; Werner, H. J.; Manby, F. R.; Knowles, P. J. Mol. Phys. 2004, 102, 2311–2321. 73. MOLPRO, version 2010.1, a package of ab initio programs, H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Schütz, and others, see http://www.molpro.net (accessed 02/27/2015). 74. Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M. WIREs Comput. Mol. Sci. 2012, 2, 242–253. 75. Dunning, T. H. J. Chem. Phys. 1989, 90, 1007–1023. 76. Weigend, F. Phys. Chem. Chem. Phys. 2002, 4, 4285–4291. 77. Weigend, F.; Köhn, A.; Hättig, C. J. Chem. Phys. 2002, 116, 3175–3183. 78. Grüning, M.; Gritsenko, O. V.; van Gisbergen, S. J. A.; Baerends, E. J. J. Chem. Phys. 2001, 114, 652–660. 79. Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865–3868. 80. Adamo, C.; Barone, V. J. Chem. Phys. 1999, 110, 6158–6170. 81. Lias, S. G., "Ionization Energy Evaluation" in NIST Chemistry WebBook, NIST Standard Reference Database Number 69, Eds. P.J. Linstrom and W.G. Mallard, National Institute of Standards and Technology, Gaithersburg MD, 20899, http://webbook.nist.gov, (retrieved February 7, 2013). 82. Heßelmann, A. J. Phys. Chem. A 2011, 115, 11321–11330. 83. Podeszwa, R.; Cencek, W.; Szalewicz, K. J. Chem. Theory Comput. 2012, 8, 1963–1969. 84. Singh, N. J.; Min, S. K.; Kim, D. Y.; Kim, K. S. J. Chem. Theory Comput. 2009, 5, 515–529. 85. Boys, S. F.; Bernardi, F. Mol. Phys. 1970, 19, 553–566. 86. Gaussian 09 Revision B.01, Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; et al. Gaussian Inc. Wallingford CT 2010. 87. Korona, T. Mol. Phys. 2013, 111, 3705–3715. 88. Moszynski, R.; Heijmen, T. G. A.; Jeziorski, B. Mol. Phys. 1996, 88, 741–758. 89. Patkowski, K.; Szalewicz, K.; Jeziorski, B. Theor. Chem. Acc. 2010, 127, 211–221. 90. Misquitta, A. J. J. Chem. Theory Comput. 2013, 9, 5313–5326. 91. Stone, A. J.; Alderton, M. Mol. Phys. 1985, 56, 1047. 92. Stone,

A.

J.

J.

Chem.

Theory

Comput.

2005,

1,

1128–1132;

stone.ch.cam.ac.uk/programs.html (accessed 02/27/2015). 93. Spackman, M. A. Chem. Phys. Lett. 2006, 418, 158–162.


GDMA

2.2.04

at

http://www–

Page 25 of 28

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


94. Yáñez, M.; Stewart, R. F.; Pople, J. A. Acta Cryst. A 1978, 34, 641–648. 95. Coppens, P.; Guru Row, T. N.; Leung, P.; Stevens, E. D.; Becker, P. J.; Yang, Y. W. Acta Cryst. A 1979, 35, 63. 96. Schnieders, M. J.; Fenn, T. D.; Pande, V. S.; Brunger, A. T. Acta Cryst. D 2009, 65, 952–965. 97. Wolff, D.; Rudd, W. G. Comput. Phys. Commun. 1999, 120, 20–32. 98. Ng, K. C.; Meath, W. J.; Allnatt, A. R. Mol. Phys. 1979, 37, 237–253. 99. Williams, T.;

Kelley, C.;

many others,

Gnuplot 4.6 patchlevel 3:

an interactive plotting program,

http://gnuplot.sourceforge.net (accessed 02/27/2015), 2013. 100. DeFusco, A.; Schofield, D. P.; Jordan, K. D. Mol. Phys. 2007, 105, 2681–2696. 101. Podeszwa, R.; Pernal, K.; Patkowski, K.; Szalewicz, K. J. Phys. Chem. Lett. 2010, 1, 550–555. 102. Kumar, R.; Wang, F.-F.; Jenness, G. R.; Jordan, K. D. J. Chem. Phys. 2010, 132, 014309. 103. Nicolini, P.; Guardia, E.; Masia, M. J. Chem. Phys. 2013, 139, 184111. 104. Grimme, S. WIREs Comput. Mol. Sci. 2011, 1, 211–228. 105. Hodges, M. P.; Stone, A. J. Mol. Phys. 2000, 98, 275–286. 106. Sanz-Garcia, A.; Wheatley, R. J. Phys. Chem. Chem. Phys. 2003, 5, 801–807. 107. Price, S. L.; Leslie, M.; Welch, G. W. A.; Habgood, M.; Price, L. S.; Karamertzanis, P. G.; Day, G. M. Phys. Chem. Chem. Phys. 2010, 12, 8478–8490. 108. Tinker software tools for molecular design. J. W. Ponder, P. Ren, R. V. Pappu, R. K. Hart, M. E. Hodgson, D. P. Cistola, C. E. Kundrot and F. M. Richards, Washington University School of Medicine, version 6.3, February 2014, at http://dasher.wustl.edu/tinker/ (accessed 02/27/2015). 109. Kim, Y. S.; Kim, S. K.; Lee, W. D. Chem. Phys. Lett. 1981, 80, 574–575. 110. Wheatley, R. J.; Price, S. L. Mol. Phys. 1990, 69, 507–533. 111. Valderrama, E.; Wheatley, R. J. J. Comput. Chem. 2003, 24, 2075–2082. 112. Bergner, A.; Dolg, M.; Küchle, W.; Stoll, H.; Preuß, H. Mol. Phys. 1993, 80, 1431–1441. 113. Dunning Jr., T. H.; Hay, P. J. in Modern Theoretical Chemistry; Ed. H. F. Schaefer III, Plenum, New York, 1977; Vol. 3, pp 1–28. 114. Charbonneau,

P.

Astrophys.

J.

Suppl.

S.

1995,

101,

http://www.hao.ucar.edu/modeling/pikaia/pikaia.php (accessed 02/27/2015).


309–334;

PIKAIA

1.2

at


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 26 of 28

115. Laury, M. L.; Wang, L.-P.; Pande, V. S.; Head-Gordon, T.; Ponder, J. W. J. Phys. Chem. B 2015, 119, 9423–9437. 116. Choi, T. H. Bull. Korean Chem. Soc. 2014, 35, 2906–2910. 117. Wang, B.; Truhlar, D. G. J. Chem. Theory Comput. 2014, 10, 4480–4487. 118. Riley, K. E.; Hobza, P. Acc. Chem. Res. 2013, 46, 927–936. 119. Kumar, P.; Bojarowski, S. A.; Jarzembska, K. N.; Domagala, S.; Vanommeslaeghe, K.; MacKerell, A. D., Jr.; Dominiak, P. M. J. Chem. Theory Comput. 2014, 10, 1652–1664. 120. Toczyłowski, R. R.; Cybulski, S. M. J. Chem. Phys. 2005, 123, 154312. 121. Stone, A. J. The Theory of Intermolecular Forces; Oxford University Press: Oxford, 2013; pp 141–145. 122. Cisneros, G. A. J. Chem. Theory Comput. 2012, 8, 5072–5080. 123. Åstrand, P. O.; Linse, P.; Karlström, G. Chem. Phys. 1995, 191, 195–202. 124. Batista, E. R.; Xantheas, S. S.; Jónsson, H. J. Chem. Phys. 1998, 109, 4546–4551. 125. Elking, D. M.; Perera, L.; Duke, R.; Darden, T.; Pedersen, L. G. J. Comput. Chem. 2011, 32, 3283–3295. 126. Loboda, O.; Ingrosso, F.; Ruiz-López, M. F.; Szalewicz, K.; Millot, C. J. Chem. Phys. 2016, 144, 034304. 127. Schäffer, R.; Jansen, G. Theor. Chem. Acc. 2012, 131, 1235. 128. Schäffer, R.; Jansen, G. Mol. Phys. 2013, 111, 2570–2584. 129. Piquemal, J. P.; Cisneros, G. A.; Reinhardt, P.; Gresh, N.; Darden, T. A. J. Chem. Phys. 2006, 124, 104101. 130. Elking, D. M.; Cisneros, G. A.; Piquemal, J.-P.; Darden, T. A.; Pedersen, L. G. J. Chem. Theory Comput. 2010, 6, 190–202. 131. Cisneros, G. A.; Elking, D.; Piquemal, J.-P.; Darden, T. A. J. Phys. Chem. A 2007, 111, 12049–12056. 132. Hodges, M. P.; Wheatley, R. J. Chem. Phys. Lett. 2000, 326, 263–268. 133. Söderhjelm, P.; Karlström, G.; Ryde, U. J. Chem. Phys. 2006, 124, 244101. 134. Cisneros, G. A.; Piquemal, J.-P.; Darden, T. A. J. Chem. Phys. 2006, 125, 184101. 135. Duke, R. E.; Starovoytov, O. N.; Piquemal, J.-P.; Cisneros, G. A. J. Chem. Theory Comput. 2014, 10, 1361–1365. 136. Stevens, W. J.; Fink, W. H. Chem. Phys. Lett. 1987, 139, 15–22. 137. Jeziorski, B.; Van Hemert, M. Mol. Phys. 1976, 31, 713–729.


Page 27 of 28

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


138. Wang, L.; Gao, J.; Bi, F.; Song, B.; Liu, C. J. Phys. Chem. A 2014, 118, 9140–9147. 139. Chaplin, M. F. (2015). Water structure science. http://www1.lsbu.ac.uk/water (accessed 02/27/2015). 140. Clark, T.; Politzer, P.; Murray, J. S. WIREs Comput. Mol. Sci. 2015, 5, 169–177. 141. Grimme, S. J. Chem. Theory Comput. 2014, 10, 4497–4514.



SAPT: Energy Decomposition Analysis Relative Energy (kcal/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

2


1.5 1 0.5 0 −0.5 −1 −1.5

−40

−20

0 20 α (deg)

40

60

Figure 10: For table of contents only.

TOC


Page 28 of 28

Toward a Physically Motivated Force Field: Hydrogen Bond

Recommend Documents