Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2832208

Formats

Article sections

Authors

Related links

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/ct900348bPMCID: PMC2832208

NIHMSID: NIHMS161432

Dennis M. Elking,^{1} G. Andrés Cisneros,^{1} Jean-Philip Piquemal,^{2} Thomas A. Darden,^{1} and Lee G. Pedersen^{1,}^{3,}^{*}

See other articles in PMC that cite the published article.

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 *Q _{lm}* 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.

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 models^{1}^{–}^{6} have been incorporated into force fields^{7}^{–}^{16} in order to account for many body effects^{17}^{–}^{20} in polar environments. In the SIBFA^{10}^{–}^{12} force field, multipoles are placed on atoms and bond barycenters in order to accurately account for the anisotropy in electrostatic interactions. AMOEBA^{13}^{–}^{16} 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 investigated^{21}^{–}^{24}.

Atomic point multipole models derived from distributed multipole expansions^{25}^{–}^{28} or fit to the electrostatic potential^{29} (ESP) have been proposed to improve the description of electrostatic interactions. However, at short range, it has been noted^{25}^{,}^{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 functions^{31}^{–}^{34} 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 densities^{35}^{–}^{37} have been used in models for liquids. Wheatley^{38}^{,}^{39} has studied Cartesian Gaussian multipole charge distributions^{40}^{,}^{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 multipoles^{38}. As an example, consider a Cartesian Gaussian dipole charge distribution ρ_{μ} with dipole moment $\overrightarrow{\mu}$, exponent α, and nuclear center $\overrightarrow{R}$ given by

$${\rho}_{\mu}\left(\overrightarrow{r}\right)\equiv \overrightarrow{\mu}\cdot {\nabla}_{R}{\left(\frac{{\alpha}^{2}}{\pi}\right)}^{\frac{3}{2}}\mathrm{exp}(-{\alpha}^{2}{\mid \overrightarrow{r}-\overrightarrow{R}\mid}^{2}),$$

(1)

where $\overrightarrow{r}$ is the field point and _{R} is the Cartesian gradient with respect to the nuclear center $\overrightarrow{R}$. The electrostatic potential ${\varphi}_{\mu}\left(\overrightarrow{r}\right)$ arising from ρ_{μ} is given by

$${\varphi}_{\mu}\left(\overrightarrow{r}\right)=\int {d}^{3}{r}^{\prime}\frac{{\rho}_{\mu}\left({\overrightarrow{r}}^{\prime}\right)}{\mid \overrightarrow{r}-{\overrightarrow{r}}^{\prime}\mid}=\overrightarrow{\mu}\cdot {\nabla}_{R}\frac{erf\left(\alpha \mid \overrightarrow{r}-\overrightarrow{R}\mid \right)}{\mid \overrightarrow{r}-\overrightarrow{R}\mid},$$

(2)

where erf(*x*) is the error function defined by

$$\mathrm{erf}\left(x\right)\equiv \frac{2}{\sqrt{\pi}}{\int}_{0}^{x}\mathrm{exp}(-{u}^{2})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.

$${\varphi}_{\mu}\left(\overrightarrow{r}\right)\approx \overrightarrow{\mu}\cdot {\nabla}_{R}\frac{1}{\mid \overrightarrow{r}-\overrightarrow{R}\mid}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\alpha \mid \overrightarrow{r}-\overrightarrow{R}\mid >>1.$$

(4)

Recently, Giese and York^{42} 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 *Q*_{lm} with the solid harmonic gradient operator^{43}^{,}^{44} *C _{lm}* () acting upon a simple Gaussian function. In the Supplementary Information (SI), a summary of the mathematical background

Methods of calculating the Gaussian charge density parameters are also being explored. Previously, we have proposed the Gaussian Electrostatic Model (GEM)^{10}^{,}^{48}^{–}^{51} 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 energy^{52}^{–}^{56} Δ*E _{self}* given by

$$\Delta {E}_{\mathit{self}}=\langle {\rho}^{\mathrm{QM}}\left(r\right)-\rho \left(r\right)\mid \frac{1}{r}\mid {\rho}^{\mathrm{QM}}\left(r\right)-\rho \left(r\right)\rangle .$$

(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 grid^{57}. 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 Δ*E _{self}*. 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-type^{58}^{,}^{59} contracted Gaussian function fit to exp(-λ*r*). For each atom, the Gaussian multipole moments *Q _{lm}* 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

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 *Q _{lm}* 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

The inter-molecular density overlap integral calculated by Gaussian multipoles can be applied to the exchange-overlap model proposed by Wheatley and Price^{61}^{,}^{62} for the exchange-repulsion energy. The first order inter-molecular exchange-repulsion energy *E _{exch}* is defined as the total

$${E}_{\mathit{exch}}\cong K{S}^{a},$$

(6)

where *a* is an empirical exponent parameter used to improved the quality of fit^{62}. The expression in eqn. 6 is commonly generalized to a pair-wise sum over atom-atom^{61}^{,}^{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 work^{48}^{,}^{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 *Q _{lm}* 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.

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.

The definition we have used for a contracted Gaussian multipole charge density is similar to the one given by Giese and York^{42}. In that work^{42}, 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 background^{43}^{–}^{47} 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 functions^{63}^{–}^{65} 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 *Q _{lm}* and the solid harmonic gradient operator

A contracted Gaussian multipole charge density $\rho (\overrightarrow{r},\overrightarrow{R})$ with moments *Q*_{lm} and nuclear center $\overrightarrow{R}$ evaluated at the point $\overrightarrow{r}$ is given by

$$\rho (\overrightarrow{r},\overrightarrow{R})\equiv \sum _{l=0}^{{l}_{\mathrm{max}}}\sum _{\mid m\mid \le l}\frac{{Q}_{lm}{C}_{lm}^{\ast}(\overrightarrow{r}-\overrightarrow{R})}{(2l-1)!!}{\rho}_{l}(\mid \overrightarrow{r}-\overrightarrow{R}\mid ;{\alpha}_{\mu}),$$

(7)

where *l _{max}* is the maximum order of the Gaussian multipoles (e.g.

$${\rho}_{l}(r;{a}_{\mu})\equiv {\left(-\frac{1}{r}\frac{d}{dr}\right)}^{l}\sum _{\mu =1}^{{N}_{c}}{c}_{\mu}{\left(\frac{{\alpha}_{\mu}^{2}}{\pi}\right)}^{3\u22152}\mathrm{exp}(-{\alpha}_{\mu}^{2}{r}^{2}),$$

(8)

where *N _{c}* is the degree of contraction. For

$$\int {d}^{3}r\rho (\overrightarrow{r},\overrightarrow{R}){C}_{lm}(\overrightarrow{r}-\overrightarrow{R})={Q}_{lm}.$$

(9)

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

$$\varphi \left(\overrightarrow{r}\right)=\int {d}^{3}{r}^{\prime}\frac{\rho (\overrightarrow{r},\overrightarrow{R})}{\mid \overrightarrow{r}-{\overrightarrow{r}}^{\prime}\mid}=\sum _{l=0}^{{l}_{\mathrm{max}}}\sum _{\mid m\mid \le l}\frac{{Q}_{lm}{C}_{lm}^{\ast}(\overrightarrow{r}-\overrightarrow{R})}{(2l-1)!!}{\varphi}_{l}(\mid \overrightarrow{r}-\overrightarrow{R}\mid ;{\alpha}_{\mu}),$$

(10)

where ϕ* _{l}* is defined by

$${\varphi}_{l}(r;{\alpha}_{\mu})\equiv {\left(-\frac{1}{r}\frac{d}{dr}\right)}^{l}\sum _{\mu =1}^{{N}_{c}}{c}_{\mu}\frac{\mathrm{erf}\left({\alpha}_{\mu}r\right)}{r},$$

(11)

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

The electrostatic interaction energy between two Gaussian multipole charge densities ${\rho}_{1}(\overrightarrow{r},{\overrightarrow{R}}_{1})$ and ${\rho}_{2}(\overrightarrow{r},{\overrightarrow{R}}_{2})$ and their nuclei Z_{1} and Z_{2} is given by

$$U=\int \int {d}^{3}r{d}^{3}{r}^{\prime}\frac{{\rho}_{1}(\overrightarrow{r},{\overrightarrow{R}}_{1}){\rho}_{2}({\overrightarrow{r}}^{\prime},{\overrightarrow{R}}_{2})}{\mid \overrightarrow{r}-{\overrightarrow{r}}^{\prime}\mid}+{Z}_{1}{\varphi}_{2}\left({\overrightarrow{R}}_{1}\right)+{Z}_{2}{\varphi}_{1}\left({\overrightarrow{R}}_{2}\right)+\frac{{Z}_{1}{Z}_{2}}{\mid {\overrightarrow{R}}_{1}-{\overrightarrow{R}}_{2}\mid},$$

(12)

where ${\varphi}_{2}\left({\overrightarrow{R}}_{1}\right)$ is the potential at ${\overrightarrow{R}}_{1}$ due to ρ_{2} and ${\varphi}_{1}\left({\overrightarrow{R}}_{2}\right)$ is defined by a similar expression. The density overlap integral *S* between two Gaussian multipole charge densities is defined by

$$S=\int {d}^{3}r{\rho}_{1}\left(\overrightarrow{r}\right){\rho}_{2}\left(\overrightarrow{r}\right).$$

(13)

The electrostatic and density overlap integrals in eqns. 12 and 13 can be expressed by^{42}

$$\int \int {d}^{3}r{d}^{3}{r}^{\prime}{\rho}_{1}\left(\overrightarrow{r}\right)O(\overrightarrow{r},{\overrightarrow{r}}^{\prime}){\rho}_{2}\left({\overrightarrow{r}}^{\prime}\right)=\sum _{{l}_{1},{m}_{1}}\sum _{{l}_{2},{m}_{2}}{Q}_{{l}_{1}{m}_{1}}^{1}{T}_{{l}_{1}{m}_{1};{l}_{2}{m}_{2}}{Q}_{{l}_{2}{m}_{2}}^{2},$$

(14)

where $O(\overrightarrow{r},{\overrightarrow{r}}^{\prime})\equiv 1\u2215\mid \overrightarrow{r}-{\overrightarrow{r}}^{\prime}\mid $ for the electrostatic integral, $O(\overrightarrow{r},{\overrightarrow{r}}^{\prime})\equiv \delta (\overrightarrow{r}-{\overrightarrow{r}}^{\prime})$ for the density overlap integral, and the interaction matrix ${T}_{{l}_{1}{m}_{1};{l}_{2}{m}_{2}}$ is given below. The electrostatic energy and density overlap integrals between two normalized contracted Gaussian monopole charge densities with unit charge are given by

$${F}_{0}\left(R\right)\equiv \sum _{\mu ,v}{c}_{\mu}^{1}{c}_{v}^{2}\frac{\mathrm{erf}\left({\alpha}_{\mu v}R\right)}{R},\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}{F}_{0}\left(R\right)\equiv \sum _{\mu ,v}{c}_{\mu}^{1}{c}_{v}^{2}{\left(\frac{\pi}{{{\alpha}_{\mu v}}^{2}}\right)}^{3\u22152}\mathrm{exp}(-{{\alpha}_{\mu v}}^{2}{R}^{2}),$$

(15)

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

$$R\equiv \mid {\overrightarrow{R}}_{1}-{\overrightarrow{R}}_{2}\mid \phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}{\alpha}_{\mu v}\equiv \frac{{\alpha}_{\mu}{\alpha}_{v}}{\sqrt{{\alpha}_{\mu}^{2}+{\alpha}_{v}^{2}}}.$$

(16)

In addition, the constants *A _{lm}* and

$${A}_{lm}\equiv \sqrt{(l+m)!(l-m)!},\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}{B}_{lm}\equiv \frac{{A}_{lm}}{(2l-1)!!},$$

(17)

for *l* ≠ 0, and *A*_{00} = *B*_{00} = 1. The scaled solid harmonic function is defined by ${R}_{lm}\left(\overrightarrow{r}\right)\equiv {C}_{lm}(\overrightarrow{r}\u2215{A}_{lm})$. The interaction matrix ${T}_{{l}_{1}{m}_{1};{l}_{2}{m}_{2}}$ from eqn. 14 is given by

$${T}_{{l}_{1}{m}_{1};{l}_{2}{m}_{2}}={B}_{{l}_{1}{m}_{1}}\sum _{l=0}^{\mathrm{min}({l}_{2},{l}_{2})}\sum _{m=-1}^{l}\frac{{(-1)}^{{l}_{2}+m}}{{A}_{lm}{B}_{lm}}{R}_{{l}_{2}-l,{m}_{2}+m}^{\ast}\left(\overrightarrow{R}\right){R}_{{l}_{1}-l,{m}_{2}-m}^{\ast}\left(\overrightarrow{R}\right){F}_{{l}_{2}+{l}_{1}-1}\left(R\right),$$

(18)

where *F _{l}*(

The atomic Gaussian multipole moments *Q*_{lm} are commonly defined^{25} with respect to a local frame of the atom ${Q}_{lm}^{\mathit{local}}$ and then rotated to the system or global frame ${Q}_{lm}^{\mathit{local}}\equiv {Q}_{lm}$ through Wigner rotation matrices ${D}_{{m}^{\prime}m}^{l}$

$${Q}_{lm}^{\mathit{global}}=\sum _{{m}^{\prime}}{D}_{{m}^{\prime}m}^{l}\left[{\mathbf{R}}^{-1}\right]{Q}_{l{m}^{\prime}}^{\mathit{local}}.$$

(19)

Recursion formulae for evaluating ${D}_{m{m}^{\prime}}^{l}$ have been given in Choi^{66} 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 neighbors^{13}^{,}^{67}^{,}^{68}.

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

$$\mathrm{exp}(-r)\cong \sum _{\mu =1}^{{N}_{c}}{d}_{\mu}\phantom{\rule{thickmathspace}{0ex}}\mathrm{exp}(-{\alpha}_{\mu}^{2}{r}^{2}).$$

(20)

For *N _{c}* = 1 − 6, the optimized exponents agree with those used in the development of the original STONG basis sets

$$\frac{{\lambda}^{3}}{8\pi}\mathrm{exp}(-\lambda r)\cong \sum _{\mu =1}^{{N}_{c}}{c}_{\mu}{\left(\frac{{\alpha}_{\mu}^{2}{\lambda}^{2}}{\pi}\right)}^{3\u22152}\mathrm{exp}(-{\alpha}_{\mu}^{2}{\lambda}^{2}{r}^{2}),$$

(21)

where ${c}_{\mu}\equiv {d}_{\mu}\u22158\pi \cdot {(\pi \u2215{\alpha}_{\mu}^{2})}^{3\u22152}$. A full list of contracted coefficients and exponents for exp(-*r*) can be found in the SI for *N _{c}* = 1 − 14.

The model for molecular charge density presented in this work consists of an effective nuclear charge *Z*_{eff} and a set of contracted Gaussian multipole moments *Q*_{lm} 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 *Z*_{eff} = Z − N_{core}, where Z is the true nuclear charge and N_{core} is the number of core electrons. The number of core electrons N_{core} 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 *Z*_{eff} 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 $\overrightarrow{r}$ is a represented as a sum over atomic Gaussian multipole charge densities given by

$${\rho}^{\mathrm{GM}}\left(\overrightarrow{r}\right)=\sum _{a}\sum _{l=0}^{{l}_{\mathrm{max}}}\sum _{\mid m\mid \le l}\frac{{Q}_{lm,a}{C}_{lm}^{\ast}(\overrightarrow{r}-{\overrightarrow{R}}_{a})}{(2l-1)!!}{\rho}_{l}(\mid \overrightarrow{r}-{\overrightarrow{R}}_{a}\mid ;{\lambda}_{a}{\alpha}_{\mu})$$

(22)

where ${\overrightarrow{R}}_{a}$ is the nuclear center of atom *a*, *Q*_{lm,a} are the atomic Gaussian multipole moments of atom *a*, λ* _{a}* is the Slater exponent on atom

$$\varphi (\overrightarrow{r};{Q}_{lm,a},{\lambda}_{a})=\sum _{a}\frac{{Z}_{eff,a}}{\mid \overrightarrow{r}-{\overrightarrow{R}}_{a}\mid}+\sum _{l=0}^{{l}_{\mathrm{max}}}\sum _{\mid m\mid \le l}\frac{{Q}_{lm,a}{C}_{lm}^{\ast}(\overrightarrow{r}-{\overrightarrow{R}}_{a})}{(2l-1)!!}{\varphi}_{i}(\mid \overrightarrow{r}-{\overrightarrow{R}}_{a}\mid ,{\lambda}_{a}{\alpha}_{\mu}),$$

(23)

where $\overrightarrow{r}$ is the field point, Z_{eff,a} is the effective nuclear charge of atom *a*, and ϕ* _{l}* is defined by eqn. 11. For each atom,

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

$$\mathcal{W}\left(\overrightarrow{R}\right)\equiv \{\begin{array}{cc}\mathrm{exp}\{-\sigma {[\mathrm{ln}{\rho}^{QM}\left(\overrightarrow{r}\right)-\mathrm{ln}{K}_{0}]}^{2}\}\hfill & {\rho}^{QM}\left(\overrightarrow{r}\right)\ge {K}_{0}\hfill \\ 1\hfill & {\rho}^{QM}\left(\overrightarrow{r}\right)\le {K}_{0}\hfill \end{array}\phantom{\}}$$

(24)

where ρ* ^{QM}* is the

Contour plot of the average RMSD error in electrostatic dimer energy (kcal/mol) as a function of parameters σ and ln*K*_{0} for the ESP weighting function. The electrostatic energies are calculated on hydrogen bonded dimers at equilibrium geometries **...**

The fitting function χ^{2} is given by

$${\chi}^{2}({Q}_{lm},\lambda )=\int {d}^{3}r\mathcal{w}\left(\overrightarrow{r}\right){[{\varphi}^{GM}({Q}_{lm},\lambda ;\overrightarrow{r})-{\varphi}^{QM}\left(\overrightarrow{r}\right)]}^{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 NWChem^{69}^{,}^{70} and optimized using a Levenberg-Marquardt non-linear least-squares fit algorithm^{71}. The *ab initio* electrostatic potential is calculated at both the B3LYP/6-31G* or HF/aug-cc-pVTZ levels using the Gaussian 03 software package^{72}. Earlier in our study, we had experimented with using rectangular grids similar to the CHELPG^{73} 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 (*N _{c}* = 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 10

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) decomposition^{74} method using a modified version of the HONDO^{75}^{,}^{76} quantum chemistry program. For the model fit to HF/aug-cc-pVTZ data, *ab initio* electrostatic energies are calculated by the Reduced Variational Space^{77}^{,}^{78} (RVS) decomposition method using the GAMESS^{79} 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 algorithm^{80}.

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.

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 *Q _{lm}* at position $\overrightarrow{R}$ are translated to the origin by the following expression

$${Q}_{lm}^{\mathit{origin}}=\sum _{{l}_{1}=0}^{{l}_{\mathrm{max}}}\sum _{{m}_{1}=-{1}_{1}}^{{l}_{1}}\sqrt{\left(\begin{array}{c}\hfill l+m\hfill \\ \hfill {l}_{1}+{m}_{1}\hfill \end{array}\right)\left(\begin{array}{c}\hfill l-m\hfill \\ \hfill {l}_{1}-{m}_{1}\hfill \end{array}\right)}{Q}_{{l}_{1}{m}_{1}}{C}_{l-{l}_{1},m-{m}_{1}}\left(\overrightarrow{R}\right),$$

(26)

where *C _{lm}* is a solid harmonic function. Note that ${Q}_{lm}^{\mathit{origin}}\ne 0$ for all values of

$${\mu}_{x}=-\sqrt{2}{Q}_{11}^{r},\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}{\mu}_{y}=-\sqrt{2}{Q}_{11}^{i},\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}{\mu}_{z}={Q}_{10}.$$

(27)

For *l* = 2, the traceless Cartesian quadrupoles ${\Theta}_{\alpha \beta}^{TL}$ are related to *Q*_{2m} by

$$\begin{array}{cc}{\Theta}_{xx}^{TL}=\frac{1}{2}(-{Q}_{20}+\sqrt{6}{Q}_{22}^{r})\hfill & \hfill {\Theta}_{xz}^{TL}=-\sqrt{\frac{3}{2}}{Q}_{21}^{r}\hfill \\ {\Theta}_{yy}^{TL}=\frac{1}{2}(-{Q}_{20}-\sqrt{6}{Q}_{22}^{r})\hfill & \hfill {\Theta}_{yz}^{TL}=-\sqrt{\frac{3}{2}}{Q}_{21}^{i}\hfill \\ {\Theta}_{zz}^{TL}={Q}_{20}\hfill & \hfill {\Theta}_{xy}^{TL}=\sqrt{\frac{3}{2}}{Q}_{22}^{r}\hfill \end{array}$$

(28)

Note the trace of the quadrupoles is zero, i.e. $\mathrm{Tr}\left({\Theta}^{TL}\right)={\Theta}_{xx}^{TL}+{\Theta}_{yy}^{TL}+{\Theta}_{zz}^{TL}=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 03^{72} and then converted to their traceless forms^{25}. For example, the expression for converting Cartesian quadrupoles Θ_{αβ} into traceless Cartesian quadrupoles ${\Theta}_{\alpha \beta}^{TL}$ is given by

$${\Theta}_{\alpha \beta}^{TL}=\frac{3}{2}{\Theta}_{\alpha \beta}-\frac{1}{2}{\delta}_{\alpha \beta}\sum _{p}{\Theta}_{pp}.$$

(29)

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

The atomic Gaussian multipole moments *Q*_{lm} 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 *Q*_{lm} and exponent parameters λ (*N _{c}* = 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 ${Q}_{lm}^{\mathit{local}}$ are plotted as a function bond length and bond angle.

The exchange-overlap model^{61}^{,}^{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 *E _{exch}* 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

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 ${Q}_{lm}^{\mathit{local}}$ 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 *N _{c}* is studied. For

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 *N _{c}* = 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 (

Electrostatic energies (kcal/mol) for equilibrium hydrogen bonded dimers (X-water). The electrostatic energies predicted by Gaussian monopoles *E*_{M}, Gaussian dipoles *E*_{DM}, and Gaussian quadrupoles *E*_{QDM} (*N*_{c} = 4) are compared to their reference B3LYP/6-31G* **...**

A similar analysis is performed at the HF/aug-cc-pVTZ level. Gaussian multipoles (*N _{c}* = 8) are fit to the ESP calculated at the HF/aug-cc-pVTZ level, while the reference

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 (*N _{c}* = 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

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.

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 dimers^{82}^{,}^{83} which represent local minima on the water-water potential energy surface. Previously, GEM^{57} 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 (*N _{c}* = 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.

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 (*N _{c}* = 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

Inter-molecular density overlap integrals (multiply by 10^{−3} e^{2}/Å^{3}) for equilibrium hydrogen bonded dimers (X-water). The overlap integrals predicted by Gaussian monopoles *S*_{M}, Gaussian dipoles *S*_{DM}, and Gaussian quadrupoles *S*_{QDM} (*N*_{c} = 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 (*N _{c}* = 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

The inter-molecular density overlap integral (10^{−3} e^{2}/Å^{3}) calculated by atomic Gaussian quadrupoles and *ab initio* is plotted as a function of H..O distance (Å) for the water – water dimer.

The permanent *molecular* multipole moments up to hexadecapole are calculated for the *atomic* Gaussian monopole, dipole, and quadrupole models (*N _{c}* = 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

Non-zero components of the molecular quadrupole moment D-Å of ammonia for Gaussian monopoles, dipoles, and quadrupoles *(N*_{c} = 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.

The atomic Gaussian multipole moments are investigated as a function of bond length *r* and bond angle θ for the case of water. Gaussian quadrupoles (*N _{c}* = 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}_{lm}^{\mathit{local}}(r,\theta )\cong {Q}_{lm}^{\mathit{local},0}+({r}_{1}-{r}_{1}^{0})\frac{\partial {Q}_{lm}^{\mathit{local},0}}{\partial {r}_{1}}+({r}_{2}-{r}_{2}^{0})\frac{\partial {Q}_{lm}^{\mathit{local},0}}{\partial {r}_{2}}+(\theta -{\theta}^{0})\frac{\partial {Q}_{lm}^{\mathit{local},0}}{\partial \theta},$$

(29)

where ${Q}_{lm}^{\mathit{local},0}$ is the atomic Gaussian multipole at the geometry optimized equilibrium structure. $\partial {Q}_{lm}^{\mathit{local},0}\u2215\partial {r}_{1}$, $\partial {Q}_{lm}^{\mathit{local},0}\u2215\partial {r}_{2}$, and $\partial {Q}_{lm}^{\mathit{local},0}\u2215\partial \theta $ are the finite difference partial derivatives of ${Q}_{lm}^{\mathit{local}}$ with respect to *r*_{1}, *r*_{2}, and θ, respectively, evaluated at the equilibrium structure.

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 Å.

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 (*N _{c}* = 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

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 ${Q}_{lm}^{\mathit{local}}$ 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 model^{61}^{,}^{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.

Click here to view.^{(170K, pdf)}

Click here to view.^{(291K, pdf)}

Click here to view.^{(531K, text)}

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.

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 *N _{c}* 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 ${Q}_{lm}^{\mathit{local}}$ as a function of bond length and bond angle are given. For the exchange-overlap model, exchange-repulsion energies, molecular pair

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.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |