PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
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 2011 January 1.
Published in final edited form as:
J Chem Theory Comput. 2010; 6(1): 190–202.
doi:  10.1021/ct900348b
PMCID: PMC2832208
NIHMSID: NIHMS161432

Gaussian Multipole Model (GMM)

Abstract

An electrostatic model based on charge density is proposed as a model for future force fields. The model is composed of a nucleus and a single Slater-type contracted Gaussian multipole charge density on each atom. The Gaussian multipoles are fit to the electrostatic potential (ESP) calculated at the B3LYP/6-31G* and HF/aug-cc-pVTZ levels of theory and tested by comparing electrostatic dimer energies, inter-molecular density overlap integrals, and permanent molecular multipole moments with their respective ab initio values. For the case of water, the atomic Gaussian multipole moments Qlm are shown to be a smooth function of internal geometry (bond length and bond angle), which can be approximated by a truncated linear Taylor series. In addition, results are given when the Gaussian multipole charge density is applied to a model for exchange-repulsion energy based on the inter-molecular density overlap.

Keywords: Gaussian multipoles, charge density, electrostatic model, multipole, overlap

1. Introduction

Force fields are routinely used to simulate biological molecules in order to study structure and function. Recently, attention has been focused on developing accurate force field models that are able to provide a more realistic account of inter-molecular interactions. For example, polarization models16 have been incorporated into force fields716 in order to account for many body effects1720 in polar environments. In the SIBFA1012 force field, multipoles are placed on atoms and bond barycenters in order to accurately account for the anisotropy in electrostatic interactions. AMOEBA1316 places multipoles on atoms and has been developed as a force field for protein simulations. In addition, electrostatic models, which employ geometry dependent charges, are being investigated2124.

Atomic point multipole models derived from distributed multipole expansions2528 or fit to the electrostatic potential29 (ESP) have been proposed to improve the description of electrostatic interactions. However, at short range, it has been noted25,30,31 that atomic point multipole electrostatic models significantly underestimate electrostatic interactions at dimer distances. This effect, called penetration error, becomes important at dimer distances where there is significant overlap of molecular charge densities. Damping functions3134 have been proposed as a short range correction to atomic point multipoles in order to account for penetration effects.

Another approach for calculating short range electrostatic interactions has been to model the electron density. For example, simple Gaussian charge densities3537 have been used in models for liquids. Wheatley38,39 has studied Cartesian Gaussian multipole charge distributions40,41, which are obtained by differentiating simple normalized Gaussian functions. It was shown that at long range, the interactions between Gaussian multipole charge densities behave asymptotically as point multipoles38. As an example, consider a Cartesian Gaussian dipole charge distribution ρμ with dipole moment μ, exponent α, and nuclear center R given by

ρμ(r)μR(α2π)32exp(α2rR2),
(1)

where r is the field point and [big down triangle, open]R is the Cartesian gradient with respect to the nuclear center R. The electrostatic potential ϕμ(r) arising from ρμ is given by

ϕμ(r)=d3rρμ(r)rr=μRerf(αrR)rR,
(2)

where erf(x) is the error function defined by

erf(x)2π0xexp(u2)du.
(3)

Note that for large x, erf(x) ≈ 1. Thus, for large exponents or large separations, the electrostatic potential arising from a Gaussian dipole behaves as a point dipole, i.e.

ϕμ(r)μR1rRαrR>>1.
(4)

Recently, Giese and York42 have derived the electrostatic energy integral, density overlap integral, and gradients of their matrix elements between contracted solid harmonic Gaussian (multipole) charge densities in the spherical tensor framework. In analogy to taking derivatives of simple normalized Gaussian functions to create Cartesian Gaussian multipoles, spherical tensor Gaussian multipoles can be constructed by contracting spherical tensor multipole moments Qlm with the solid harmonic gradient operator43,44 Clm ([big down triangle, open]) acting upon a simple Gaussian function. In the Supplementary Information (SI), a summary of the mathematical background4347 used in this work, including definitions and theorems for the solid harmonic function Clm(x,y,z) and the solid harmonic gradient operator Clm ([big down triangle, open]), is provided.

Methods of calculating the Gaussian charge density parameters are also being explored. Previously, we have proposed the Gaussian Electrostatic Model (GEM)10,4851 as an electrostatic model based on Gaussian charge density. The GEM density ρ(r) is represented by a linear expansion of conventional auxiliary Gaussian basis sets. The density coefficients are fit to the ab initio density ρQM(r) through a least-squares fit to the error in self-interaction electrostatic energy5256 ΔEself given by

ΔEself=ρQM(r)ρ(r)1rρQM(r)ρ(r).
(5)

GEM has been shown to accurately reproduce inter-molecular electrostatic interaction energies at both short and long range distances. In addition, it was also shown that a GEM type model can be constructed by fitting to the electrostatic potential (ESP) numerically on a grid57. When fitting to ESP, it was found that fewer Gaussian basis functions were needed to reproduce ab initio inter-molecular electrostatic energies as compared to fitting to the error in self-interaction electrostatic energy ΔEself. In that work, we have also given some preliminary results for fitting a single s-type Gaussian charge function on each atom to the ESP. It was found that a simple model consisting of a single negative Gaussian charge distribution and positive nucleus on each atom is able to account for a large fraction of the penetration error and represents a significant improvement over atomic point charge models for short range electrostatic interactions.

In the present study, we propose a natural extension to the single (uncontracted) s-type Gaussian charge model by generalizing to a model based on a single diffuse contracted Gaussian multipole charge density on each atom. The radial part of the Gaussian multipole charge density is a Slater-type58,59 contracted Gaussian function fit to exp(-λr). For each atom, the Gaussian multipole moments Qlm and a single Slater-type exponent parameter λ are fit to the electrostatic potential (ESP) surrounding the molecule calculated at the B3LYP/6-31G* or HF/aug-cc-pVTZ levels of theory. The model is tested by comparing electrostatic dimer energies, inter-molecular density overlap integrals, and permanent molecular moments with their reference ab initio values. In addition, a significantly improved method of fitting the Gaussian multipoles to the ESP numerically on a grid is adopted. In particular, angular grids available in most quantum chemistry codes are used, and a smooth weighting function, similar to the one proposed by Hu60 et al., is used to filter out points near the nuclear centers.

In this work, we will show that a single diffuse contracted Gaussian multipole (shell) on each atom is capable of reproducing ab initio inter-molecular electrostatic energies and density overlap integrals on hydrogen bonded dimers at equilibrium geometries. For inter-molecular electrostatic energies, the accuracy attained by Gaussian multipoles is shown to be approximately 0.1 kcal/mol, which is comparable to that of our original GEM model. Since the GEM model is represented by an auxiliary Gaussian basis set consisting of many uncontracted Gaussian shells (multipoles) on each atom, the proposed single Gaussian multipole model is expected to be an efficient model suitable for developing a force field for molecular dynamics simulations. In addition, since there are fewer fitted parameters used in Gaussian multipoles, the fit is more stable. Based on this property, the atomic Gaussian multipoles Qlm are shown to be a continuous function of internal geometry (bond lengths, bond angles, etc). More specifically, it is shown that the atomic Gaussian multipole moments Qlm can be approximated as a truncated linear Taylor series in bond length and bond angle for the case of water.

The inter-molecular density overlap integral calculated by Gaussian multipoles can be applied to the exchange-overlap model proposed by Wheatley and Price61,62 for the exchange-repulsion energy. The first order inter-molecular exchange-repulsion energy Eexch is defined as the total ab initio dimer energy minus the inter-molecular electrostatic energy calculated using the frozen monomer wave functions. The ab initio exchange-repulsion energy Eexch can be modeled by fitting a proportionality constant K to the inter-molecular density overlap integral S by

