fields combining GEM with polarizable force fields. GEM* combines the Coulomb and exchange-repulsion terms from GEM with the polarization, van der Waals (modified) and bonded terms from AMOEBA. We have also introduced a multi-scale implementation denoted S/G-1 that combines GEM and SIBFA demonstrating the advantage of the incorporation of the electron density force field when applied to the description of metals in biomolecular systems. Overall, the GEM framework and its hybrid approaches offer a new strategy for molecular dynamics with enhanced accuracy at a cost that will enable the possibility of adequate sampling.

Introduction

The use of classical potentials for simulations of chemical and biochemical systems with molecular dynamics has been a field of intense research. Currently, it is possible to simulate systems with millions of atoms and millisecond time scales [START_REF] Schulten | Biomolecular modeling in the era of petascale computing[END_REF]; [START_REF] Shaw | Atomic-Level Characterization of the Structural Dynamics of Proteins[END_REF]]. With exa-scale computing, i.e. 10 18 floating point operations per second (FLOPs), on the horizon it is necessary to evaluate the performance of the current potentials. Indeed, long-time biomolecular simulations have revealed some issues already. For example Raval et al. carried out a study on 24 proteins (both homology models and experimental structures) used in recent CASP competitions involving at least 100µs MD simulations [START_REF] Raval | Refinement of protein structure homology models via long, all-atom molecular dynamics simulations[END_REF]]. For most systems, the structures drifted away from the native state, even when starting from the experimental structure. Although only two conventional force fields were employed, the authors concluded that this is most likely a limitation of the available point-charge force fields. As simulations on these and longer scales grow more widespread with improvements in computing power, node inter-connect, and graphical processing unit (GPU) hardware [START_REF] Stone | Accelerating molecular modeling applications with graphics processors[END_REF]], the accuracy of these classical potentials will be further tested.

In this context, there has been a recent impetus to develop more accurate force fields. One of the main thrusts has been to improve the description of the bonded interactions by including anharmonicity and of the non-bonded interactions by introducing ex-

Introduction 3
plicit polarization and a better description of the charge anisotropy terms. Several force fields that employ distributed multipoles and use explicit polarization (or are QM based using simpliefied MO schemes) have been proposed including AMOEBA, SIBFA, EFP, X-Pol, mDC and NEMO among others [START_REF] Hermida-Ramón | Inter-and intramolecular potential for the n-formylglycinamidewater system. a comparison between theoretical modeling and empirical force fields[END_REF]; [START_REF] Gresh | Anisotropic, polarizable molecular mechanics studies of inter-, intra-molecular interactios, and ligand-macromolecule complexes. a bottom-up strategy[END_REF]; [START_REF] Ponder | Current status of the amoeba polarizable force field[END_REF]; [START_REF] Day | An effective fragment method for modeling solvent effects in quantum mechanical calculations[END_REF]; [START_REF] Xie | The design of a next generation force field: The x-pol potential[END_REF]; [START_REF] Xie | X-pol potential: An electronic structure-based force field for molecular dynamics simulation of a solvated protein in water[END_REF]; [START_REF] Mills | Polarisable multipolar electrostatics from the machine learning method kriging: an application to alanine[END_REF]; [START_REF] Popelier | A generic force field based on quantum chemical topology[END_REF]; [START_REF] Giese | A variational linear-scaling framework to build practical, efficient next-generation orbital-based quantum force fields[END_REF]; [START_REF] Babin | Development of a first principles water potential with flexible monomers. ii: Trimer potential energy surface, third virial coefficient, and small clusters[END_REF]; [START_REF] Giese | Parametrization of an orbital-based linear-scaling quantum force field for noncovalent interactions[END_REF]]. The use of distributed multipoles results in an improved description of the charge density anisotropy and provides more accurate electrostatic interactions [START_REF] Stone | The theory of intermolecular forces[END_REF]; [START_REF] Price | Toward more accurate model intermolecular potentials for organic molecules[END_REF]; [START_REF] Popelier | Atoms in Molecules: An Introduction[END_REF]; [START_REF] Kosov | Atomic partitioning of molecular electrostatic potentials[END_REF]; Popelier et al. (2001a); Popelier and Kosov (2001); [START_REF] Mcdaniel | First-principles many-body force fields from the gas phase to liquid: A universal approach[END_REF]]. However, distributed multipoles suffer from one drawback since they cannot describe the overlap of charge density as two molecules get close to each other. This is known as the penetration effect [START_REF] Stone | The theory of intermolecular forces[END_REF]; [START_REF] Freitag | Evaluation of charge penetration between distributed multipolar expansions[END_REF]]. It is possible to reduce the penetration error by employing empirical damping functions [START_REF] Kairys | Evaluation of the charge penetration energy between non-orthogonal molecular orbitals using the spherical gaussian overlap approximation[END_REF]; [START_REF] Freitag | Evaluation of charge penetration between distributed multipolar expansions[END_REF]; [START_REF] Piquemal | Improved formulas for the calculation of the electrostatic contribution to the intermolecular interaction energy from multipolar expansion of the electronic distribution[END_REF]; [START_REF] Cisneros | Simple formulas for improved point-charge electrostatics in classical force fields and hybrid quantum mechanical/molecular mechanical embedding[END_REF]; [START_REF] Wang | Including charge penetration effects in molecular modeling[END_REF]; [START_REF] Stone | Electrostatic damping functions and the penetration energy[END_REF]]. It is also possible to include this effect via the use of neural networks [START_REF] Handley | Potential energy surfaces fitted by artificial neural networks[END_REF]].

