Nature of Halide–Water Interactions: Insights from Many-Body

ACS2GO © 2019. ← → → ←. loading. To add this web app to the home screen open the browser option menu and tap on Add to homescreen...
0 downloads 0 Views 1MB Size
Subscriber access provided by ALBRIGHT COLLEGE

Quantum Electronic Structure

On the Nature of Halide-Water Interactions: Insights from Many-Body Representations and Density Functional Theory Brandon B. Bizzarro, Colin K. Egan, and Francesco Paesani J. Chem. Theory Comput., Just Accepted Manuscript • Publication Date (Web): 26 Mar 2019 Downloaded from http://pubs.acs.org on March 26, 2019

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

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

Page 1 of 39 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

On the Nature of Halide-Water Interactions: Insights from Many-Body Representations and Density Functional Theory Brandon B. Bizzarro,†,§ Colin K. Egan,∗,†,§ and Francesco Paesani∗,†,‡,¶ †Department of Chemistry and Biochemistry, University of California San Diego, La Jolla, California 92093, United States ‡Materials Science and Engineering, University of California San Diego, La Jolla, California 92093, United States ¶San Diego Supercomputer Center, University of California San Diego, La Jolla, California 92093, United States §Contributed equally to this work E-mail: [email protected]; [email protected]

Abstract Interaction energies of halide-water dimers, X− (H2 O), and trimers, X− (H2 O)2 , with X = F, Cl, Br, and I, are investigated using various many-body models and exchange correlation functionals selected across the hierarchy of density functional theory (DFT) approximations. Analysis of the results obtained with the many-body models demonstrates the need to capture important close-range interactions in the regime of large intermolecular orbital overlap, such as charge transfer and charge penetration. Failure to reproduce these effects can lead to large deviations relative to reference data calculated at the coupled cluster level of theory. Decompositions of interaction energies

1

ACS Paragon Plus Environment

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

carried out with the absolutely localized molecular orbital energy decomposition analysis (ALMO-EDA) method demonstrate that permanent and inductive electrostatic energies are accurately reproduced by all classes of XC functionals (from generalized gradient corrected (GGA) to hybrid and range-separated hybrid functionals), while significant variance is found for charge transfer energies predicted by different XC functionals. Since GGA and hybrid XC functionals predict the most and least attractive charge transfer energies, respectively, the large variance is likely due to the delocalization error. In this scenario, the hybrid XC functionals are then expected to provide the most accurate charge transfer energies. The sum of Pauli repulsion and dispersion energies are the most varied among the XC functionals, but it is found that a correspondence between the interaction energy and the ALMO-EDA total frozen energy may be used to determine accurate estimates for these contributions.

1

Introduction

Molecular dynamics (MD) simulations have become an invaluable tool, potentially allowing for detailed investigation of chemical phenomena that are experimentally inaccessible. 1–9 Presently, the major computational bottleneck associated with MD simulations is handling the electronic degrees of freedom 10 that determine the potential energy landscape of the molecular system of interest. 11 Electronic wave function methods and, in particular, configuration interaction and its approximations, such as the coupled cluster hierarchy, 11,12 offer “chemical” accuracy but scale very poorly with system size, 13 and are therefore computationally intractable in MD simulations. 14,15 One of the most efficient alternatives to wave function methods is the use of analytical potential energy functions (PEFs), 14,15 commonly knows as force fields (FFs), which effectively integrate out the electronic degrees of freedom. The ideal construction of a FF entails generating a lossless representation of the potential energy surface (PES) associated with the molecular system of interest, which is then used to evolve in time the dynamics of the

2

ACS Paragon Plus Environment

Page 2 of 39

Page 3 of 39 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

remaining nuclear degrees of freedom. Unfortunately, the construction of accurate PESs suffers from numerous complications due to the difficulties in describing many-electron systems. For example, eliminating the electronic degrees of freedom requires that FFs be specifically tailored to particular molecular systems to achieve sufficient accuracy, and that each atom be assigned a unique identity that is preserved during the MD simulation. Even when remaining within the above-mentioned restrictions, FF construction is not trivial since there does not exist a unique functional representation that accurately describes the target PES. 14 The most common approach to the development of FFs is to express the total potential energy of a molecular system as a sum of contributions arising from different types of interactions. Intramolecular interactions are often partitioned into terms involving stretching, bending, and dihedral coordinates. Intermolecular interactions tend to be much more difficult to represent, especially in configurations involving large inter-monomer orbital overlap (IMOO) at close range. Intermolecular interactions are usually represented as a sum of individual terms describing different types of physical interactions: permanent electrostatics (ELEC), polarization (POL), charge transfer (CT), London dispersion (DISP), and Pauli repulsion (PAULI). Most common FFs use fixed, atom-centered point charges interacting via Coulomb’s law to represent the ELEC term and relatively simple and computationallyefficient representations of the combined PAULI and DISP terms, such as Lennard-Jones and Buckingham potentials, hereafter defined as P+D terms. 14 These approximations are known to become less accurate at short intermolecular distances, where point charges fail to reproduce effects due to the diffusivity of electron clouds (e.g., charge penetration, CP, effects), and the DISP term cannot be defined unambiguously, being a component of the electron correlation energy. With advances in both MD algorithms and computational power, it has become more common for a FF to include a POL term, usually through the use of atom-centered point-dipole polarizabilities 16 or Drude oscillators. 17 An explicit representation of CT is in general omitted, although attempts have been made to developing effective functional forms representing CT effects. 18–23

3

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 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 advent of explicit many-body PEFs for water, 24–34 which are rigorously derived from many-body expansions of the interaction energies 35 calculated at the coupled cluster level of theory with single, double, and iterative triple excitations, CCSD(T), the current “gold standard” for molecular interactions, 36 has boosted the prospects for realistic MD simulations of aqueous systems. 37 An efficient and accurate approach to constructing explicit manybody PEFs is represented by the MB-nrg methodology, which has so far been applied to pure water systems through the MB-pol PEF 32–34,38 as well as aqueous systems containing halide and alkali-metal ions. 39–44 The MB-nrg PEFs integrate classical representations of the ELEC, DISP, and POL terms, which correctly reproduce interactions at large molecular separations, with data-driven two-body (2B) and three-body (3B) terms that smoothly turn on as monomers approach each other. These 2B and 3B terms are derived from large sets of CCSD(T) reference data and expressed in terms of permutationally invariant polynomials (PIPs) 45 that effectively represent all non-classical interactions (e.g., PAULI and CT terms) and correct for inadequacies associated with a purely classical representation of electrostatic interactions (e.g., CP effects) at close range. In this context, the PIPs can thus be viewed as corrections to purely classical representations of molecular interactions, which effectively describe non-classical interactions resulting from IMOO (i.e., PAULI and CP contributions) and electron delocalization (i.e., CT contributions). 14,46 Density functional theory (DFT) is another common approach to representing molecular interactions in MD simulations. Although, unlike FFs, density functionals are not constructed to reproduce individual types of interactions, several energy decomposition approximations (EDAs) have been proposed which attempt to determine individual contributions to interaction energies calculated with a given density functional. 47 Since there is no unique way to decompose interaction energies, particularly in the regime of large IMOO, 48–50 one must interpret any decomposition as being the consequence of the definitions built into the particular EDA scheme used. In spite of the unavoidable ambiguity in the definitions of individual interaction terms, there are a number of properties that can be satisfied by these

4

ACS Paragon Plus Environment

Page 4 of 39

Page 5 of 39 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

definitions which effectively bound their variability, including: 1. The terms in the decomposition exactly sum to the total interaction energy 2. Individual terms are calculated variationally 3. Individual terms have meaningful basis set limits 4. Rigorously known asymptotic limits of individual terms are recovered in a smooth way The 2nd generation Absolutely Localized Molecular Orbital Energy Decomposition Analysis (ALMO-EDA) method 48 follows this approach, 51 defining terms in order to satisfy a number of “desirable qualities” (including those listed above). In order to provide fundamental insights into the nature of ion–water interactions and assess the ability of different models commonly used in MD simulations to describe these interactions, we report here a systematic analysis of halide–water interactions as predicted by three prototypical implicit and explicit many-body PEFs and various XC functionals selected across the Jacob’s ladder of DFT approximations. 52,53 This analysis sheds light on both strengths and weaknesses of different many-body PEFs in representing individual contributions to halide–water interactions. In addition, ALMO-EDA calculations carried out with a hierarchy of XC functionals, from generalized gradient corrected (GGA) to hybrid and range-separated hybrid functionals, provide the opportunity to conduct detailed comparisons between different XC functionals, which allows for assessing the relative accuracy of different XC approximations to the description of halide–water interactions.

2 2.1

Theoretical and Computational Methodology Electronic structure calculations

All X− (H2 O) dimers and X− (H2 O)2 trimers were optimized at the RI-MP2 54,55 level of theory with the aug-cc-pVTZ 56 basis set. Calculations involving Br− and I− used effective core potentials with 10 and 28 electrons, respectively. 57 A threshold of 1.0 × 10−6 a.u., a step size 5

ACS Paragon Plus Environment

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