EexchKSa,
(6)

where a is an empirical exponent parameter used to improved the quality of fit62. The expression in eqn. 6 is commonly generalized to a pair-wise sum over atom-atom61,62 contributions between the two monomers. One of the challenges of developing parameters and applying the exchange-overlap model is first finding an accurate and convenient representation of the molecular charge density. In a previous work48,50, we have applied the exchange-overlap model in eqn. 6 to the GEM molecular charge density. A molecular pair K parameter was fit to the ab initio exchange-repulsion energy for the water-water dimer over several randomly oriented water-water geometries. Although fitting molecular pair K parameters over atomic pair parameters may not be ideal for constructing a general force field, we are interested in studying the effects of applying an anisotropic charge density to the exchange-overlap model. In the present study, a single molecular pair K parameter (eqn. 6) is fit to the ab initio exchange energies using the inter-molecular density overlap integrals calculated by atomic Gaussian monopoles, dipoles, and quadrupoles for small molecule hydrogen bonded dimers. Including anisotropy in the description of charge density is shown to make a significant improvement in reproducing ab initio exchange-repulsion energies.

In the following section, a Gaussian multipole charge density is defined and expressions for the electrostatic potential, electrostatic energy, and density overlap integral are given. Details on fitting Gaussian multipoles to the electrostatic potential are provided along with a discussion on how the model is tested with ab initio inter-molecular electrostatic energies, inter-molecular density overlap integrals, and permanent molecular multipole moments. This is followed by a brief discussion on how the inter-molecular density overlap integrals calculated by Gaussian multipoles are applied to the exchange-overlap model. In the next section, results for Gaussian multipoles are presented. Inter-molecular electrostatic energies, inter-molecular density overlap integrals, and permanent molecular multipole moments calculated by Gaussian multipoles are compared with their respective ab initio values. The geometry dependence of atomic Gaussian multipole moments Qlm is presented for the case of water. Results are given when the Gaussian multipole inter-molecular density overlap integrals are applied to the exchange-overlap model. Lastly, the results are summarized and future applications are discussed in the Conclusions.

2. Methods

In this section, a definition of a contracted Gaussian multipole charge density is given along with expressions for the electrostatic potential, electrostatic energy, and density overlap integrals. This is followed by a brief discussion of Slater-type contracted Gaussian functions and a description of how the Gaussian multipoles are fit to the ESP. In the next sub-section, computational details on the calculation of ab initio electrostatic energies and inter-molecular density overlap integrals are discussed. The calculation of permanent molecular multipole moments are described. This section concludes with computational details of how the exchange-overlap parameters are fit.

2.1 Gaussian Multipoles

The definition we have used for a contracted Gaussian multipole charge density is similar to the one given by Giese and York42. In that work42, expressions for efficiently calculating electrostatic and density overlap matrix element integrals between real regular solid harmonic contracted Gaussian functions are given along with gradients of their matrix elements. In order to test the model for Gaussian multipole charge density, we have implemented similar expressions using complex solid harmonic Gaussian (multipole) functions, which are given below. A derivation of the following results along with a summary of necessary mathematical background4347 is given in the SI for the interested readers. Many of the theorems quoted in the SI have been used in the evaluation of integrals between solid harmonic Gaussian basis functions6365 for quantum chemistry calculations.

In the Introduction, it was mentioned that a spherical tensor Gaussian multipole charge distribution can be defined in terms of multipole moments Qlm and the solid harmonic gradient operator Clm ([big down triangle, open]) acting upon a normalized Gaussian function. As shown in the SI, the solid harmonic gradient operator is especially useful in deriving integral quantities such as electrostatic energies and density overlap integrals. In the following discussion, the final results for electrostatic energy and overlap integral between two Gaussian multipole charge densities are given after all derivative operations have been evaluated.

A contracted Gaussian multipole charge density ρ(r,R) with moments Qlm and nuclear center R evaluated at the point r is given by

ρ(r,R)l=0lmaxmlQlmClm(rR)(2l1)!!ρl(rR;αμ),
(7)

where lmax is the maximum order of the Gaussian multipoles (e.g. lmax = 0 for monopoles, 1 for dipoles, etc.), Clm is the complex conjugate of a solid harmonic function, and ρl is a derivative of a contracted Gaussian charge density defined by

ρl(r;aμ)(1rddr)lμ=1Nccμ(αμ2π)32exp(αμ2r2),
(8)

where Nc is the degree of contraction. For l = 0, the density ρ0 is normalized to unity, (Σμcμ=1). multipole moments of the charge density ρ(r,R) with respect to the center R are the coefficients Qlm, i.e.

d3rρ(r,R)Clm(rR)=Qlm.
(9)

The electrostatic potential ϕ arising from ρ in eqn. 7 is given by

ϕ(r)=d3rρ(r,R)rr=l=0lmaxmlQlmClm(rR)(2l1)!!ϕl(rR;αμ),
(10)

where ϕl is defined by

ϕl(r;αμ)(1rddr)lμ=1Nccμerf(αμr)r,
(11)

and erf(x) is the error function defined in eqn. 3.

The electrostatic interaction energy between two Gaussian multipole charge densities ρ1(r,R1) and ρ2(r,R2) and their nuclei Z1 and Z2 is given by

U=d3rd3rρ1(r,R1)ρ2(r,R2)rr+Z1ϕ2(R1)+Z2ϕ1(R2)+Z1Z2R1R2,
(12)

where ϕ2(R1) is the potential at R1 due to ρ2 and ϕ1(R2) is defined by a similar expression. The density overlap integral S between two Gaussian multipole charge densities is defined by

S=d3rρ1(r)ρ2(r).
(13)

The electrostatic and density overlap integrals in eqns. 12 and 13 can be expressed by42

d3rd3rρ1(r)O(r,r)ρ2(r)=l1,m1l2,m2Ql1m11Tl1m1;l2m2Ql2m22,
(14)

where O(r,r)1rr for the electrostatic integral, O(r,r)δ(rr) for the density overlap integral, and the interaction matrix Tl1m1;l2m2 is given below. The electrostatic energy and density overlap integrals between two normalized contracted Gaussian monopole charge densities with unit charge are given by

F0(R)μ,vcμ1cv2erf(αμvR)R,F0(R)μ,vcμ1cv2(παμv2)32exp(αμv2R2),
(15)

where R is the distance between nuclear centers and αμv is the Gaussian product exponent defined by

RR1R2αμvαμαvαμ2+αv2.
(16)

In addition, the constants Alm and Blm are defined by

Alm(l+m)!(lm)!,BlmAlm(2l1)!!,
(17)

for l ≠ 0, and A00 = B00 = 1. The scaled solid harmonic function is defined by Rlm(r)Clm(rAlm). The interaction matrix Tl1m1;l2m2 from eqn. 14 is given by

Tl1m1;l2m2=Bl1m1l=0min(l2,l2)m=1l(1)l2+mAlmBlmRl2l,m2+m(R)Rl1l,m2m(R)Fl2+l11(R),
(18)

where Fl(R) [equivalent] 2l (d/dR2)1 F0(R). Finally, the point multipole results for electrostatic potential and electrostatic energy can be found by considering an uncontracted Gaussian multipole (Nc = 1), taking the large exponent limit (α→∞), and noting erf (x) → 1 as x → ∞.

The atomic Gaussian multipole moments Qlm are commonly defined25 with respect to a local frame of the atom Qlmlocal and then rotated to the system or global frame QlmlocalQlm through Wigner rotation matrices Dmml

Qlmglobal=mDmml[R1]Qlmlocal.
(19)