Another possibility to avoid the charge penetration and anisotropy shortcomings is to use a continuous description of the molecular charge density. Several methods that describe the electronic distribution explicitly have been proposed [START_REF] Wheatley | Gaussian multipole functions for describing molecular charge distributions[END_REF]; [START_REF] Gavezzotti | Calculation of intermolecular interaction energies by direct numerical integration over electron densities i. electrostatic and polarization energies in molecular crystals[END_REF]; [START_REF] Eckhardt | Computer simulations and analysis of structural and energetic features of some crystalline energetic materials[END_REF]; [START_REF] Volkov | Calculation of electrostatic interaction energies in molecular dimers from atomic multipole moments obtained by different methods of electron density partitioning[END_REF]; [START_REF] Coppens | The interplay between experiment and theory in charge-density analysis[END_REF]; [START_REF] Paricaud | From dimer to condensed phases at extreme conditions: Accurate predictions of the properties of water by a gaussian charge polarizable model[END_REF]]. We have introduced the Gaussian Electrostatic Model (GEM) [Cisneros et al. (2005a); Piquemal et al. (2006a); [START_REF] Cisneros | Generalization of the gaussian electrostatic model: Extension to arbitrary angular momentum, distributed multipoles, and speedup with reciprocal space methods[END_REF]]. GEM uses density fitting (DF) techniques [Boys and Shavit (1959); [START_REF] Dunlap | On first-row diatomic molecules and local density models[END_REF]; [START_REF] Köster | Scalable evaluation of polarization energy and associated forces in polarizable molecular dynamics: I. toward massively parallel direct space computations[END_REF]] to reproduce the molecular electronic density using Hermite Gaussian auxiliary basis sets (ABSs). These fitted densities are employed to calculate each intermolecular component as obtained from energydecomposition (EDA) prodecures. The reason for the use of EDA methods for the parametrization of GEM is that it enables the separation of each of the components of the intermolecular inter-actions. There are a variety of EDA approaches that can be employed including symmetry-adapted perturbation theory (SAPT), Kitaura-Morokuma (KM), restricted variational space (RVS), constrained space orbital variations (CSOV), to name a few [START_REF] Eisenschitz | Perturbation theory[END_REF]; Hirshfelder (1967a,b); [START_REF] Murrel | Intermolecular forces in the region of small orbital overlap[END_REF]; [START_REF] Kitaura | A new energy decomposition scheme for molecular interactions within the hartree-fock approximation[END_REF]; [START_REF] Bagus | A new analysis of charge transfer and polarization for ligand-metal bonding: model studies for al4co and al4nh3[END_REF]; [START_REF] Stevens | Frozen fragment reduced variational space analysis of hydrogen bonding interactions. applications to the water dimer[END_REF]; [START_REF] Jeziorski | Perturbation theory approach to intermolecular potential energy surfaces of van der waals complexes[END_REF]; [START_REF] Glendening | Natural energy decomposition analysis: An energy partitioning procedure for molecular interactions with application to weak hydrogen bonding strong ionic, and moderate donor-acceptor interactions[END_REF]; [START_REF] Glendening | Natural energy decomposition analysis: explicit evaluation of electrostatic and polarization effects with application to aqueous clusters of alkali metal cations and neutrals[END_REF]; [START_REF] Mo | Energy decomposition analysis of intermolecular interactions using a block-localized wave function approach[END_REF]; [START_REF] Heßelmann | Dft-sapt with density fitting: a new efficient method to study intermolecular interaction energies[END_REF]; [START_REF] Piquemal | A csov study of the difference between hf and dft intermolecular interaction energy values: the importance of the charge transfer contribution[END_REF]; [START_REF] Khaliullin | An efficient self-consistent field method for large systems of weakly interacting components[END_REF]; [START_REF] Wu | Density-based energy decomposition analysis for intermolecular interactions with variationally determined intermediate state energies[END_REF]; [START_REF] Lu | Directional dependence of hydrogen bonds: A density-based energy decomposition analysis and its implications on force field development[END_REF].

In this contribution we present the theory behind the GEM method and recent advances and results on the application of two hybrid GEM potentials. In section 1.3 we provide a brief review of the analytical and numerical density fitting methods and it's implementation including the methods employed to control numerical instabilities. This is followed by a review of the procedure to obtain distributed site multipoles from the fitted Hermite coefficients in section 1.4. Section 1.5 describes the extension of reciprocal space methods for continuous densities. Section 1.6 describes the complete form for GEM and a novel hybrid force field, GEM*, which combines term from GEM and AMOEBA for MD simulations. Finally, section 1.7 describes the implementation and initial applications of a multi-scale program that combines GEM and SIBFA.

Density fitting methods

The use of ABSs for density fitting is a field of intense study. This method relies on the use of auxiliary basis functions (ABS), generally Gaussians, to expand the molecular electron density

ρ(r) = k c k Λ(r).
(1.1)

For GEM the ABSs consist of Hermite Gaussians, Λ tuv (r). The expansion coefficients c k for the approximate density ρ may be obtained by minimizing Eq. 1.2 using some metric Ô [START_REF] Dunlap | On first-row diatomic molecules and local density models[END_REF]; [START_REF] Eichkorn | Auxiliary basis sets to approximate coulomb potentials[END_REF]; [START_REF] Köster | Efficient recursive computation of molecular integrals for density functional methods[END_REF]; [START_REF] Köster | Scalable evaluation of polarization energy and associated forces in polarizable molecular dynamics: I. toward massively parallel direct space computations[END_REF]].

Density fitting methods 5

E self =< ρ(r) -ρ(r)| Ô|ρ(r) -ρ(r) > (1.2)
Several operators Ô can be employed including the overlap operator Ô = 1, the Coulomb operator Ô = 1/r or the damped Coulomb operator Ô = erf c(βr)/r [START_REF] Jung | Auxiliary basis expansions for large scale electronic structure calculations[END_REF]]. The minimization of Eq. 1.2 with respect to the expansion coefficients c k leads to a linear system of equations:

∂E self ∂c l = - µ,ν P µν < µν| Ô|l > + k c k < k| Ô|l > (1.3)
The solution of Eq. 1.3 requires the inversion of a the ABS matrix G =< k| Ô|l >. In principle this matrix should be positive definite and symmetric. In practice however, this matrix is almost singular and therefore the diagonalization to obtain its inverse must be done with care. To this end we have explored analytical and numerical procedures to obtain G and G -1 .

Analytical Fitting

The analytical fitting procedure involves the explicit evaluation of all the matrix elements of G and its subsequent inversion, which is achieved by diagonalization. We have implemented several methods for the diagonalization step. Initially we employed singular value decomposition (SVD) [START_REF] Press | Numerical recipes in fortran77; the art of scientific computing[END_REF]] by setting the inverse of the eigenvalue to zero if it is below a certain cutoff. However, this method produces undesirable numerical instabilities (noise) when the number of basis functions starts to grow as we and others have discussed previously [Cisneros et al. (2005b); [START_REF] Podeszwa | Density-fitting method in symmetry-adapted perturbation theory based on kohnsham description of monomers[END_REF]].

In the current implementation we employ the Tikhonov regularization method [START_REF] Press | Numerical recipes in fortran77; the art of scientific computing[END_REF]]. This approach is similar to the constrained density fitting algorithm of Misquitta and Stone [START_REF] Misquitta | Distributed polarizabilities obtained using a constrained density-fitting algorithm[END_REF]]. Here, the redundant basis set contributions are penalized by minimizing E self + λ k x 2 k , resulting in a more stable diagonalization procedure [Cisneros et al. (2006a)].

For problematic systems we also implemented the damped Coulomb operator Ô = erf c(βr)/r proposed by [START_REF] Jung | Auxiliary basis expansions for large scale electronic structure calculations[END_REF]] to attenuate the near singular behavior due to long range interactions present in G. In our studies we observed that noise is still a problem for some systems. The noise in the fit is known to arise due to the attempt of the ABSs to fit the density at the nuclear cores [Cisneros et al. (2005b[Cisneros et al. ( , 2006a))].

The use of Hermite Gaussians with angular moment greater than 0 requires the rotation of the fitting coefficients. In GEM this is addressed in a similar manner to multipolar force fields by defining a global molecular frame and a reference (local) site frame [START_REF] Toukmaji | Efficient PME-based approach to fixed and induced dipolar interactions[END_REF]; [START_REF] Sagui | Towards an accurate representation of electrostatics in classical force fields: Efficient implementation of multipolar interactions in biomolecular simulations[END_REF]; Cisneros et al. (2006a)]. The use of Hermite Gaussians provides a straightforward solution to the rotation since they are defined by partial derivatives of a spherical Gaussian which can be taken either with respect to the local (reference) frame or with respect to the global coordinates [Cisneros et al. (2006a)]. Moreover, the rotation frames are the same for the distributed multipoles.

Numerical fitting

As mentioned above, the numerical instabilities in the fit arise from the attempt to fit the nuclear cores. Thus, if the density at the cores is discarded then the fit should become more stable. This can be achieved by using numerical grids to evaluate a given molecular property and discarding points at and near the core. This can be achieved by minimizing the following fitting function:

χ 2 = i W (r i ) (y(r i ) -ỹ(r i , c k )) 2 , (1.4)
where y(r i ) denotes the ab initio molecular property of interest at point i and ỹ(r i , c k ) is the same property evaluated with the kth ABS element at the same point on the grid. Finally, W (r) is the weighting function for the point on the grid, which can be defined in several ways [START_REF] Bayly | A well-behaved electrostatic potential base method using charge restraints for deriving atomic charges: the resp method[END_REF]; [START_REF] Hu | Fitting molecular electrostatic potentials from quantum mechanical calculations[END_REF]]. Hu et al have proposed a weighting function that provides a smooth cutoff near the cores and at long distance to avoid any discontinuities [START_REF] Hu | Fitting molecular electrostatic potentials from quantum mechanical calculations[END_REF]]:

W (r i ) = exp[-σ(log ρ promol (r) -log ρ ref ) 2 ], (1.5)
where ρ promol is a reference promolecular atomic density, and σ and ρ ref are adjustable parameters. It has been shown that the surface for σ and ρ ref is relatively flat [START_REF] Hu | Fitting molecular electrostatic potentials from quantum mechanical calculations[END_REF]; [START_REF] Elking | Gaussian multipole model (gmm)[END_REF]]. We have implemented a modified version of this weighting function previously [START_REF] Elking | Gaussian multipole model (gmm)[END_REF]]. The main differences between the original Hu et al. weight and our implementation are the re-optimization of the promolecular atomic electron densities at the MP2/aug-cc-pVQZ level and the values for σ and ρ ref which correspond to 0.42 and -7.0 respectively. The minimization of Equation 1.5 leads to a linear system of equations that can be expressed as: c -c 0 = -H 0,-1 g 0 . As was the case for the analytic DF, we employ Tikhonov regularization for the inversion of the Hessian that arises for the linear-least-squares procedure. In our initial implementation of numerical fitting we explored different molecular properties including electronic density, molecular electrostatic potential (mESP), and the three components of the electric field [START_REF] Cisneros | Numerical fitting of molecular properties to hermite gaussians[END_REF]]. All the properties were gridded on rectangular grids. Subsequently we showed that the use of spherical molecular grids based on the scheme proposed by Becke [START_REF] Becke | A multicenter numerical integration scheme for polyatomic molecules[END_REF]] significantly reduce the number of fitting points [START_REF] Elking | Gaussian multipole model (gmm)[END_REF]; [START_REF] Cisneros | Application of gaussian electrostatic model (gem) distributed multipoles in the amoeba force field[END_REF]].

Distributed multipoles

In this subsection we present the methodology to obtain Cartesian point multipoles from the Hermite coefficients obtained in the fitting procedure. In all our work we have purposefully employed ABSs with a maximum angular momentum of 2, which results in distributed multipoles only up to quadrupoles. This ensures that the distributed multipoles obtained can be directly employed in the AMOEBA force field. However, higher order multipoles can be obtained if an ABS with higher angular momentum is used.