Figure 1: Schematic representation of the geometry used in the scans along the RX − O distance carried out for each of the four X− (H2 O) dimers, with X = F, Cl, Br, and I. of 1.0 × 10−6 a.u., and a gradient precision of 1.0 × 10−8 a.u. were used in all optimizations. Scans of each halide ion along the direction of the hydrogen bonded O-H bond of the water molecule were performed by optimizing the water monomer at the RI-MP2/aug-cc-pVTZ level of theory and initially placing the ion 2 ˚ A away from the oxygen atom, with an O-HX− angle of 180o (Fig. 1). Each ion was then displaced along the direction of the hydrogen bonded O-H bond (RX − O ) up to 7 ˚ A from the oxygen, in 0.1 ˚ A increments. Reference interaction energies for all X− (H2 O) dimers and X− (H2 O)2 trimers were calculated at the explicitly correlated coupled cluster level of theory, CCSD(T)-F12b, 58,59 in the complete basis set (CBS) limit that was achieved via a two-point extrapolation 60 between the values obtained with the aug-cc-pVTZ and aug-cc-pVQZ basis sets, 56 using effective core potentials with 10 and 28 electrons 57 for calculations of clusters containing Br− and I− , respectively. All CCSD(T)-F12b and RI-MP2 calculations were performed using MOLPRO, version 2015.1. 61

2.2

The AMOEBA force field

The Atomic Multipole Optimized Energetics for Biomolecular Applications (AMOEBA) is an implicit many-body (i.e., polarizable) FF commonly used in MD simulations. 16,62–71 The AMOEBA functional form can be expressed as

AM OEBA AM OEBA OEBA OEBA VIN = VELEC + VPAM + VPAM OL +D T

6

ACS Paragon Plus Environment

(1)

Page 7 of 39 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

AM OEBA VELEC is modeled using point-multipoles (through the quadrupole) centered on each atom

site, which were calculated via distributed multipole analysis from MP2/aug-cc-pVTZ enerOEBA gies. 62 VPAM is treated explicitly by including self-consistent induced dipoles at atomOL OEBA centered polarizable sites. In order to avoid overpolarization in the close-range, VPAM OL

adopts a Thole-type electrostatic damping scheme, with exponential-3 damping, 62,72 while AM OEBA none of the terms in VELEC are damped. OEBA The remaining term in Eq. 1, VPAM , represented by a buffered 14-7 potential, 73 +D

describes repulsive interactions at close range, which mainly correspond to Pauli repulsion contributions within common energy decomposition schemes, and London dispersion,

VP +D

7 

 1+γ −2 = ij ρ7ij + γ  7   7  1+δ 1+γ 1+δ = ij − 2ij ρij + δ ρ7ij + γ ρij + δ 

1+δ ρij + δ

(2)

with ρij = |Ri − Rj |/Rij . Here, |Ri − Rj | is the distance between atoms i and j, and Rij and ij are the minimum energy distance and the depth of the potential well between atoms i and j, respectively. δ and γ are buffering constants, and serve to damp the divergence of the repulsive R−14 term at short distances. 73 Both δ and γ were originally determined from fits to noble gas data. 73 Values of Rij and ij for homonuclear dimers were fitted to MP2/aug-ccpVTZ energies, and were refined using experimental data. Values for heteronuclear dimers were then determined using the following combination rules: 62

Rij = ij =

(Rii )3 + (Rjj )3 (Rii )2 + (Rjj )2 4ii jj 1/2

1/2

(ii + jj )2

7

ACS Paragon Plus Environment

(3) (4)

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

Page 8 of 39

The TTM-nrg PEF

Within the TTM-nrg PEF, 41,74 water–water interactions are represented by the MB-pol PEF, 32–34 while, similarly to the AMOEBA FF, ion–water interactions are described by the following classical terms: T T M −nrg T T M −nrg T M −nrg T M −nrg T T M −nrg VIN = VELEC + VPTOL + VPTAU + VDISP T LI

(5)

The TTM-nrg model employs an extended Thole-type model 72 describing ion–water perT T M −nrg T M −nrg manent electrostatic (VELEC ) and polarization (VPTOL ) interactions. To guarantee

full compatibility with the MB-pol PEF for water, the TTM-nrg functional form employs smeared induced dipoles (centered at the atom sites) and smeared, geometry-dependent point charges within the exponential-4 smearing scheme introduced in Ref. 75). The charges for each halide ion are set to −1 and the geometry-dependent atom-centered charges are defined as in the MB-pol model to reproduce the reference water dipole moment surface introduced in Ref. 76. The ion dipole polarizabilites were calculated at the CCSD(T) level of theory with the t-aug-cc-pV5Z basis set for F− and Cl− and the t-aug-cc-pV5Z-PP basis set for Br− and I− , following the procedure outlined in Ref. 77. T M −nrg T T M −nrg In the TTM-nrg model, (VPTAU ) and (VDISP ) are represented as sums of pairwiseLI

interactions between the ion and the oxygen (O) and the two hydrogen (H1 and H2 ) atoms T M −nrg of a water molecule. The repulsive terms, VPTAU , describing Pauli repulsion between the LI

ion and the water molecule is expressed as a sum of Born-Mayer functions, T M −nrg VPTAU = AOX ebOX ROX + AH1 X ebH1 X RH1 X + AH2 X ebH2 X RH2 X LI

(6)

where Aαβ and bαβ are fitted parameters for each pair of atom types (α and β), and Rij is T T M −nrg the distance between atoms i and j. VDISP is represented by a sum of Tang-Toennies

8

ACS Paragon Plus Environment

Page 9 of 39 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

dispersion terms f6T T (δαβ ), 78 T T M −nrg VDISP = −f6T T (δXO )

C6,XH C6,XH C6,XO − f6T T (δXH1 ) 6 1 − f6T T (δXH2 ) 6 2 6 RXO RXH1 RXH2

(7)

where C6,αβ is the dispersion coefficient for a given pair of atom types (α and β), and δαβ is a damping parameter, which is set to be equal to the corresponding parameter, bαβ in the corresponding Born-Mayer function (Eq 6). 74 All C6,αβ coefficients were calculated using the exchange-hole dipole model (XDM), 79 as implemented in the Postg software. 80,81

2.4

The MB-nrg PEF