Recursion formulae for evaluating Dmml have been given in Choi66 et al. For each atom, the Cartesian transformation matrix R between the local and global frames is defined with respect to the relative positions of the atom and its neighbors13,67,68.

2.2 Slater-Type Contracted Gaussian functions

The contraction coefficients dμ and exponents αμ are fit to a simple Slater function exp(-r) over all space for Nc = 1 − 14

exp(r)μ=1Ncdμexp(αμ2r2).
(20)

For Nc = 1 − 6, the optimized exponents agree with those used in the development of the original STONG basis sets58,59. The contraction coefficients dμ and exponents αμ fit to exp(-r) are used to find the corresponding coefficients and exponents for exp(-λr) by a scaling argument. The final expression for a normalized Slater-type charge density with unit charge and exponent λ is given by

λ38πexp(λr)μ=1Nccμ(αμ2λ2π)32exp(αμ2λ2r2),
(21)

where cμdμ8π(παμ2)32. A full list of contracted coefficients and exponents for exp(-r) can be found in the SI for Nc = 1 − 14.

2.3 Non-linear Fit to Potential

The model for molecular charge density presented in this work consists of an effective nuclear charge Zeff and a set of contracted Gaussian multipole moments Qlm with a single diffuse Slater exponential parameter λ centered on each atom. Only the valence charge density is modeled, and the core electron density near the nuclear centers is neglected by using screened nuclear charges Zeff = Z − Ncore, where Z is the true nuclear charge and Ncore is the number of core electrons. The number of core electrons Ncore is taken to be 0 for hydrogen, 2 for the first row elements, and 10 for the second row elements. Thus, the screened nuclear charges Zeff are set to 1.0 for H, 4.0 for C, 5.0 for N, 6.0 for O, 7.0 for F, and 7.0 for Cl. Initially, we experimented with using the true nuclear charges, e.g. Z = 8 for O. However, when the true nuclear charges are used, the inter-molecular density overlap integrals at equilibrium dimer distances are consistently overestimated by 10–15%. If the true nuclear charges are used, the model is forced to account for both the core and valence electron density with a single diffuse Slater-type Gaussian function. By only modeling the valence charge density, we had found that using effective screened nuclear charges gave significantly smaller errors when comparing to the ab initio inter-molecular electrostatic energy and density overlap integral.

The model for Gaussian multipole molecular charge density ρGM evaluated at the point r is a represented as a sum over atomic Gaussian multipole charge densities given by

ρGM(r)=al=0lmaxmlQlm,aClm(rRa)(2l1)!!ρl(rRa;λaαμ)
(22)

where Ra is the nuclear center of atom a, Qlm,a are the atomic Gaussian multipole moments of atom a, λa is the Slater exponent on atom a, and ρl is defined by eqn. 8. The electrostatic potential due to the effective nuclear charges and Gaussian multipole charge density can be found from eqn. 10 as

ϕ(r;Qlm,a,λa)=aZeff,arRa+l=0lmaxmlQlm,aClm(rRa)(2l1)!!ϕi(rRa,λaαμ),
(23)

where r is the field point, Zeff,a is the effective nuclear charge of atom a, and ϕl is defined by eqn. 11. For each atom, Qlm,a and λa are treated as optimizable parameters and fit to the ab initio electrostatic potential ϕQM surrounding the molecule. Gaussian quadrupoles are defined by lmax = 2, i.e. a Gaussian monopole, dipole, and quadrupole with the same atomic exponent parameter λa are placed on a given atom a. Similarly, Gaussian dipoles are defined by lmax = 1, i.e. a Gaussian monopole and dipole with the same atomic exponent parameter λa are placed on a given atom a.

Points near the nuclei are filtered out or discarded by using a weighting function w(r). The weighting function w(r) used in this work is a modified version of a weighting function taken from Hu et al.60 In this present study, w(r) is a sigmoid function of the form

W(R){exp{σ[lnρQM(r)lnK0]2}ρQM(r)K01ρQM(r)K0}
(24)

where ρQM is the ab initio electron density. The weighting function w(r) is small for regions of high electron density, while w(r)=1 for regions of low electron density. The adjustable parameters σ and lnK0 control the curvature and location of the sigmoid function, respectively. In the limit of large σ, w(r) becomes a step function. For Gaussian multipoles, the σ and lnK0 parameters are set to 0.3 and −6.0. The σ and lnK0 parameters are selected by performing a 2D scan of parameters and observing the average error in electrostatic dimer energy of hydrogen bonded dimers at equilibrium geometries. A plot of the average RMSD error in electrostatic dimer energy is given in Fig. 1 for the case of contracted (Nc = 4) Gaussian quadrupoles fit to the ESP calculated at the B3LYP/6-31G* level. The surface is flat and a range of values perform equally well, e.g. (σ, lnK0) = (0.2, −7), (0.3, −6), (0.4, −5), etc. Similar parameter scans were performed with Gaussian monopoles, dipoles, and quadrupoles at various degrees of contraction Nc, and the same set of parameters were found to be local minima. In Hu et al.60, a weighting function for fitting atomic point charges to the ab initio ESP is proposed. There, it was shown that the point charges are stable with respect to varying conformations. The main differences between the weighting function ww(r) used in this work and the one proposed by Hu60 et al. are that the ab initio electron density ρQM is used, rather than an empirical model for ρQM, and points far away from the molecule are given a weight of 1.0. We have also fit atomic point multipoles to the ESP using the σ and lnK0 weight parameters of 0.8 and −9.0 given by Hu60 et al.

Figure 1
Contour plot of the average RMSD error in electrostatic dimer energy (kcal/mol) as a function of parameters σ and lnK0 for the ESP weighting function. The electrostatic energies are calculated on hydrogen bonded dimers at equilibrium geometries ...

The fitting function χ2 is given by

χ2(Qlm,λ)=d3rw(r)[ϕGM(Qlm,λ;r)ϕQM(r)]2,
(25)

where ϕGM and ϕQM are the electrostatic potentials (ESP) calculated by Gaussian multipoles (eqn. 23) and ab initio, respectively. χ2 is approximated numerically on a coarse grain molecular grid taken from a modified version of NWChem69,70 and optimized using a Levenberg-Marquardt non-linear least-squares fit algorithm71. The ab initio electrostatic potential is calculated at both the B3LYP/6-31G* or HF/aug-cc-pVTZ levels using the Gaussian 03 software package72. Earlier in our study, we had experimented with using rectangular grids similar to the CHELPG73 type grids used in optimizing atomic point charges. A relatively fine grid spacing of 0.05 Å was used and points within 1 Å of any nuclei were discarded. For uncontracted Gaussian multipole (Nc = 1), this procedure gave Gaussian multipole parameters which predicted electrostatic energies in approximate agreement to the results presented here using the molecular grids with a smooth weighting function. However, for water, the number of rectangular grid points needed was on the order of 106 – 107, which can be compared to 103 – 104 grid points used in the coarse grain molecular grids.

2.4 Ab initio Dimer Energy and Molecular Density Overlap Test

The model is tested by comparing inter-molecular electrostatic dimer energies and density overlap integrals on equilibrium dimer geometries of various molecules hydrogen bonded to water. The geometries of the dimers are optimized at both the B3LYP/6-31G* or the HF/aug-cc-pVTZ levels, while keeping the monomers rigid in their respective monomer-optimized geometries. For the model fit to B3LYP/6-31G* data, the model is tested by comparing to ab initio electrostatic energies calculated by the Constrained Space Orbital Variation (CSOV) decomposition74 method using a modified version of the HONDO75,76 quantum chemistry program. For the model fit to HF/aug-cc-pVTZ data, ab initio electrostatic energies are calculated by the Reduced Variational Space77,78 (RVS) decomposition method using the GAMESS79 quantum chemistry program. In addition, we have developed code to calculate ab initio inter-molecular density overlap integrals from the ab initio density matrix using the McMurchie-Davidson algorithm80.

