Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
J Chem Theory Comput. Author manuscript; available in PMC 2010 November 23.
Published in final edited form as:
J Chem Theory Comput. 2009 January 1; 5(11): 3022–3031.
doi:  10.1021/ct9004189
PMCID: PMC2990227

Polarizable intermolecular potentials for water and benzene interacting with halide and metal ions


A complete derivation of polarizable intermolecular potentials based on high-level, gas-phase quantum-mechanical calculations is proposed. The importance of appreciable accuracy together with inherent simplicity represents a significant endeavor when enhancement of existing force fields for biological systems is sought. Toward this end, symmetry-adapted perturbation theory (SAPT) can provide an expansion of the total interaction energy into physically meaningful e.g. electrostatic, induction and van der Waals terms. Each contribution can be readily compared with its counterpart in classical force fields. Since the complexity of the different intermolecular terms cannot be fully embraced using a minimalist description, it is necessary to resort to polyvalent expressions capable of encapsulating overlooked contributions from the quantum-mechanical expansion. This choice results in consistent force field components that reflect the underlying physical principles of the phenomena. This simplified potential energy function is detailed and definitive guidelines are drawn. As a proof of concept, the methodology is illustrated through a series of test cases that include the interaction of water and benzene with halide and metal ions. In each case considered, the total energy is reproduced accurately over a range of biologically relevant distances.

1 Introduction

The development of increasingly scalable software and larger access to massively parallel computer architectures allow longer numerical simulations to be performed on continuously more complex molecular assemblies. This, in turn, has brought the investigation of biologically relevant systems over realistic time scales within reach, e.g. atomistic simulations of a membrane for hundred of microseconds1 or folding of small proteins.2? In the latter examples, the use of a pairwise additive force field, e.g. Amber,3 Charmm,4 Gromos5 or Opls-AA,6 provides in general a reasonable reproduction of experimental data as long as polarization effects are not dominant.

When, on the contrary, induction effects can no longer be neglected, or vary significantly in the course of the simulation, non-polarizable potential energy functions rapidly attain the inherent limits of their validity. Two recent studies7,8 endeavored to describe the folding of a small protein, the human Pin1 WW domain, by means of the Charmm macromolecular force field. Whereas this potential energy function was able to preserve the β-sheet-like native structure over the 200 ns time scale, it failed to fold within 10 µs the protein chain from an extended motif into the native conformation. It was demonstrated that this failure is due to force-field accuracy rather than insufficient sampling — whether or not neglect of polarization phenomena is responsible for protein misfolding remains at this stage a matter of debate, which would require additional investigation to find a definitive answer. Metalloproteins constitute another class of biologically relevant systems where induction effects need to be accounted for. In this case, a metal ion polarizes the environment, which, in turn, alters the binding hierarchy9 of neighboring chelating agents. In particular,10,11 when using a standard non-polarizable potential energy function, the stability of the active site is perturbed on account of the exaggerated interaction of the ion with water, compared to that with either the vicinal carbonyl and carboxylate groups. Ion channels also represent important biological systems, in which small ions permeate through narrow pores, strongly polarizing the walls of the latter,12 and modulating as a result the ionic selectivity.13 Formation of crystals can also be governed by polarization effects.1417 When decomposing the total interaction energy into electrostatic, induction and van der Waals contributions, the importance of the induction energy suggests that it may play a predominant role in the stability of the crystal and, hence, should not be ignored.

The need for taking into account polarization phenomena explicitly in molecular-dynamics simulations constitutes an ongoing effort that can be traced back to the eighties.18 One of the grand challenges for polarizable force fields is the derivation of functions that are suitable for biological simulations. Towards this end, a number of routes have been explored to include polarization effects explicitly. A first route is based on Drude oscillators,19 which describe electronic induction by the movement of a fictitious mass bonded to polarizable atoms by means of a stiff spring. Another approach, referred to as fluctuating charges,2023 consists in changing the charge of the atoms along the simulation. These charges are adjusted according to the principle of partial equalization of orbital electronegativities.24 A third path, namely multipole expansions,21,2530 uses distributed polarizabilities to model induction energies following a self-consistent procedure to determine the induced dipoles. These methods have been probed in a variety of applications ranging from liquid simulations3134 to simulations of DNA.35 They have been combined with molecular-dynamics simulations of increasing complexity, such as that of polarizable membranes — see for instance Ref. 36. In the latter reference, Drude oscillators model induction effects and allow a better agreement to be reached between simulation and experiment for the reproduction of the dipole potential that arises at the water-lipid interface.

The main thrust of the present work is to build a polarizable force field that reconciles simplicity and accuracy. The chosen approach relies on the central idea37 of partitioning the interaction energy into physically meaningful contributions and to determine the associated parameters using quantum-mechanical methods. A logical road map is followed for the construction of a consistent potential energy function.9,3843 This consistency imposes a de novo derivation of gas-phase atomic charges. The induction contribution is then derived from ab initio calculations and modified by ad hoc damping functions. Last, van der Waals parameters are determined using an appropriate mathematical expression that guarantees the faithful description of the overall intermolecular interactions at play. Despite the precision of such a description, a number of key elements are still missing for a faithful reproduction of non-bonded interaction energies between molecules over a wide range of distances. In the reproduction of the total quantum-mechanical interaction energy, the rudimentary description needed for large-scale simulations of biological systems is bound to failure if high accuracy is sought, as several terms, like non-multipolar contributions are not described explicitly. Choices have, therefore, to be made to obtain a model of sufficient accuracy, while remaining concise. Appropriate functions ought to be determined to include implicitly the missing contributions. In this study, a complete set of definitive guidelines is proposed to derive compact polarizable intermolecular potentials. The approach is illustrated through a series of test cases that include the interaction of water and benzene with halide and metal ions.

2 Theoretical underpinnings

The quantum-mechanical (QM) calculations reported here rely on symmetry-adapted perturbation theory44,45 (SAPT), which provides a formal, rigorous framework for the expansion of the total interaction energy into meaningful contributions utilized as a reference. In particular, use is made in this work of the SAPT2 expansion, henceforth called SAPT for simplicity, which accounts for correlation effects at the second-order Møller-Plesset46 (MP2) level of theory. Within SAPT, the interaction energy is decomposed into the individual contributions shown in Eq. 1. In the latter expression, the first four terms stand for the electrostatic (Uele), induction (Uind), exchange (Uexc) and dispersion (Udis) components. In addition, the SAPT expansion embraces two terms corresponding to the coupling between exchange and induction (Uexc-ind), and between exchange and dispersion (Uexc-dis). Finally, the last term (δHF) accounts for a collection of higher-order induction and exchange-induction terms. Altogether, the contributions are written


The stated energy components cannot, when taken separately, be studied in a straightforward fashion with a classical decomposition.47 For a comparison on a one-to-one basis, certain individual contributions need to be summed, as has been done previously by several authors.9,4851 It is then possible to cast the description into three potential energy functions, namely, for electrostatic, induction and van der Waals contributions.

The electrostatic component (Δ𝒰eleSAPT) corresponds to the first-order contribution to the electrostatic energy. The pure induction term given by SAPT ought to be added to the exchange-induction term, as well as to higher-order induction and exchange-induction terms to encompass the full induction potential (Δ𝒰indSAPT). The last contribution, referred to as van der Waals potential (Δ𝒰vdWSAPT), consists of the sum of exchange and dispersion terms.

As was put forth previously, a molecular-mechanical expansion can mimic the above summed QM energies. This includes an electrostatic potential that relies upon atomic charges, an induction term associated to an appropriate damping function and an ad hoc van der Waals function.9

The first order contribution, of electrostatic nature, is modeled by a set of atom-centered point charges and is a rudimentary approach to describe the molecular electrostatic potential. While this description is simple, it has proven52 to reproduce under most circumstances the target quantity with an appreciable accuracy. Determination of atomic charges is achieved53 by means of regular grids surrounding the molecule of interest, over which the QM electrostatic potential is mapped. Penetration effects, which represent an important, short-range component of the electrostatic interaction, are evidently absent from the classical, Coulomb potential.54 The incompleteness of the model is rooted in the computation of the molecular electrostatic potential outside the sphere of convergence.55 To a large extent, these effects stem from the overlap of electron densities of interacting compounds, which cannot be accounted for by a localized atomic description of the electrostatic potential.56 Ad hoc functions5760 have been shown to mimic the respective short-range interactions. Here, use will only be made of an ansatz inferred from the formulation of Freitag et al. by Piquemal and coworkers60 — henceforth referred to as Cisneros et al. It holds the advantage of being derived explicitly for atomic charges and, compared to other functions, it offers a simpler expression and parametrization. Charges are, therefore, modified as stated in equation (3) in Ref. 60.

As detailed above, different routes can be chosen for the description of induction phenomena in the classical simulations of biological systems. Here, polarization is modeled by means of distributed polarizabilities61,62 over subsets of atoms. Previous articles9,43,63 have demonstrated that the use of charge-flow polarizabilities between chemically bonded atoms supplemented by isotropic dipole polarizabilities on heavy atoms is able to recover the anisotropy of the molecular polarizability with both simplicity and effectiveness. Polarizabilities were determined by means of a fitting procedure41 that relies on three-dimensional maps of the induction energy evaluated quantum-mechanically. As was asserted previously,9,43,63 the classical expansion of polarization ought to be damped in order to reproduce the destabilizing effect of exchange-induction effects occurring at short intermolecular distances. In the present work, the Jensen et al.64 and the Tang and Toennies65,66 (TT) functions have been chosen. In the former expression, the distances in the interaction tensor are modified by a scalar parameter, a, that appears in equation (16) in Ref. 64. In the latter expression, truncation of the series is done at the sixth order, as suggested by Millot and Stone,67 and subsequently rewritten at the third order, as proposed by Meredith and Stone.68 The final expressions for the induction energy utilized herein are:




where Qta and ΔQta are, respectively, the permanent moment and the induced moment of rank t located at site a of molecule A, Ttuab is the element of the electrostatic tensor for the interaction of the multipole moments of ranks t and u, found, respectively, at sites a and b, and f3(β;rab) is the truncated Tang and Toennies damping function, in which β is a scalar parameter. αttaa is the atomic polarizability, which describes the change in the multipole moment of rank t at site a resulting from the t′–th derivative of the potential created by all other moments at site a′. The above expressions are also detailed in equation (21) and equation (22) in Ref. 50.

The last term of the classical intermolecular potential is the van der Waals contribution, which is modeled in the present work by three alternate formulations depending on the r distance. The first one is the standard 6–12 Lennard-Jones (LJ) expression,69 viz. Δ𝒰vdWMM/LJ=ε[(σ/r)12(σ/r)6], where the parameters are ε and σ correspond respectively to the depth of the potential well and the distance at which the interparticle potential is zero. The second one, used among others in the Amoeba70 force field, is the Halgren71 potential, viz. Δ𝒰vdWMM/Halgren=ε[1.07R*/(r+0.07R*)]7×[1.12R*7/(r7+0.12R*7)2] , where R* is the minimum-energy distance. The third one is the exp–6 Buckingham function,72 viz. Δ𝒰vdWMM/Buckingham=εexp(A1r)(A2/r)6, where A1 and A2 are shape parameters.