In the MB-nrg PEF, the interactions between individual ions and water molecules are by explicit terms rigorously derived from many-body expansions of the corresponding interaction energies calculated at the CCSD(T)-F12b level of theory. 39,40,42 Briefly, the halide–water MB-nrg PEFs adopt the same functional form and parameterization as the corresponding TTM-nrg PEFs to describe interaction terms associated with permanent electrostatics (i.e., M B−nrg T T M −nrg B−nrg T M −nrg (VELEC = VELEC ), polarization (i.e., VPMOL = VPTOL ), and dispersion energy M B−nrg T T M −nrg 41,74 (i.e., VDISP = VDISP ). However, in the MB-nrg PEFs, these classical terms are M B−nrg M B−nrg supplemented with explicit 2B (V2B ) and 3B (V3B ) terms that are expressed by

permutationally invariant polynomials 45 in variables that are functions of the distances between the ion and the six sites of an MB-pol water molecule, and smoothly reduce to zero beyond a cutoff oxygen-ion distance defined as the distance at which the interactions can accurately be represented by the remaining classical terms. The coefficients of both 2B and 3B permutationally invariant polynomials were optimized using Tikhonov regression (also known as ridge regression) 82 to reproduce reference 2B and 3B energies calculated at the CCSD(T)-F12b level of theory 58,59 in the CBS limit. 39,40 The 2B and 3B MB-nrg permutationally invariant polynomials effectively recover quantummechanical effects in ion–water interactions (e.g., charge transfer and penetration) which

9

ACS Paragon Plus Environment

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

cannot be accounted for by purely classical expressions as those employed in the AMOEBA FF and TTM-nrg PEFs. Specific details on the development of the MB-nrg PEFs can be found in the original references. 39,40

2.5

ALMO-EDA calculations

For a given base functional, E base [·], ALMO-EDA decomposes the interaction energy calculated by that functional, VIN T , into a sum of frozen, polarization, and charge transfer energies, VF RZ , VP OL , and VCT , respectively, 48

VIN T = VF RZ + VP OL + VCT

(8)

Additionally, VF RZ can be further decomposed into a sum of permanent electrostatic, Pauli repulsion, and dispersion energies, VELEC , VDISP , and VP AU LI , respectively, 48,49

VF RZ = VELEC + VDISP + VP AU LI

(9)

Since a detailed description of the ALMO-EDA method is reported in the original references, 48,49,51,83 we only summarize the main aspects of the method here. The three terms on the right-hand side of Eq. 8 are total-system interaction energy contributions, with each term corresponding to its own total-system electron density, each with different degrees of relaxation (optimization). The least-relaxed (highest energy) density, the frozen density, ρF RZ (r), results from the antisymmetric product of the (fully-relaxed) isolated monomer occupied orbitals, superimposed into the complex geometry. The frozen density is “frozen” in the sense that monomer densities are unchanged relative to their optimized, isolated densities, apart from the antisymmetrization. The frozen energy, VF RZ is taken to be the difference from subtracting the sum of the isolated monomer total energies, E1B , from the total energy of

10

ACS Paragon Plus Environment

Page 11 of 39 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

ρF RZ (r) calculated with the base functional. VF RZ = E base [ρF RZ (r)] − E1B = EFtotRZ − E1B

(10)

The decomposition of the frozen energy (Eq. 9) is described in references 48 and 49, but we point out here that the calculation of the dispersion energy, VDISP , requires the use of a “dispersion-free” functional, specific to the base functional, which ideally would be identical to the base functional apart from lacking any dispersion interaction. Therefore, in the interest of fair comparisons between XC functionals, only the sum VP +D = VP AU LI + VDISP is considered in the analyses presented in Section 3. Additionally, we only consider the ALMO-EDA 2 “classical” electrostatic energy 48,49 so that decompositions of XC functionals can be compared with those obtained with many-body models. The polarization energy, VP OL , corresponds to the total-system electron density that has been relaxed (with respect to the base functional) such that each monomer polarizes each other monomer (self-consistently) under the constraint that inter-monomer charge transfer does not occur. In an idealized configuration for which it is possible to uniquely partition the complete set of frozen orbitals between the monomers (i.e. in a configuration with no IMOO) this polarization would correspond to orbital-relaxation resulting from mixing each monomer’s (frozen) occupied orbitals with the same monomer’s (frozen) virtual orbitals in the field of the surrounding monomers. In order to handle real systems, ALMO-EDA 2 determines polarization subspaces for each monomer with the use of fragment electric field response functions, which ensures that polarization in the asymptotic limit is reproduced in a smooth way. 46,48,51 VP OL is then computed by subtracting E1B and VF RZ from the total energy of the polarized density, ρP OL (r) calculated with the base functional. VP OL = E base [ρP OL (r)] − E1B − VF RZ

(11)

The charge transfer energy corresponds to the fully-relaxed total-system electron density, 11

ACS Paragon Plus Environment

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

ρSCF (r), that one would obtain from an unconstrained DFT calculation with the base functional. This electron density can be thought of as resulting from mixing each monomer’s (polarized) occupied orbitals with its own virtual orbitals, and with the occupied and virtual orbitals of the surrounding monomers. 46,48,51 The total interaction energy, VIN T , results from subtracting E1B from the total energy of ρSCF (r) calculated by the base functional, and the charge transfer energy, VCT , is taken to be the remainder from subtracting VF RZ and VP OL from VIN T , tot VIN T = min E base [ρSCF (r)] − E1B = ESCF − E1B

(12)

ρSCF (r)

VCT = VIN T − VF RZ − VP OL

(13)

The ALMO-EDA 2 method 48 was used in all energy decompositions that were carried out using the following XC functionals: BLYP, 84,85 PBE, 86 revPBE, 87 SCAN, 88 TPSS, 53 B97M-rV, 89 B3LYP, 90 PBE0, 91 revPBE0, M06-2X, 92 ωB97X, 93 ωB97X-D, 94 and ωB97M-V. 95 Excluding B97M-rV, ωB97X, ωB97X-D, and ωB97M-V, all DFT energies were calculated with and without the D3(0) empirical dispersion correction. 96 All ALMO-EDA 2 calculations were carried out with Q-Chem and used the aug-cc-pVQZ basis set for all dimers, 56 except for the I− (H2 O) dimer for which an effective core potential for 28 electrons was also used for I− . 57

3 3.1

Results Many-body halide–water models

In the AMOEBA model and both TTM-nrg and MB-nrg PEFs, the point-multipoles and atomic polarizabilities, which determine the ELEC and POL terms, respectively, were defined independently of any fitting to reference interaction energies. Additionally, in the case of both TTM-nrg and MB-nrg PEFs, the C6 coefficients, which determine the DISP term, 12

ACS Paragon Plus Environment

Page 13 of 39 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

were obtained from XDM calculations, independently of any fitting to reference interaction energies (see Section 2), while in the AMOEBA model, the dispersion energy is included in the fitting. Consequently, by construction, the remaining energy contributions (VREM ) in all three models should recover the difference between the reference interaction energies (VREF ) and the sum of these pre-determined terms according to

VREM = VREF − VELEC − VP OL − (VDISP )

(14)

where the parenthesis around VDISP indicates that this term is pre-determined only in the TTM-nrg and MB-nrg PEFs. In the AMOEBA model and TTM-nrg PEF, VREM is expressed as sums of buffered 14-7 and Born-Mayer potentials between pairs of atoms i and j, respectively, each containing two fitting parameters: ij and ρij for each buffered 14-7 potentials in the AMOEBA model, and Aij and bij for each Born-Mayer TTM-nrg (see Section 2). Comparing Eqs. 2 and 6, one can see that the two pairs of fitting parameters have the same effect on their parent potentials, i.e., ij and Aij scale the energy axis, while ρij and bij scale the distance axis. Therefore, given the different functional form, the AMOEBA model and TTM-nrg PEF can differ in their accuracy in representing VREM , despite having the same degree of flexibility. On the other hand, VREM in the MB-nrg PEF is represented by a permutationally invariant polynomial in variables that are functions of the distances between the ion and the six sites of an MB-pol water molecule (see Section 2). It follows from Eqs. 8, 9, and 14 that VREM contains VP AU LI and VCT as well as charge penetration effects that cannot be captured by classical electrostatic models based on Tholetype expressions. As mentioned above, in the case of the AMOEBA model, VREM also contains VDISP . Although the AMOEBA buffered 14-7 and the TTM-nrg Born-Mayer potentials were designed to be limited in scope—modeling VP +D , and VP AU LI , respectively— both potentials were fitted to the total VREM , resulting in the contamination of both potentials with additional close-range interactions that cannot be described by classical expressions

13

ACS Paragon Plus Environment

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

(e.g., charge transfer and penetration). This contamination likely imposes a greater burden on the AMOEBA model because of the additional need to simultaneously represent VDISP with a function of equivalent flexibility. In order to assess the ability of the AMOEBA model and both TTM-nrg and MB-nrg PEFs to describe halide–water interactions, Fig. 2 shows comparisons with CCSD(T)-F12b interaction energies calculated for radial scans obtained by moving each halide ion along the direction of the hydrogen bonded O-H bond of the water molecule, as described in Section 2.1. In all cases, the MB-nrg PEFs quantitatively reproduce the CCSD(T)-F12b data. The TTM-nrg curves are quite accurate in the attractive regions of the scans involving Cl− , Br− and I− , but display deviations of ∼1.3 kcal/mol in the minimum-energy region of the F− (H2 O) scan. Importantly, the TTM-nrg PEFs also show relatively large deviations in the repulsive regions of each scan, consistently being more repulsive than the corresponding CCSD(T)-F12b curves. Since the only difference between MB-nrg and TTM-nrg PEFs is in the representation of VREM , the different accuracy exhibited by the two models directly reflects the ability of the permutationally invariant polynomials adopted by the MB-nrg PEF to represent VP AU LI as well as to effectively capture non-classical contributions (e.g., charge transfer and charge penetration). On the other hand, the AMOEBA model has difficulties in reproducing the CCSD(T)F12b reference values in all four scans, predicting the onset of repulsive interactions at too large intermolecular separations, and often yielding minima that are too attractive. Additionally, the deviations from the CCSD(T)-F12b curves appear to be larger for smaller ions, with the description of the F− (H2 O) scan being the most inaccurate. These deficiencies are likely due to the inability of the buffered 14-7 potential to accurately capture P+D, CT, and CP contributions.

14

ACS Paragon Plus Environment

Page 14 of 39

20 15

5

CCSD(T)-F12b

CCSD(T)-F12b

AMOEBA

AMOEBA

TTM-nrg

TTM-nrg

MB-nrg

MB-nrg

15 10

0

5

-5

0

-10

-5

-15 -20

-10

-25 -30

20





a) F (H2O) 2.0

2.5

3.0

3.5

4.0

4.5

20 15 10

5.0

5.5

b) Cl (H2O) 6.0

2.0

2.5

3.0

3.5

4.0

4.5

5.0

5.5

-15 6.0

CCSD(T)-F12b

CCSD(T)-F12b

AMOEBA

AMOEBA

TTM-nrg

TTM-nrg

MB-nrg

MB-nrg

20 15 10

5

5

0

0

-5

-5

-10

Energy (kcal/mol)

Energy (kcal/mol)

10

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

Journal of Chemical Theory and Computation

Energy (kcal/mol)

Page 15 of 39

-10 −



c) Br (H2O)

-15 2.0

2.5

3.0

3.5

4.0 RX-O (Å)

4.5

5.0

5.5

d) I (H2O) 6.0

2.0

2.5

3.0

3.5

4.0 RX-O (Å)

4.5

5.0

5.5

-15 6.0

Figure 2: Comparisons between CCSD(T)-F12b (open circles) and AMOEBA (red), TTMnrg (green), and MB-nrg (blue) interaction energies calculated for radial scans of X− (H2 O) dimers, with X = F (a), Cl (b), Br (c), and I (d), in which each halide ion is displaced along the direction of the hydrogen bonded O–H bond of the water molecule.

3.2

2B interaction energies

The accuracy of AMOEBA, TTM-nrg and MB-nrg is further assessed through comparisons with various XC functionals in Fig. 3, which shows deviations relative to CCSD(T)F12b/CBS interaction energies calculated for each the four X− (H2 O) dimers in the corresponding minimum-energy geometries obtained at the RI-MP2/aug-cc-pVTZ(-PP) level of theory. As a reference, also shown in Fig. 3 is the ±1 kcal/mol threshold (dashed line) that is generally used to define “chemical accuracy”. Similar comparisons for other halide–water dimers are reported in the Supporting Information. General trends can be identified from the analysis of Fig. 3 which provide guidance in 15

ACS Paragon Plus Environment



a) F (H2O) (ref: -32.79 kcal/mol)