For the water-water dimer, the model is tested by calculating inter-molecular electrostatic energies and density overlap integrals on non-equilibrium dimer geometries. Several water-water dimer geometries are generated by rigidly translating one water molecule with respect to the other in the direction of the inter-molecular H..O hydrogen bond. The inter-molecular electrostatic energies and density overlap integrals calculated by Gaussian multipoles are plotted as function of H..O distance and compared to their ab initio values. In addition, 100 water-water dimer geometries are generated in random orientations, while the relative center of masses lie between 2.5 and 5.0 Å. Scatter plots of the inter-molecular electrostatic energy and density overlap integral are presented comparing the results calculated by Gaussian quadrupoles with their ab initio values.

2.5 Molecular Multipole Moments

The Gaussian multipoles are further tested by comparing permanent molecular dipoles (l = 1), quadrupoles (l = 2), octapoles (l = 3), and hexadecapoles (l = 4) with their ab initio molecular multipoles. The atomic Gaussian multipole moments Qlm at position R are translated to the origin by the following expression25

Qlmorigin=l1=0lmaxm1=11l1(l+ml1+m1)(lml1m1)Ql1m1Cll1,mm1(R),
(26)

where Clm is a solid harmonic function. Note that Qlmorigin0 for all values of l and m. For example, an atomic Gaussian dipole (with respect to its atomic position) contributes to the total molecular dipole, quadrupole, octapole,.. (with respect to the origin). The translated atomic Gaussian multipoles at the origin Qlmorigin are summed to give the total spherical tensor molecular moment. Expressions for converting real spherical tensor multipoles into traceless Cartesian multipoles can be found in Özdoğan81. Below, the results for converting complex spherical tensor multipoles Qlm [equivalent] Qrlm + iQilm into their traceless Cartesian forms are given. For l = 1, the Cartesian dipole μα is related to Q1m by

μx=2Q11r,μy=2Q11i,μz=Q10.
(27)

For l = 2, the traceless Cartesian quadrupoles ΘαβTL are related to Q2m by

ΘxxTL=12(Q20+6Q22r)ΘxzTL=32Q21rΘyyTL=12(Q206Q22r)ΘyzTL=32Q21iΘzzTL=Q20ΘxyTL=32Q22r
(28)

Note the trace of the quadrupoles is zero, i.e. Tr(ΘTL)=ΘxxTL+ΘyyTL+ΘzzTL=0. For the conversion formulae of traceless Cartesian octapoles (l = 3) and hexadecapoles (l = 4) from their complex spherical tensor moments, see the SI. The ab initio molecular multipole moments are calculated by Gaussian 0372 and then converted to their traceless forms25. For example, the expression for converting Cartesian quadrupoles Θαβ into traceless Cartesian quadrupoles ΘαβTL is given by

ΘαβTL=32Θαβ12δαβpΘpp.
(29)

Similar expressions for converting Cartesian octapoles and hexadecapoles into their traceless forms are given in the SI.

2.6 Gaussian Multipole Geometry Dependence

The atomic Gaussian multipole moments Qlm in the local frame (eqn. 19) are calculated as a function of bond length and bond angle for a water molecule. The geometry of water is optimized at the B3LYP/6-31G* level. Several geometries of water are found by performing two separate one dimensional scans of perturbing one of the bond lengths and the bond angle away from equilibrium in increments of 0.1 Å and 1°, respectively. Atomic Gaussian quadrupoles Qlm and exponent parameters λ (Nc = 4) are fit to the B3LYP/6-31G* ESP for the optimized water geometry. For each perturbed geometry, new atomic Gaussian quadrupoles are fit the B3LYP/6-31G* ESP calculated for that geometry, while the exponent parameters λ are kept at their geometry optimized values. The atomic Gaussian multipoles in the local frame Qlmlocal are plotted as a function bond length and bond angle.

2.7 Overlap-Exchange Model

The exchange-overlap model61,62 is tested using the model for Gaussian multipole charge density for hydrogen bonded dimer pairs. For a given dimer, several geometries are generated by translating one of the monomers in increments of 0.1 Å along the axis defined by the two atoms forming the hydrogen bond. For each dimer geometry, the B3LYP/6-31G* exchange energy Eexch is calculated through CSOV decomposition. Only dimer geometries for which the total dimer energy is within +5 kcal/mol of the total minimum dimer energy are kept. Typically, this entails 40 – 60 dimer geometries, of which the exchange energy lies between 0 and 30 kcal/mol. For example, the exchange overlap model for the water-water dimer is fit to 63 geometries, whose exchange energies lie between 0 and 27 kcal/mol. For each dimer pair, a single molecular pair K parameter is fit to the ab initio exchange energy (eqn. 6) using the inter-molecular density overlap integrals S calculated by Gaussian multipoles with a = 0.95.

3. Results

