CationâÏ Interactions: Accurate Intermolecular Potential from

Subscriber access provided by Otterbein University

Article

Cation-# Interactions: Accurate Intermolecular Potential from Symmetry-Adapted Perturbation Theory Kay Ansorg, Maxim Tafipolsky, and Bernd Engels J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/jp403578r • Publication Date (Web): 07 Aug 2013 Downloaded from http://pubs.acs.org on August 13, 2013

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

The Journal of Physical Chemistry B 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 23

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

The Journal of Physical Chemistry

Cation-π Interactions: Accurate Intermolecular Potential from Symmetry-Adapted Perturbation Theory Kay Ansorg, Maxim Tafipolsky*, and Bernd Engels* July 26, 2013 Institut für Physikalische und Theoretische Chemie, Universität Würzburg, Campus Hubland Nord, Emil-Fischer-Strasse 42, D-97074 Würzburg, Germany

Principal Corresponding Author * Prof. Dr. Bernd Engels, 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 9313185394, E-mail: [email protected] 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] The authors declare no competing financial interest.

Abstract Symmetry-adapted perturbation theory (SAPT) is used to decompose the total intermolecular interaction energy between the ammonium cation and a benzene molecule into four physically-motivated individual contributions: electrostatics, exchange, dispersion and induction. Based on this rigorous decomposition it is shown unambiguously that both the electrostatic and the induction energy components contribute almost equally to the attractive forces stabilizing the dimer with a nonnegligible contribution coming from the dispersion term. A polarizable potential model for the interaction of ammonium cation with benzene is parametrized by fitting these four energy components separately using the functional forms of the AMOEBA force field augmented with the missing charge penetration energy term calculated as a sum over pairwise electrostatic energies between spherical atoms. It is shown that the proposed model is able to produce accurate intermolecular interaction energies as compared to ab initio results thus avoiding error compensation to a large extent.

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

Page 2 of 23

Keywords: Force field, AMOEBA, parameterization, genetic algorithm, charge penetration

1

Introduction

Understanding the forces acting between atoms and molecules is crucial for elucidation the organization and as a consequence the physicochemical properties of matter being it in gas phase or in condensed phase (in solids and in solutions). As has been stressed in a recent book of Stone, 1 all of the important intermolecular forces “arise ultimately from the electrostatic interaction between the particles comprising the two molecules”. Importantly, at short intermolecular distances, the electron exchange, a purely quantum mechanical effect, should be taken into account, that is, the total wavefunction should be properly antisymmetrised. This particular energy contribution arises from the fermionic character of the electrons. At the same time, the tendency to give intermolecular forces different names such as hydrogen bonds (conventional or nonconventional), salt bridge, hydrophobic (and solvophobic in general), ion−π, XH−π (X=C, N, O) or π−π interactions, etc. is dictated more by tradition (adopted in (bio)chemistry) than by strict physical criteria. Consider two interacting molecules (when two monomers form a dimer) and apply the theory of intermolecular interactions. From a quantum-mechanical perspective, one can use either a supermolecular approach, where the interaction energy is calculated simply as a difference between the dimer energy and the sum of energies of two monomers, or a perturbation approach, where the interaction energy is given as a sum of the first-, second- and higher-order terms of perturbation. Recent developments within the Symmetry-Adapted intermolecular Perturbation Theory (SAPT) 2 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. 3–8 The efficient implementations 4,7,9 of this theory allow 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 as reference data. We should also note, that it is more appropriate to discuss the forces underlying the intermolecular interactions using physically sound concepts developed, for example, within the SAPT framework, which are at the same time already quite familiar to most chemists and are widely used in many popular force fields, namely electrostatics, exchange, induction and dispersion. Of course, we should understand that only the total energy is quantum-mechanical observable and all partitions into different components are ad hoc. Most force fields represent the intermolecular interaction energy as a sum of various components such as electrostatics, exchange-repulsion, dispersion and induction (if polarizability is taken into account explicitly or implicitly). These force fields are usually parametrized using experimental (thermodynamic) data and sometimes including ab initio calculations as reference. An impressive number of attempts exist in the literature aimed at parametrizing the intermolecular potentials based on ab initio data including intermolecular perturbation theory. 10–38 The use of the intermolecular perturbation theory would be an ideal basis to parametrize a force field utilizing the ab initio data exclusively. We should stress, however, that a direct comparison is hampered by the fact that the meaning of various components obtained from the


Page 3 of 23

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


perturbation theory is different from that used in force fields. What is usually called the electrostatic energy within force fields (described by point charges, bond dipoles or even distributed multipoles) actually represents its long-range part only. The full electrostatics must include the charge penetration effects as well. These short-range interactions contribute significantly when two molecules are so close to each other that their charge distributions overlap. 22,28,39–43 In force fields, this particular short-range contribution is implicitly absorbed in the so-called van der Waals (vdW) term which also includes exchange-repulsion and dispersion parts of the total intermolecular energy (see, for example, 44,45 ). Interestingly, one of the earliest attempts to parameterize the main terms in a force field separately (including charge penetration) based on the data from the intermolecular perturbation theory was made almost 30 years ago by Sokalski et al. 46 The aim of this work is twofold. First, the symmetry-adapted perturbation theory is used to shed more light on a particular noncovalent interaction called “cation−π interaction”, where a cationic species interacts with an aromatic system. This attractive interaction contributes to the stabilization of biomolecules and plays an important role in molecular recognition processes (for some recent reviews on this subject, see 47–50 ). Second, employing the ammonium cation−benzene dimer as a model system, we develop a polarizable potential capable of reproducing the results from SAPT calculations. As in our previous work, 28 the main purpose here is to be able to describe the repulsion-attraction forces in a broad range of intermolecular distances and orientations in a balanced way thus avoiding error cancellation to a large extent. As before, this will be achieved by a separate fit of four individual contributions of the interaction energy using the functional forms of the AMOEBA 51 force field. We then compare our newly developed polarizable potential with the original AMOEBA model and also with highly accurate calculations. We note, that the AMOEBA force field is a polarizable one capable of including non-pairwise-additive many body effects in the force field via induced atomic dipoles. It has been shown previously that both the electrostatic term and the induction (polarization) term are mainly responsible for large attractive forces between a cation and an aromatic system. 52–56 The advantages of using a polarizable potential to study cation−π interactions have been recently demonstrated by Orabi and Lamoureux 57 and by Zheng et al. 58 These findings motivated us to use a polarizable potential for our parametrization.

2

Computational Details

The reference data for parametrization were calculated in the framework of SAPT as a sum of individual first (electrostatic) and second (induction and dispersion) order interaction terms accompanied by their respective exchange counterparts. The intramolecular correlation effects of the individual SAPT terms are described by static and time-dependent density functional theory (DFT-SAPT). 4,9 In all our DFT-SAPT calculations the dimer-centered basis sets were used, i.e. all atomic basis functions of both monomers were employed in both monomer and dimer calculations. To speed up the calculations, we employed the density-fitting variant 59 of DFT-SAPT as implemented in the MOLPRO 60,61 program package using the augmented correlation consistent aug-cc-pVXZ (X=T, Q) basis sets of Dunning 62 (abbreviated as avtz and avqz, respectively) along with an appropriate combination of cc-pV(X+1)Z JK and aug-cc-pVXZ MP2-fitting basis



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 23