b) Cl (H2O) (ref: -15.58 kcal/mol)

11.0

Page 16 of 39

2.0 1.5

3.0

1.0

1.5

0.5

0.0

0.0 -0.5

-1.5

-1.0 -3.0 −



1.5

d) I (H2O) (ref: -10.99 kcal/mol)

c) Br (H2O) (ref: -13.35 kcal/mol)

2.5 2.0 1.5

1.0

1.0

0.5

0.5

SS

AN

TP

-1.5

SC

D3 )

PB

E(-

D3 rev

(-D

E(-

PB

BL YP

(-D 3) (-D 3) B9 7M -rV B3 LY P(D3 PB ) E0 (-D rev 3) PB E0 (-D M0 3) 62 X(D3 ) ωB 97 X-D ωB 97 MV MB -nr g TT Mnrg AM OE BA

-1.0 )

-0.5

-1.0 3)

0.0

-0.5

P(D3 ) PB E(D3 rev ) PB E(D3 TP ) SS (-D 3) SC AN (-D 3) B9 7M -rV B3 LY P(D3 PB ) E0 (-D rev 3) PB E0 (-D M0 3) 62 X(D3 ) ωB 97 X-D ωB 97 MV MB -nr g TT Mnrg AM OE BA

0.0

Emodel - ECCSD(T) (kcal / mol) Emodel - ECCSD(T) (kcal / mol)

4.5

BL Y

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

Emodel - ECCSD(T) (kcal / mol) Emodel - ECCSD(T) (kcal / mol)

Journal of Chemical Theory and Computation

Figure 3: Deviations from the CCSD(T)-F12b interaction energies of the X− (H2 O) dimers, with X = F (a), Cl (b), Br (c), and I (d), in the corresponding RI-MP2/aug-cc-pVTZ(-PP) minimum-energy configurations calculated with various XC functionals as well as AMOEBA, TTM-nrg, and MB-nrg models. assessing the accuracy of different models in describing halide–water interactions. Among the GGA functionals considered in this study, BLYP and revPBE without D3 corrections significantly underestimate the interaction strength in all halide–water dimers, with deviations from the CCSD(T)-F12b reference values ranging from +1.1 kcal/mol for I− (H2 O) to +3.0 kcal/mol for F− (H2 O). Inclusion of the D3 corrections improves the accuracy of both XC functionals, although the interaction energies predicted for F− (H2 O) still remain outside the boundaries of “chemical accuracy”. Interestingly, PBE, the other GGA functional included in this analysis, accurately reproduces the CCSD(T)-F12b interaction energy for the fluoride–water dimer, which seems to be fortuitous, considering that its accuracy significantly deteriorates for the other three halide–water dimers. Contrary to BLYP and revPBE, the addition of the D3 correction has a deleterious effect on the accuracy of the PBE functional, with the deviations from the CCSD(T) reference data increasing progressively from ∼0.3 kcal/mol for F− (H2 O) to ∼1.0 kcal/mol for I− (H2 O). Among the three meta-GGA functionals, SCAN systematically overestimates the interaction strengths in all halide–water dimers, with the D3 correction being always relatively 16

ACS Paragon Plus Environment

Page 17 of 39 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

small (∼0.1 kcal/mol). An opposite trend is predicted by TPSS, another meta-GGA functional, which underestimates the interaction strengths in all halide–water dimers, although the deviations from the CCSD(T)-F12b reference values are always within 1 kcal/mol. Including the D3 correction has relatively large effects on the accuracy of the TPSS functional, leading to overbound dimers, with the exception of F− (H2 O). The last meta-GGA functional considered in this study, B97M-rV, performs significantly better, on average, than both SCAN(-D3) and TPSS(-D3), with its accuracy improving from F− (H2 O) to I− (H2 O). In general, all hybrid GGA functionals considered in this study (B3LYP, PBE0, and revPBE0) exhibit relatively higher accuracy when compared to their parent GGA models, although they follow similar trends. Both B3LYP and revPBE0 underestimate the interaction strength in all halide–water dimers, with the associated deviations from the CCSD(T)F12b reference values being between ∼1.5 kcal/mol and ∼0.5 kcal/mol, respectively, and decreasing significantly after inclusion of the D3 corrections. On the other side, PBE0 systematically overestimates the interaction strengths in all halide–water dimers and inclusion of the D3 corrections increases the deviations from the CCSD(T)-F12b value (up to −1.5 kcal/mol for the fluoride–water dimer). The meta-GGA and hybrid M06-2X functional is not able to correctly reproduce the interaction energy in the F− (H2 O) dimer, deviating by more than −2.0 kcal/mol from the CCSD(T)-F12b reference value, although its accuracy noticeably improves for the other halide–water dimers. Interestingly, the D3 dispersion correction has little effect on the performance of M06-2X. This can be traced back to the ability of M06-2X to capture “mediumrange dispersion-like” interactions. 92 The range-separated hybrid, meta-GGA functional, ωB97M-V, on average, provides the closest agreement with the CCSD(T)-F12b reference values among all XC functionals considered in this analysis. In the same family, ωB97X-D, is only marginally less accurate than ωB97M-V overall, providing the closest agreement with the CCSD(T)-F12b reference values for F− (H2 O) and I− (H2 O) dimers. AMOEBA and TTM-nrg predict interaction energies that are within “chemical accuracy”

17

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

for all dimers except F− (H2 O), which is significantly underbound in both models. These deviations can be attributed to the inability of purely classical induction schemes adopted by both AMOEBA and TTM-nrg to correctly account for quantum-mechanical effects such as charge transfer and penetration, and Pauli repulsion, which are emphasized in F− (H2 O) by the covalent-like nature of the interactions at short and medium ranges. Fig. 3 shows that the MB-nrg PEF provides the most accurate description of the interaction energies of all X− (H2 O) dimers among all XC functionals and many-body models analyzed in this study, with the largest deviation being +0.18 kcal/mol for F− (H2 O). Considering that both TTMnrg and MB-nrg adopt the same description of all classical contributions to the halide–water interactions, the different performance of these two many-body PEFs must be related to their different ability to represent VREM in Eq. 14.

3.3

3B interaction energies

Deviations of 3B interaction energies calculated by each method relative to CCSD(T)-F12b

3.0 2.5

1.5





b) Cl (H2O)2 (ref: 0.97 kcal/mol)

a) F (H2O)2 (ref: 4.56 kcal/mol)

1.0

2.0 1.5

0.5

1.0 0.5

0.0

0.0

-0.5

-0.5

1.5

Emodel - ECCSD(T) (kcal / mol) Emodel - ECCSD(T) (kcal / mol)

for 3B interaction energies of halide–water global minimum energy trimers are shown in

1.5





d) I (H2O)2 (ref: 0.23 kcal/mol)

c) Br (H2O)2 (ref: 0.57 kcal/mol)

BA OE

AM

g

nrg M-

TT

V

MB

-nr

D

M97

ωB

X

X97

ωB

E0

62

M0

rev PB

P

E0

PB

-rV

LY

B3

AN

7M

SC

B9

E

SS TP

E

PB rev

YP

PB

BL

MV MB -nr g TT Mnrg AM OE BA

97 XD

97

ωB

ωB

E0

62 X

PB

M0

rev

PB

LY

B3

7M

SC

B9

TP

PB rev

PB

P

-0.5

E0

-0.5

-rV

0.0

AN

0.0

E

0.5

SS

0.5

E

1.0

YP

1.0

BL