Briefly, we have expanded on the work by Challacombe et al., who have shown that Hermite Gaussians have a simple relation to elements of the Cartesian multipole tensor [START_REF] Challacombe | Modern developments in Hartree-Fock theory: Fast methods for computing the Coulomb matrix[END_REF]]. Once the Hermite coefficients have been determined, they may be employed to calculate point multipoles centered at the ex-pansion sites. Thus, if hc tuv represents the coefficient of a Hermite Gaussian of order Λ tuv , then if this Hermite is normalized we have:

hc 000 Λ 0 dr = hc 000 (1.6)
This guarantees that higher order multipole integrals will integrate to integer numbers, for example, for the dipole integral in the z direction, d z :

hc 001 zΛ 001 dxdydz = hc 001 z ∂ ∂S z Λ 0 dz = -hc 001 z ∂ ∂z Λ 0 dz = hc 001 (1.7)
For quadrupole and higher order integrals the same relationships hold, although different cases need to be considered (see ref. [Cisneros et al. (2006a)]). In practice, following Stone's definition [START_REF] Stone | The theory of intermolecular forces[END_REF]], we have used traceless quadrupoles. Furthermore, the use of GEM distributed multipoles (GEM-DM) for multipolar force fields provides a straightforward way to determine the penetration error in the site-site Coulomb interaction energy due to the connection with the GEM Hermites. Thus, this connection provides a natural way to generate damping functions to lessen the penetration error [START_REF] Piquemal | Improved formulas for the calculation of the electrostatic contribution to the intermolecular interaction energy from multipolar expansion of the electronic distribution[END_REF]; [START_REF] Cisneros | Simple formulas for improved point-charge electrostatics in classical force fields and hybrid quantum mechanical/molecular mechanical embedding[END_REF]].

A further advantage of this approach to distributed multipoles is that, unlike some conventional multipole expansions [START_REF] Stone | Distributed multipole analysis: Stability for large basis sets[END_REF]; Popelier et al. (2001b)], the (spherical) multipole expansion obtained from Hermite Gaussians in this way is intrinsically finite of order t + u + v (i.e. the highest angular momentum in the ABS) as shown in [Cisneros et al. (2006a)], similar to the multipoles obtained by Volkov and Coppens [START_REF] Volkov | Calculation of electrostatic interaction energies in molecular dimers from atomic multipole moments obtained by different methods of electron density partitioning[END_REF]].

Reciprocal space methods for integral evaluation

The use of molecular densities results in the need to compute a large number of two center integrals for the intermolecular interaction. A significant computational speedup can be achieved by using reciprocal space methods based on Ewald sums. In this way, the integrals are calculated in direct or reciprocal space depending on the exponent of the Gaussian Hermites.

Here we describe how the Ewald formalism can be extended to take into account the Gaussian distribution. Let U denote a unit cell whose edges are given by the vectors a 1 , a 2 and a 3 . An idealized infinite crystal can be generated by all periodic translations n = n 1 a 1 + n 2 a 2 + n 3 a 3 for all integer triples (n 1 , n 2 , n 3 ), with n 1 , n 2 , n 3 not all zero. Now, consider a collection of N normalized spherical Gaussian charge densities ρ 1 . . . ρ N (e.g. GEM-0) centered at {R 1 . . . R N } ∈ U with exponents α i , i.e. ρ i (r) = q i (α i /π) 3/2 exp(-α i (r -R i ) 2 ), and let q 1 + . . . + q N = 0. Note that N need not be limited only to atomic positions, e.g., GEM-0 includes sites on the oxygen lone pairs and the bisector line between the two hydrogens [Piquemal et al. (2006a)]. The Coulomb energy of the central unit cell within a large spherical crystal, due to the interactions of the Gaussian charge distributions with each other and all periodic images within the crystal can be calculated using Ewald methods.

In particular, to determine the reciprocal part in the Ewald sum it is necessary to grid the Gaussian densities. However, this can become intractable for Gaussian functions with large exponents (compact Gaussians). In the initial implementation the charge densities were classified into compact or diffuse Hermite Gaussians based on a given Ewald exponent β. Therefore, if the exponent of a given Hermite was above the cutoff it was considered compact, and diffuse (α i < β) otherwise. With this, the contributions involving diffuse Hermites can be calculated in reciprocal space exclusively [START_REF] Cisneros | Generalization of the gaussian electrostatic model: Extension to arbitrary angular momentum, distributed multipoles, and speedup with reciprocal space methods[END_REF]].

This was later improved by the realization that the Ewald exponent, β may be different for each pair ij [START_REF] Darden | Dual bases in crystallographic computing[END_REF]]. Thus, β is chosen to be infinite for ij pairs where at least one of the Gaussians is diffuse. In this way, all pairs that involve diffuse Hermites are evaluated in reciprocal space. That is, given θ > 0 a Gaussian distribution q i ρ i is classified as compact (i ∈ c) if α i ≥ 2θ and diffuse (i ∈ d) otherwise. Subsequently, for i, j ∈ {c}, select β so that 1/θ = 1/α i + 1/α j + 1/β, otherwise β = ∞. In the case of GEM, the fitted densities are expanded in a linear combination of Hermite Gaussians Λ tuv (r, α, R). Thus, the charge distribution is given by

ρ i (r, α, R) = L l=1
tuv c i,l,tuv Λ tuv (r, α l , R i ), where c i,l,tuv are the Hermite coefficients and L denotes the different ABS exponents on center i. With this, the Ewald expression becomes