sets of Weigend et al. 63,64 To produce accurate results, the asymptotic correction to the exchange-correlation potential should be applied. We used the gradient-regulated asymptotic correction approach of Grüning et al., 65 to correct the PBE0 66,67 hybrid exchange and PW91 correlation functionals by employing the calculated ionization potentials of 0.3423 (exp: 0.3397 68 ) and 0.9764 Hartree for benzene and the ammonium cation, respectively, and the HOMO energy of −0.2679 (benzene) and −0.8365 (ammonium cation) Hartree obtained using PBE0 26 functional. The corresponding shift parameters were approximated by the difference between the HOMO energy and the (negative) ionization potential of each monomer. As no higher than second-order terms are currently implemented in SAPT, third- and higher-order induction and exchange-induction contributions were estimated using the supermolecular approach at the Hartree-Fock level 69 (counterpoise corrected for basis set superposition error) and added to the sum of first- and second-order DFT-SAPT energy contributions. Accurate evaluation of the dispersion (and its exchange counterpart) contribution is the most timeconsuming and difficult part of the SAPT calculation (see, e.g. 70,71 and references therein). This particular term usually requires the use of extended basis sets. 72 Therefore, we utilized a larger augmented correlation consistent aug-cc-pVQZ (for C and N) and cc-pVQZ (for H) basis sets to check the convergence of this particular second-order energy component. We found that for the present purposes, the use of a smaller aug-cc-pVTZ basis set is justified, which only insignificantly underestimates the dispersion contribution. Moreover, ignoring the diffuse functions on hydrogens (i.e. using the cc-pVTZ basis set) only insignificantly underestimates the dispersion contribution for this system as well. The total interaction energy calculated within the SAPT framework, truncated at second order, is usually written as the sum of individual components (the superscripts 1 or 2 in parentheses refer to the order of the perturbation correction):

(1)

(1)

(2)

(2)

(2)

(2)

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

(1)

The δ(HF) term (mainly due to higher-order induction and exchange-induction) is calculated as the difference between the counterpoise corrected Hartree-Fock (HF) total dimer interaction energy and the sum of the uncorrelated SAPT contributions through second order in the intermolecular interaction operator:

(1)

(1)

(2)

(2)

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

(2)

We then fit our analytical model potential based on the AMOEBA functional forms to the corresponding terms of the total DFT-SAPT interaction energies separately as done in our previous work. 28 As a counterpart of the vdW term, we combine three terms (exchange, dispersion and exchange-dispersion) together:

(1)

(2)

(2)

SAPT EvdW = Eexch + Edisp + Eexch−disp

(3)

whereas the sum of induction, exchange-induction and higher-order (δ(HF)) contributions is used to represent the induction energy to be compared with the induction energy from a polarizable model:


Page 5 of 23

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)

(2)

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

(4)

We remark here, that from a symmetry-adapted perturbation theory perspective, the usual pairwise additive vdW term (based on Lennard-Jones, Born-Mayer, Buckingham or Morse potentials) used in force fields mixes up the short(1)

(2)

(2)

range first-order exchange (Eexch ) and the long-range second-order dispersion and exchange-dispersion (Edisp + Eexch−disp ) contributions (Eq. 3). 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, 73 and utilizing the aug-ccpVTZ basis set for the C and N atoms and the cc-pVTZ basis set for the H atoms 62 (abbreviated as “avtz*”) using the Gaussian program package. 74 For all dimer calculations, the monomers were kept rigid at the geometries optimized at the SCS-RI-MP2/TZVPP level of theory. In AMOEBA, 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 (GDMA) program 75 (for more details, see 51 ). As in our previous work, the electrostatic contribution due to charge penetration is included explicitly 41 in order to account for short-range quantum effects that are not accounted for by the classical multipolar expansion valid only at long range. We estimate 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, see 28 ). The expansion-contraction parameters for the aromatic carbon (1.0) and hydrogen (1.4) atoms were taken from our previous work, whereas the values for the N (2.5) and H (1.5) atoms of the ammonium cation were optimized to (1)

reproduce the first-order electrostatic term (Eelst ) for a number of cation−benzene dimer configurations. To describe the pairwise additive vdW interactions, we have adopted the buffered 14-7 functional form 76 used in the AMOEBA force field (for details, see 51 ). It is, however, necessary to adjust the vdW parameters appropriately, 41 since the electrostatic part of the force field is modified by adding the short-range interactions as described above. We therefore reparametrize the values of the atom size (in Å) and homoatomic well depth (in kcal/mol), whereas a reduction factor for hydrogen atoms was kept unchanged. All parameter refinements utilized the same genetic algorithm PIKAIA 77 as done before 28 with the root-mean-square deviation (rmsd) between the SAPT reference data and the corresponding part of the potential used as a fitness function. Induction is the interaction of an induced dipole on one fragment with the permanent dipole on another fragment, expressed in terms of the dipole polarizability. In AMOEBA, the molecular polarizability is expressed as a sum over atomic isotropic dipole polarizabilities. Iterating the dipole−induced dipole interaction to self-consistency captures many body effects. Isotropic polarizabilities for carbon (1.750 Å3 ) and hydrogen (0.696 Å3 ) in benzene were adopted from the AMOEBA force field, 78 whereas the values of 1.0 Å3 for nitrogen (slightly less than the recommended value of 1.073 51 ) and 0.15 Å3 (recommended value: 0.496 Å3 , see the amoeba09.prm parameter file) for hydrogen in the ammonium cation were employed. It was deemed neccesary to use approximately three times smaller value for the hydrogen isotropic polarizability



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 23

Figure 1: Monodentate (left), bidentate (middle) and tridentate (right) configurations of the ammonium cation−benzene dimer. 0 −2 −4 −6 −8 −10 −12 −14 −16 −18

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

CCSD(T) DFT−SAPT

2.5

3

3.5

4

4.5

5

5.5

6

CCSD(T) DFT−SAPT

2.5

3

R (Å)

3.5

4

4.5

5

5.5

R (Å)

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

CCSD(T) DFT−SAPT

2.5

3

3.5

4

4.5

5

5.5

R (Å)

Figure 2: Total interaction energy (in kcal/mol) for monodentate (left), bidentate (middle) and tridentate (right) configurations of the ammonium cation−benzene dimer as a function of the distance between the nitrogen atom and the center of the benzene ring (in Angstrom). in order to reproduce the calculated static dipole−dipole molecular polarizability tensor for the ammonium cation. The default value of 0.39 for the damping factor was adopted.

3

Results and Discussion

Total reference interaction energy. In Figure 2 we compare the DFT-SAPT total interaction energies and the CCSD(T) supermolecular interaction energies for monodentate (one hydrogen atom points to the center of the benzene ring), bidentate (two hydrogens point to the benzene ring), and tridentate (three hydrogens point to the benzene ring) relative orientations of the cation and the benzene ring as a function of the intermolecular distance (see Figure 1). As can be judged from this Figure 2 (numerical values are given as Supplementary Data), the agreement between two methods (CCSD(T) and DFT-SAPT) is good, which makes us confident of using the DFT-SAPT energies as the reference data for the parametrization. We note, however, that the DFT-SAPT energies start to deviate from the CCSD(T) results at shorter distances. Excluding the diffuse functions on hydrogens in the DFT-SAPT calculations (i.e. using the cc-pVTZ insted of the aug-cc-pVTZ basis set as is done in our CCSD(T) calculations) decreases the deviation only marginally (less than 2%). Our total interaction energy for the most attractive monodentate orientation of the ammonium cation (at the bottom of the curve shown in Fig. 2left) calculated at the CCSD(T)/avtz* level of theory (−17.5 kcal/mol) is larger than that


Page 7 of 23

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


published recently by Sherrill et al. 56 , (−16.4 kcal/mol) based on the CCSD(T)/6-311++G(2d,2p) level theory. This is probably due to the different basis sets used in both calculations (aug-cc-pVTZ(C,N)/cc-pVTZ(H) vs. 6-311++G(2d,2p)). Our value corresponds to the intermolecular distance (the distance between the nitrogen atom and the centroid of the benzene ring) of around 3.0 Å, whereas they give the value of 3.1 Å. We also confirm that the interaction energy is almost independent of the orientation of the cation relative to the benzene ring, but slight preference exists for both the monodentate and bidentate orientations (compare three potential energy curves in Figure 2). Our interaction energy at the intermolecular distance of around 2.9 Å for the bidentate orientation (−17.8 kcal/mol) is somewhat more negative than the value of −17.4 kcal/mol calculated by Sa et al. 79 at the MP2/6-31G(d,p) level of theory and the value of −17.6 kcal/mol obtained recently by Orabi and Lamoureux 57 at the MP2/6-311++G(d,p) level of theory. All of them are significantly more negative than the value of −16.0 kcal/mol calculated by Reddy and Sastry 80 at the same level of theory (MP2/6-311++G(d,p)). This discrepancy is probably due to differently optimized geometries used for the calculations of binding energies. Orabi and Lamoureux 57 have optimized the dimer geometry at the MP2/6-311++G(d,p) level of theory and obtained the value of 2.92 Å for the intermolecular distance, whereas Reddy and Sastry 80 used the preoptimized geometry at the B3LYP/6-31G(d,p) level and arrived at a somewhat (0.02 Å) larger intermolecular distance of 2.94 Å. Dehez et al. 81 have calculated the values of −16.1 kcal/mol (monodentate, at 3.1 Å) and −16.5 kcal/mol (bidentate, at 3.0 Å) at the MP2/6-311++G(d,p) level of theory (corrected for BSSE), which are in accord with their values obtained from SAPT(MP2)/6-311++G(d,p) calculations. To investigate the possible influence of the geometry relaxation on the intermolecular energy, we have fully optimized the ammonium-benzene dimer at the MP2/aug-cc-pVTZ level of theory followed by a single-point calculation at the CCSD(T)/avtz* level. The obtained stationary point was confirmed to be the true minimum by frequency calculations and is believed to be the global minimum on the MP2/aug-cc-pVTZ potential energy surface. Our value of −18.7 kcal/mol (CCSD(T)/avtz*, corrected for BSSE) is in accord with our DFT-SAPT value of −19.0 kcal/mol and with a previous study of Kim et al. 53 (−18.6 kcal/mol) for a fully optimized dimer geometry at the MP2/aug-cc-pVDZ level of theory (corrected for BSSE). We note, that the dimer geometry optimized by Kim et al. 53 is very similar to ours with the distance between the nitrogen atom and the center of the benzene ring of 2.88 Å (see Supplementary Data). Wang et al. 82 reported the value of −17.5 kcal/mol for the geometry (optimized at the PBE/6-31++G(d,p) level) very close to the bidentate orientation (at 2.96 Å) based on their SAPT(MP2)/6-31++G(d,p) calculations. We also mention here an alternative energy decomposition analysis by Mo et al. 83 based on a block-localized wavefunction and multistate density functional theory using HF/6-311G(d,p), where they have found the value of −18.6 kcal/mol for the total binding energy for the fully optimized ammonium-benzene dimer with the dispersion contribution evaluated at the MP2/6-311G(d,p) level of theory. Yet another energy decomposition scheme (reduced variational space self-consistent field method) used by Aschi et al. 84 gave the value of −17.9 kcal/mol for the bidentate dimer orientation (at 3.16 Å) geometry-optimized at the HF/6-311+G(d,p) level of theory. Surprisingly, Singh et al. 72 have estimated the binding energy of −21.4 kcal/mol for the bidentate dimer orientation (at 2.94 Å with a slight offset of 0.06 Å) at the complete basis set limit (CCSD(T)/CBS).



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 23

Table 1: DFT-SAPT energy components for the most attractive monodentate, bidentate, tridentate and fully optimized ammonium cation−benzene configurations (in kcal/mol) monodentate bidentate tridentate optimized (1)

Eelst (1) Eexch (2) Eind

(2)

Eexch−ind (2) Edisp (2)

Eexch−disp δ(HF) SAPT Etot

−11.32 10.52 −13.48 4.08 −5.78

−12.68 12.68 −14.74 5.25 −6.47

−12.01 10.99 −13.32 4.85 −5.98

−13.68 14.77 −16.40 5.91 −7.05

0.71 −2.34 −17.61

0.85 −2.81 −17.92

0.78 −1.57 −16.27

0.95 −3.51 −19.01

The origin of such a large discrepancy with our calculations as well as with the results from others is still unclear. Table 1 compares the decomposition of the DFT-SAPT total interaction energy into physically-motivated contributions (see Eq. 1) calculated at the most attractive intermolecular distance (at the bottom of each curve) for monodentate (3.0 Å), bidentate (2.9 Å) and tridentate (2.9 Å) configurations as well as for a fully optimized dimer geometry. The data presented in Table 1 unambiguously indicate that the most important attractive contributions to the interaction forces (1)

come from the electrostatic (Eelst ) and induction (Eq. 4) energy terms, which are almost equally crucial in stabilizing the dimer. We should also note that for the most attractive monodentate and bidentate cation−benzene configurations (in the region of minima of the potential energy curves, see Figure 2) approximately 75-85% of the electrostatic energy comes from the multipolar (long-range) part, whereas the remaining 15-25% originates from the short-range (charge penetration) contribution (see Figure 3). This short-range term even dominates for intermolecular distances smaller than 3 Å, where (2)

(2)

the multipole expansion is not valid anymore. The dispersion (Edisp + Eexch−disp ) contributes around 18 % to the overall (1)

(2)

(2)

(2)

(2)

attractive forces (Eelst + Eind + Eexch−ind + δ(HF) + Edisp + Eexch−disp ), which is not negligible. Training data set. Our training data set consists of homodimers only, namely, we have generated a number of representative configurations of the benzene dimer (155) and the ammonium cation dimer (54). These are used to refine the vdW parameters as well as the expansion-contraction parameters (called “kappa” in 28 ) to reproduce the exchange(1)

dispersion (see Eq. 3) and the electrostatic (Eelst ) terms, respectively, calculated with DFT-SAPT. Validation data set. We have generated a large number of different dimer configurations including both on- and off-axis face-on orientations as well as a number of side-on ones (overall around 350 configurations) and used them for validation purposes. The on-axis face-on orientation corresponds to the geometry where the nitrogen atom of the cation is placed on the line perpendicular to the benzene ring and passing through its center (as in mono- bi- and tridentate orientations above). Despite the fact that the side-on configurations are found to be energetically much less favorable then the face-on ones, 56 it was important to include them in order to take various geometries, which might result due to constraints arising, for example, from the protein environment, into account. Electrostatic energy. The long-range part of the electrostatic energy described by Coulomb interactions between atomic multipoles (up to and including quadrupoles) starts to deviate (underestimate) the reference (first-order) electro-


Page 9 of 23