3 Computational details

All molecular geometries were optimized using the Gaussian03 suite of programs73 at the MP2(Full)/6–311++G(2d,2p) level of approximation. Electronic properties were computed at the MP2(Full)/Sadlej level of approximation, considering that the Sadlej7477 basis set provides a very good compromise between number of Gaussian contractions and accuracy.78 All individual contributions to the interaction energy were determined with the SAPT2008 program of Jeziorski et al.,44 interfaced to the ATMOL integral and self-consistent-field package.79 The QM contributions to the total interaction energy serve as a basis of comparison for the parametrization of the classical model.

The atomic charges of water and benzene were derived employing the Opep package,42 using 2,375 and 8,070 grid points, respectively. As shown in Table 1, the molecular electrostatic potential of an isolated water molecule is reproduced with a mean error of about 48 % when use is made of atom-centered charges only. As water is one of the most important components in computer simulations of condensed phases, in particular in biological systems, an appreciable error in the description of its electrostatic potential is not acceptable. An improvement of the rudimentary three-point (3-p) charge model is, therefore, proposed. The simplest way to improve the reproduction of the potential consists in adding a fictitious site along the bisector of the molecule, in the spirit of the TIP4P80,81 model.

Table 1
Models of net atomic charges and regenerated multipole moments for water and benzene at the MP2(Full)/Sadlej//MP2(Full)/6-311++G(2d,2p) level of approximation a

As can be seen in Figure 1, a minimum in the root-mean-square deviation between the QM and the point-charge derived potentials can be found when moving the fictitious site. At this minimum, the additional charge is located 0.276 Å below the oxygen atom (see “Water (4-p)” in Table 1). The resulting enhanced agreement with the target electrostatic potential can be largely explained by the position of the supplementary charge, closer to the center of the molecule. This result is consistent with the nearly isotropic molecular electrostatic potential of water.54 In the case of benzene, as was shown in a recent article,82 the use of atom-centered charges is sufficient to model the electrostatic potential with an acceptable accuracy. As can be seen in Table 1, a comparable agreement between the QM potential and that regenerated from point charges is also found here.

Figure 1
(a) Root-mean-square deviation (RMSD) between the electrostatic potentials determined quantum-mechanically and the point charge model expressed as a function of the position of the dummy atom for water along the C2 axis. (b) Corresponding mean error (Δε) ...

Penetration of electron clouds results in a short-range contribution, which, for completeness and for attaining the desired accuracy in the reproduction of the intermolecular interaction energy, ought to be included in the electrostatic model. Towards this end, a novel parametrization of the function of Cisneros et al. was probed, wherein the term Ωij of equation (4) in Ref. 60 had been modified. Though conceptually simple, the expression proposed to account for penetration effects lacks generality to be of routine practical use for the modeling of large biomolecular assemblies. For this reason, no analytical correction of penetration effects will be introduced explicitly in the intermolecular potential, but will instead be incorporated in the van der Waals potential.

The induction energy was mapped employing the Opep package,42 with grids of 3,905 and 3,192 points, for water and benzene, respectively (see Ref. 9 and 43 for details). The models of distributed polarizabilities utilized herein rely on a multipole expansion of the molecular polarizability reduced to a combination of zeroth-order charge-flow (α00,00) and first-order isotropic dipole (α1κ,1κ) polarizabilities. For clarity, it should be noted that charge-flow polarizabilities are not equivalent to the popular fluctuating charges,2023 which are based on a partial equalization of orbital electronegativities and, hence, are very different in spirit. The complete set of parameters fitted simultaneously to quantum-mechanically determined induction energy maps53 have been shown to reproduce with an appreciable accuracy the average molecular polarizability alongside with its intrinsic anisotropy. The distributed polarizabilities are reported in Table 2 for water and benzene and in Table 3 for all the ions considered as polarizable in this work.

Table 2
Models of distributed polarizabilities and regenerated molecular polarizabilities of water and benzene at the MP2(Full)/Sadlej//MP2(Full)/6–311++G(2d,2p) level of approximation a
Table 3
Dipole polarizability of ions computed at the MP2(Full)/Sadlej level of approximation

Damping parameters are determined with the objective to reproduce the induction energy supplied by SAPT. The first damping function considered is that of Jensen et al., the parameters of which can be found for a number of atom types. For the Tang and Toennies damping function, parameters are given for pair interactions. In some cases, no damping is necessary, due to the importance of the non-multipolar contribution to the induction energy, as has been documented by Torheyden and Jansen.49 This contribution to the induction cannot be extracted from the SAPT expansion and, hence, cannot be modeled explicitly in the model proposed here. This missing information can still be added in an appropriately chosen form of the van der Waals potential, compatible with the exponential, repulsive nature of non-multipolar contributions.

The last contribution to the total energy is the van der Waals potential, which, in the SAPT expansion is considered as being the sum of exchange and dispersion terms. As outlined previously, both the electrostatic and the induction models adopted here lack a number of contributions that cannot be accounted for to satisfy the stringent criteria of simplicity and tractability. Since the overall model is not sufficiently precise at short distances and when dealing with appreciable non-multipolar effects, determination of the van der Waals parameters was not carried out based on the exact contribution supplied by the SAPT expansion, but using a slightly modified term, which is the total energy (Δ𝒰totSAPT) minus the electrostatics and induction components of the molecular-mechanical potential. This statement implies that the optimized van der Waals parameters inherently depend on the electrostatic and induction models that have been used here. In other words, modification of either the electrostatic or the polarization term necessarily imposes that the van der Waals potential be adjusted. Derivation of these parameters was achieved employing a Levenberg-Marquart algorithm.85 The range of chemically and biologically relevant distances spans from a minimum separation, where the energy is roughly equal to Δ𝒰tot,minSAPT10kcal/mol – 10 kcal/mol, to a quasi-infinite one. The van der Waals contribution was determined based on three alternate functions, i.e. Lennard-Jones,69 Halgren71 and Buckingham.72 The exponential nature of the repulsive van der Waals contribution makes the Buckingham exp-6 potential a well-suited candidate to model van der Waals contribution and short-range interactions, such as non-multipolar terms.

4 Results and Discussion

4.1 Interaction of water with metal cations

Classical models derived to describe the interaction of cations and anions with water molecules are to a large extent suboptimal. This can be ascribed to three main reasons, recently highlighted in a study covering the interaction of water with a divalent calcium cation.9 First, a significant error arises on account of the rudimentary representation of the molecular electrostatic potential for water by means of a three-point (3-p) charge model, which misreproduces the target quantity by up to 50 %. Second, the description of the exchange-induction term of the SAPT expansion using the damping function devised by Jensen et al., while constituting a judicious choice, appears to suffer from a number of shortcomings. Last, the most questionable aspect of the proposed polarizable model lies in the use of the 6–12 Lennard-Jones potential, traditionally utilized in macromolecular force fields. As was shown, however, this potential systematically fails to reproduce in a consistent fashion the exchange and dispersion contributions to the intermolecular energy computed within the SAPT framework. This shortcoming has prompted several authors to turn to more appropriate functions.6,27,8688 To address the above issues, the interaction of a series of metal ions with water has been investigated (see Figure 2). For conciseness, only the water-sodium interaction will be detailed here (see Figure 3); results for all the systems are provided in Supporting Information.

Figure 2
(a) Interaction of a water molecule with a cation along the C2 axis. (b) Interaction of a water molecule with an anion along the O–H bond. (c) Interaction of a water molecule with another water molecule defining a hydrogen bond. (d) Interaction ...
Figure 3
Interaction of a monovalent sodium ion with water. The reference (SAPT) curve is shown as a solid line. Three models for the electrostatic contribution are represented through either no penetration with a three-point (3-p) charge model (dotted line), ...

As can be observed in Figure 3, the electrostatic contribution is equally well reproduced using a 4-p or the 3-p charge model. This result is not surprising given that the cation approaches water along the C2 axis of the latter. Error in the reproduction of the molecular electrostatic potential in the direction of the bisector, on the order of 5 % at the equilibrium intermolecular distance, is equivalent for the two models. This, unfortunately, is no longer true when the molecular electrostatic potential is measured quantum-mechanically along the O–H bond. In this event, the error reaches 15 % for the 3-p charge model, but only 5 % for the 4-p one. The electrostatic term modified by a penetration function is also displayed in Figure 3. This modification suggested by Cisneros et al. improves the accuracy of the model at short distances, even though, for water-cation dimers, the effect of electron-cloud penetration is nearly negligible at distances around the intermolecular equilibrium separation (around 2.2 Å). As underlined in the previous section, this correction to the electrostatic potential will not be introduced explicitly, but will, instead, be embedded in the relevant representation of the van der Waals potential, thereby simplifying the model, while remaining physically consistent.

The undamped induction energy of the classical model overestimates the quantity supplied by the SAPT expansion.9 The damping functions considered here are that of Jensen et al. and that of Tang and Toennies. The latter function is to be favored over the former, as is reflected in the better behavior of the corresponding damped potential. All the parameters determined for these damping functions are gathered in Table 4 and Table 5.

Table 4
Values of the Jensen et al. damping parameter (a) for each atom type
Table 5
Values of the damping parameter (β) in the Tang and Toennies damping function. Each value corresponds to a pair of interaction

The negligible errors in the reproduction of the electrostatic and induction contributions in water-cation complexes makes the total energy (Δ𝒰totSAPT) minus Δ𝒰ele4p and Δ𝒰ind4p/TT very similar for the SAPT van der Waals contribution (Δ𝒰vdWSAPT). Employing a traditional 6–12 LJ function, the repulsive component of the latter results in a systematic failure to match the target van der Waals potential — the associated error is magnified when the induction effects dominate over the Coulomb and dispersion components. Conversely, if use is made of a Halgren or a Buckingham function, reproduction of this potential proves to be extremely accurate. Although depicted in Figure 3, the Halgren function will not be used hereafter, because the exponential part of the exp–6 Buckingham function constitutes a more physically sound framework, capable of absorbing the errors due to penetration effects. Interestingly enough, non-multipolar contributions to the induction energy appear to be negligible in water-cation interactions, and can, thus, be safely omitted. All the parameters of the aforementioned three van der Waals functions are given in Table 6. These parameters are the raw coefficients for each interaction in the absence of combination rules.

Table 6
van der Waals parameters for H2O–X systems. The parameters are obtained using a model of point charges and an induction model of distributed polarizabilities damped by the function of Tang and Toennies. (1) ε and An (n = 1, 2) correspond ...