E(ρ {N } ) = 1 2 n ′ N i=1 li∈c tiuivi c i,li,tiuivi N j=1 lj ∈c tj uj vj (-1) (tj +uj +vj ) c j,lj,tj uj vj × ∂ ∂R ijx ti+tj ∂ ∂R ijy ui+uj ∂ ∂R ijz vi+vj ×    erf c(θ 1/2 |R ij -n|) -erf c(µ 1/2 lilj |R ij -n|) |R ij -n|    + 1 2πV m =0 1 m 2 exp(-π 2 m 2 /2θ) l1∈c S l1 (m) ×exp(-π 2 m 2 /2θ) l2∈c S l2 (-m) (1.8) + 1 2πV m =0 1 m 2 N (l1,l2) / ∈c×c exp(-π 2 m 2 /α l1 ) ×exp(-π 2 m 2 /α l2 )S l1 (m)S l2 (-m) - π 2V N i=1 N j=1 l1∈c l2∈c c i,l1,000 c j,l2,000 1 θ - 1 α l1 - 1 α l2 - N i=1 E self (ρ i ) + 2πD 2 3V + ε(K),
where the first term corresponds to the direct part of the Ewald sum, the second and third terms to the reciprocal part, R ij = R i -R j , the structure factors S l (m) involve derivatives of the Fourier exponential with respect to the Hermite centers, E self (ρ i ) is the correction due to the self energy of each Hermite interacting with its replicate, the term involving the unit cell dipole D = q 1 R 1 +. . .+ q N R N is the surface term, ε(K) denotes a quantity that converges to 0 as K -→ ∞, m denotes the reciprocal lattice vectors, and 1/µ lilj = 1/α li + 1/α lj [START_REF] Cisneros | Generalization of the gaussian electrostatic model: Extension to arbitrary angular momentum, distributed multipoles, and speedup with reciprocal space methods[END_REF]; [START_REF] Darden | Dual bases in crystallographic computing[END_REF]]..

Since the ABSs include Hermites with l > 0, the direct space contributions can be efficiently evaluated by using the McMurchie-Davidson (MD) recursion [START_REF] Mcmurchie | One-and two-electron integrals over cartesian gaussian functions[END_REF]]. This recursion has been used to calculate the required erfc and higher derivatives for multipole interactions [START_REF] Sagui | Towards an accurate representation of electrostatics in classical force fields: Efficient implementation of multipolar interactions in biomolecular simulations[END_REF]]. This approach was also employed for the Hermite Gaussians [START_REF] Cisneros | Generalization of the gaussian electrostatic model: Extension to arbitrary angular momentum, distributed multipoles, and speedup with reciprocal space methods[END_REF]], where it was shown that the MD recursion is applicable to other types of operators besides 1/r. For the reciprocal sums three methods were implemented: full [START_REF] Ewald | Die Berechnung optischer und elektrostatischer Gitterpotentiale[END_REF]], sPME [START_REF] Essmann | A smooth particle mesh Ewald method[END_REF]] and FFP [START_REF] York | The Fast Fourier Poisson (FFP) method for calculationg Ewald sums[END_REF]].

The latter two methods rely on the use of fast fourier transforms to approximate the structure factors that arise in the reciprocal term, which results in the efficient evaluation of this term and has been shown to scale as O(N logN ) for sPME [START_REF] Essmann | A smooth particle mesh Ewald method[END_REF]].

The GEM and GEM* force fields

The initial implementation of the full GEM potential involved the use of spherical type Hermites only, resulting in what was termed GEM-0 [START_REF] Cisneros | Generalization of the gaussian electrostatic model: Extension to arbitrary angular momentum, distributed multipoles, and speedup with reciprocal space methods[END_REF]]. This initial parametrization included the terms described below.

The GEM functional form

The idea for GEM is to employ the fitted Hermite Gaussians to evaluate each term in

E T otal = E GEM Coulomb + E GEM exch-rep +E GEM polarization + E GEM charge-transf er , (1.9)
where the Coulomb term is given by

E GEM Coulomb = A>B Z A Z B r AB + A>B Z A ρB (r B ) r AB dr + A>B Z B ρA (r A ) r AB dr + A>B ρA (r A )ρ B (r B ) r AB dr.
(1.10)

The exchange-repulsion term is calculated by means of the charge density overlap following the Wheatley-Price overlap model [START_REF] Wheatley | An overlap model for estimating the anisotropy of repulsion[END_REF]; [START_REF] Domene | Overlap-model and ab initio cluster calculations of ion properties in distorted environments[END_REF]]:

E GEM exch-rep = K A>B ρA (r A )ρ B (r B ). (1.11)
The polarization term is approximated by the use of dipole polarizabilities, which yield a very good results for the polarization energies (if the electric fields are not large) [START_REF] Böttcher | Theory of Electric Polarization[END_REF]]. To this end, the electric fields are calculated with the fitted densities and interacted with distributed dipolar polarizabilidies with the Garmer and Steven's approach [START_REF] Day | An effective fragment method for modeling solvent effects in quantum mechanical calculations[END_REF]]

E GEM polarization = 1 2 xyz j ∆µ(i)(γE 0 (j)), (1.12)
where ∆µ(i) = α(i) xyz j E(∆µ(i)) + (γE 0 (j)), and γ is a scaling factor for the permanent electric fields [START_REF] Cisneros | Generalization of the gaussian electrostatic model: Extension to arbitrary angular momentum, distributed multipoles, and speedup with reciprocal space methods[END_REF]].

Finally, the charge transfer term is evaluated using the semiempirical formalism implemented in the SIBFA force field [START_REF] Gresh | Sibfa[END_REF]; [START_REF] Piquemal | Toward a separate reproduction of the contributions to the hartree-fock and dft intermolecular energies by polarizable molecular mechanics with the sibfa potential[END_REF]]:

E GEM charge-transf er = 2C Lα (I * αβ ) 2 ∆E * αβ , (1.13)
where C = 3.5 is a constant parametrized to reproduce the value of the charge-transfer energy (obtained with CSOV) for the canonical water dimer at equilibrium distance. I * αβ is a function

The GEM and GEM* force fields 13

of the overlap between the localized molecular orbital (LMO) for the donor lone pair and antibonding LMO of the acceptor, as well as the electrostatic potential on site A arising from all other sites obtained with the GEM densities. ∆E * αβ is a difference between the ionization potential of A and the electron affinity of the acceptor site.