0 −5 −10 −15

E(1)elst (DFT−SAPT) Multipoles (AMOEBA)

−20 −25

2.5

3

3.5

4 4.5 R (Å)

5

5.5

6

Figure 3: Long-range multipole electrostatics against the DFT-SAPT reference first-order electrostatic energy (in kcal/mol) for the monodentate orientation of the ammonium cation−benzene dimer as a function of the intermolecular distance (in Angstrom)

ESAPT − EFF (kcal/mol)

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


−5 −10 −15 −20 −25 −25 −20 −15 −10 −5 ESAPT (kcal/mol)

0

4 3 2 1 0 −1 −2 −3 −4 −5 −25 −20 −15 −10 −5 ESAPT (kcal/mol)

0

Figure 4: Scatter plots of electrostatic energies (left) and their differences (right) calculated for 350 ammonium cation−benzene dimer configurations (in kcal/mol). static contribution at around 3.5 Å as is evident from Figure 3, where the data for the monodentate orientation as a function of the intermolecular distance is shown. Whereas the refined kappa value for the nitrogen atom (0.93) is very close to 1.0, the optimized kappa value for the hydrogen atom in the ammonium cation (3.5) deviates significantly from 1.0 indicating large contraction of its electron density. This is meaningful taking into account that the cation bears an overall positive charge and the fact that the nitrogen atom is more electronegative than the hydrogen atom. In Figure 4 (1)

we compare the reference data (ESAPT = Eelst ) with the sum of the long-range (multipolar) and the optimized short-range contributions obtained from our force field (EFF ). As can be seen, the (abs.) deviation for the majority of configurations does not exceed 1.0 kcal/mol, which is quite encouraging. Induction energy. We (see Table 1) and others 56,57,81 have shown that the induction term contributes significantly to the stabilization of the ammonium cation−benzene dimer. As will be shown, the AMOEBA polarizable model is capable



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

−5

−5

−10

−10

−15

−15

Page 10 of 23

0 −5 −10

−20 −25

−20

DFT−SAPT Modified AMOEBA

−15


−25 2.5

3

3.5

4

4.5

5

5.5

6

2.5

R (Å)

3

3.5

4

4.5

5

5.5

R (Å)

−20


2.5

3

3.5

4

4.5

5

5.5

R (Å)

Figure 5: Induction energy (in kcal/mol) for monodentate (left), bidentate (middle) and tridentate (right) configurations of the ammonium cation−benzene dimer as a function of the distance between the nitrogen atom and the center of the benzene ring (in Angstrom). of reproducing the reference induction contribution (Eq. 4) with only one parameter adjusted to reproduce the static molecular polarizability tensor calculated in this work at the MP2/aug-cc-pVTZ level. Our calculated value of 1.35 Å3 for the ammonium cation is in excellent agreement with some earlier estimates by Voisin et al. 85 (1.36 Å3 ) and by Dehez et al. 81 (1.34 Å3 ). To reproduce our MP2/aug-cc-pVTZ value, it was deemed neccesary to use a significantly reduced value for the isotropic polarizability of the H(−N) atom (0.15 Å3 ). The use of the recommended value of 0.496 Å3 86 results in overestimation of the ammonium cation polarizability by 70 % (2.28 Å3 ). We found that this recommended value is more suitable for the hydrogen atom in neutral molecules such as methane. Our findings are in accord with the recent studies of Elking et al. 87 and Kim and Rhee, 88 where it has been found that a much smaller values for the isotropic polarizabilities of the nitrogen and hydrogen atoms in the ammonium cation should be used (despite the different damping factors employed) to reproduce the calculated cation polarizability tensor (at the B3LYP/cc-pVTZ or B3LYP/aug-cc-pVTZ levels). This actually makes sense since due to positive charge of the cation, the polarizability of the hydrogen atom should be much reduced as compared to the hydrogen atom in neutral species such as methane. We have also carried out additional calculations (at the MP2/aug-cc-pVTZ level) of the molecular polarizability tensors of methyl-substituted ammonium ion (up to tetramethylammonium) and found out that the recommended value (0.496 Å3 ) for the isotropic polarizability of hydrogens attached to the carbon atom is probably still too large to fairly reproduce the reference ab initio molecular polarizabilities (data not shown). Additional evidence for a much smaller value needed for hydrogens comes from a comparison between the DFT-SAPT reference values and their force field counterparts for all 300 cation-benzene configurations considered. Moreover, we have calculated the induction energy for a number of cation-cation dimers (33 configurations) and found that the agreement is significantly improved (rmsd was three times smaller) when the new value (0.15 Å3 ) was used. In Figure 5, we compare the DFT-SAPT reference induction energies with the values obtained from the AMOEBA polarizable model estimated by iterating the induced dipoles to self-consistency as implemented in TINKER using the modified isotropic polarizabilities of the nitrogen and hydrogen atoms in the ammonium cation (labeled as “Modified AMOEBA”). The total polarization energy (both intra- and intermolecular) estimated with the AMOEBA potential is used for comparison with the reference values calculated according to Eq. 4.


Page 11 of 23

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


As can be judged from a good agreement with the reference values, the AMOEBA polarizable potential is sufficiently accurate for reliable treating these cation−π interactions. We should stress, however, that the hybrid approach used in the present study based on second-order intermolecular perturbation theory (see Eqs. 1 and 4) to evaluate the induction contribution is only an approximation to the true induction component of the total energy, where the third and higher-order terms are estimated at the Hartree-Fock level, i.e. neglecting intramolecular electron correlation (for a more thorough discussion, see 89–91 and references therein). One further note is in order here. The DFT-SAPT implementation in the MOLPRO program we are using is based on the so-called S 2 (or single-exchange of electrons) approximation for the (2)

(2)

evaluation of the second-order Eexch−ind and Eexch−disp terms, whereas multiple electron exchange (S ∞ - approximation) (1)

is taken into account in the calculation of the first-order Eexch term both at the HF-SAPT and DFT-SAPT levels. We (1)

always use the Eexch (S ∞ ) term in the evaluation of both the δ(HF) correction for the higher-order induction/exchangeinduction terms (Eq. 4) and for the exchange plus dispersion component (Eq. 3). Thus, the discrepancy seen between the DFT-SAPT and CCSD(T) curves at shorter intermolecular distances (see Fig. 2) might be due to underestimation of the exchange-induction component in DFT-SAPT. We, therefore, anticipate that the attractive forces due to the induction contribution (see Table 1 and Fig. 5) are somewhat exaggerated by the present DFT-SAPT calculations. As Schäffer and Jansen 91 have recently shown, the underestimation introduced by the S 2 -approximation is more pronounced for the secondorder exchange-induction component than for the first-order exchange energy. They developed an approach to include multiple electron exchange for the second-order exchange terms as well. With their implementation (not yet available to us), it would be indeed very interesting and informative to see how large this underestimation of the second-order exchange (2)

(2)

(2)

(2)