In this section, results for Gaussian monopoles, dipoles, and quadrupoles are presented. Recall that `Gaussian quadrupoles' refers to a model in which Gaussian monopoles, dipoles, and quadrupoles with the same atomic exponent parameter λ are placed on each atom. In the next few sub-sections, results for inter-molecular electrostatic energy, density overlap integral, and permanent molecular multipole moment calculated by Gaussian multipoles are compared with their respective ab initio values. For the case of water, the atomic local frame Gaussian multipole moments Qlmlocal are plotted as a function of bond length r and bond angle θ. Lastly, results are given when Gaussian multipoles are applied the exchange-overlap model.

The dependence of inter-molecular electrostatic energy and density overlap integral on the degree of Slater-type contraction Nc is studied. For Nc = 1, the model for charge density is a single Gaussian function. In the limit of large Nc, the model for charge density is equivalent to using a Slater function exp(−λr). For Gaussian multipoles fit to the B3LYP/6-31G* ESP, the results for energy and density overlap favored Nc = 4. However, for the Gaussian multipoles fit to the HF/aug-cc-pVTZ ESP, the errors in inter-molecular electrostatic energy and density overlap decreased for larger values of Nc. The results indicate the optimal degree of contraction Nc depends on the size of the ab initio basis set. The smaller 6-31G* basis set favors a smaller degree of contraction, while the larger aug-cc-pVTZ basis set prefers a larger degree of contraction. The values Nc = 4 and Nc = 8 are chosen for the Gaussian multipoles fit to the B3LYP/6-31G* ESP and the HF/aug-cc-pVTZ ESP, respectively. For more details on the Nc dependence, see the SI.

3.1 Electrostatic Energy

Electrostatic dimer energies for several molecules hydrogen bonded to water are calculated at their equilibrium geometries. In Table I, electrostatic dimer energies are given for Gaussian monopoles, dipoles, and quadrupoles with Nc = 4, which are fit to the ESP calculated at the B3LYP/6-31G* level. The Gaussian multipole electrostatic dimer energies are compared with their reference B3LYP/6-31G* values. The root mean square deviation (RMSD) errors in electrostatic dimer energy are 0.568, 0.567, and 0.094 kcal/mol for Gaussian monopoles, dipoles, and quadrupoles, respectively. On average, the errors for Gaussian monopoles and dipoles (Nc = 4) are quite similar, while a significant improvement is gained for Gaussian quadrupoles. As a representative example, the electrostatic dimer energies of the water–methanol(1) dimer are −8.751, −8.103, and −8.558 kcal/mol for Gaussian monopoles, dipoles, and quadrupoles, respectively. These numbers can be compared to the reference ab initio electrostatic energy for the water–methanol(1) dimer of −8.524 kcal/mol. Gaussian quadrupoles are found to be particularly important in predicting the electrostatic dimer energies of organic halides with water. For example, the ab initio electrostatic energy for the water–CH3Cl(1) dimer is −0.372 kcal/mol. This result can be compared to the electrostatic dimer energies predicted by Gaussian monopoles, dipoles, and quadrupoles of −1.070, −0.719, and −0.362 kcal/mol, respectively.

Table I
Electrostatic energies (kcal/mol) for equilibrium hydrogen bonded dimers (X-water). The electrostatic energies predicted by Gaussian monopoles EM, Gaussian dipoles EDM, and Gaussian quadrupoles EQDM (Nc = 4) are compared to their reference B3LYP/6-31G* ...

A similar analysis is performed at the HF/aug-cc-pVTZ level. Gaussian multipoles (Nc = 8) are fit to the ESP calculated at the HF/aug-cc-pVTZ level, while the reference ab initio electrostatic energies are calculated at the same level of theory using the RVS decomposition method. Due to computational limitations, 11 of the original 25 hydrogen bonded dimers are studied at this level. The dimers with the smallest monomers are chosen (water, ammonia, methanol, CH3F, and CH2F2). As expected, a significant improvement is found by increasing the multipole order from Gaussian monopoles to Gaussian quadrupoles. The RMSD errors in electrostatic dimer energy are 0.885, 0.366, and 0.133 kcal/mol for Gaussian monopoles, dipoles, and quadrupoles, respectively. As an example, the electrostatic energies for the water–methanol(1) dimer predicted by Gaussian monopoles, dipoles, and quadrupoles (Nc = 8) are −9.233, −8.187, and −8.847 kcal/mol, respectively. These results can be compared to the HF/aug-cc-pVTZ electrostatic energy of −8.753 kcal/mol. For more individual results, see the SI.

The results for inter-molecular electrostatic energy given above are calculated on equilibrium dimer geometries. For non-equilibrium dimer geometries, the inter-molecular electrostatic energy are calculated for the water-water dimer. The electrostatic energies calculated by Gaussian quadrupoles (Nc = 4) and ESP fitted atomic point quadrupoles are plotted for various hydrogen bond distances H..O in Fig. 2 for the water-water dimer and compared to their reference B3LYP/6-31G* values. As in the case of Gaussian quadrupoles, we call `atomic point quadrupoles' as a model in which ESP fitted atomic point monopoles, dipoles, and quadrupoles are placed on each atom. The optimized equilibrium dimer H..O distance is found to be 1.94 Å at the B3LYP/6-31G* level. At the equilibrium dimer separation, Gaussian quadrupoles predict an electrostatic energy of −8.195 kcal/mol, atomic point quadrupoles predict −6.422 kcal/mol, while the reference ab initio result is −8.235 kcal/mol. The under-estimation of the atomic point quadrupole energy is an example of the penetration error25,30,31 for atomic point multipoles. At long range H..O distances, beyond 2.3 Å, both Gaussian quadrupoles and atomic point quadrupoles accurately reproduce the ab initio electrostatic energy. At short range H..O distances less than 1.64 Å, the Gaussian quadrupole electrostatic energy begins to slowly deviate from its reference B3LYP/6-31G* result. For example at 1.54 Å, the B3LYP/6-31G* ab initio electrostatic energy is −20.11 kcal/mol, which can be compared to the result of −21.54 kcal/mol calculated by Gaussian quadrupoles and −11.74 kcal/mol calculated by atomic point quadrupoles. In Fig. 3, the inter-molecular electrostatic energies calculated by Gaussian quadrupoles are compared with their ab initio reference values for several randomly oriented water-water dimers. The electrostatic energies calculated by Gaussian quadrupoles agree with their ab initio reference values for energies ranging from −10 kcal/mol to +5 kcal/mol. Additional scatter plots of inter-molecular electrostatic energy can be found in the SI for randomly oriented hydrogen bonded dimers.

Figure 2
The electrostatic energy (kcal/mol) calculated by atomic point quadrupoles, Gaussian quadrupoles, and ab initio is plotted as a function of H..O distance for the water – water dimer.
Figure 3
The inter-molecular electrostatic energy (kcal/mol) calculated by Gaussian quadrupoles (y-axis) and ab initio (x-axis) is plotted for randomly oriented water – water dimer geometries.

In order to compare Gaussian multipoles with the Gaussian Electrostatic Model (GEM)57, electrostatic dimer energies are calculated on 10 water dimers82,83 which represent local minima on the water-water potential energy surface. Previously, GEM57 was fit to the B3LYP/6-31G* ESP using the A1 and P1 auxiliary Gaussian basis sets (ABS). Two GEM models for water were developed. In a 3 pt. GEM water model, ABS's have been placed on atomic centers only. A second GEM water model has been developed by placing ABS's on both the atomic centers and bond midpoints resulting in a 5 pt. GEM water model. The average absolute errors in the electrostatic energy for the 3 pt. GEM water model are 0.11 and 0.12 kcal/mol for the A1 and P1 basis sets, respectively. Including basis functions on bond midpoints results in an improved fit for GEM. The errors in the 5 pt. GEM water fit to the ESP are 0.06 and 0.04 kcal/mol for the A1 and P1 basis sets, respectively. The average errors in electrostatic dimer energy for GEM can be compared to the error of 0.06 kcal/mol for Gaussian quadrupoles (Nc = 4). The GEM 5 pt. water model with the A1 and P1 basis sets have a total of 110 and 213 primitive Gaussian functions, respectively. The total number of uncontracted basis functions in the GEM water models can be compared to the total number of contracted Gaussian quadrupoles of 27.

3.2 Molecular Density Overlap Integral

The model for Gaussian multipoles is further tested by comparing inter-molecular density overlap integrals with their ab initio values. In Table II, the inter-molecular density overlap integrals are given for the hydrogen bonded dimers at their equilibrium geometries. The inter-molecular density overlap integrals for Gaussian monopoles, dipoles, and quadrupoles (Nc = 4) are compared with their respective B3LYP/6-31G* values. On average, there is a significant improvement in going up in multipole order from Gaussian monopoles to Gaussian quadrupoles. The RMSD errors in inter-molecular density overlap integrals are 2.571, 0.752, and 0.195 10−3 e23 for Gaussian monopoles, dipoles, and quadrupoles, respectively. A representative example is the water-methylamine complex at equilibrium, which has an ab initio value of 15.89 10−3 e23. This value can be compared to the inter-molecular density overlap calculated by Gaussian monopoles, dipoles, and quadrupoles of 19.78, 17.05, and 15.56 10−3 e23, respectively. Similar trends can be found for the Gaussian multipoles fit to the HF/aug-cc-pVTZ ESP. At the HF/aug-cc-pVTZ level, the RMSD errors in inter-molecular density overlap integral are 2.668, 0.502, and 0.268 10−3 e23 for Gaussian monopoles, dipoles, and quadrupoles (Nc = 8), respectively. For individual results on HF/aug-cc-pVTZ Gaussian multipoles, see the SI.

Table II
Inter-molecular density overlap integrals (multiply by 10−3 e23) for equilibrium hydrogen bonded dimers (X-water). The overlap integrals predicted by Gaussian monopoles SM, Gaussian dipoles SDM, and Gaussian quadrupoles SQDM (Nc = 4) ...