In our initial implementations of GEM-0 and GEM we have not introduced an explicit term for the dispersion interactions. This is because these force fields have been originally parametrized using the CSOV method at the DFT level, which, by definition, does not include a dispersion contribution. However, it is possible to include this term in a similar way to the SIBFA potential [START_REF] Gresh | Sibfa[END_REF]; [START_REF] Piquemal | Toward a separate reproduction of the contributions to the hartree-fock and dft intermolecular energies by polarizable molecular mechanics with the sibfa potential[END_REF]].

GEM*: molecular dynamics with fitted densities

After our implementation of GEM-0, we extended the Coulomb and exchange-repulsion terms to enable the use or arbitrary angular momentum Hermites [START_REF] Cisneros | Generalization of the gaussian electrostatic model: Extension to arbitrary angular momentum, distributed multipoles, and speedup with reciprocal space methods[END_REF]]. However, both implementations only enabled energy calculations. In order to carry out MD simulations it is necessary to evaluate the associated forces efficiently. Until recently, it was impractical to do this since the analitical form of the force for the charge-transfer term was unavailable.

To enable the performance of MD simulations, a hybrid force field called GEM* was developed. GEM* combines the Coulomb and exchange-repulsion terms from GEM with the polarization, van der Waals (modified) and bonded terms from AMOEBA. The functional form for GEM* is thus:

E T otal = E GEM Coulomb + E GEM exch-rep + E AMOEBA polarization +E AMOEBA V dW + E AMOEBA bonded .
(1.14)

The Coulomb and exchange-repulsion terms for GEM* are evaluated with the same expressions as for GEM (Eqns. 1.10 and 1.11). Since GEM* includes an explicit term for exchange, it was necessary to modify the original van der Waals function implemented in AMOEBA. In this case, we have modified the buffered Halgren function (modHalgren) by removing the repulsive term as follows:

E modHalgren = -ǫ ij 1.07R * ij (R ij + 0.07R * ij ) 7
.

(1.15)

The polarization and bonded terms are the same as those in the original AMOEBA implementation [START_REF] Ren | Polarizable atomic multipole water model for molecular mechanics simulation[END_REF]]. The only difference in the polarization is that the permanent electric fields for the calculation of the induced dipoles are calcualted with the distributed multipoles obtained from the fitted Hermites for consistency between the Coulomb and polarization terms [START_REF] Cisneros | Application of gaussian electrostatic model (gem) distributed multipoles in the amoeba force field[END_REF]].

The initial implementation of GEM* was tested by fitting parameters for a water potential [START_REF] Duke | Gem*: A molecular electronic density-based force field for classical molecular dynamics simulations[END_REF]]. These parameters were compared to reference ab initio values for total intermolecular interactions corrected for basis set superposition error via the counterpoise correction . The reference data was calculated at the MP2(full)/aug-cc-pVTZ level to match the original AMOEBA parametrization [Ren andPonder (2003, 2002); [START_REF] Ponder | Current status of the amoeba polarizable force field[END_REF]; [START_REF] Ren | Polarizable atomic multipolebased molecular mechanics for organic molecules[END_REF]]. The molecular density used to obtain the fitting coefficients for GEM* was calculated at the same level of theory as above for a water molecule at the AMOEBA equilibrium geometry.

Three parametrizations were investigated, termed models 1-3 in the discussion below. The difference among the three models involves the use of different ABSs, A1 or A2 [START_REF] Andzelm | Density functional gaussian-typeorbital approach to molecular geometries, vibrations and reaction energies[END_REF]; [START_REF] Godbout | DGauss Version 2[END_REF]], and/or the dataset of water oligomers used for the parametrization. Model 1 was fitted using the A2 ABS to reproduce intermolecular interaction energies for the canonical water dimer (see Figure 1.1.), several random dimers, and selected water clusters from [START_REF] Temelso | Benchmark structures and binding energies of small water clusters with anharmonicity corrections[END_REF]]. Models 2 and 3 were parametrized to reproduce intermolecular energies for the canonical water dimer only using the A2 (model 2) and A1 (model 3) ABSs. All calculations for GEM* were performed with a modified pmemd version in the AMOEBA suite of programs [START_REF] Case | The amber biomolecular simulation programs[END_REF]].

The GEM and GEM* force fields 15

The comparison of the QM reference for energies and forces calculated with GEM* models 2 and 3 showed that both models reproduce the total intermolecular interactions well. However, model 3 deviated in the forces due to the limited accuracy because of the small number of Gaussians in the ABS employed for model 3. In addition, both models were observed to produce significant errors in the intermolecular energies for random dimers and binding energies for different oligomers. These errors are due to the improper description of these two models to describe interactions between H atoms since both models were fitted only to the canonical water dimer.

Conversely, all results for model 1 showed good agreement with the QM references for dimers as well as larger clusters. The results show that a better parametrization can be obtained once a slightly larger data set that included different dimer orientations was considered. Recently Babin et al. have developed a novel water model parametrized only from QM data [START_REF] Babin | Toward a universal water model: First principles simulations from the dimer to the liquid phase[END_REF][START_REF] Babin | Development of a first principles water potential with flexible monomers: Dimer potential energy surface, vrt spectrum, and second virial coefficient[END_REF]] using results from 40000 dimers calculated at the CCSD/CBS level.

The performance of GEM* was tested by performing 100 MD steps in the NVE ensemble with a series of water boxes of increasing size (216, 512, 1024, 2048 and 4096 molecules). All MD calculations were done on a single Xeon X5550 CPU with 12 GB of memory at 300 K with an 8 Å cutoff for Van der Waals interactions, using the Beeman integrator, a 1 fs time step and a dipole tolerance for the SCF convergence of 10 -6 . The calculation of the polarization with the induced dipoles was performed using the PME method with a B-spline order of 5 and a grid-size of 24. During the parametrization and testing, it was realized that the overlap integrals for both sets of ABSs employed (A1 and A2) tend to 0 at distances greater than 5 Å. Therefore, although it is also possible to perform the Overlap integrals in reciprocal space as described in [START_REF] Cisneros | Generalization of the gaussian electrostatic model: Extension to arbitrary angular momentum, distributed multipoles, and speedup with reciprocal space methods[END_REF]; [START_REF] Darden | Dual bases in crystallographic computing[END_REF]], the exchange integrals were evaluated only in direct space with the same cutoff as the van der Waals interactions (8 Å), or with a reduced cutoff of 6 Å.