Emodel - ECCSD(T) (kcal / mol) Emodel - ECCSD(T) (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

Page 18 of 39

Figure 4: Deviations from the CCSD(T)-F12b 3B energies of the X− (H2 O)2 trimers, with X = F (a), Cl (b), Br (c), and I (d), in the corresponding RI-MP2/aug-cc-pVTZ(-PP) minimumenergy configurations calculated with various XC functionals as well as AMOEBA, TTM-nrg, and MB-nrg models.

18

ACS Paragon Plus Environment

Page 19 of 39 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

Fig. 4. Similar comparisons for other halide–water trimers are reported in the Supporting Information. All XC functionals display deviations within 1 kcal/mol for each trimer, such that the least accurate 3B energies are more repulsive than CCSD(T)-F12b. Overall, revPBE0, ωB97X-D, and ωB97M-V provide the most accurate 3B energies, while PBE is consistently the least accurate, displaying deviations greater than 0.5 kcal/mol for each trimer. TTM-nrg deviates from CCSD(T)-F12b by +1.2 kcal/mol for F− (H2 O)2 , but remains within 0.25 kcal/mol of our reference 3B energies for all other trimers. Similarly, AMOEBA has trouble with F− (H2 O)2 , yielding a 3B energy ∼3.0 kcal/mol more repulsive than that of CCSD(T)-F12b. It should be noted that the 3B energy of both AMOEBA and TTM-nrg are due entirely to 3B POL. Among both XC functionals and many-body models included in Fig. 4, MB-nrg provides the closest agreement with the CCSD(T)-F12b reference values, displaying deviations no greater than 0.1 kcal/mol for each trimer. Considering the MBnrg shares the same VELEC , VP OL , and VDISP as TTM-nrg, this level of accuracy provides further evidence for the ability of the 3B PIP adopted by MB-nrg to quantitatively reproduce non-classical 3B effects at close range.

3.4

Energy decomposition analysis

To provide fundamental insights into the ability of both many-body PEFs and DFT models to describe halide–water interactions, the ALMO-EDA method is used to decompose the interaction energies of all X− (H2 O) dimers, calculated by each XC functional, into individual contributions defined in Section 2.5. Additionally, ALMO-EDA analysis within families of XC functionals offers some perspective regarding the construction of a DFT model. For example, PBE and revPBE differ from each other by the value of a single parameter in their exchange functionals (specifically in their exchange enhancement factors), 87 with the revPBE exchange energy being more sensitive to the electron density gradient, and having a lower bound in the large-gradient limit, effectively making revPBE more non-local. 97 Consequently, 19

ACS Paragon Plus Environment

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

Figure 5: Energy decompositions for each minimum energy X− (H2 O) dimer, with X = F (a), Cl (b), Br (c), and I (d), with numerical energy values (in kcal/mol) shown on the right-hand side of each bar.

20

ACS Paragon Plus Environment

Page 20 of 39

Page 21 of 39 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

differences in energy decompositions between PBE and revPBE directly reflect the nature of this single parameter. Moreover, the hybrid functionals, PBE0 and revPBE0 differ from their parent functionals in that their respective GGA exchange term is scaled by 0.75, and Hartree-Fock exchange (scaled by 0.25) is added. Therefore, comparing decompositions of PBE with those of PBE0, and revPBE with revPBE0 will reveal the nature of taking GGA to “GGA0.” Fig. 5, shows that all DFT models tend to provide similar values for the ELEC and POL terms for each halide–water dimer, with the largest standard deviation (among all XC functionals) being 0.3 kcal/mol for both contributions to the interaction energy of the F− (H2 O) dimer. Consensus among DFT models belonging to all classes suggests that ELEC and POL energies of halide–water interactions are weakly-dependent on non-local effects associated with including the kinetic energy density and/or Hartree-Fock exchange in the XC functional. The average values of the charge transfer energy predicted by each class of XC functionals tend to become less negative from the GGA to the hybrid functionals. This progression can be explained by considering the delocalization error, 98 which is mostly due to the selfinteraction error, which is particularly significant for the GGA functionals. This finding is consistent with earlier work studying intermolecular energy decompositions. 99–101 Importantly, all three GGA functionals predict similar values for the CT term, with a standard deviation of at most 0.2 kcal/mol in the case of the F− (H2 O) dimer. The largest disagreement among all functionals considered in this study is found for the P+D term, and as with all other contributions to the total interaction energy, the values for this term vary the most for the F− (H2 O) dimer. Notably, differences between PBE and revPBE almost entirely manifest in their P+D energies, with the revPBE P+D term being at least 1 kcal/mol more repulsive than that of PBE for each minimum energy dimer. BLYP and revPBE differ in their P+D terms by at most 0.4 kcal/mol for all dimers. On average, P+D energies from hybrid functionals are less-repulsive than those from the GGAs, but

21

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

B3LYP and revPBE0 P+D terms are consistently the most repulsive among the hybrids (note that B3LYP uses BLYP exchange scaled by a factor of 0.72 90 ). It is interesting to note that P+D energies calculated with PBE0 and revPBE0 are closer in magnitude than those calculated with PBE and revPBE by a factor of ∼0.75, reflecting the scaling of the PBE and revPBE exchange functionals in PBE0 and revPBE0, respectively. Since, as mentioned above, the disagreement in the CT terms predicted by different classes of XC functionals can be traced back to the delocalization error, it seems reasonable to assume that, among all classes, the hybrid functionals provide the most accurate value -3 −



a) F (H2O)

-4

BLYP PBE revPBE SCAN TPSS B97M-rV B3LYP PBE0 revPBE0 M06-2X ωB97X ωB97X-D ωB97M-V Best Fit Hybrid Fit

4

2

1.63 kcal/mol

1.05 kcal/mol

0

-2 -3.0 -3

-2.0

-1.0

0.0

1.0

2.0

3.0

-5

-5.96 kcal/mol

-6 -6.40 kcal/mol

Total Frozen Energy (kcal/mol)

b) Cl (H2O)

-7

4.0

-1.0



-0.5

0.0

0.5

1.0

1.5

2.0

-8 -3



c) Br (H2O)

d) I (H2O)

-4

-4

-5 -5.28 kcal/mol

-5

-5.84 kcal/mol

-6

-5.60 kcal/mol

-6

-6.22 kcal/mol

Total Frozen Energy (kcal/mol)

Total Frozen Energy (kcal/mol)

6

Total Frozen 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

Page 22 of 39

-7 -7

-1.0

-0.5 0.0 0.5 1.0 1.5 Interaction Energy Error Relative to CCSD(T)-F12b (kcal/mol)

-1.0

-0.5 0.0 0.5 1.0 Interaction Energy Error Relative to CCSD(T)-F12b (kcal/mol)

1.5

Figure 6: Total frozen energy from ALMO-EDA plotted against interaction energy deviations (relative to CCSD(T)-F12b) for X− (H2 O) dimers with X = a) F, b) Cl, c) Br, d) I, for each XC functional, demonstrating the correlation between the two quantities. The best fit lines with respect to all functionals are shown as red solid lines, and the best fit lines with respect to all hybrid functionals are shown as black solid lines. The y-intercepts for each best fit line are also shown, giving an estimate of the “most accurate” total frozen energies for each dimer.

22

ACS Paragon Plus Environment

Page 23 of 39 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

for this contribution. Therefore, since both ELEC and POL contributions are of similar magnitude among all XC functionals, it appears that, within the ALMO-EDA scheme, the P+D term is the most uncertain contribution to the interaction energies of the halide– water dimers. Given the similarity between Eqs. 8 and 9, it seems likely that the frozen energy (P+D + ELEC) should be correlated with the total interaction energy. To test this hypothesis, Fig. 6 shows correlation plots between the total frozen energies predicted by the different functionals for each X− (H2 O) dimer and the corresponding deviations from the CCSD(T)-F12b total interaction energies. A noticeable correlation indeed exists between these two quantities among all functionals, which becomes even more apparent if only the hybrid functionals are considered. From linear fits to the data shown in Fig. 6, it is thus possible to derive an estimate for the “most accurate” total frozen energies associated with each halide–water dimer. Fig. 5 also shows energy decompositions calculated by the three many-body models studied here. Comparing decompositions of XC functionals with those of many-body models term-by-term demonstrates significant differences between the two approaches. Keeping in mind that low-energy halide–water dimers have small intermolecular separations, these differences are attributed to the failure of simple functional forms adopted by AMOEBA and TTM-nrg in capturing complicated close-range interactions. In contrast, the MB-nrg models accurately reproduce CCSD(T)-F12b interaction energies for all four halide–water systems, demonstrating the ability of PIPs to recover all these close-range effects. In order to further emphasize that the differences between many-body models and DFT energy decomposition terms is due to close-range effects such as CT and CP, energy decompositions are calculated for all many-body models as well as ωB97M-V along the scans described in Sections 2.1 and 3.1. ELEC and POL energies calculated along the F− (H2 O) and I− (H2 O) scans are shown in Fig. 7. It should be noted that both ELEC and POL terms are common between MB-nrg and TTM-nrg. These scans show that AMOEBA and MB-nrg converge to the ωB97M-V value for ELEC

23

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

-10



0



a) F (H2O) ELEC

b) I (H2O) ELEC -10

Energy (kcal/mol)

-20

-20

-30

-30

-40

-40

-50 ωB97M-V ELEC MB-nrg ELEC AMOEBA ELEC

-60

ωB97M-V ELEC MB-nrg ELEC AMOEBA ELEC

-50

Energy (kcal/mol)

0

-60

-70 -70 -80 2.0 0

2.5

3.0

3.5

4.0

4.5



2.5

3.0

3.5

4.0

4.5

0



c) F (H2O) POL

5.0

d) I (H2O) POL -10

-20 -20 -30

-30

ωB97M-V POL MB-nrg POL AMOEBA POL

ωB97M-V POL MB-nrg POL AMOEBA POL

-50

-40

2.0

-40

Energy (kcal/mol)

-10 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

Page 24 of 39

2.5

3.0

3.5

4.0

4.5

2.5

3.0

RX-O (Å)

3.5 RX-O (Å)

4.0

4.5

-60 5.0

Figure 7: ELEC and POL energies calculated with each of the three many-body models analyzed in this study and the ωB97M-V functional along the scans described in Sections 2.1 and 3.1 for the F− (H2 O) (left) and I− (H2 O) (right) dimers. at large distances but deviate significantly at close range. These discrepencies in the ELEC contribution, where the values produced by both MB-nrg and AMOEBA are too repulsive relative to ωB97M-V, can be explained by attractive close-range effects, like CP, which are absent from the electrostatic terms adopted by many-body models but are included in the ALMO-EDA ELEC energy. It should be noted that the deviations between MB-nrg and AMOEBA ELEC energies below 2.5–3.0 ˚ A are due to Thole damping which is absent in the AMOEBA model. Comparing the POL energies along these scans, deviations of many-body models from ωB97M-V again occur at close range. However, it appears that these differences are an artifact of the Thole damping implemented in both MB-nrg and AMOEBA to avoid the polarization catastrophe, which is a consequence of using fixed atomic

24

ACS Paragon Plus Environment

Page 25 of 39 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