The results given above are for inter-molecular density overlap integrals calculated at equilibrium dimer distances. For non-equilibrium dimer geometries, the inter-molecular density overlap integrals are compared with their ab initio values for the water-water dimer. In Fig. 4, the inter-molecular density overlap integrals calculated by B3LYP/6-31G* Gaussian quadrupoles (Nc = 4) and the B3LYP/6-31G* reference values are plotted as a function of H..O distance. The inter-molecular overlap integrals predicted by Gaussian quadrupoles agree with their ab initio values for H..O distances ranging from 2.3 Å to 1.7 Å. For shorter separations, the inter-molecular overlap integrals predicted by Gaussian quadrupoles begin to overestimate the ab initio result. In Fig. 5, the inter-molecular density overlap integral calculated by Gaussian quadrupoles are compared with their respective ab initio values for several randomly oriented water-water geometries. The inter-molecular overlaps integrals range from 0 to 30.0 10−3 e23. There is a small over-estimation (~6 %) of the inter-molecular density overlap integrals calculated by Gaussian quadrupoles as compared with their ab initio values. In the SI, additional scatter plots of inter-molecular density overlap integral can be found for other randomly oriented hydrogen bonded dimers.

Figure 4
The inter-molecular density overlap integral (10−3 e23) calculated by atomic Gaussian quadrupoles and ab initio is plotted as a function of H..O distance (Å) for the water – water dimer.
Figure 5
The inter-molecular density overlap integral (10−3 e23) calculated by Gaussian quadrupoles (y-axis) and ab initio (x-axis) is plotted for randomly oriented water – water dimer geometries.

3.3 Molecular Multipole Moments

The permanent molecular multipole moments up to hexadecapole are calculated for the atomic Gaussian monopole, dipole, and quadrupole models (Nc = 4) and compared with their reference B3LYP/6-31G* values. As an example, the non-zero components of the molecular quadrupole for ammonia are given in Table III. A significant improvement is found by increasing the atomic Gaussian multipole order from Gaussian monopoles to Gaussian quadrupoles. For example, the ab initio Qxx component of the molecular quadrupole is 1.3140 D-Å, which can be compared to the value predicted by Gaussian monopoles, dipoles, and quadrupoles of 0.9211, 1.2701, and 1.3139 D-Å, respectively.

Table III
Non-zero components of the molecular quadrupole moment D-Å of ammonia for Gaussian monopoles, dipoles, and quadrupoles (Nc = 4) fit to B3LYP/6-31G* ESP.

The results for molecular multipole moments given for ammonia are representative of the other 14 monomers studied in this work. The RMSD errors in molecular multipole moment up to hexadecapole are averaged over all the molecules and presented in Table IV. As expected, there is a significant decrease in average RMSD error for increasing atomic Gaussian multipole order. For example, the average RMSD errors in molecular hexadecapole are 3.402 DÅ3 for atomic Gaussian monopoles, 0.969 for Gaussian dipoles, and 0.130 for Gaussian quadrupoles. For more individual results, see the SI.

Table IV
Average RMSD error (Δ) in Traceless Molecular Dipole, Quadrupole, Octapole, and Hexadecapole moment over 14 molecules predicted by Gaussian monopoles (M), dipoles (D), and quadrupoles (Q) (Nc = 4) fit to B3LYP/6-31G* ESP.

3.4 Gaussian Multipole Geometry Dependence

The atomic Gaussian multipole moments are investigated as a function of bond length r and bond angle θ for the case of water. Gaussian quadrupoles (Nc = 4) are calculated at the B3LYP/6-31G* level for several geometries of water obtained by perturbing either the O-H1 bond length or the H1-O-H2 bond angle. The atomic Gaussian multipole moments in the local frame (eqn. 19) are converted to their traceless Cartesian moments (eqns. 27 and 28). In Fig. 6, the atomic Gaussian monopole moment q on oxygen is plotted when the H1-O bond length is varied. The Gaussian monopole charge q is a smooth function of bond length, which can be approximated as a straight line for bond lengths between 0.92 and 1.00 Å (the equilibrium bond length is 0.9684 Å). In addition, the zz component of the local frame traceless Cartesian quadrupole moment Θzzlocal on the H1 hydrogen is plotted as a function of H1-O-H2 bond angle in Fig. 7. Θzzlocal can be approximated as a straight line for bond angles between 95° and 111° (the equilibrium bond angle is 103.66°). In both Figs. 6 and and7,7, the straight lines are determined by the value of the local frame atomic multipoles Qlmlocal at the equilibrium geometry and the finite difference derivative of Qlmlocal with respect to the bond length or bond angle at the equilibrium geometry. The results given in Figs. 6 and and77 are examples, and the other atomic Gaussian moments Qlmlocal follow similar trends (see the SI for more examples). For small internal geometry perturbations, the above results suggest that the local frame atomic Gaussian multipole moments can be approximated by a truncated linear Taylor Series as a function of the two bond lengths r1 (H1-O) and r2 (H2-O) and the bond angle θ (H1-O-H2) as

Qlmlocal(r,θ)Qlmlocal,0+(r1r10)Qlmlocal,0r1+(r2r20)Qlmlocal,0r2+(θθ0)Qlmlocal,0θ,
(29)

where Qlmlocal,0 is the atomic Gaussian multipole at the geometry optimized equilibrium structure. Qlmlocal,0r1, Qlmlocal,0r2, and Qlmlocal,0θ are the finite difference partial derivatives of Qlmlocal with respect to r1, r2, and θ, respectively, evaluated at the equilibrium structure.