When adding the different terms of the potential energy function, the total interaction energy is reproduced accurately, hence, confirming the applicability of the model (see Table 7 and Table 8 for geometrical and energetic properties of the minimum-energy complex). Not too unexpectedly, calculations based on a LJ potential lead to a severely flawed reproduction of the total interaction energy.

Table 7
Values of the interaction distance for all the minimum-energy complexes. The reference quantum-mechanical values are given in bold
Table 8
Values of the interaction energy for all the minimum-energy complexes. The reference quantum-mechanical values are given in bold

4.2 The water dimer and interaction of water with anions

In this section, the first four halide ions have been considered to interact with water along one O–H bond of the latter, which corresponds to the most favorable approach of these chemical species. In addition, a water dimer has been examined in the geometry provided in Ref. 89 (see also Figure 2). The results obtained for these two classes of complexes show similar trends. For conciseness, only the interaction in the water dimer will be detailed here (see Figure 4); results for all the systems are provided in Supporting Information.

Figure 4
Interaction of two water molecules. The reference (SAPT) curve is shown as a solid line. Three models for the electrostatic contribution are represented through either a three-point (3-p) model with no penetration (dotted line), a four-point (4-p) model ...

The importance of the 4-p charge model ought to be emphasized for the reproduction of the electrostatic potential. Even at distances where penetration effects are known to be negligible, the 3-p charge model fails to reproduce the electrostatic contribution, whereas the 4-p charge model describes it successfully. Inclusion of a penetration correction yields an overall agreement between the classical model and the QM reference, irrespective of the intermolecular distance explored. Penetration effects for water-anion systems, as well as for the water dimer are crucial, as they account for up to 50 % of the total electrostatic contribution at very short distances. They are, nonetheless, noticeably smaller (approximately 8 %) at an intermolecular distance corresponding to the minimum-energy complex. Accordingly, it is reasonable to expect that penetration effects can be transferred in an effective fashion into the van der Waals potential, as discussed previously.

When penetration effects are large, non-multipolar contributions are anticipated to be equally substantial. In this event, the induction term is burdened by overwhelming non-multipolar stabilizing effects, which counterbalance the destabilizing effect of the exchange-induction component. These antagonistic phenomena preclude the use of a damping function for the present complexes. As can be seen in Table 5, no damping function is used, except for fluoride, which is the smallest anion with the least non-multipolar effects. The missing non-multipolar contribution to the induction energy will, thus, be incorporated into the Buckingham potential.

The conclusions reached in light of the calculations based on the water-cation dimers appear to hold also for the present complexes. Furthermore, the assumption that an exp–6 Buckingham potential, unlike a 6–12 LJ potential, can compensate for the incompleteness of the classical model by embracing its overlooked, sizable contributions, proves to be perfectly legitimate and reasonable (see Table 7 and Table 8).

In the particular instance of negatively charged compounds, the SAPT expansion can, under certain circumstances, fail to converge.90 For this reason, the results (not shown here) supplied by an SAPT and an RVS91 expansion, using the Gamess92 suite of programs, have been compared. If intramolecular correlation is ignored, the RVS and the SAPT expansions yield the same energies for all the dimers considered here, which provides a safeguard, in as much as convergence of the SAPT reference calculation is concerned. Moreover, the recent work of Kim et al.93 on halide anions interacting with π-systems further suggests convergence of the SAPT expansion for such dimers.

4.3 Interaction of benzene with metal cations

The same series of metal cations utilized in the above study of water-cation dimers was considered for the interaction with benzene along its C6 axis (see Figure 2). This approach corresponds to the most favorable interaction pattern94 for a π-electron cloud interacting with a positively charged species. For conciseness, only the benzene-sodium interaction will be displayed in Figure 5; results for all systems investigated (viz. Li+, K+, Mg2+ and Ca2+) are provided in Supporting Information.

Figure 5
Interaction of a monovalent sodium ion with benzene. The reference (SAPT) curve is shown as a solid line. The electrostatic contribution obtained with a twelve-point (12-p) model with no penetration is represented as a dotted line. Two models for the ...

A twelve-point (12-p) charge model derived from the molecular electrostatic potential tends to overestimate the interaction in cation-π systems, despite the appropriate reproduction of the molecular quadrupole moment of benzene and, more generally, its full molecular electrostatic potential. This overestimation is more prone to be manifested with ions bearing a highly-localized charge and, on the contrary, is expected to be less perceptible for a diffuse charge, as in the case of ammonium,43 let alone guanidinium. The results suggest that the interaction of cations with π-systems is not strictly speaking an attractive quadrupole-charge one.94 A case in point is the benzene-lithium dimer, where the electrostatic energy changes its slope at very short distances, which, in all likelihood, is due to the positively charged lobes of the quadrupole moment interacting with the cation. The resulting destabilizing effect, which largely counterbalances the penetration contribution, is not accounted for in the electrostatic potential and is transferred to a short-distance, repulsive van der Waals term.

As noted above for the cation-water complexes, the undamped induction energy overestimates the values obtained from an SAPT analysis, but inclusion of the Tang and Toennies damping formalism reproduces the SAPT induction profile very nicely. Similarly, the 6–12 LJ potential departs markedly from the SAPT van der Waals contribution, which is always better reproduced when Halgren or Buckingham functions are used. Conclusions drawn on the basis of the water-cation complexes are shown to hold for all other π-cation interactions reported herein (van der Waals parameters are gathered in Table 9 and the geometrical and energetic properties of the minimum-energy complex are given in Table 7 and Table 8).

Table 9
van der Waals parameters for benzene–Xn+ systems. The parameters are obtained using a model of point charges on atoms with no penetration and an induction model of distributed polarizabilities damped by the function provided by Tang and Toennies. ...