terms is for a cation-benzene dimer (i.e. Eexch−ind (S ∞ ) vs. Eexch−ind (S 2 ) and Eexch−disp (S ∞ ) vs. Eexch−disp (S 2 )). In addition, the evaluation of the third-order terms at the DFT-SAPT level will definitely contribute to our understanding of the validity of the δ(HF) correction term used routinely at present. These issues should be clarified before implementing a more sophisticated model to describe the induction energy within a force field (see, for example, 32,33 ). van der Waals energy. In the majority of commonly used force fields, the so-called “repulsion” part of the vdW term absorbs implicitly the attractive short-range electrostatic contribution due to charge penetration. Since we explicitly include this contribution in the electrostatic part of our force field (for more details, see 28 ), the parameters of the vdW term have to be readjusted appropriately. As already mentioned above, we adopted the functional form for the vdW energy (repulsion-dispersion) used in the AMOEBA force field and reoptimized the two parameters for each atom class (atom size, r, and well depth, ε) based on the sum of the three separate contributions calculated using SAPT, namely exchange-repulsion, dispersion and exchange-dispersion (Eq. 3), denoted by ESAPT . Table 2 compares the old and the new values for the carbon, nitrogen and hydrogen atoms obtained by fitting the exchange-dispersion contributions calculated with DFT-SAPT for homodimers. We found that the refinement of the so-called reduction parameters brings about an additional improvement of the fit. These parameters are used to shift repulsion-dispersion hydrogen sites inward from the H atom to the neighbor heavy atom (C or N in our case). A much larger reduction of ca. 20% of the N−H bond length found here is probably due to the larger inward shift of the electron density of hydrogens towards the nitrogen atom.



Table 2: van der Waals parameters for the AMOEBA atom class AMOEBA r, Å ε, kcal/mol reduction N (“Ammonium Ion N+”) 3.71 0.1050 H (“Ammonium Ion H4N+”) 2.48 0.0115 0.90 C (“Benzene C”) 3.80 0.0890 H (“Benzene HC”) 2.98 0.0260 0.92


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

50 40 30 20 10 0

Page 12 of 23

and FF-SAPT force field. FF-SAPT r, Å ε, kcal/mol reduction 3.41(3) 0.0568(6) 3.18(3) 0.0143(1) 0.78(1) 4.23(4) 0.0412(4) 3.18(3) 0.0143(1) 0.91(1)

4 2 0 −2 −4 −6 −8

0 10 20 30 40 50 60 ESAPT (kcal/mol)

0 10 20 30 40 50 60 ESAPT (kcal/mol)

Figure 6: Scatter plots of van der Waals energies (left) and their differences (right) calculated for 350 ammonium cation−benzene dimer configurations (in kcal/mol). The values given in parentheses are uncertainties in the fitted parameters estimated as the change in the parameter that increases the rms error by 25%. 44 We should stress that the difference between the two sets of parameters is due to the fact that the original AMOEBA parameters absorb implicitly the attractive component of the interaction energy due to charge penetration effects, and, therefore, the corresponding term of the force field is better described as the repulsion-dispersion-penetration term. 44 We validate our new vdW parameters by reproducing the DFT-SAPT reference data on the benzene-cation dimers over a wide range of energy values. The agreement is fairly good as can be judged from Figure 6. The Figure also shows that with our new potential even the energies in a very “repulsive” region of the ammonium cation−benzene dimer configurations are accurately modeled. Finally, in Figure 7 is shown the comparison between the total CCSD(T) and DFT-SAPT interaction energy and its counterparts calculated with our newly developed potential (labeled as “FF-SAPT”) and with the original AMOEBA force field. Overall, our refined AMOEBA model reproduces more closely the reference data than the original AMOEBA intermolecular potential. We also note that our new potential is able to treat various configurations of the dimer in a more balanced way than the original AMOEBA potential does. This is particularly evident for bidentate orientations of the ammonium cation with respect to the benzene molecule. To judge the global performance of our new force field, in Figure 8 we compare the total intermolecular energies with those obtained with the original AMOEBA force field. It can be seen that the overall performance is improved compared to the reference values with the reduced rmsd. We have also tested the performance of our newly developed force field by reproducing the stationary points located


−15 −18

2.6

2.8

3

CCSD(T) DFT−SAPT FF−SAPT AMOEBA

0

3.2

−5 −10

−12 −5

0

Binding Energy (kcal/mol)

0


−12


5


−15 −18 2.6

2.8

3

3.2

−10

−15


−8 −12

−5

−16 2.6

2.8

3

3.2

−10

−15

−15 −20 2.5

3

3.5

4

4.5

5

5.5

6

−20 2.5

R (Å)

3

3.5

4

4.5

5

5.5

R (Å)

2.5

3

3.5

4

4.5

5

5.5

R (Å)

Figure 7: Total interaction energy (in kcal/mol) for monodentate (left), bidentate (middle) and tridentate (right) configurations of the ammonium cation−benzene dimer as a function of the distance between the nitrogen atom and the center of the benzene ring (in Angstrom). The inset shows the region near the bottom of the curves at a larger scale.


EFF (kcal/mol)

50 40 30 20 10

5 0 −5

−10

0

−10

−15

−20 −20−10 0 10 20 30 40 50 ESAPT (kcal/mol)

−20−10 0 10 20 30 40 50 ESAPT (kcal/mol)