Timings for all the tested systems are shown in Figure 1.2.. For comparison we performed the calculations for all cases with the Coulomb integrals evaluated completely in direct space (all Gaus-sians set as compact), or by placing two Gaussians in the diffuse set and employing sPME or FFP with two different exchange cutoffs. As discussed previously, the evaluation of the Hermite Gaussians in reciprocal space requires significantly larger grids and B-spline orders [START_REF] Cisneros | Generalization of the gaussian electrostatic model: Extension to arbitrary angular momentum, distributed multipoles, and speedup with reciprocal space methods[END_REF]; [START_REF] Darden | Dual bases in crystallographic computing[END_REF]]. For the three smallest boxes the calculations are faster when all the integrals are evaluated in direct space. This is due to the large overhead for the FFTs due to the fine grids required for accurate evaluation of the energies and forces in reciprocal space. As the system grows in size, the calculation becomes faster by using the compact/diffuse density split method with PME. The smallest water box comprising 216 molecules takes 100 s for the evaluation of the 100 energy/force calculations including the evaluation of Coulomb and overlap integrals for 15120 basis functions. That is, our code is able to evaluate all Coulomb and overlap integrals for 15120 basis functions for a single step in 1 s. For the 4096 box, comprising 286720 basis functions, when the exchange cutoff is reduced to 6 Åthe time is 2363 s. Moreover, after an initial optimization of the code, this time is reduced to 789 s. This is only 4 times slower than the same 4096 system calculated with AMOEBA on the same CPU. We expect more gains in performance as the code is further optimized.

Combining SIBFA and GEM: S/G-1

In this initial implementation the direct coupling between GEM and SIBFA has been only performed at the induction level (polarization and charge transfer energies) in the spirit of QM/MM techniques [START_REF] Chaudret | Towards improved treatment of metal cations in polarizable molecular mechanics using the hybrid gaussian electrostatics / distributed multipoles gem/sibfa approach[END_REF]]. Indeed, the GEM equations for exchange-repulsion involve overlap integrals between densities of both interacting fragments. Thus a mixed S/G-1 scheme is not possible for this term since the overlap between a GEM density and SIBFA's multipoles would be zero. Therefore, in the present S/G-1 implementation, the electrostatic, exchange-repulsion and dispersion energies are computed at the sole SIBFA level and include SIBFA's short range corrections. In a forthcoming work the full multiscale implementation including full Gaussian electrostatics first order energy will be reported.

Combining SIBFA and GEM: S/G-1 17

Therefore in the first implementation, GEM is only used to compute the second-order polarization and charge-transfer contributions between the cation and its bound ligands. Finally, the dispersion equations are the same for both methods, as they do not depend on electric fields or potentials. As a proof of concept, we limited ourselves to a Hartree-Fock level parametrization (no dispersion) of the method.

Briefly, within S/G-1, the evaluation of the polarization and charge-transfer energies use the same formalism. The differences between the two levels of computations are linked to the level of calculation of electric fields and potentials (i.e. using electronic densities for GEM or using distributed multipoles for SIBFA). Indeed, differences between GEM and SIBFA electric fields and electrostatic potentials can arise at short distances since GEM is identical to an ab initio field. Both procedures converge to the same solutions upon increasing the distances when the multipolar approximation starts to be valid as GEM densities act as a continuous electrostatic model. Therefore, in a similar spirit as in QM/MM approaches, specific fragments can be defined so as to be handled with either the GEM density or the SIBFA multipoles. For the first implementation, the use of GEM densities were restricted to the metal cations whereas the rest of molecules were described using SIBFA.

In order to try to include the different previously discussed physical effects within a MM scheme, we show here some results focusing on the polarization contribution in the case of the Ca(II)-H2O complex.

Figure 1.3. displays four curves, namely the reference ab initio CSOV polarization contribution, the undamped full GEM polarization energy, the full GEM + damping approach, and results obtained upon computing the polarization energy obtained with the exact ab initio undamped field values extracted from a quantum mechanical Gaussian 09 computation. The damping procedure is identical to the one used by SIBFA and is detailed in the technical appendix in ref. [START_REF] Chaudret | Towards improved treatment of metal cations in polarizable molecular mechanics using the hybrid gaussian electrostatics / distributed multipoles gem/sibfa approach[END_REF]].

It is important to point out that the GEM fields alone, in spite of their quasi-perfect match with their ab initio counterparts, do not provide a good reproduction of ab initio results at short-range: a damping of the fields is required to gain accuracy at very short distances. Such conclusions are confirmed as our results obtained from undamped exact ab initio fields (i.e. computed using the original molecular orbitals) are basically identical to the undamped GEM results (see Figure I). This clearly shows that inclusion of short-range quantum effects, inherently present within QM and GEM fields, is not sufficient to reproduce the true polarization energy. This is because the final ab initio CSOV polarization energy embodies both penetration and exchange-polarization effects. The first quantity is present in GEM (as in QM), however, the exchange-polarization arises from the required orthogonalization of molecular orbitals of both Ca(II) and H2O fragments within the constrained self-consistent field procedure. Therefore, since GEM does not include this repulsive effect, the computed polarization energy is overestimated. A straightforward solution to the problem is to apply the exact same field damping procedure that is used for the SIBFA polarization contribution. As can be sene from Figure I, the GEM+damping approach accurately reproduces the CSOV reference by selectively including the different effects.