5 Conclusion

The present study addresses the shortcomings of a simplified polarizable potential energy function and how these shortcomings can be overcome to model with an appreciable accuracy intermolecular interactions. It is desirable to have access to a simple polarizable potential energy function that satisfies the criteria of tractability and precision. This necessarily demands that important choices be made, inasmuch as the contributions to be included are concerned. In other words, it is pivotal that all the underlying physics of the interactions be properly introduced to reflect the various contributions at play.

With this objective in mind, a series of guidelines is proposed. These guidelines can be summarized as follows. The electrostatic contribution to the total energy is described by a truncation of the molecular electrostatic potential at the monopole level, which has proven in many instances to constitute a very reasonable approximation.52 As indicated by Swart et al.,95 this description leads to small errors for non-symmetric systems, but can fail in the case of small, symmetric molecules. Such is the case for the water molecule, where fictitious sites are needed to enhance the reproduction of the true molecular electrostatic potential. Accordingly, a TIP4P-like80,81 potential was proposed for water. In addition to the first-order contribution to the total energy, penetration effects are present at short distances in all interacting systems. This term can be modeled with an exponential function and can be introduced through the exp–6 Buckingham potential. The total interaction energy is further decomposed into an induction component, which, at the classical level, relies on the conjunction of charge-flow polarizabilities between covalently bonded atoms and isotropic dipole polarizabilities determined from QM induction energy maps.

The resulting polarizable models are simple, yet sufficiently accurate for recovering the anisotropy of the molecular polarizability. In most circumstances, the damping function proposed by Tang and Toennies recovers the exchange-induction contribution. This particular function can be expanded up to third order only, which greatly simplifies its calculation. Furthermore, when the non-multipolar component to the induction energy cannot be neglected, a classical multipole expansion systematically underestimates the polarization part of the interaction energy. Should this situation arise, it is necessary to remove the damping function from the model and encapsulate this stabilizing effect through another suitable function. The remaining terms of the QM expansion that are described explicitly can be incorporated in a single contribution that involves the exp–6 Buckingham potential. The latter is highly appropriate because it encompasses each and every missing term that can be expressed by way of an exponential, repulsive representation.

Extension of the present work to an all-atom description for the van der Waals potential is underway to provide a better three-dimensional description of the short-range potential. The very encouraging results reported herein suggest that the proposed methodology for deriving polarizable force fields on the basis of a small set of prototypical interacting chemical species can be readily extended to more complex molecular systems.

Supplementary Material


The authors thank János G. Ángyán, Georg Jansen and Krzysztof Szalewicz for stimulating discussions. The authors are indebted to the Centro de Supercomputación de Catalunya, Barcelona, Spain and to the Centre Informatique National de l’Enseignement Supérieur, Montpellier, France, for provision of generous amounts of computational time.


Supporting Information Available

Intermolecular energy profiles for water and benzene interacting with halide and metal ions. This material is available free of charge via the Internet at