polarizabilities, an approximation which fails to describe the finite extent of the electron density. 102 DFT models do not suffer from the polarization catastrophe given that no such approximation is made. Although Fig. 5 shows that MB-nrg differs greatly from ωB97M-V in its ELEC and POL energies at close range, in addition to the lack of an explicit CT term, Fig. 2 demonstrates that the MB-nrg PIPs provide QCs for all close-range effects, allowing for a lossless representation of CCSD(T)-F12b interaction energies. Furthermore, Fig. 3 shows that MB-nrg is able to accurately reproduce these effects more consistently than most common XC functionals. Statistics for DFT energy decompositions of X− (H2 O)2 clusters can be found in the Supporting Information. In general, the contribution energy values within and between all classes agree for all ion configurations. For each trimer, the total 3B interaction energy is predominantly composed of POL and charge transfer. The 3B energy of AMOEBA and TTM-nrg consists entirely of 3B POL, while the MB-nrg 3B energy consists of 3B PIP QCs as well as TTM-nrg POL, so their energy decompositions follow directly from their interaction energies shown in Fig. 4.

4

Conclusions

We have analyzed the ability of prototypical implicit and explicit many-body models as well as various XC functionals selected across the hierarchy of DFT approximations to describe the interactions between halide ions and water through the decomposition of interaction energies of halide–water dimers and trimers into their fundamental physical contributions. The energy decompositions of XC functionals have been carried out within the ALMO-EDA scheme and the accuracy of the different models has been assessed through comparisons with CCSD(T)-F12b reference interaction energies. Our analysis indicates that many-body models that do not directly attempt to capture close-range effects in the regime of large inter-monomer orbital overlap are unable to accurately reproduce CCSD(T)-F12b interaction

25

ACS Paragon Plus Environment

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

energies. In addition, while XC functionals belonging to different classes of DFT approximations tend to provide similar descriptions of permanent electrostatics and polarization energies, they differ significantly in their ability to represent Pauli repulsion, dispersion, and charge transfer contributions. The large variability in representing charger transfer by different XC functionals is largely due to the delocalization error, which is particularly large in GGA functionals. Finally, we have identified a seemingly-meaningful correlation between the total interaction energy and total frozen energy among all XC functionals, which allows for estimating the magnitude of the energy term representing the sum of Pauli repulsion and London dispersion from energy decompositions carried out with various XC functionals simultaneously.

Supporting Information Available Cartesian coordinates of low-lying isomers of X− (H2 O) and X− (H2 O)2 complexes, with X = F, Cl, Br, and I, along with additional analyses of 2B and 3B interactions.

This material

is available free of charge via the Internet at http://pubs.acs.org/.

Acknowledgement This research was supported by the National Science Foundation through grant CHE-1453204. All calculations used resources of the Extreme Science and Engineering Discovery Environment (XSEDE), 103 which is supported by the National Science Foundation through grant ACI-1053575, under allocation TG-CHE110009, the High Performance Computing Modernization Program (HPCMP) through grant FA9550-16-1-0327 by the Air Force Office of Scientific Research, and the Triton Shared Computing Cluster (TSCC) at the San Diego Supercomputer Center.

26

ACS Paragon Plus Environment

Page 27 of 39 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

References (1) van Gunsteren, W. F.; Berendsen, H. J. C. Computer Simulation of Molecular Dynamics: Methodology, Applications, and Perspectives in Chemistry. Angew. Chem., Int. Ed. Engl. 1990, 29, 992–1023. (2) Warshel, A. Computer Simulations of Enzyme Catalysis: Methods, Progress, and Insights. Annu. Rev. Biophys. Biomol. Struct. 2003, 32, 425–443. (3) Garcia-Viloca, M.; Gao, J.; Karplus, M.; Truhlar, D. G. How Enzymes Work: Analysis by Modern Rate Theory and Computer Simulations. Science 2004, 303, 186–195. (4) Yamakov, V.; Wolf, D.; Phillpot, S. R.; Mukherjee, A. K.; Gleiter, H. DeformationMechanism Map for Nanocrystalline Metals by Molecular-Dynamics Simulation. Nat. Mater. 2004, 3, 43. (5) Nielsen, S. O.; Lopez, C. F.; Srinivas, G.; Klein, M. L. Coarse Grain Models and the Computer Simulation of Soft Materials. J. Phys.: Condens. Matter 2004, 16, R481. (6) Markvoort, A. J.; Hilbers, P. A. J.; Nedea, S. V. Molecular Dynamics Study of the Influence of Wall-Gas Interactions on Heat Flow in Nanochannels. Phys. Rev. E 2005, 71, 066702. (7) Karplus, M.; Kuriyan, J. Molecular Dynamics and Protein Function. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 6679–6685. (8) Voth, G. A. Computer Simulation of Proton Solvation and Transport in Aqueous and Biomolecular Systems. Acc. Chem. Res. 2006, 39, 143–150. (9) De Vivo, M.; Masetti, M.; Bottegoni, G.; Cavalli, A. Role of Molecular Dynamics and Related Methods in Drug Discovery. J. Med. Chem. 2016, 59, 4035–4061. (10) Tuckerman, M. Statistical Mechanics: Theory and Molecular Simulation; Oxford university press, 2010. 27

ACS Paragon Plus Environment

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

(11) Szabo, A.; Ostlund, N. S. Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory; Courier Corporation, 2012. (12) Bartlett, R. J.; Musial, M. Coupled-Cluster Theory in Quantum Chemistry. Rev. Mod. Phys. 2007, 79, 291. (13) Head-Gordon, M.; Artacho, E. Chemistry on the Computer. Phys. Today 2008, 61, 58. (14) Mao, Y.; Demerdash, O.; Head-Gordon, M.; Head-Gordon, T. Assessing Ion–Water Interactions in the AMOEBA Force Field Using Energy Decomposition Analysis of Electronic Structure Calculations. J. Chem. Theory Comput. 2016, 12, 5422–5437. (15) Schatz, G. C. The Analytical Representation of Electronic Potential-Energy Surfaces. Rev. Mod. Phys. 1989, 61, 669. (16) Ponder, J. W.; Wu, C.; Ren, P.; Pande, V. S.; Chodera, J. D.; Schnieders, M. J.; Haque, I.; Mobley, D. L.; Lambrecht, D. S.; DiStasio Jr, R. A.; ; Head-Gordon, M.; Clark, G. N. I.; Johnson, M. E.; Head-Gordon, T. Current Status of the AMOEBA Polarizable Force Field. J. Phys. Chem. B 2010, 114, 2549–2564. (17) Anisimov, V. M.; Lamoureux, G.; Vorobyov, I. V.; Huang, N.; Roux, B.; MacKerell, A. D. Determination of Electrostatic Parameters for a Polarizable Force Field Based on the Classical Drude Oscillator. J. Chem. Theory Comput. 2005, 1, 153–168. (18) Millot, C.; Soetens, J.-C.; Martins Costa, M. T.; Hodges, M. P.; Stone, A. J. Revised Anisotropic Site Potentials for the Water Dimer and Calculated Properties. J. Phys. Chem. A 1998, 102, 754–770. (19) Van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard, W. A. ReaxFF: a Reactive Force Field for Hydrocarbons. J. Phys. Chem. A 2001, 105, 9396–9409.

28

ACS Paragon Plus Environment

Page 28 of 39

Page 29 of 39 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

(20) Chelli, R.; Pagliai, M.; Procacci, P.; Cardini, G.; Schettino, V. Polarization Response of Water and Methanol Investigated by a Polarizable Force Field and Density Functional Theory Calculations: Implications for Charge Transfer. J. Chem. Phys. 2005, 122, 074504. (21) Youngs, T. G. A.; Hardacre, C. Application of Static Charge Transfer Within an IonicLiquid Force Field and Its Effect on Structure and Dynamics. ChemPhysChem 2008, 9, 1548–1558. (22) Lee, A. J.; Rick, S. W. The Effects of Charge Transfer on the Properties of Liquid Water. J. Chem. Phys. 2011, 134, 184507. (23) Soniat, M.; Rick, S. W. The Effects of Charge Transfer on the Aqueous Solvation of Ions. J. Chem. Phys. 2012, 137, 044511. (24) Bukowski, R.; Szalewicz, K.; Groenenboom, G. C.; van der Avoird, A. Predictions of the Properties of Water from First Principles. Science 2007, 315, 1249–1252. (25) Bukowski, R.; Szalewicz, K.; Groenenboom, G. C.; van der Avoird, A. Polarizable Interaction Potential for Water from Coupled Cluster Calculations. I. Analysis of Dimer Potential Energy Surface. J. Chem. Phys. 2008, 128, 094313. (26) Bukowski, R.; Szalewicz, K.; Groenenboom, G. C.; van der Avoird, A. Polarizable Interaction Potential for Water from Coupled Cluster Calculations. II. Applications to Dimer Spectra, Virial Coefficients, and Simulations of Liquid Water. J. Chem. Phys. 2008, 128, 094314. (27) Wang, Y.; Shepler, B. C.; Braams, B. J.; Bowman, J. M. Full-Dimensional, Ab Initio Potential Energy and Dipole Moment Surfaces for Water. J. Chem. Phys. 2009, 131, 054511.

29

ACS Paragon Plus Environment

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