50 EFF (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 13 of 23

40 30 20 10

5 0 −5

−10

0

−10

−15

−20 −20−10 0 10 20 30 40 50 ESAPT (kcal/mol)

−20−10 0 10 20 30 40 50 ESAPT (kcal/mol)

Figure 8: Scatter plots of total intermolecular energies and their differences (top: AMOEBA, bottom: FF-SAPT) calculated for 350 ammonium cation−benzene dimer configurations (in 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 14 of 23

using the MP2/aug-cc-pVTZ calculations as reference. Three stationary points corresponding to on-axis mono-, bi- and tridentate orientations, which are found to be saddle points on the MP2/aug-cc-pVTZ energy surface, were located with our new force field, whereas the original AMOEBA force field was able to locate only two stationary points corresponding to bi- and tridentate configurations, displaying a clear preference for tridentate configurations. Moreover, with the new force field, the optimized geometry of the dimer was very close to that found at the MP2/aug-cc-pVTZ level, whereas this was not the case with the original AMOEBA force field. We believe that this failure is not just due to the overestimation of the induction component of the total energy mentioned above (we checked this by modifying only the atomic polarizabilities in the original AMOEBA force field), but is a direct consequence of a more balanced treatment of various terms in the new force field. At the same time, our new force field underestimates somewhat the interaction energy for the on-axis monodentate configuration of the dimer (see Figure 7). To exploit the benefit of using physically-motivated contributions calculated within the SAPT framework even further, it would be advantageous to fit the exchange and dispersion parts separately. This will allow, in addition, the use of more sophisticated (and physically motivated) functional forms than those already implemented in the widely used force fields (including AMOEBA). 45 Taking into account the importance of this particular cation−π interactions for biochemistry, we suggest to include the ammonium-benzene dimer in the well-established database 92 of highly accurate intermolecular energies as an example of an induction dominated system for benchmarking.

4

Conclusions

In this work we investigate the use of ab initio high level quantum-mechanical method based on symmetry-adapted intermolecular perturbation theory (SAPT) to study the interactions in vacuum between an organic cation (modeled by the ammonium ion) and a prototypical aromatic system (modeled by a benzene molecule). The so-called cation−π interactions, which have long been considered a challenge for molecular modeling, can be modeled quite accurately and reliably by applying a physics-based potential. Motivated by the need for a force field that can capture the right physics of intermolecular interactions, we incorporated an extra pairwise-additive energy term, which describes the short-range contribution to the electrostatic energy due to the charge penetration, into the AMOEBA polarizable force field. Subsequently, we fitted the three physically-motivated contributions, namely electrostatic, exchange-repulsion and dispersion, to their SAPT counterparts. Combined with the reoptimized parameters for the van der Waals term, we obtained a force field, which avoids to a large extent the error cancellation between different terms, and hence delivers a considerably improved description of the ammonium−benzene interactions. Equipped with the advanced intermolecular perturbation theory, it is now very timely to rethink a number of ad hoc concepts such as “π−π” or “ion−π interactions” invoked to describe the forces that drive association between molecules (see, for example 93 and references therein).


Page 15 of 23

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


Acknowledgements We are grateful to the Deutsche Forschungsgemeinschaft (GRK1221, SFB 630) and the Volkswagen Stiftung for financial support.

Supporting Information Available: The Cartesian coordinates of the mono-, bi-, tridentate and a fully optimized configurations of the ammonium cation−benzene dimers together with the corresponding total imtermolecular interaction energies calculated by the CCSD(T) and DFTSAPT methods. The induction contributions calculated both within the DFT-SAPT framework and with the force field are also given. This information is available free of charge via the Internet at http://pubs.acs.org.

References [1] Stone, A. J. The Theory of Intermolecular Forces; Oxford University Press: Oxford, 2013. [2] Jeziorski, B.; Moszynski, R.; Szalewicz, K. Perturbation-Theory Approach to Intermolecular Potential-Energy Surfaces of Van-der-Waals Complexes. Chem. Rev. 1994, 94, 1887–1930. [3] Heßelmann, A.; Jansen, G. The Helium Dimer Potential from a Combined Density Functional Theory and SymmetryAdapted Perturbation Theory Approach Using an Exact Exchange-Correlation Potential. Phys. Chem. Chem. Phys. 2003, 5, 5010–5014. [4] Heßelmann, A.; Jansen, G.; Schutz, M. Density-Functional Theory-Symmetry-Adapted Intermolecular Perturbation Theory with Density Fitting: A New Efficient Method to Study Intermolecular Interaction Energies. J. Chem. Phys. 2005, 122, 014103. [5] Bukowski, R.; Podeszwa, R.; Szalewicz, K. Efficient Calculation of Coupled Kohn-Sham Dynamic Susceptibility Functions and Dispersion Energies with Density Fitting. Chem. Phys. Lett. 2005, 414, 111–116. [6] Podeszwa, R.; Bukowski, R.; Szalewicz, K. Density-Fitting Method in Symmetry-Adapted Perturbation Theory Based on Kohn-Sham Description of Monomers. J. Chem. Theory Comput. 2006, 2, 400–412. [7] Hohenstein, E. G.; Sherrill, C. D. Density Fitting of Intramonomer Correlation Effects in Symmetry-Adapted Perturbation Theory. J. Chem. Phys. 2010, 133, 014101. [8] Tekin, A.; Jansen, G. How Accurate Is the Density Functional Theory Combined With Symmetry-Adapted Perturbation Theory Approach for CH-π and π-π Interactions? A Comparison To Supermolecular Calculations for the Acetylene-Benzene Dimer. Phys. Chem. Chem. Phys. 2007, 9, 1680–1687.



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 23

[9] Misquitta, A. J.; Szalewicz, K. Symmetry-Adapted Perturbation-Theory Calculations of Intermolecular Forces Employing Density-Functional Description of Monomers. J. Chem. Phys. 2005, 122, 214109. [10] Mooij, W. T. M.; van Duijneveldt, F. B.; van Duijneveldt-van de Rijdt, J. G. C. M.; van Eijck, B. P. Transferable Ab Initio Intermolecular Potentials. 1. Derivation from Methanol Dimer and Trimer Calculations. J. Phys. Chem. A 1999, 103, 9872–9882. [11] Špirko, V.; Engvist, O.; Soldán, P.; Selzle, H. L.; Schlag, E. W.; Hobza, P. Structure and Vibrational Dynamics of the Benzene Dimer. J. Chem. Phys. 1999, 111, 572–582. [12] Mitchell, J. B. O.; Price, S. L. A Systematic Nonempirical Method of Deriving Model Intermolecular Potentials for Organic Molecules: Application To Amides. J. Phys. Chem. A 2000, 104, 10958–10971. [13] Hloucha, M.; Sum, A. K.; Sandler, S. I. Computer simulation of acetonitrile and methanol with ab initio-based pair potentials. J. Chem. Phys. 2000, 113, 5401–5406. [14] Engkvist, O.; Åstrand, P. O.; Karlström, G. Accurate Intermolecular Potentials Obtained from Molecular Wave Functions: Bridging the Gap between Quantum Chemistry and Molecular Simulations. Chem. Rev. 2000, 100, 4087– 4108. [15] Donchev, A. G.; Ozrin, V. D.; Subbotin, M. V.; Tarasov, O. V.; Tarasov, V. I. A Quantum Mechanical Polarizable Force Field for Biomolecular Interactions. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 7829–7834. [16] Torheyden, M.; Jansen, G. A New Potential Energy Surface for the Water Dimer Obtained from Separate Fits of Ab Initio Electrostatic, Induction, Dispersion and Exchange Energy Contributions. Mol. Phys. 2006, 104, 2101–2138. [17] Podeszwa, R.; Bukowski, R.; Szalewicz, K. Potential Energy Surface For the Benzene Dimer and Perturbational Analysis of π-π Interactions. J. Phys. Chem. A 2006, 110, 10345–10354. [18] Li, X.; Volkov, A. V.; Szalewicz, K.; Coppens, P. Interaction Energies between Glycopeptide Antibiotics and Substrates in Complexes Determined by X-ray Crystallography: Application of a Theoretical Databank of Aspherical Atoms and a Symmetry-Adapted Perturbation Theory-Based Set of Interatomic Potentials. Acta Cryst. D 2006, 62, 639–647. [19] Hellmann, R.; Bich, E.; Vogel, E. Ab Initio Intermolecular Potential Energy Surface and Second Pressure Virial Coefficients of Methane. J. Chem. Phys. 2008, 128, 214303. [20] Szalewicz, K.; Leforestier, C.; van der Avoird, A. Towards the Complete Understanding of Water by a First-Principles Computational Approach. Chem. Phys. Lett. 2009, 482, 1–14. [21] Archambault, F.; Chipot, C.; Soteras, I.; Javier Luque, F.; Schulten, K.; Dehez, F. Polarizable Intermolecular Potentials for Water and Benzene Interacting with Halide and Metal Ions. J. Chem. Theory Comput. 2009, 5, 3022–3031.


Page 17 of 23

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


[22] Wang, B.; Truhlar, D. G. Including Charge Penetration Effects in Molecular Modeling. J. Chem. Theory Comput. 2010, 6, 3330–3342. [23] Totton, T. S.; Misquitta, A. J.; Kraft, M. A First Principles Development of a General Anisotropic Potential for Polycyclic Aromatic Hydrocarbons. J. Chem. Theory Comput. 2010, 6, 683–695. [24] 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. [25] Boese, A. D.; Forbert, H.; Masia, M.; Tekin, A.; Marx, D.; Jansen, G. Constructing Simple Yet Accurate Potentials For Describing the Solvation of HCl/Water Clusters in Bulk Helium and Nanodroplets. Phys. Chem. Chem. Phys. 2011, 13, 14550–14564. [26] Leforestier, C.; Tekin, A.; Jansen, G.; Herman, M. First Principles Potential for the Acetylene Dimer and Refinement by Fitting to Experiments. J. Chem. Phys. 2011, 135, 234306. [27] Taylor, D. E.; Rob, F.; Rice, B. M.; Podeszwa, R.; Szalewicz, K. A Molecular Dynamics Study of 1,1-diamino-2,2dinitroethylene (FOX-7) Crystal Using a Symmetry Adapted Perturbation Theory-Based Intermolecular Force Field. Phys. Chem. Chem. Phys. 2011, 13, 16629–16636. [28] Tafipolsky, M.; Engels, B. Accurate Intermolecular Potentials with Physically Grounded Electrostatics. J. Chem. Theory Comput. 2011, 7, 1791–1803. [29] McDaniel, J. G.; Schmidt, J. R. Physically-Motivated Force Fields from Symmetry-Adapted Perturbation Theory. J. Phys. Chem. A 2013, 117, 2053–2066. [30] Yokogawa, D. Development of the Isotropic Site-Site Potential for Exchange Repulsion Energy and Combination with the Isotropic Site-Site Potential for Electrostatic Part. J. Chem. Phys. 2012, 137, 204101. [31] Wang, L.-P.; Chen, J.; Van Voorhis, T. Systematic Parametrization of Polarizable Force Fields from Quantum Chemistry Data. J. Chem. Theory Comput. 2013, 9, 452–460. [32] Misquitta, A. J.; Stone, A. J. Accurate Induction Energies for Small Organic Molecules: 1. Theory. J. Chem. Theory Comput. 2008, 4, 7–18. [33] Misquitta, A. J.; Stone, A. J.; Price, S. L. Accurate Induction Energies for Small Organic Molecules. 2. Development and Testing of Distributed Polarizability Models against SAPT(DFT) Energies. J. Chem. Theory Comput. 2008, 4, 19–32. [34] Holt, A.; Boström, J.; Karlström, G.; Lindh, R. A NEMO Potential that Includes the Dipole-Quadrupole and Quadrupole-Quadrupole Polarizability. 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 18 of 23

[35] Arnautova, Y. A.; Jagielska, A.; Pillardy, J.; Scheraga, H. A. Derivation of a New Force Field for Crystal-Structure Prediction Using Global Optimization: Nonbonded Potential Parameters for Hydrocarbons and Alcohols. J. Phys. Chem. B 2003, 107, 7143–7154. [36] Bordner, A. J.; Cavasotto, C. N.; Abagyan, R. A. Direct Derivation of van der Waals Force Field Parameters from Quantum Mechanical Interaction Energies. J. Phys. Chem. B 2003, 107, 9601–9609. [37] Hayes, J. M.; Greer, J. C.; Morton-Blake, D. A. A Force-Field Description of Short-Range Repulsions for High Density Alkane Molecular Dynamics, Simulations. J. Comput. Chem. 2004, 25, 1953–1966. [38] Tsuzuki, S.; Uchimaru, T.; Tanabe, K.; Kuwajima, S. Refinement of Nonbonding Interaction Potential Parameters for Methane on the Basis of the Pair Potential Obtained by MP3/6-311G(3d,3p)-Level Ab-Initio Molecular-Orbital Calculations - the Anisotropy of H/H Interaction. J. Phys. Chem. 1994, 98, 1830–1833. [39] Wheatley, R. J.; Mitchell, J. B. O. Gaussian Multipoles in Practice - Electrostatic Energies for Intermolecular Potentials. J. Comput. Chem. 1994, 15, 1187–1198. [40] Kairys, V.; Jensen, J. H. Evaluation of the Charge Penetration Energy between Non-orthogonal Molecular Orbitals Using the Spherical Gaussian Overlap Approximation. Chem. Phys. Lett. 1999, 315, 140–144. [41] Spackman, M. A. The Use of the Promolecular Charge Density to Approximate the Penetration Contribution to Intermolecular Electrostatic Energies. Chem. Phys. Lett. 2006, 418, 158–162. [42] Donchev, A. G. Ab Initio Quantum Force Field for Simulations of Nanostructures. Phys. Rev. B 2006, 74, 235401. [43] Rob, F.; Podeszwa, R.; Szalewicz, K. Electrostatic Interaction Energies with Overlap Effects from a Localized Approach. Chem. Phys. Lett. 2007, 445, 315–320. [44] Fraschini, E.; Stone, A. J. H - H Model Potential for Exchange-Repulsion Energy of Methane Dimer. J. Comput. Chem. 1998, 19, 847–857. [45] Stone, A. J.; Misquitta, A. J. Atom-Atom Potentials from Ab Initio Calculations. Int. Rev. Phys. Chem. 2007, 26, 193–222. [46] Sokalski, W. A.; Lowrey, A. H.; Roszak, S.; Lewchenko, V.; Blaisdell, J.; Hariharan, P. C.; Kaufman, J. J. Nonempirical Atom-Atom Potentials for Main Components of Intermolecular Interaction Energy. J. Comput. Chem. 1986, 7, 693– 700. [47] Ma, J. C.; Dougherty, D. A. The Cation-π Interaction. Chem. Rev. 1997, 97, 1303–1324. [48] Salonen, L. M.; Ellermann, M.; Diederich, F. Aromatic Rings in Chemical and Biological Recognition: Energetics and Structures. Angew. Chem. Int. Ed. 2011, 50, 4808–4842.


Page 19 of 23

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


[49] Kim, K. S.; Tarakeshwar, P.; Lee, J. Y. Molecular Clusters of π-Systems: Theoretical Studies of Structures, Spectra, and Origin of Interaction Energies. Chem. Rev. 2000, 100, 4145–4186. [50] Dougherty, D. A. The Cation-π Interaction. Acc. Chem. Res. 2013, 46, 885–893. [51] 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., Current Status of the AMOEBA Polarizable Force Field. J. Phys. Chem. B 2010, 114, 2549–2564. [52] Cubero, E.; Luque, F. J.; Orozco, M. Is Polarization Important in Cation-π Interactions? Proc. Natl. Acad. Sci. U.S.A. 1998, 95, 5976–5980. [53] Kim, D.; Hu, S.; Tarakeshwar, P.; Kim, K. S.; Lisy, J. M. Cation-π Interactions: A Theoretical Investigation of the Interaction of Metallic and Organic Cations with Alkenes, Arenes, and Heteroarenes. J. Phys. Chem. A 2003, 107, 1228–1238. [54] Pratuangdejkul, J.; Jaudon, P.; Ducrocq, C.; Nosoongnoen, W.; Guerin, G. A.; Conti, M.; Loric, S.; Launay, J. M.; Manivet, P. Cation-π Interactions in Serotonin: Conformational, Electronic Distribution, and Energy Decomposition Analysis. J. Chem. Theory Comput. 2006, 2, 746–760. [55] Soteras, I.; Orozco, M.; Luque, F. J. Induction Effects in Metal Cation-Benzene Complexes. Phys. Chem. Chem. Phys. 2008, 10, 2616–2624. [56] Marshall, M. S.; Steele, R. P.; Thanthiriwatte, K. S.; Sherrill, C. D. Potential Energy Curves for Cation-π Interactions: Off-Axis Configurations Are Also Attractive. J. Phys. Chem. A 2009, 113, 13628–13632. [57] Orabi, E. A.; Lamoureux, G. Cation-π and π-π Interactions in Aqueous Solution Studied Using Polarizable Potential Models. J. Chem. Theory Comput. 2012, 8, 182–193. [58] Zheng, X.; Wu, C.; Ponder, J. W.; Marshall, G. R. Molecular Dynamics of β-Hairpin Models of Epigenetic Recognition Motifs. J. Am. Chem. Soc. 2012, 134, 15970–15978. [59] Polly, R.; Werner, H. J.; Manby, F. R.; Knowles, P. J. Fast Hartree-Fock Theory Using Local Density Fitting Approximations. Mol. Phys. 2004, 102, 2311–2321. [60] 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. [61] Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M. Molpro: a General-Purpose Quantum Chemistry Program Package. WIREs Comput. Mol. Sci. 2012, 2, 242–253. [62] Dunning, T. H. Gaussian-Basis Sets for Use in Correlated Molecular Calculations. 1. The Atoms Boron Through Neon and Hydrogen. J. Chem. Phys. 1989, 90, 1007–1023.



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 23

[63] Weigend, F. A Fully Direct RI-HF Algorithm: Implementation, Optimised Auxiliary Basis Sets, Demonstration of Accuracy and Efficiency. Phys. Chem. Chem. Phys. 2002, 4, 4285–4291. [64] Weigend, F.; Köhn, A.; Hättig, C. Efficient Use of the Correlation Consistent Basis Sets in Resolution of the Identity MP2 Calculations. J. Chem. Phys. 2002, 116, 3175–3183. [65] Grüning, M.; Gritsenko, O. V.; van Gisbergen, S. J. A.; Baerends, E. J. Shape Corrections to Exchange-Correlation Potentials by Gradient-Regulated Seamless Connection of Model Potentials for Inner and Outer Region. J. Chem. Phys. 2001, 114, 652–660. [66] Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865–3868. [67] Adamo, C.; Barone, V. Toward Reliable Density Functional Methods Without Adjustable Parameters: The PBE0 Model. J. Chem. Phys. 1999, 110, 6158–6170. [68] 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). [69] Moszynski, R.; Heijmen, T. G. A.; Jeziorski, B. Symmetry-Adapted Perturbation Theory for the Calculation of Hartree-Fock Interaction Energies. Mol. Phys. 1996, 88, 741–758. [70] Heßelmann, A. Comparison of Intermolecular Interaction Energies from SAPT and DFT Including Empirical Dispersion Contributions. J. Phys. Chem. A 2011, 115, 11321–11330. [71] Podeszwa, R.; Cencek, W.; Szalewicz, K. Efficient Calculations of Dispersion Energies for Nanoscale Systems from Coupled Density Response Functions. J. Chem. Theory Comput. 2012, 8, 1963–1969. [72] Singh, N. J.; Min, S. K.; Kim, D. Y.; Kim, K. S. Comprehensive Energy Analysis for Various Types of π-Interaction. J. Chem. Theory Comput. 2009, 5, 515–529. [73] Boys, S. F.; Bernardi, F. Calculation of Small Molecular Interactions by Differences of Separate Total Energies - Some Procedures with Reduced Errors. Mol. Phys. 1970, 19, 553–566. [74] Gaussian 09 Revision B.01. Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; et al. Gaussian Inc. Wallingford CT 2010. [75] Stone, A. J. Distributed Multipole Analysis: Stability for Large Basis Sets. J. Chem. Theory Comput. 2005, 1, 1128–1132; GDMA 2.2.04 at http://www–stone.ch.cam.ac.uk/programs.html. [76] Halgren, T. A. 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.


Page 21 of 23

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


[77] Charbonneau, P. Genetic Algorithms in Astronomy and Astrophysics. Astrophys. J. Suppl. S. 1995, 101, 309–334; PIKAIA 1.2 at http://www.hao.ucar.edu/modeling/pikaia/pikaia.php. [78] Ren, P.-Y.; Wu, C.-J.; Ponder, J. W. Polarizable Atomic Multipole-Based Molecular Mechanics for Organic Molecules. J. Chem. Theory Comput. 2011, 7, 3143–3161. [79] Sa, R. J.; Zhu, W. L.; Shen, J. H.; Gong, Z.; Cheng, J. G.; Chen, K. X.; Jiang, H. L. How Does Ammonium Dynamically Interact with Benzene in Aqueous Media? A First Principle Study Using the Car-Parrinello Molecular Dynamics Method. J. Phys. Chem. B 2006, 110, 5094–5098. [80] Reddy, A. S.; Sastry, G. N. Cation [M = H+, Li+, Na+, K+, Ca2+, Mg2+, NH4+, and NMe4+] Interactions with the Aromatic Motifs of Naturally Occurring Amino Acids: A Theoretical Study. J. Phys. Chem. A 2005, 109, 8893–8903. [81] Dehez, F.; Angyan, J. G.; Gutierrez, I. S.; Luque, F. J.; Schulten, K.; Chipot, C. Modeling Induction Phenomena in Intermolecular Interactions with an Ab Initio Force Field. J. Chem. Theory Comput. 2007, 3, 1914–1926. [82] Wang, Z.-X.; Zhang, J.-C.; Cao, W.-L. Study on the Cation-π Interactions between Ammonium Ion and Aromatic π Systems. Chin. J. Chem. 2006, 24, 1523–1530. [83] Mo, Y.; Bao, P.; Gao, J. L. Energy Decomposition Analysis Based on a Block-Localized Wavefunction and Multistate Density Functional Theory. Phys. Chem. Chem. Phys. 2011, 13, 6760–6775. [84] Aschi, M.; Mazza, F.; Di Nola, A. Cation-π Interactions between Ammonium Ion and Aromatic Rings: an Energy Decomposition Study. THEOCHEM 2002, 587, 177–188. [85] Voisin, C.; Cartier, A.; Rivail, J. L. Computation of Accurate Electronic Molecular Polarizabilities. J. Phys. Chem. 1992, 96, 7966–7971. [86] 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 5.1, February 2010, at http://dasher.wustl.edu/tinker/. [87] Elking, D.; Darden, T.; Woods, R. J. Gaussian Induced Dipole Polarization Model. J. Comput. Chem. 2007, 28, 1261–1274. [88] Kim, H. W.; Rhee, Y. M. Molecule-Specific Determination of Atomic Polarizabilities with the Polarizable Atomic Multipole Model. J. Comput. Chem. 2012, 33, 1662–1672. [89] Patkowski, K.; Szalewicz, K.; Jeziorski, B. Orbital Relaxation and the Third-Order Induction Energy in SymmetryAdapted Perturbation Theory. Theor. Chem. Acc. 2010, 127, 211–221. [90] Lao, K. U.; Herbert, J. M. Breakdown of the Single-Exchange Approximation in Third-Order Symmetry-Adapted Perturbation Theory. J. Phys. Chem. A 2012, 116, 3042–3047.



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 23

[91] Schäffer, R.; Jansen, G. Intermolecular Exchange-Induction Energies Without Overlap Expansion. Theor. Chem. Acc. 2012, 131, 1235. [92] Hobza, P. Calculations on Noncovalent Interactions and Databases of Benchmark Interaction Energies. Acc. Chem. Res. 2012, 45, 663–672. [93] Martinez, C. R.; Iverson, B. L. Rethinking the Term ”pi-stacking”. Chem. Sci. 2012, 3, 2191–2201.


Page 23 of 23

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 OF CONTENTS GRAPHIC

R

ACS Paragon Plus Environment

CationâÏ Interactions: Accurate Intermolecular Potential from

Recommend Documents

CationâÏ Interactions: Accurate Intermolecular Potential from