Figure 6
Atomic Gaussian monopole moment q (e) on oxygen in a water molecule as a function of H1-O bond length. The B3LYP/6-31G* equilibrium H1-O bond length is 0.9684 Å.
Figure 7
Atomic local frame Gaussian quadrupole moment Θzzlocal (eÅ2) on the `H1' hydrogen in a water molecule as a function of H1-O-H2 bond angle. The B3LYP/6-31G* equilibrium H1-O-H2 bond angle at is 103.66°.

3.5 Exchange-Overlap Model

The exchange-overlap model is fit to exchange-repulsion energies calculated at the B3LYP/6-31G* level through CSOV decomposition as described in section 2.7. The inter-molecular density overlap integrals are calculated from Gaussian quadrupoles (Nc = 4) fit to the B3LYP/6-31G* ESP. For the water-water dimer, the RMSD error of fit is 0.350 kcal/mol for Gaussian quadrupoles, over a range of exchange energies from 0.0 to 27.0 kcal/mol. At the equilibrium water-water dimer distance, the exchange energy calculated by the exchange-overlap model is 5.449 kcal/mol, compared to the ab initio result of 5.362 kcal/mol. In Fig. 8, the water-water exchange energy calculated by the model is compared with the ab initio value and plotted as a function of H..O distance. A similar analysis is performed on the other hydrogen bonded dimers. The RMSD fit errors are averaged over all the dimers and given by 0.764, 0.379, and 0.275 kcal/mol for Gaussian monopoles, dipoles, and quadrupoles, respectively. This result indicates that including anisotropy in the model for charge density makes a significant improvement when applying the exchange-overlap model. At equilibrium distances, the exchange energies lie between 0.1 and 10 kcal/mol. The RMSD in exchange energy at the equilibrium dimer distance averaged over the hydrogen bonded dimers are 1.782, 0.262, and 0.276 kcal/mol for Gaussian monopoles, dipoles, and quadrupoles, respectively. For Gaussian quadrupoles, the molecular pair K parameters range from 0.5786 kcal/mol (103 Å3/e2)0.95 for the water-CH3F dimer to 0.7499 for the water-ammonia dimer. The average value of the molecular pair K parameter over all the dimers is 0.6607. Because of the variability of the K parameters, a single molecular pair K parameter may not be sufficient for larger molecules or when a large amount of configuration space is sampled. For more individual results, including exchange parameters K, RMSD errors, and exchange dimer energies, see the SI.

Figure 8
The exchange energy calculated by Gaussian quadrupoles and ab initio is plotted as a function of H..O distance for the water–water dimer.

4. Conclusion

We have proposed a model based on contracted Gaussian multipole charge density. The atomic Gaussian multipoles are fit to the ab initio electrostatic potential and are shown to reproduce ab initio electrostatic dimer energies, inter-molecular density overlap integrals, and permanent molecular multipole moments. For the case of water, the local frame atomic Gaussian multipole moments Qlmlocal are shown to be a smooth function of bond length r and bond angle θ, which can be approximated as a truncated linear Taylor series. In a follow up work, we will present analytic atomic force expressions for geometry dependent Gaussian multipoles and show that geometry dependent electrostatic models are capable of reproducing ab initio electrostatic atomic forces. In addition, the inter-molecular density overlap integrals calculated by Gaussian multipoles has been applied to a model61,62 for exchange-repulsion energy based on inter-molecular density overlap integral. A molecular pair K parameter is fit to the ab initio exchange-repulsion energy for hydrogen bonded dimers. A significant improvement is found in going from Gaussian monopoles to Gaussian quadrupoles, indicating that including anisotropy in the description of atomic charge density is important. Though the preliminary results of applying the exchange-overlap model to the Gaussian multipole charge density are encouraging, a more extensive investigation would be useful, possibly by studying atomic pair K parameters fit to a larger molecular data set. We plan to further study the exchange-overlap model using Gaussian multipoles, and we hope to propose a general set of transferable exchange-overlap parameters in the near future.

Supplementary Material

1_si_001

2_si_002

3_si_003

5. Acknowledgments

This research was supported in part by the Intramural Research Program of the NIH and NIEHS (Z01 ES9043010-23). LGP acknowledges support from NSF FRG DMR 0804549 and from NIH L06350. We would like to thank the reviewers for their helpful comments, which improved the manuscript.

Footnotes

There is a three part Supplementary Information (SI) available. In part I of the SI, additional mathematical details and background are given. Elementary properties of spherical, solid, and scaled solid harmonic functions are given along with some theorems for the solid harmonic function and the solid harmonic gradient operator. Expressions for electrostatic energy and density overlap integral are derived for complex Gaussian multipoles along with the Cartesian gradients of their matrix elements. In addition, formulae for converting complex spherical tensor multipole moments into their traceless Cartesian forms are given. In part II of the SI, additional results for Gaussian multipoles are presented. The dependence on the degree of Slater-type contraction Nc is discussed. Additional tables and figures of electrostatic dimer energies and inter-molecular density overlap integrals are presented. For the case of water, additional plots of Qlmlocal as a function of bond length and bond angle are given. For the exchange-overlap model, exchange-repulsion energies, molecular pair K parameters, and RMSD of fits are provided for B3LYP/6-31G* Gaussian multipoles. In part III of the SI, the Gaussian multipole parameters, the Slater-type contraction coefficients/exponents, and the equilibrium dimer geometries are given. This information is available free of charge via the Internet at http://pubs.acs.org/.

REFERENCES

1. Thole BT. Chem. Phys. 1981;59:341–350.
2. Elking DM, Darden T, Woods RJ. J. Comput. Chem. 2007;28:1261–1274. [PMC free article] [PubMed]
3. Lamoureux G, MacKerell AD, Roux B. J. Chem. Phys. 2003;119(10):5185–5197.
4. Rick SW, Stuart SJ, Berne BJ. J. Chem. Phys. 1994;101(7):6141–6156.
5. Stern HA, Kaminski GA, Banks JL, Zhou R, Berne BJ, Friesner R. J. Phys. Chem. 1999;103:4730–4737.
6. Piquemal J-P, Chelli R, Procacci P, Gresh N. J. Phys. Chem. A. 2007;111:8170–8176. [PubMed]
7. Patel S, Brooks CL., III J. Comput. Chem. 2004;25:1–16. [PubMed]
8. Patel S, MacKerell AD, Brooks CL., III J. Comput. Chem. 2004;25:1504–1514. [PubMed]
9. Kaminski GA, Stern HA, Berne BJ, Friesner R. J. Phys. Chem. A. 2004;108:621–627.
10. Gresh N, Cisneros GA, Darden TA, Piquemal J-P. J. Chem. Theory Comput. 2007;3:1960–1986. [PMC free article] [PubMed]
11. Gresh N. J. Comput. Chem. 1995;16:856–882.
12. Piquemal J-P, Chevreau H, Gresh N. J. Chem. Theory Comput. 2007;3:824–837.
13. Ren P, Ponder JW. J. Phys. Chem. B. 2003;107:5933–5947.
14. Ren P, Ponder JW. J. Comput. Chem. 2002;23:1497–1506. [PubMed]
15. Grossfield A, Ren P, Ponder JW. J. Am. Chem. Soc. 2003;125:15671–15682. [PubMed]
16. Piquemal J-P, Perera L, Cisneros GA, Ren P, Pedersen LG, Darden TA. J. Chem. Phys. 2006;125:054511. [PubMed]
17. Iuchi S, Izvekov S, Voth GA. J. Chem. Phys. 2007;126:124505–12. [PubMed]
18. Paesani F, Iuchi S, Voth GA. J. Chem. Phys. 2007;127(7):074506–15. [PubMed]
19. Wang F, Jordan KD. J. Chem. Phys. 2002;116:6973–6981.
20. Jiang H, Jordan KD, Taylor CE. J. Phys. Chem. B. 2007;111:6486–6492. [PubMed]
21. Koch U, Popelier PLA, Stone AJ. Chem. Phys. Lett. 1995;2:253–260.
22. Koch U, Stone AJ. J. Chem. Soc., Faraday Trans. 1996;92(10):1701–1708.
23. Cho K, Kang YK, No KT, Scheraga HA. J. Phys. Chem. B. 2001;105:3624–3634.
24. Mankoo P, Keyes T. J. Chem. Phys. 2008;129:034504–9. [PubMed]
25. Stone AJ. The Theory of Intermolecular Forces. Oxford University Press; Oxford, UK: 2000.
26. Stone AJ. Chem. Phys. Lett. 1981;83:233–239.
27. Vigné-Maeder F, Claverie P. J. Chem. Phys. 1998;88:4934–4948.
28. Náray-Szabó G, Ferenczy GG. Chem. Rev. 1995;95(4):829–847.
29. Ángyán JG, Chipot C, Dehez F, Hättig C, Jansen G, Millot C. J. Comput. Chem. 2003;24:997–1008. [PubMed]
30. Qian W, Krimm S. J. Phys. Chem. A. 2005;109:5608–5618. [PubMed]
31. Freitag MA, Gordon MS, Jensen JH, Stevens WJ. J. Chem. Phys. 2000;112(17):7300–7306.
32. Piquemal J-P, Gresh N, Giessner-Prettre C. J. Phys. Chem. A. 2003;107:10353–10359. [PubMed]
33. Cisneros GA, Tholander S. Na-Im, Elking DM, Darden TA, Parisel O, Piquemal J-P. Int. J. Quant. Chem. 2008;108:1905–1912. [PMC free article] [PubMed]
34. Slipchenko LV, Gordon MS. Mol. Phys. 2009;107:999–1016.
35. Chelli R, Righini R, Califano S, Procacci P. J. Mol. Liq. 2002;96–97:87–100.
36. Masia M, Probst M, Rey R. J. Chem. Phys. 2005;123:164505–13. [PubMed]
37. Paricaud P, Předota M, Chialvo A, Cummings P. J. Chem. Phys. 2005;122:244511–14. [PubMed]
38. Wheatley RJ. Mol. Phys. 1993;79:597–610.
39. Wheatley RJ, Mitchell JBO. J. Comp. Chem. 1994;15:1187–1198.
40. Hall GG. Adv. at. molec. Phys. 1985;20:41–63.
41. Martin D, Hall GG. Theor. chim. Acta. 1981;59:281–290.
42. Giese TJ, York DM. J. Chem. Phys. 2008;128(6):064104–6. [PMC free article] [PubMed]
43. Hobson EW. The Theory of Spherical and Ellipsoidal Harmonics. Chelsea; New York, NY: 1955. p. 93.
44. Bayman BF. J. Math. Phys. 1978;19:2558–2562.
45. Chakrabarti S, Dewangan DP. J. Phys. B: At. Mol. Opt. Phys. 1995;28:L769–774.
46. Helgaker T, Jorgensen P, Olsen J. Molecular Electronic-Structure Theory. Wiley; Chichester, UK: 2004. pp. 337–424.
47. Arken GB. Mathematical Methods for Physicists. 5th ed. Academic Press; San Diego, CA: 2000. pp. 693–765.
48. Cisneros GA, Piquemal J-P, Darden TA. J. Chem. Phys. 2005;123:044109–10. [PMC free article] [PubMed]
49. Cisneros GA, Piquemal J-P, Darden TA. J. Chem. Phys. 2006;125(18):184101–16. [PMC free article] [PubMed]
50. Piquemal J-P, Cisneros GA, Reinhardt P, Gresh N, Darden TA. J. Chem. Phys. 2006;124(10):104101–12. [PMC free article] [PubMed]
51. Cisneros GA, Darden TA, Gresh N, Reinhardt P, Parisel O, Pilmé J, Piquemal J-P. Design of Next Generation Force Fields from Ab Initio Computations: Beyond Point Charges. In: York DM, Lee T-S, editors. Electrostatics Multi-scale Quantum Models for Biocatalysis: Modern Techniques and Applications, for the Book Series: Challenges and Advances in Computational Chemistry and Physics. Vol. 7. Springer Verlag; New Amsterdam, The Netherlands: 2009. pp. 137–172.
52. Dunlap BI, Connolly WD, Sabin JR. J. Chem. Phys. 1979;71:4993–4999.
53. Köster AM. J. Chem. Phys. 1996;104:4114–4124.
54. Köster AM. J. Chem. Phys. 2003;118:9943–9951.
55. Eichkorn K, Treutler O, Ohm H, Haser M, Ahlrichs R. Chem. Phys. Lett. 1995;240:283–289.
56. Jung Y, Sodt A, Gill P. M. W. Gill, Head-Gordon M. Proc. Natl. Acad. Sci. U.S.A. 2005;102:6692–6697. [PubMed]
57. Cisneros GA, Elking DM, Piquemal J-P, Darden TA. J. Phys. Chem. A. 2007;111:12049–12056. [PMC free article] [PubMed]
58. Hehre WJ, Stewart RF, Pople JA. J. Chem. Phys. 1969;51(6):2657–2664.
59. Stewart RF. J. Chem. Phys. 1970;52(1):431–438.
60. Hu H, Lu Z, Yang W. J. Chem. Theory Comput. 2007;3:1004–1013.
61. Wheatley RJ, Price SL. Mol. Phys. 1990;69:507–533.
62. Mitchell JBO, Price SL. J. Phys. Chem. A. 2000;104:10958–10971.
63. Dunlap BI. Phys. Rev. A. 1990;42(3):1127–1137. [PubMed]
64. Dunlap BI. Int. J. Quant. Chem. 2001;81:373–383.
65. Dunlap BI. J. Chem. Phys. 2003;118(3):1036–1043.
66. Choi CH, Ivanic J, Gordon MS, Ruedenberg K. J. Chem. Phys. 1999;11:8825–8831.
67. Toukmaji A, Sagui C, Board JA, Darden TA. J. Chem. Phys. 2000;113:10913–10927.
68. Sagui C, Pedersen LG, Darden TA. J. Chem. Phys. 2004;120:73–87. [PubMed]
69. Bylaska EJ, de ong WA, Govind N, Kowalski K, Straatsma TP, Valiev M, Wang D, Apra E, Windus TL, Hammond J, Nichols P, Hirata S, Hackler MT, Zhao Y, Fan P-D, Harrison RJ, Dupuis M, Smith DMA, Nieplocha J, Tipparaju V, Krishnan M, Wu Q, Van Voorhis T, Auer AA, Nooijen M, Brown E, Cisneros G, Fann GI, Fruchtl H, Garza J, Hirao K, Kendall R, Nichols JA, Tsemekhman K, Wolinski K, Anchell J, Bernholdt D, Borowski P, Clark T, Clerc D, Dachsel H, Deegan M, Dyall K, Elwood D, Glendening E, Gutowski M, Hess A, Jaffe J, Johnson B, Ju J, Kobayashi R, Kutteh R, Lin Z, Littlefield R, Long X, Meng B, Nakajima T, Niu S, Pollack L, Rosing M, Sandrone G, Stave M, Taylor H, Thomas G, van Lenthe J, Wong A, Zhang Z. NWChem, A Computational Chemistry Package for Parallel Computers, Version 5.1. Pacific Northwest National Laboratory; Richland, WA: 2007. A modified version.
70. Kendall RA, Aprà E, Bernholdt DE, Bylaska EJ, Dupuis M, Fann GI, Harrison RJ, Ju J, Nichols JA, Nieplocha J, Straatsma TP, Windus TL, Wong AT. High Performance Computational Chemistry: an Overview of NWChem a Distributed Parallel Application. Computer Phys. Comm. 2000;128:260–283.
71. Press WH, Flannery BP, Teukolsky SA, Vetterling WT. Numerical Recipes in C: The Art of Scientific Computing. 2nd ed. Cambridge University Press; Cambridge: 1992. p. 683.
72. Frisch MJ, Trucks GW, Schlegel HB, Scuseria GE, Robb MA, Cheeseman 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, Revision C.02. Gaussian, Inc.; Wallingford CT: 2004.
73. Breneman CM, Wiberg KB. J. Comp. Chem. 1990;11:361–373.
74. Bagus PS, Hermann K, Bauschlicher CW., Jr. J. Chem. Phys. 1984;80:4378–4386.
75. Piquemal J-P, Marquez A, Parisel O, Giessner-Prettre C. J. Comput. Chem. 2005;26:1052–1062. [PubMed]
76. Dupuis M, Marquez A, Davidson ER. HONDO95.3. Indiana University; Bloomington, IN: 1995.
77. Stevens WJ, Fink WH. Chem. Phys. Lett. 1987;139(1):15–22.
78. Chen W, Gordon MS. J. Phys. Chem. 1996;100:14316–14328.
79. Schmidt MW, Baldridge KK, Boatz JA, Elbert ST, Gordon MS, Jensen JJ, Koseki S, Matsunaga N, Nguyen KA, Su S, Windus TL, Dupuis M, Montgomery JA. J. Comput. Chem. 1993;4:1347.
80. McMurchie LE, Davidson ER. J. Comp. Phys. 1978;26:218–231.
81. Özdoğan T. J. Math. Chem. 2006;42(2):201–214.
82. Tschumper GS, Leininger ML, Hoffman BC, Valeev EF, Schaffer HF, III, Quack M. J. Chem. Phys. 2002;116:690–701.
83. van Duijneveldt-van, de Rijdt JGCM, Mooij WTM, Duijneveldt FB. Phys. Chem. Chem. Phys. 2003;5:1169–1180.