(28) Wang, Y.; Huang, X.; Shepler, B. C.; Braams, B. J.; Bowman, J. M. Flexible, Ab Initio Potential, and Dipole Moment Surfaces for Water. I. Tests and Applications for Clusters up to the 22-mer. J. Chem. Phys. 2011, 134, 094509. (29) Wang, Y.; Bowman, J. M. Ab Initio Potential and Dipole Moment Surfaces for Water. II. Local-Monomer Calculations of the Infrared Spectra of Water Clusters. J. Chem. Phys. 2011, 134, 154510. (30) Babin, V.; Medders, G. R.; Paesani, F. Toward a Universal Water Model: First Principles Simulations from the Dimer to the Liquid Phase. J. Phys. Chem. Lett. 2012, 3, 3765–3769. (31) Medders, G. R.; Babin, V.; Paesani, F. A Critical Assessment of Two-Body and ThreeBody Interactions in Water. J. Chem. Theory Comput. 2013, 9, 1103–1114. (32) Babin, V.; Leforestier, C.; Paesani, F. Development of a “First Principles” Water Potential with Flexible Monomers: Dimer Potential Energy Surface, VRT Spectrum, and Second Virial Coefficient. J. Chem. Theory Comput. 2013, 9, 5395–5403. (33) Babin, V.; Medders, G. R.; Paesani, F. Development of a “First Principles” Water Potential with Flexible Monomers. II: Trimer Potential Energy Surface, Third Virial Coefficient, and Small Clusters. J. Chem. Theory Comput. 2014, 10, 1599–1607. (34) Medders, G. R.; Babin, V.; Paesani, F. Development of a “First-Principles” Water Potential with Flexible Monomers. III. Liquid Phase Properties. J. Chem. Theory Comput. 2014, 10, 2906–2910. (35) Hankins, D.; Moskowitz, J. W.; Stillinger, F. H. Water Molecule Interactions. J. Chem. Phys. 1970, 53, 4544–4554. (36) Rezac, J.; Hobza, P. Benchmark Calculations of Interaction Energies in Noncovalent Complexes and Their Applications. Chem. Rev. 2016, 116, 5038–5071. 30

ACS Paragon Plus Environment

Page 30 of 39

Page 31 of 39 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

(37) Cisneros, G. A.; Wikfeldt, K. T.; Ojam¨ae, L.; Lu, J.; Xu, Y.; Torabifard, H.; Bart´ok, A. P.; Cs´anyi, G.; Molinero, V.; Paesani, F. Modeling Molecular Interactions in Water: From Pairwise to Many-Body Potential Energy Functions. Chem. Rev. 2016, 116, 7501–7528. (38) Paesani, F. Getting the Right Answers for the Right Reasons: Toward Predictive Molecular Simulations of Water with Many-Body Potential Energy Functions. Acc. Chem. Res. 2016, 49, 1844–1851. (39) Bajaj, P.; G¨otz, A. W.; Paesani, F. Toward Chemical Accuracy in the Description of Ion-Water Interactions Through Many-Body Representations. I. Halide-Water Dimer Potential Energy Surfaces. J. Chem. Theory Comput. 2016, 12, 2698–2705. (40) Bajaj, P.; Wang, X.-G.; Jr., T. C.; Paesani, F. Vibrational Spectra of Halide-Water Dimers: Insights on Ion Hydration from Full-Dimensional Quantum Calculations on Many-Body Potential Energy Surfaces. J. Chem. Phys. 2018, 148, 102321. (41) Riera, M.; G¨otz, A. W.; Paesani, F. The i-TTM model for Ab Initio-Based Ion–Water Interaction Potentials. II. Alkali Metal Ion–Water Potential Energy Functions. Phys. Chem. Chem. Phys. 2016, 18, 30334–30343. (42) Riera, M.; Mardirossian, N.; Bajaj, P.; G¨otz, A. W.; Paesani, F. Toward Chemical Accuracy in the Description of Ion–Water Interactions Through Many-Body Representations. Alkali-Water Dimer Potential Energy Surfaces. J. Chem. Phys. 2017, 147, 161715. (43) Riera, M.; Brown, S. E.; Paesani, F. Isomeric Equilibria, Nuclear Quantum Effects, and Vibrational Spectra of M+ (H2 O)n=1−3 Clusters, with M = Li, Na, K, Rb, and Cs, Through Many-Body Representations. J. Phys. Chem. A 2018, 122, 5811–5821. (44) Zhuang, D.; Riera, M.; Schenter, G. K.; Fulton, J. L.; Paesani, F. Many-Body Effects

31

ACS Paragon Plus Environment

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

Determine the Local Hydration Structure of Cs+ in Solution. J. Phys. Chem. Lett. 2019, 10, 406–412. (45) Braams, B. J.; Bowman, J. M. Permutationally Invariant Potential Energy Surfaces in High Dimensionality. Int. Rev. Phys. Chem. 2009, 28, 577–606. (46) Mao, Y.; Ge, Q.; Horn, P. R.; Head-Gordon, M. On the Computational Characterization of Charge-Transfer Effects in Noncovalently Bound Molecular Complexes. J. Chem. Theory Comput. 2018, 14, 2401–2417. (47) Phipps, M. J. S.; Fox, T.; Tautermann, C. S.; Skylaris, C.-K. Energy Decomposition Analysis Approaches and Their Evaluation on Prototypical Protein–Drug Interaction Patterns. Chem. Soc. Rev. 2015, 44, 3177–3211. (48) Horn, P. R.; Mao, Y.; Head-Gordon, M. Probing Non-Covalent Interactions with a Second Generation Energy Decomposition Analysis Using Absolutely Localized Molecular Orbitals. Phys. Chem. Chem. Phys. 2016, 18, 23067–23079. (49) Horn, P. R.; Mao, Y.; Head-Gordon, M. Defining the Contributions of Permanent Electrostatics, Pauli Repulsion, and Dispersion in Density Functional Theory Calculations of Intermolecular Interaction Energies. J. Chem. Phys. 2016, 144, 114107. (50) Horn, P. R.; Head-Gordon, M. Alternative Definitions of the Frozen Energy in Energy Decomposition Analysis of Density Functional Theory Calculations. J. Chem. Phys. 2016, 144, 084118. (51) Horn, P. R.; Head-Gordon, M. Polarization Contributions to Intermolecular Interactions Revisited with Fragment Electric-Field Response Functions. J. Chem. Phys. 2015, 143, 114111. (52) Perdew, J. P.; Schmidt, K. Jacob’s Ladder of Density Functional Approximations for the Exchange-Correlation Energy. AIP Conf. Proc. 2001; pp 1–20. 32

ACS Paragon Plus Environment

Page 32 of 39

Page 33 of 39 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

(53) Tao, J.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Climbing the Density Functional Ladder: Nonempirical Meta–Generalized Gradient Approximation Designed for Molecules and Solids. Phys. Rev. Lett. 2003, 91, 146401. (54) Feyereisen, M.; Fitzgerald, G.; Komornicki, A. Use of Approximate Integrals in Ab Initio Theory. An Application in MP2 Energy Calculations. Chem. Phys. Lett. 1993, 208, 359–363. (55) Weigend, F.; H¨aser, M. RI-MP2: First Derivatives and Global Consistency. Theor. Chem. Acc. 1997, 97, 331–340. (56) Dunning, T. H. Gaussian Basis Sets for Use in Correlated Molecular Calculations. I. The Atoms Boron Through Neon and Hydrogen. J. Chem. Phys. 1989, 90, 1007–1023. (57) Peterson, K. A.; Figgen, D.; Goll, E.; Stoll, H.; Dolg, M. Systematically Convergent Basis Sets with Relativistic Pseudopotentials. II. Small-Core Pseudopotentials and Correlation Consistent Basis Sets for the Post-d Group 1618 Elements. J. Chem. Phys. 2003, 119, 11113–11123. (58) Adler, T. B.; Knizia, G.; Werner, H. J. A Simple and Efficient CCSD(T)-F12 Approximation. J. Chem. Phys. 2007, 127, 221106. (59) Knizia, G.; Adler, T. B.; Werner, H. J. Simplified CCSD(T)-F12 Methods: Theory and Benchmarks. J. Chem. Phys. 2009, 130, 054104. (60) Hill, J. G.; Peterson, K. A.; Knizia, G.; Werner, H.-J. Extrapolating MP2 and CCSD Explicitly Correlated Correlation Energies to the Complete Basis Set Limit with First and Second Row Correlation Consistent Basis Sets. J. Chem. Phys. 2009, 131, 194105. (61) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Sch¨ utz, M.; Celani, P.; Gy¨orffy, W.; Kats, D.; Korona, T.; Lindh, R.; Mitrushenkov, A.; Rauhut, G.; Shamasundar, K. R.; Adler, T. B.; Amos, R. D.; Bernhardsson, A.; Berning, A.; 33

ACS Paragon Plus Environment

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