The initial implementation of the S/G-1 method has been developed to describe metal cations in a ligand environment. To this end, S/G-1 has been parametrized to model Zn(II) and Hg(II) cations. The Zn(II) parametrization was performed on a series of representative mono-ligand complexes and subsequently employed to calculate the polarization and charge-transfer energies for a series of multi-valent complexes as shown in Figure 1.4.. S/G-1 matches the anti-cooperative behavior of the polarization and charge-transfer (at the RVS level) energies. For example, for [Zn(H 2 O) 4/2 ] 2+ the charge-transfer energy from S/G-1(RVS) corresponds to -26.8(-28.7) kcal/mol. Correspondingly, the polarization energy for S/G-1 is -103.8 kcal/mol compared to -101.7 kcal/mol for its RVS counterpart.

As a test of the applicability of S/G-1 to large systems, we applied this method to a Zn(II) dependent system recently studied using the SIBFA procedure. To this end, the Zn(II)-alcohol dehydrogenase (ADH) active site [START_REF] De Courcy | Energy analysis of zn polycoordination in a metalloprotein environment and of the role of a neighboring aromatic residue. what is the impact of polarization?[END_REF]] was calculated, GEM was used to model the Zn(II) cation and the remaining system was treated with SIBFA. As can be seen in Figure 1.4., S/G-1 successfully reproduced the RVS values for a complicated hetero-polyligated complex [START_REF] Chaudret | Towards improved treatment of metal cations in polarizable molecular mechanics using the hybrid gaussian electrostatics / distributed multipoles gem/sibfa approach[END_REF]].

To further test the performance of S/G-1 to model heavy metals where relativistic and correlation effects are important , a calibration for Hg(II) was performed. The cation polarization energy requires two components. The first arises from the cation dipolar polarizability and depends on the electrostatic field that the cation is subjected to. The second component arises from the quadrupolar polarizability and depends upon the field gradient [START_REF] Buckingham | multipolar expansion[END_REF]]. The magnitude of the second component was found to be important in the case of some metal cations, such as Cu(I) and Hg(II), and this component had to be explicitly formulated in SIBFA [START_REF] Gresh | Modeling copper(i) complexes: Sibfa molecular mechanics versus ab initio energetics and geometrical arrangements[END_REF]]. However, although the reference quantum dipolar polarizability for Hg(II) can be easily obtained, its quadrupolar polarizability could not be derived by QM calculations using a small core pseudopotential. Therefore we resorted to theo available Cu(I) [START_REF] Buckingham | multipolar expansion[END_REF]] value as a starting point. The values of the polarization energy for Hg(II) were obtained from RVS calculations on the [Hg(H 2 O) 2 ] 2+ complex, with Hg(II) equidistant from the two water molecules. In this complex, the field undergone by Hg(II) is zero, but the gradient is non-zero. Subsequently, the Cu(I) quadrupolar polarizability was employed as a starting point for the Hg(II) and scaled to match the polarization energy of the complex.

Figure 1.5. shows the second-order polarization and chargetransfer energies and the total intermolecular interaction energies for a di-aquo mercury complex at different distances. Except at very short distances, the error between SIBFA/GEM with respect to the ab initio calculations are very small. The S/G-1 chargetransfer values are very close to one another. The main differences arise from the polarization contribution. Moreover, the use of the GEM densities for the calculation of the second-order components results in better agreement at short range with respect to RVS than the original SIBFA method. The close agreement found for the monoligated Hg(II) complex used for the parametrization is conserved in the polyligated complexes. This shows that the nonadditivity of both the polarization and charge-transfer components can be reproduced with S/G-1 for Hg(II) as it has been demonstrated for Zn(II).

Conclusion and Perspective

The present contribution reports the development of the ab initio GEM force field that uses Hermite Gaussian electrostatics to include short-range quantum mechanical effects into molecular mechanics. It also details the GEM* and S/G-1 approaches that blend together GEM with the AMOEBA and SIBFA polarizable potentials, both of which rely on distributed multipoles. As the capability to achieve high accuracy on the interaction energies is demonstrated, first application to molecular dynamics have been detailed as the potentialities of the GEM densities to accurately treat difficult systems such as metalloproteins have been exposed. In term of perspectives, all models should be able to be used directly in MD simulations as all the required gradients were recently coded. It should open the possibility of large-scale molecular dynamics using Gaussian Hermite functions as the models will benefit from recent advances in algorithmic and in hybrid MPI/OPEN-MP parallelism that use new scalable strategies with gains going from 2 to 3 orders in magnitude in time within the present framework [Lipparini et al. (2014)]. Overall, efforts will be devoted to propose a scalable integrate methodology incorporating both distributed Multipoles and Hermite Gaussian densities in popular packages such as Amber and Tinker.

Figure 1

 1 Figure 1.1. Total intermolecular interaction energy for the canonical water dimer calculated with the three GEM* models compared to ab initio and AMOEBA

Figure 1

 1 Figure1.3. Ab initio polarization energies (kcal/mol) for the Ca(II)-H2O complex, computed using the CSOV procedure (blue), polarization energy computed using: i) distributed polarizabilities+ab-initio fields (grey),ii) distributed polarizabilities +GEM with (green) and iii) without fields damping (red).

Figure 1

 1 Figure 1.4. Polarization (A), charge-transfer (B) and second order induction (C) (pol+c-t) energies in poly-ligated Zn(II) complexes calculated with RVS and S/G-1. The RVS charge-transfer and induction energies are BSSE corrected. The complexes are as follows: 1: [Zn(CH3S)3] -, 2: [Zn(CH3S)4] 2-, 3: [Zn(imidazole)3] 2+ , 4: [Zn(imidazole)4] 2+ , 5: [Zn(H2O)6] 2+ , 6: [Zn(H2O) 5/1 ] 2+ , 7: [Zn(H2O) 4/2 ] 2+ , 8: cluster model for alcohol dehydrogenase active site

Acknowledgments

This work was supported by Wayne State University. Computing time from Wayne State's C&IT is gratefully acknowledged. This work was also supported by the French state funds managed by CALSIMLAB and the ANR within the Investissements dAvenir program under referenceANR-11-IDEX-0004-02