1. Yin Y, Arkhipov A, Schulten K. Structure. 2009;17:882–892. [PMC free article] [PubMed]
2. Lei H, Duan Y. J. Mol. Biol. 2007;370:196–206. [PMC free article] [PubMed]
3. Cornell WD, Cieplak P, Bayly CI, Gould IR, Merz KM, Jr, Ferguson DM, Spellmeyer DC, Fox T, Caldwell JC, Kollman PA. J. Am. Chem. Soc. 1995;117:5179–5197.
4. MacKerell> AD., Jr et al. J. Phys. Chem. B. 1998;102:3586–3616. [PubMed]
5. Oostenbrink C, Villa A, Mark AE, van Gunsteren WF. J. Comput. Chem. 2004;25:1656–1676. [PubMed]
6. Kaminski GA, Stern HA, Berne BJ, Friesner RA. J. Phys. Chem. A. 2003;108:621–627.
7. Freddolino PL, Liu F, Gruebele M, Schulten K. Biophys. J. 2008;94:75–77. [PubMed]
8. Freddolino PL, Park S, Roux B, Schulten K. Biophys. J. 2009;96:3772–3780. [PubMed]
9. Dehez F, Archambault F, Soteras I, Luque FJ, Chipot C. Mol. Phys. 2008;106:1685–1696.
10. Shiratori Y, Nakagawa S. J. Comp. Chem. 1991;12:717–730.
11. Project E, Nachliel E, Gutman M. J. Comp. Chem. 2008;29:1163–1169. [PubMed]
12. Bucher D, Guidoni L, Maurer P, Rothlisberger U. J. Chem. Theory Comput. 2009 In press. [PubMed]
13. Allen T, Andersen O, Roux B. Biophys. Chem. 2006;124:251–267. [PubMed]
14. Fowler PW, Stone AJ. J. Phys. Chem. 1987;91:509–511.
15. Gavezzotti A. J. Phys. Chem. B. 2002;106:4145–4154.
16. Gavezzotti A. J. Phys. Chem. B. 2003;107:2344–2353.
17. Welch GWA, Karamertzanis PG, Misquitta AJ, Stone AJ, Price SL. J. Chem. Theory Comput. 2008;4:522–532. [PubMed]
18. Barnes P, Finney JL, Nicholas JD, Quinn JE. Nature. 1979;282:459–464.
19. Lamoureux G, Roux B. J. Chem. Phys. 2003;119:3025–3039.
20. Rappe AK, Goddard WA. J. Phys. Chem. 1991;95:3358–3363.
21. Stern HA, Kaminski GA, Banks JL, Zhou R, Berne BJ, Friesner RA. J. Phys. Chem. B. 1999;103:4730–4737.
22. Patel S, Brooks CL. J. Comput. Chem. 2003;25:1–16. [PubMed]
23. Patel S, Mackerell AD, Jr, Brooks CL. J. Comput. Chem. 2004;25:1504–1514. [PubMed]
24. Sanderson RT. Chemical Bonds and Bond Energy. New York: Academic Press; 1976. pp. 1–15.
25. Le Sueur CR, Stone AJ. Mol. Phys. 1994;83:293–307.
26. Engkvist O, Å strand P-O, Karlström G. Chem. Rev. 2000;100:4087–4108. [PubMed]
27. Ren P, Ponder JW. J. Comput. Chem. 2002;23:1497–1506. [PubMed]
28. Hagberg D, Brdarski S, Karlström G. J. Phys. Chem. B. 2005;109:4111–4117. [PubMed]
29. Gresh N, Cisneros GA, Darden TA, Piquemal J-P. J. Chem. Theory Comput. 2007;3:1960–1986. [PMC free article] [PubMed]
30. Darley MG, Popelier PLA. J. Phys. Chem. A. 2008;112:12954–12965. [PubMed]
31. Berendsen HJC, Grigera JR, Straatsma TP. J. Phys. Chem. 1987;91:6269–6271.
32. Masella M, Flament JP. J. Chem. Phys. 1997;107:9105–9116.
33. Zhang Q, Zhang X, Yu L, Zhao D-X. J. Mol. Liq. 2009;145:58–66.
34. Zhang Q, Zhang X, Zhao D-X. J. Mol. Liq. 2009;145:67–81.
35. Babin V, Baucom J, Darden TA, Sagui C. J. Phys. Chem. 2006;110:11571–11581. [PubMed]
36. Harder E, MacKerell AD, Jr, Roux B. J. Am. Chem. Soc. 2009;131:2760–2761. [PMC free article] [PubMed]
37. Claverie P. In: Intermolecular Interactions: From Diatomics to Biopolymers. Pullman B, editor. Wiley Interscience; 1978. pp. 69–305.
38. Celebi N, Ángyán JG, Dehez F, Millot C, Chipot C. J. Chem. Phys. 2000;112:2709–2717.
39. Dehez F, Soetens J-C, Chipot C, Ángyán JG, Millot C. J. Phys. Chem. A. 2000;104:1293–1303.
40. Chipot C, Dehez F, Ángyán JG, Millot C, Orozco M, Luque FJ. J. Phys. Chem. A. 2001;105:11505–11514.
41. Dehez F, Chipot C, Millot C, Ángyán JG. Chem. Phys. Lett. 2001;338:180–188.
42. Ángyán JG, Chipot C, Dehez F, Hättig C, Jansen G, Millot C. J. Comput. Chem. 2003;24:997–1008. [PubMed]
43. Dehez F, Ángyán JG, Soteras I, Luque FJ, Schulten K, Chipot C. J. Chem. Theory Comput. 2007;3:1914–1926. [PubMed]
44. Jeziorski B, Moszynski R, Szalewicz K. Chem. Rev. 1994;94:1887–1930.
45. Bukowski R, Cencek W, Jankowski P, Jeziorska M, Jeziorski B, Kucharski SA, Lotrich VF, Misquitta AJ, Moszyński R, Patkowski K, Podeszwa R, Rybak S, Szalewicz K, Williams HL, Wheatley RJ, Wormer PES, Zuchowski PS. SAPT2008: An ab initio program for many-body symmetry-adapted perturbation theory calculations of intermolecular interaction energies. Sequential and parallel Versions. Newark: Department of Physics and Astronomy, University of Delaware; DE 19716 and Department of Chemistry, University of Warsaw: ul. Pasteura 1, 02-093 Warsaw, Poland, 2008.
46. Møller C, Plesset MS. Phys. Rev. 1934;46:618–622.
47. Burcl R, Chalasiński G, Bukowski R. J. Chem. Phys. 1995;103:1498–1507.
48. Szalewicz K, Patkowski K, Jeziorski B. Structure & Bonding. 2005;116:43–117.
49. Torheyden M, Jansen G. Mol. Phys. 2006;104:2101–2138.
50. Misquitta AJ, Stone AJ. J. Chem. Theory Comput. 2008;4:7–18. [PubMed]
51. Berka K, Laskowski R, Riley KE, Hobza P, Vondrasek J. J. Chem. Theory Comput. 2009;5:982–992. [PubMed]
52. Chipot C, Ángyán JG. New J. Chem. 2005;29:411–420.
53. Ángyán JG, Chipot C. Int. J. Quant. Chem. 1994;52:17–37.
54. Vigné-Maeder F, Claverie P. J. Chem. Phys. 1988;88:4934–4948.
55. Stone AJ. In: Theoretical models of chemical bonding. Maksić H, editor. Vol. 4. Berlin: Springer–Verlag; 1991. pp. 103–131.
56. Stone AJ. The theory of intermolecular forces. Oxford: Clarendon Press; 1996.
57. Freitag MA, Gordon MS, Jensen JH, Stevens WJ. J. Chem. Phys. 2000;112:7300–7306.
58. Gordon MS, Freitag MA, Bandyopadhyay P, Jensen JH, Kairys V, Stevens WJ. J. Phys. Chem. A. 2001;105:293–307.
59. Piquemal J-P, Gresh N, Giessner-Prettre C. J. Phys. Chem. A. 2003;107:10353–10359. [PubMed]
60. Cisneros GA, Tholander SN-I, Parisel O, Darden TA, Elking D, Perera L, Piquemal J-P. Int. J. Quant. Chem. 2008;108:1905–1912. [PMC free article] [PubMed]
61. Applequist J, Carl JR, Fung K-K. J. Am. Chem. Soc. 1972;94:2952–2960.
62. Stone AJ, Alderton M. Mol. Phys. 1985;56:1047–1064.
63. Soteras I, Curutchet C, Bidon-Chanal A, Dehez F, Ángyán JG, Orozco M, Chipot C, Luque FJ. J. Chem. Theory Comput. 2007;3:1901–1913. [PubMed]
64. Jensen L, Åstrand P-O, Osted A, Kongsted J, Mikkelsen KV. J. Chem. Phys. 2002;116:4001–4010.
65. Tang KT, Toennies JP. J. Chem. Phys. 1984;80:3726–3741.
66. Tang KT, Toennies JP. J. Chem. Phys. 1977;66:1496–1506.
67. Millot C, Stone AJ. Mol. Phys. 1992;77:439–462.
68. Meredith AW, Stone AJ. J. Phys. Chem. A. 1998;102:434–445.
69. Jones J. Proc. R. Soc. London, Ser. A. 1924;106:441–462.
70. Ren P, Ponder JW. J. Phys. Chem. B. 2003;107:5933–5947.
71. Halgren TA. J. Am. Chem. soc. 1992;114:7827–7843.
72. Buckingham RA, Hamilton J, Massey HSW. Proc. R. Soc. A. 1941;179:103–122.
73. Frisch MJ, Trucks GW, Schlegel HB, Scuseria GE, Robb MA, Cheese-man JR, Montgomery JA, Jr, Vreven T, Kudin KN, Burant JC, Millam JM, Iyengar SS, Tomasi J, Barone V, Mennucci B, Cossi M, Scalmani G, Rega N, Petersson GA, Nakatsuji H, Hada M, Ehara M, Toyota K, Fukuda R, Hasegawa J, Ishida M, Nakajima T, Honda Y, Kitao O, Nakai H, Klene M, Li X, Knox JE, Hratchian HP, Cross JB, Bakken V, Adamo C, Jaramillo J, Gomperts R, Stratmann RE, Yazyev O, Austin AJ, Cammi R, Pomelli C, Ochterski JW, Ayala PY, Morokuma K, Voth GA, Salvador P, Dannenberg JJ, Zakrzewski VG, Dapprich S, Daniels AD, Strain MC, Farkas O, Malick DK, Rabuck AD, Raghavachari K, Foresman JB, Ortiz JV, Cui Q, Baboul AG, Clifford S, Cioslowski J, Stefanov BB, Liu G, Liashenko A, Piskorz P, Komaromi I, Martin RL, Fox DJ, Keith T, Al-Laham MA, Peng CY, Nanayakkara A, Challacombe M, Gill PMW, Johnson B, Chen W, Wong MW, Gonzalez C, Pople JA. Gaussian 03. Wallingford, CT: Gaussian, Inc.; 2004. Revision C.02.
74. Sadlej AJ. Theor. Chim. Acta. 1992;79:123–140.
75. Sadlej AJ. Theor. Chim. Acta. 1992;81:45–63.
76. Sadlej AJ. Theor. Chim. Acta. 1992;81:339–354.
77. Sadlej AJ. Collec. Czech. Chem. Commun. 1988;53:1995–2016.
78. Misquitta AJ, Stone AJ, Price SL. J. Chem. Theory Comput. 2008;4:19–32. [PubMed]
79. Saunders VR, Guest MF. ATMOL Program Package. Great Britain: SERC Daresbury Laboratory, Daresbury; 1976.
80. Jorgensen WL, Chandrasekhar J, Madura JD, Impey RW, Klein ML. J. Chem. Phys. 1983;79:926–935.
81. Colonna F, Angyan JG, Tapia O. Chem. Phys. Lett. 1990;172:55–61.
82. Dzyabchenko AV. Russ. J. Phys. Chem. A. 2008;82:758–766.
83. Coker H. J. Phys. Chem. 1976;80:2084–2091.
84. Fajans K. J. Phys. Chem. 1970;74:3407–3410.
85. Marquart D. SIAM J. Appl. Math. 1963;11:431–441.
86. Gresh N, Garmer DR. J. Comput. Chem. 1995;17:1481–1495.
87. Brdarski S, Karlström G. J. Phys. Chem. A. 1998;102:8182–8192.
88. Lago NF, Larrañaga FH, Albertí M. Eur. Phys. J. D. 2009 In press.
89. Mas EM, Szalewicz K. J. Chem. Phys. 1996;104:7606–7614.
90. Szalewicz K. Private communication
91. Chen W, Gordon MS. J. Phys. Chem. 1996;100:14316–14328.
92. Schmidt MW, Baldridge KK, Boatz JA, Elbert ST, Gordon MS, Jensen JH, Koseki S, Matsunaga N, Nguyen KA, Su S, Windus TL, Dupuis M, Montgomery JA. J. Comput. Chem. 1993;14:1347–1363.
93. Kim D, Tarakeshwar P, Kim KS. J. Phys. Chem. A. 2004;108:1250–1258.
94. Ma JC, Dougherty DA. Chem. Rev. 1997;97:1303–1324. [PubMed]
95. Swart M, van Duijnen PT, Snijders JG. J. Comput. Chem. 2000;22:79–88.