Cooper, D. L.; Deegan, M. J.; Dobbyn, A. J.; Eckert, F.; Goll, E.; Hampel, C.; Hesselmann, A.; Hetzer, G.; Hrenar, T.; Jansen, G.; K¨oppl, C.; Liu, Y.; Lloyd, A. W.; Mata, R. A.; May, A. J.; McNicholas, S. J.; Meyer, W.; Mura, M. E.; Nicklass, A.; O’Neill, D. P.; Palmieri, P.; Peng, D.; Pfl¨ uger, K.; Pitzer, R.; Reiher, M.; Shiozaki, T.; Stoll, H.; Stone, A. J.; Tarroni, R.; Thorsteinsson, T.; Wang, M. MOLPRO, a Package of Ab Initio Programs, version 2015.1 ed.; 2015; see http://www.molpro.net. (62) Ren, P.; Ponder, J. W. Polarizable Atomic Multipole Water Model for Molecular Mechanics Simulation. J. Phys. Chem. B 2003, 107, 5933–5947. (63) Grossfield, A.; Ren, P.; Ponder, J. W. Ion Solvation Thermodynamics from Simulation with a Polarizable Force Field. J. Am. Chem. Soc. 2003, 125, 15671–15682. (64) Ren, P.; Ponder, J. W. Temperature and Pressure Dependence of the AMOEBA Water Model. J. Phys. Chem. B 2004, 108, 13427–13437. (65) Piquemal, J.-P.; Perera, L.; Cisneros, G. A.; Ren, P.; Pedersen, L. G.; Darden, T. A. Towards Accurate Solvation Dynamics of Divalent Cations in Water Using the Polarizable Amoeba Force Field: From Energetics to Structure. J. Chem. Phys. 2006, 125, 054511. (66) Wu, J. C.; Piquemal, J.-P.; Chaudret, R.; Reinhardt, P.; Ren, P. Polarizable Molecular Dynamics Simulation of Zn (II) in Water Using the AMOEBA Force Field. J. Chem. Theory Comput. 2010, 6, 2059–2070. (67) Wu, J. C.; Chattree, G.; Ren, P. Automation of AMOEBA Polarizable Force Field Parameterization for Small Molecules. Theor. Chem. Acc. 2012, 131, 1138. (68) Marjolin, A.; Gourlaouen, C.; Clavagu´era, C.; Ren, P. Y.; Wu, J. C.; Gresh, N.; Dognon, J.-P.; Piquemal, J.-P. Toward Accurate Solvation Dynamics of Lanthanides and Actinides in Water Using Polarizable Force Fields: From Gas-Phase Energetics to Hydration Free Energies. Theor. Chem. Acc. 2012, 131, 1198. 34

ACS Paragon Plus Environment

Page 34 of 39

Page 35 of 39 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

(69) Xiang, J. Y.; Ponder, J. W. A Valence Bond Model for Aqueous Cu (II) and Zn (II) Ions in the AMOEBA Polarizable Force Field. J. Comput. Chem. 2013, 34, 739–749. (70) Wang, L.-P.; Head-Gordon, T.; Ponder, J. W.; Ren, P.; Chodera, J. D.; Eastman, P. K.; Martinez, T. J.; Pande, V. S. Systematic Improvement of a Classical Molecular Model of Water. J. Phys. Chem. B 2013, 117, 9956–9972. (71) Shi, Y.; Xia, Z.; Zhang, J.; Best, R.; Wu, C.; Ponder, J. W.; Ren, P. Polarizable Atomic Multipole-Based AMOEBA Force Field for Proteins. J. Chem. Theory Comput. 2013, 9, 4046–4063. (72) Thole, B. T. Molecular Polarizabilities Calculated with a Modified Dipole Interaction. Chem. Phys. 1981, 59, 341–350. (73) Halgren, T. A. The Representation of van der Waals (vdW) Interactions in Molecular Mechanics Force Fields: Potential Form, Combination Rules, and vdW Parameters. J. Am. Chem. Soc. 1992, 114, 7827–7843. (74) Arismendi-Arrieta, D. J.; Riera, M.; Bajaj, P.; Prosmiti, R.; Paesani, F. i-TTM Model for Ab Initio-Based Ion–Water Interaction Potentials. 1. Halide–Water Potential Energy Functions. J. Phys. Chem. B 2015, 120, 1822–1832. (75) Burnham, C. J.; Anick, D. J.; Mankoo, P. K.; Reiter, G. F. The Vibrational Proton Potential in Bulk Liquid Water and Ice. J. Chem. Phys. 2008, 128, 154519. (76) Partridge, H.; Schwenke, D. W. The Determination of an Accurate Isotope Dependent Potential Energy Surface for Water from Extensive Ab Initio Calculations and Experimental Data. J. Chem. Phys. 1997, 106, 4618–4639. (77) Woon, D. E.; Dunning, T. H. Gaussian Basis Sets for Use in Correlated Molecular Calculations. IV. Calculation of Static Electrical Response Properties. J. Chem. Phys. 1994, 100, 2975–2988. 35

ACS Paragon Plus Environment

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

(78) Tang, K. T.; Toennies, J. P. An Improved Simple Model for the van der Waals Potential Based on Universal Damping Functions for the Dispersion Coefficients. J. Chem. Phys. 1984, 80, 3726–3741. (79) Becke, A. D.; Johnson, E. R. A Density-Functional Model of the Dispersion Interaction. J. Chem. Phys. 2005, 123, 154101. (80) Kannemann, F. O.; Becke, A. D. Van der Waals Interactions in Density-Functional Theory: Intermolecular Complexes. J. Chem. Theory Comput. 2010, 6, 1081–1088. (81) Otero-De-La-Roza, A.; Johnson, E. R. Non-Covalent Interactions and Thermochemistry Using XDM-Corrected Hybrid and Range-Separated Hybrid Density Functionals. J. Chem. Phys. 2013, 138, 204109. (82) Tikhonov, A. N. Solution of Incorrectly Formulated Problems and the Regularization Method. Soviet Math. Dokl. 1963; pp 1035–1038. (83) Khaliullin, R. Z.; Cobar, E. A.; Lochan, R. C.; Bell, A. T.; Head-Gordon, M. Unravelling the origin of intermolecular interactions using absolutely localized molecular orbitals. The Journal of Physical Chemistry A 2007, 111, 8753–8765. (84) Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A 1988, 38, 3098. (85) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B 1988, 37, 785. (86) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865. (87) Zhang, Y.; Yang, W. Comment on “Generalized Gradient Approximation Made Simple”. Phys. Rev. Lett. 1998, 80, 890.

36

ACS Paragon Plus Environment

Page 36 of 39

Page 37 of 39 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

(88) Sun, J.; Ruzsinszky, A.; Perdew, J. P. Strongly Constrained and Appropriately Normed Semilocal Density Functional. Phys. Rev. Lett. 2015, 115, 036402. (89) Mardirossian, N.; Ruiz Pestana, L.; Womack, J. C.; Skylaris, C.-K.; Head-Gordon, T.; Head-Gordon, M. Use of the rVV10 Nonlocal Correlation Functional in the B97M-V Density Functional: Defining B97M-rV and Related Functionals. J. Phys. Chem. Lett. 2016, 8, 35–40. (90) Becke, A. D. Density-Functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648–5652. (91) Adamo, C.; Barone, V. Toward Reliable Density Functional Methods Without Adjustable Parameters: The PBE0 Model. J. Chem. Phys. 1999, 110, 6158–6170. (92) Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-Class Functionals and 12 Other Functionals. Theor. Chem. Acc. 2008, 120, 215–241. (93) Chai, J.-D.; Head-Gordon, M. Systematic Optimization of Long-Range Corrected Hybrid Density Functionals. J. Chem. Phys. 2008, 128, 084106. (94) Chai, J.-D.; Head-Gordon, M. Long-Range Corrected Hybrid Density Functionals with Damped Atom–Atom Dispersion Corrections. Phys. Chem. Chem. Phys. 2008, 10, 6615–6620. (95) Mardirossian, N.; Head-Gordon, M. Mapping the Genome of Meta-Generalized Gradient Approximation Density Functionals: The Search for B97M-V. J. Chem. Phys. 2015, 142, 074111. (96) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Ini-

37

ACS Paragon Plus Environment

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

tio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (97) Perdew, J. P.; Burke, K.; Ernzerhof, M. Perdew, Burke, and Ernzerhof Reply. Phys. Rev. Lett. 1998, 80, 891. (98) Hait, D.; Head-Gordon, M. Delocalization Errors in Density Functional Theory Are Essentially Quadratic in Fractional Occupation Number. J. Phys. Chem. Lett. 2018, 9, 6280–6288. (99) Piquemal, J.-P.; Marquez, A.; Parisel, O.; Giessner-Prettre, C. A CSOV study of the difference between HF and DFT intermolecular interaction energy values: The importance of the charge transfer contribution. J. Comput. Chem. 2005, 26, 1052– 1062. (100) Wu, Q.; Ayers, P. W.; Zhang, Y. Density-based energy decomposition analysis for intermolecular interactions with variationally determined intermediate state energies. J. Chem. Phys. 2009, 131, 164112. (101) Wang, F.-F.; Jenness, G.; Al-Saidi, W.; Jordan, K. Assessment of the performance of common density functional methods for describing the interaction energies of (H 2 O) 6 clusters. J. Chem. Phys. 2010, 132, 134303. (102) Stone, A. The Theory of Intermolecular Forces; OUP Oxford, 2013. (103) Towns, J.; Cockerill, T.; Dahan, M.; Foster, I.; Gaither, K.; Grimshaw, A.; Hazlewood, V.; Lathrop, S.; Lifka, D.; Peterson, G. D.; Roskies, R.; Scott, J. R.; WilkinsDiehr, N. XSEDE: Accelerating Scientific Discovery. Comput. Sci. Eng. 2014, 16, 62–74.

38

ACS Paragon Plus Environment

Page 38 of 39

Page 39 of 39 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

TOC Figure

39

ACS Paragon Plus Environment