Multipolar Ewald methods, 1: theory, accuracy, and performance.

The Ewald, Particle Mesh Ewald (PME), and Fast Fourier–Poisson (FFP) methods are developed for systems composed of spherical multipole moment expansions. A unified set of equations is derived that takes advantage of a spherical tensor gradient operator formalism in both real space and reciprocal space to allow extension to arbitrary multipole order. The implementation of these methods into a novel linear-scaling modified “divide-and-conquer” (mDC) quantum mechanical force field is discussed. The evaluation times and relative force errors are compared between the three methods, as a function of multipole expansion order. Timings and errors are also compared within the context of the quantum mechanical force field, which encounters primary errors related to the quality of reproducing electrostatic forces for a given density matrix and secondary errors resulting from the propagation of the approximate electrostatics into the self-consistent field procedure, which yields a converged, variational, but nonetheless approximate density matrix. Condensed-phase simulations of an mDC water model are performed with the multipolar PME method and compared to an electrostatic cutoff method, which is shown to artificially increase the density of water and heat of vaporization relative to full electrostatic treatment.

Introduction

It is well-known that the molecular modeling of many biological processes requires a rigorous treatment of the long-ranged electrostatic interactions.[1−5] Furthermore, recent years have seen the development of many promising next-generation force fields and fast ab initio methods that endeavor to generalize the electrostatic interactions to higher-order atomic multipoles.[6] Specific examples include the following: polarizable force fields, such as AMOEBA;[7−12] density-based force fields, such as GEM[13−16] and S/G-1;[17] force fields that model the perturbative or many-body expansion of the energy;[18−21] hybrid quantum mechanical/molecular mechanical (QM/MM) methods;[22] and molecular orbital-based quantum mechanical force fields (QMFFs), such as the embedded fragment model,[23] X-Pol,[24−29] and the closely related modified “divide-and-conquer” (mDC) methods.[30−32] The use of higher-order atomic multipoles in these examples, and many other models, offer the promise of improved accuracy, but at a larger computational cost.[33−35] In order to elucidate the strengths and weaknesses of proposed models, it is necessary to make comparison to experiment, which often requires their application within molecular simulations under periodic boundary conditions. At the same time, these applications provide empirical measures of the preliminary model’s true computational cost. Therefore, fast electrostatic algorithms must be developed that are equipped to handle generalized charge densities, so that new models can be applied and tested.

The Ewald,[36] Particle Mesh Ewald[37−40] (PME), and Fast Fourier–Poisson[41] (FFP) methods are three algorithms used to evaluate the long-ranged electrostatic interactions of periodic systems. These methods were originally designed to accommodate systems consisting of point charges; therefore, modifications to their original formulations are required to apply them to models employing multipolar charge densities. In the present work, we derive a unified set of equations for the Ewald, PME, and FFP methods using a spherical tensor gradient operator formalism that extend these methods to arbitrary multipole order. Expressions for the energies, forces, and the generalized multipolar potentials required to incorporate the electrostatics into the QMFF self-consistent field (SCF) procedure are provided. Previous works have extended these methods within the framework of Cartesian point multipoles and Cartesian Gaussians,[42−47] as opposed to the solid harmonic multipoles considered here. Other related works were not formulated to arbitrary multipole order[48,49] or led to a formulism that did not prove to be computationally efficient.[50] The efficiency of our formulation is demonstrated through QMFF molecular dynamics applications described below and extended in Part 2 of this series.[51] During preparation of this manuscript, Simmonett[52] has independently reported an extension of PME for use with solid harmonic multipoles that shares many characteristics with the PME method presented here; however, our analysis extends significantly beyond ref (52) by comparing the accuracy and performance of Ewald, PME, and FFP methods as stand-alone electrostatic methods and upon implementation within a linear-scaling QMFF.

The paper is organized as follows. Section 2 describes the generalized charge density; derives equations for the generalized Ewald, PME, and FFP methods; and discusses net charge corrections, error analysis, and the integration of the methods into the mDC QMFF. Section 3 compares the accuracy and computational cost of the methods, as a function of multipolar order and within the context of condensed-phase molecular simulations of water.

Methods

The Charge Density

Our goal is to compute the electrostatic interaction of a neutral charge density ρ(r) with itself and its periodic images:[53]where the unit cell defining the periodicity is described by the lattice vectors a1, a2, a3 and reciprocal-space vectors a1*, a2*, a3*; and the periodic images are replicated by integer lattice translations n = n1a1 + n2a2 + n3a3. The energy (given by eq 1) is finite only when the charge density is neutral; however, a widely used approximation for the treatment of charged systems will be discussed in Section 2.6. Furthermore, we adopt the standard convention of removing the infinite self-energy of point charges (and point multipoles) whenever those energies may appear. Upon replicating all integer translations of the aperiodic density, ρ(r), the edges of a lattice with the same shape, orientation, and dimensions can be redrawn about any origin; and the periodic density filling each cell will contain exactly one instance of ρ(r) that appears to have been “wrapped” into the cell’s interior. The apparent wrapping of ρ(r) to the cell boundary is an illusion formed by the effluence of density emanating from the other translations, upon making ρ(r) periodic. The aperiodic density is not required to be confined within a primary unit cell to achieve this effect. Nor is it a requirement for the evaluation of eq 1, because the electrostatic potential of ∑ρ(r + n) is periodic.

The methods presented in this work are solutions to eq 1 for a density composed of atom-centered point multipole expansions. There are various ways of defining what is meant by a multipole expansion. The two most common ways being Cartesian multipoles[44] or those based on spherical harmonics.[52,54,55] Our definition of a multipole moment is the inner product of a density with a real-valued regular solid harmonic:[56]An auxiliary basis of charge density χ(r), satisfyingis thus well-suited to describe an atomic multipole expansion. The definition of C(r) is closely related to the scaled regular solid harmonics frequently encountered in fast multipole methods.[54,57−59] In brief, the real-valued scaled regular solid harmonics Rc/s(r) are the real Rc(r) = ReR(r) and imaginary Rs(r) = ImR(r) components of the complex-valued scaled regular solid harmonics:[58]whereis an associated Legendre polynomial. Because the complex-valued harmonics are symmetry related R*(r) = (−1)R(r), the real-valued harmonics are fully described by the set of non-negative m values. The notation is simplified, when appropriate, through the introduction of a Greek subscript whose sign merely acts to symbolize the cosine/sine designation[57]The real-valued regular solid harmonicsdiffer only by factorsthat are chosen to reproduce Racah’s normalization

In the interest of writing the expression for a point multipole, we begin by deducing the form of a Gaussian multipole function from eqs 3 and 9:whereis a Gaussian monopole, ∇ = {d/dX,d/dY,d/dZ} are gradients with respect to the Gaussian center, and C(∇) is a spherical tensor gradient operator[60] (STGO). A STGO is constructed by replacing the Cartesian coordinate arguments of the solid harmonic with their corresponding Cartesian derivative operators. The STGO obeys several mathematical properties.[61,62] The present work specifically makes use of the STGO chain-rule identity,which follows from the homogeneity of solid harmonics C(ar) = aC(r). The application of eq 12 to the last line of eq 10 establishes the relationship between χ(r – R) and a Gaussian monopole. A Dirac delta function is a Gaussian monopole in the limit of infinite exponent:and a point multipole is this limit applied to eq 10:Thus, a sum of atom-centered point multipole expansions is

Plane Waves

Methods for solving eq 1 avoid the infinite sum of explicit lattice translations by approximating a periodic density with a finite number of plane wave basis functions ⟨r|k⟩ ≡ e. The description of these methods is greatly simplified using Dirac notation, which we briefly summarize here. In this notation, the spatial representation of a function is equivalently written ⟨r|f⟩ = f(r); a complex conjugate is ⟨f|r⟩ = f*(r); the inner product of two functions ⟨f|g⟩ = ∫f*(r)g(r) d3r integrates all space or, if both functions are periodic, the volume of the unit cell; and we reserve |k⟩ to index the plane wave basis functions. This index is meant to be the function’s “angular wave number”, which is defined as k = 2π(k1a1* + k2a2* + k3a3*). From these definitions, one can show that the basis is orthogonal ⟨k|k′⟩ = δV, where V = a1 · a2 × a3 is the unit-cell volume; they obey a trivial addition theorem ⟨r + r′|k⟩ = ⟨r|k⟩⟨r′|k⟩; they are eigenfunctions of the gradient operator ∇⟨r|k⟩ = ik⟨r|k⟩; and the electrostatic potential of a plane wave basis functionbecomes particularly amenable to computation. Finally, it is easy to show that the plane waves are eigenfunctions of the STGO,via direct application of eq 12 with u = ikT·r.

Plane waves can be used as a basis to represent real, periodic functions:where c = a + ib. Complex numbers have been introduced to compactly express two independent linear least-square fit problems: the representation of the even (cosine) and odd (sine) character of the density. One can show that the two least-square fit solutions for a and b are equivalent to having chosento minimize the sum of squared errors ⟨Δ|Δ⟩, defined byIn other words, eq 19 is identical to having determined a and b by solving the two (uncoupled) equationsandrespectively. By integrating over the volume of a unit cell, the first line of eq 19 effectively integrates one instance of ρ(r) that has been wrapped into the cell’s interior. The second line of eq 19 exploits the periodicity of the plane waves by “unwrapping” the density and modifying the integration limits accordingly.

One can deduce that cos(−x) = cos(x) and sin(−x) = −sin(x) imply the symmetries a– = a, b– = −b, and c– = c*. As a consequence of these symmetries,and, therefore, the sum over allk in eq 18 naturally acts to cancel the imaginary components from each term. That is, although each term is complex, the other half of the sum adds its conjugate. The appearance of “Re” in our description of Ewald, PME, and FFP is to remind the reader that the imaginary numbers vanish and only the real component of each term thus needs to be computed. Equation 18 is often referred to as a “complex-to-real” reverse (or inverse) Fourier transform and is typically implemented within computer software to utilize only a subset of symmetry-unique c* values that are provided as input.

Weighted least-squares fits can be performed to generalize the plane wave expansion described above. These fits produce coefficientsthat minimize the sum of squared errors ⟨Δ|Ô|Δ⟩ weighted by a linear Hermitian operator Ô. In the general case, the components of c store the result of a fit that may have required a coupled (as opposed to independent) solution to a and b. For the specific case ⟨r|Ô|r′⟩ = |r – r′|–1, eq 21 is called an “electrostatic fit”. The coefficients of an electrostatic fit reduce to eq 19, because the plane waves are eigenfunctions of the Coulomb operator (eq 16). In other words, eq 19 are the plane wave expansion coefficients that best reproduce the electrostatic potential of the periodic density.

Multipolar Ewald

One could, in principle, directly project ρ(r) into the periodic basis and compute the interaction from the plane wave representation; in practice, however, this is intractable because an infinite number of plane waves would be required to model δ(r). Therefore, the Ewald method abandons the direct solution of eq 1, preferring instead to decompose the periodic density into smooth and discontinuous components:whereby the smooth model density composed of Gaussiansis well-approximated from a linear least-squares fit to reasonably few plane waves:Furthermore, the model density was chosen to reproduce the long-range electrostatic potential of the point multipoles to ensure that only short-ranged corrections are required. In our notation, ⟨k|ρ̃⟩ and ⟨k|χ⟩ = e–( are Fourier coefficients of the Gaussian density and Gaussian monopole, respectively. The summation over k in eq 24 is a reverse Fourier transform, and S is a “structure factor”.

The Ewald energy is the interaction of ρ(r) with the electrostatic potential of the plane wave projected periodic Gaussian density ϕ̃(r) (see eq 32) upon correcting the short-ranged differences between the point and Gaussian potentials Δϕ(r) (see eq 27):where p̃ = ∂⟨ρ|ϕ̃⟩/∂q and Δp = ∂⟨ρ|Δϕ⟩/∂q are the corresponding “multipolar potentials”. The real-space corrections to the potential, energy, and multipolar potential are given as follows:whereThese corrections are assumed to be sufficiently short-ranged, such that only nearby minimum image separations R and |r – R| need to be considered. An efficient algorithm for computing eq 30 and its gradients, particularly for low-order harmonics, is described in ref (56). The application of that algorithm requires the following definition of the “auxiliary vector”:where F(x) is the Boys function.[63]

The reciprocal-space Ewald potential, energy, multipolar potential, and contribution to the gradient are readily obtained upon noting ⟨δ|k⟩ = 1.The extension of Ewald to atom-centered point multipole expansions presented here differs from the standard point-charge Ewald method only by the inclusion of STGOs in the structure factor and real-space correction.

Multipolar Particle Mesh Ewald

A linear-scaling implementation of the real-space Ewald corrections requires only short-range cutoffs and a sufficiently large Gaussian exponent; however, by fixing the cutoff—and, hence, the exponent—the computational complexity of the reciprocal-space Ewald potential then scales with O(N2), becauseIn the unlikely event that the particles happened to be positioned such that they formed a uniformly spaced grid, however, then the expression for S would become a sum over grid points (see eq 48) and the contribution of all grid points to all S could be computed O(N log N) using a Fast Fourier Transform (FFT). The novelty of PME is its coercion of the particle positions to exploit this otherwise unlikely scenario to speed the calculation of the reciprocal-space potential and energy.

the number of plane waves must increase with the size of the system, to maintain consistent resolution of the periodic density, and

the calculation of each S requires a summation over the number of particles.

PME approximates the particle positions with a linear combination of predefined uniformly spaced FFT grid points R,that can be used to compute Fourier coefficients from numerical quadrature. The Gaussian quadrature of plane waves is equivalent to Simpson’s integration rule performed on a regular grid:where f is a discrete Fourier transform (DFT):t is the index of the grid (t = (t1, t2, t3)) and R is the location of a FFT grid point:N1, N2, and N3 are the number of FFT grid points in each lattice direction; N = N1N2N3.

θ is a Cardinal B-spline weight evaluated about a particle position:In our notation, the argument of M(u) is meant to be periodically wrapped into the range (−N/2,N/2].

M is a Cardinal B-spline function that holds recursion properties for efficient evaluation of their values and gradients.[38]The B-spline weights are real ⟨r|θ⟩ = ⟨θ|r⟩, have even symmetry ⟨−r|θ⟩ = ⟨r|θ⟩, sum to one (∑⟨R – r|θ⟩ = 1 ∀r), and are nonzero only for those grid points near the particle. The number of grid points contributing to the approximate position of a particle is determined from the chosen B-spline order without regard to the FFT grid spacing.

Following the work of Schoenberg,[64−66] previous descriptions[38] of PME have used Cardinal B-splines as a linear basis to construct “exponential Euler splines.” In that point of view, the structure factors are computed from approximate plane waves that are spline-evaluated at the actual particle positions. We below choose an alternate perspective by applying eq 36 to obtain the same outcome, which we then relate to ref (38).

The charge density is composed of atom-centered functions; therefore, it is convenient to re-express eq 36 aswhich is deduced from the invariance of the relative particle and grid positions upon their simultaneous reflection about the origin and subsequent translation by +r; that is, R → r – R and R → r – R. An atom-centered Gaussian is thus approximated by a linear combination of Gaussians centered about the FFT grid points:However, this has the unintended consequence of making the approximate Gaussian appear artificially diffuse as the B-spline order is increased. Fortunately, if the B-spline is of sufficient order to be accurately integrated by the discrete Fourier transform, then the Fourier coefficients of the approximate Gaussianare those of the true Gaussian,scaled by the B-spline DFT coefficients θ. Note that θ = θ*, because the B-splines are even functions. Repeating this procedure for , or any function |f⟩, yields analogous results. In other words, the Fourier coefficients of an approximately positioned function should be scaled by θ–1.

Let us continue on a brief aside that relates eq 36 to the exponential Euler splines used in ref (38) by considering its application to ⟨k|k′⟩ = δV in the manner above using either of the relationsThe result ∫⟨k|r⟩⟨r̃|k′⟩ d3r = θ⟨k|k′⟩ implies the need to effectively renormalize an “approximate plane wave” upon grid interpolation; that is, ⟨r|k⟩ ≈ ⟨r̃|k⟩θ–1, which are the exponential splines. The relationship between eq 36 and the exponential Euler splines is the interpretation of whether it is the approximate plane wave or the Fourier coefficient of the approximately positioned function that is renormalized:The subtle distinction between these interpretations has no practical consequence on the PME method other than, perhaps, how the mathematical terms are grouped together; eq 47 scales the structure factor by θ–1, whereas ref (38) absorbs it into the definition of S.

The periodic representation of the model Gaussian density (eq 23) composed of the approximate Gaussians (eq 42) is, upon scaling the Fourier coefficients,where the PME structure factor,is the forward DFT ofwhose evaluation avoids O(N2) operations by the locality of θ(r).

The PME reciprocal-space potential is the electrostatic potential of eq 47,and yields the following PME reciprocal-space energy, multipolar potential, and contribution to the gradient:The extension of PME to atom-centered point multipole expansions presented here differs from the standard point charge PME method only by the inclusion of STGOs acting upon the B-spline in eqs 49–51 and, like Ewald, the real-space correction (eqs 27 and 28).

Although the STGO has many remarkable properties[61] when acting upon a spherical function, a regular or irregular solid harmonic, or a product of the these functions, the B-spline weights cannot take advantage of these properties. Instead, one must treat the STGO as a linear combination of gradient operators,whose coefficients are exactly those that relate Cc/s(r) to Cartesian monomials:[67]whereand oc = 0 and os = 1. In the above notation, m = |μ|, C(∇) = Cc(∇), and C(∇) = Cs(∇).

Multipolar Fast Fourier–Poisson Method

Fast Fourier–Poisson (FFP), like Ewald, abandons the direct solution to eq 1 by introducing a Gaussian density whose potential can be determined from a plane wave expansion.[41] Unlike the Ewald method, FFP evaluates the plane wave expansion coefficients from numerical quadrature to take advantage of FFTs. The electrostatic potential of the plane-wave-resolved Gaussian density is computed in Fourier space, evaluated on a regular grid by means of a reverse FFT, and then used to numerically integrate the Coulomb self-energy of the Gaussian density. Therefore, the FFP real-space corrections account for the short-ranged differences between the Gaussian and point interaction energies:The plane wave representation of the periodic model density composed of the Gaussian multipoles iswhereThe evaluation of ⟨R|ρ̃⟩ avoids O(N2) operations by evaluating χ(R – R) only for minimum image distances within a suitable cutoff. The real-space energy corrections, ΔT, vanish when the Gaussians do not overlapwhere s = 1 − δδδδ). Therefore, eq 60 is well-approximated by computing the integrals for the short-ranged, minimum-image R separations only. Explicit expressions for the real-space energy correction, multipolar potentials, and gradient contribution are analogous to those shown in eqs 28–31 upon substituting ζ with ζ/2, which arises from the Gaussian product theorem.

The FFP reciprocal-space energy, multipolar potential, and contribution to the gradient are as follows:The gradient of the Gaussian basis,is efficiently computed using the following derivative properties:[58]When negative m values are encountered, the reader is implicitly instructed to apply the symmetry property Rc/s(r) = ±(−1)Rc/s(r), which follows directly from R*(r) = (−1)Rc/s(r), where the plus/minus sign (±) corresponds to the cosine/sine designation.

Correction for Charged Systems

Let ⟨r|ĵ|r′⟩ = ∑|r – r′ + n|–1; then, we can rewrite eq 1.If the net charge is nonzero (∫ρ(r) d3r = Q), then the system is nonphysical and this definition of the energy is infinite. Although it is possible to neutralize the density in arbitrarily different ways to produce different energies, a particularly convenient choice is made by introducing a uniform background density,[68]and ∑ρ̅(r + n) = Q/V, everywhere. With this choice of neutralization, the energy isHowever, 1/2 = 1/2⟨ρ|ϕ̃⟩, as written in the previous sections, by virtue of excluding the k = 0 term in the Ewald summation, because it is charge-neutral. Similarly, 1/2 vanishes, because the exclusion of the k = 0 term produces a potential composed of functions that are each orthogonal to a constant over the range of the unit cell. The last term in eq 70 vanishes, by symmetry, for all nonzero multipoles and is well-approximated by an integral over all space if ζ is sufficiently large to extinguish the correction potential within the bounds of the unit cell (in a minimum image context). As a result of these properties, the energy of a charged system isAn analogous result is written for FFP by replacing ζ with ζ/2.

For completeness, we note that it has been pointed out[69,70] that eq 1 is absolutely convergent only if the total dipole moment of the system μ is zero; however, it is otherwise only conditionally convergent. In other words, the asymptotic value of eq 1 is dependent on the shape of the supercells used to replicate the system. A “dipole surface” term 2π|μ|2/[(2ε + 1)V] can be applied to mimic an infinite spherical sample of a cubic cell; however, this energy is usually ignored in Ewald implementations,[38] because its application to charged systems has a dependence on an arbitrary origin and it produces discontinuous energy changes when charged residues are wrapped into the primary unit cell.[71] By ignoring this term, we are said to employ “tinfoil boundary conditions”, because it is equivalent to embedding the crystal inside a perfect conductor (ε → ∞).

Error Analysis

Figure 1 compares the Ewald energy (eq 26) to a brute force evaluation of eq 1, as a function of replicated unit cell length. The fundamental unit cell is a 10 Å cube containing three-point multipoles:where R = (1 Å)ẑ. The unit cell is replicated by extending the sides of the cell to form increasingly larger cubes. The various lines in Figure 1 indicate the nonzero angular moment in eq 72. The “relative energy error” is defined aswhere E is the Ewald energy evaluated with |kmax| = 256 and E is eq 1 after replication to the specified size. The relative error is bounded by ≥10–15, because of our use of double-precision arithmetic. The loss of one or two additional digits of precision should be expected, because of the accumulation of round-off errors.

Figure 1

Comparison between the Ewald energy to the energy computed from explicit replication of a unit cell, as a function of replicated cell length for systems composed of the indicated multipole orders.

Figure 2–4 make use of an ad hoc three-site water model to analyze the relative force errors and the real-space and reciprocal-space timings of Ewald, PME, and FFP, as a function of atomic multipole order. These three figures are the analysis of a box of 1024 waters whose bulk density has been equilibrated with TIP4P-Ew at 298 K. The 3-site multipolar water models used to perform the analysis are systematically constructed from a 5-site charge-only water model. The 3- and 5-site water models use the DFTB3[72] gas-phase structure of water (these coordinates are superimposed onto the TIP4P-Ew waters): ROH = 0.957143 Å, RHH = 1.572691 Å, ∠HOH = 110.4812° (the 5-site model includes two additional charges located at the O–H bond midpoints). The atom charges (a.u.) (qO5-site = −0.12404869708, qH5-site = 0.66106645662) and the bond charges (qB5-site = −0.59904210808) were computed from the isolated water DFTB3 density matrix; that is, qO5-site is the one-center atomic orbital (AO) product contribution to the Mulliken charge of O, and qB5-site is the Mulliken bond charge between O and H. The multipole moments of the 3-site model are constructed from a Mulliken-like partitioning and subsequent solid harmonic translation of the 5-site model bond charges:The intrawater interactions are excluded when analyzing the relative force errors with these models.

Figure 2

Decomposition of the intermolecular electrostatic forces using a box of 1024 waters. Lmax indexes the ad hoc 3-site water model being considered, as distinguished by the multipole expansion order on the oxygen.

Figure 4

Reciprocal-space evaluation timings (left column), relative force errors (RFEs) (middle column), and the Gaussian exponents (right column) chosen to minimize the RFEs of Ewald (top row), PME (middle row), and FFP (bottom row), using the box of 1024 waters. Lmax indexes the ad hoc 3-site water model being used.

In order to interpret the results in Figure 4, it is useful to understand how the multipole order influences the magnitude of the forces. For this purpose, Figure 2 decomposes the intermolecular electrostatic forces within the 1024 water box by the percent contributed from the charge–charge interactions, the non-“charge−charge” interactions, and those interactions that involve a particular angular momentum L. Let F be the 3 N vector of atomic forces for a system of waters composed of the ad hoc 3-site water model using an order L multipole expansion. The “percent contribution” of the charges, the multipoles, and a particular multipole L to the forces are respectively chosen to be

The relative force errors (RFEs) in Figure 4 are calculated using the expressionwhere Fref, is the vector of Ewald forces evaluated with |kmax| = 128 and whose real-space corrections include all N2 minimum image interactions. The Gaussian exponent has been nonlinearly optimized to minimize each RFE, and their values are shown in Figure 4 using the convention β = ζ1/2.

The timings reported in Figures 3, 4, and 6 were performed on an Intel Xeon E5520 2.27 GHz workstation consisting of a total of 8 cores, and the software was rudimentarily parallelized with OpenMP to make use of all cores.

Figure 3

Wall clock time required to evaluate the real-space corrections (eq 29) for a box of 1024 waters using different real-space cutoffs. Lmax indexes the ad hoc 3-site water model being used. The inset values are timing ratios T(Lmax)/T(0), that is, the slow-down relative to the charge-only model.

Figure 6

SCF and SCF + gradient timings using the mDC quantum force field for a series of water boxes. “L = 0” indicates that only atomic charges were used to compute electrostatics, whereas “L = 2” indicates that up to quadrupoles were used on the heavy atoms.

Implementation within the mDC Quantum Force Field

The analyses in Figures 5, 6, and 7, and the data shown in Table 1, were produced from condensed-phase calculations using the mDC linear scaling quantum force field.[31] In brief, the quantum mechanical treatment of the entire system is replaced by a series of quantum calculations for each molecule. Although the molecular orbitals of each fragment are not directly coupled, the subsystems remain coupled through the interactions of their electron densities and empirical potentials. Specifically, the mDC method performs DFTB3/3OB semiempirical calculations[72] of each molecule while subjecting them to an effective chemical potential arising from the intermolecular interactions. This idea is conceptually similar to “density functional embedding theory”;[73] however, mDC replaces the rigorous evaluation of the embedding potential with computational tractable empirical approximations that are tuned for high accuracy. The standard DFTB3 semiempirical model evaluates the electrostatic interactions from Mulliken charges that are updated until self-consistency is reached. In the mDC formalism described below, the standard DFTB3 treatment continues to be used for those atoms within a common fragment; however, we concoct an auxiliary set of atomic multipoles from the density matrix to improve the intermolecular interactions. The multipolar potentials arising from those interactions enter the fragment Fock matrices as an effective external chemical potential, and the multipoles are updated at each step of the SCF procedure. Let E be the DFTB3 energy for molecule A consisting of atoms located at R; then, the mDC energy iswhere N is the number of electrons in fragment A, q is an atomic multipole moment determined from the DFTB3 density matrix, and p = ∂Einter/∂q is the multipolar potential describing the effective external potential experienced by atom a. With this convention, the σ-spin Fock matrix used to construct the orbitals of molecule A is given asAs explained in ref (31), the atomic charges are computed from a biased Mulliken partition of the density matrix, and the higher-order multipole moments (up to quadrupoles for heavy elements) are computed from the single-center AO–product components of the density matrix. The intermolecular interaction energy Einter consists of Lennard-Jones potentials and the electrostatic interaction of the atomic multipoles. The electrostatic energy and multipolar potentials are evaluated using either the Ewald, PME, or FFP methods described in the previous sections. Because the multipole moments are dependent on the density matrix and the Fock matrix is dependent on the interactions with the other molecules, the mDC energy must be optimized until a mutual convergence is met for all molecules.

Figure 5

Conservation of the mDC total energy for 1 million steps in a NVE simulation of N,N-dimethylglycine, using the multipolar PME method with 6th order B-splines.

Figure 7

RFEs observed with the mDC linear-scaling quantum force field using various methods of computing the electrostatic interactions. The top portion of the figure shows a box of 1024 waters, whereas the bottom of the figure shows a supercell of N,N-dimethylglycine containing 4608 atoms. Errors are reported relative to FFP evaluated with 3 points/Å and a minimum image (“min img.”) short-range correction. The label “9 Å” indicates the range used for the short-range real-space corrections and, for FFP, the evaluation of the model Gaussian density at the FFT grid points. Dotted lines indicate the ability of the electrostatics protocol to reproduce the electrostatic component of the force, using the atomic multipoles obtained from the reference FFP calculation. The dashed lines are mDC RFEs, using cutoff-based electrostatics without any Ewald treatment.

Table 1

Effect of Electrostatic Cutoff on the Condensed-Phase Properties of Water at 298 K

	density, ρ  (g/cm3)	heat of vaporization, ΔHvap (kcal/mol)
expta	0.9970	10.51
PMEb	0.9969	10.62
9 Å cutoff	1.0110	10.87

Data taken from refs (86) and (87).

Sixth-order B-spline interpolation, grid spacing of 1 point/Å, and short-range cutoff of 9 Å.

Figure 5 illustrates the quality of the mDC forces by demonstrating conservation of the total energy in a NVE simulation of an N,N-dimethylglycine (DMG) crystal at 225 K. The crystal consists of 288 DMG residues (4608 atoms) constructed from supercell replication of the experimentally determined X-ray structure[74] and was simulated at the experimental density with a locally modified development version of PMEMD[75] that incorporates the mDC method. The simulation was performed using a 0.5 fs time step for 1 million steps, and the electrostatics were computed with sixth-order B-spline PME, a 9 Å real-space cutoff, and a 1 pt/Å FFT grid density.

Figure 6 demonstrates linear scaling of the mDC method for a series of water boxes. In addition, timings were recorded for an ad hoc mDC model that limits the intermolecular interactions to monopoles for the sole purpose of determining the relative cost of incorporating atomic quadrupoles. The PME calculations were performed with sixth-order B-splines, a real-space cutoff of 9 Å, and an FFT grid density of 1 point/Å. The FFP calculations were performed with an FFT grid density of 1 point/Å and a cutoff of 9 Å in both the real- and reciprocal-space calculations, that is, the evaluation of the model Gaussian density at the FFT grid points.

Unlike the ad hoc water models used in Figure 4, mDC self-consistently models each molecule’s polarization in response to its environment (the other molecules). Because the polarization is induced from the environment’s electrostatic potential, the electrostatic protocol indirectly introduces atomic force errors by altering the electronic polarization. Within the mDC method, the change in polarization is produced from the difference in the multipolar potential contributions to the Fock matrix. Figure 7 decomposes the mDC force errors to quantify the extent to which the protocols indirectly effect the forces through their difference in polarizations. First, a mDC calculation is performed using a reference electrostatic protocol to produce a SCF converged energy (EmDC,ref) and charge density (qref). The dotted lines in Figure 7 are the RFEs computed from the vectors of intermolecular electrostatic forces evaluated with qrefwhere Eelec,model and Eelec,ref are the model and reference intermolecular electrostatic energies, respectively. A second mDC calculation is then performed using a model electrostatic protocol throughout the SCF procedure to produce a converged energy EmDC,model and charge density qmodel. The solid lines in Figure 7 are the RFEs computed fromwhich are evaluated with different electrostatic protocols and charge densities. The dashed lines in Figure 7 are computed from eqs 83 and 84, where EmDC,model is evaluated using a molecule-based cutoff electrostatic protocol; that is, no treatment for periodicity was applied beyond employing the minimum image convention within the cutoff radius. The electrostatic interactions between the waters were smoothly switched off from 8 Å to 9 Å, based on the O–O separation. All reference calculations were performed with FFP using a 3 point/Å grid density and a full minimum image treatment of the short-range corrections and reciprocal-space evaluation.

The PME and electrostatic cutoff methods are compared in Table 1, which examines how these methods effect the density and heat of vaporization of liquid water at 298 K. These properties were computed from 8-ns simulations of 512 waters in the NPT ensemble using the Monte Carlo barostat,[75] and the Andersen thermostat,[76] as implemented in PMEMD. A time step of 1 fs was used in conjunction with SHAKE[77,78] to constrain the internal structure of the waters to the isolated DFTB3 geometry. The parameters of the mDC water model used to generate these properties are described in Part 2 of this series.[51]

Results and Discussion

Comparison to Brute Force Replication

Figure 1 establishes the correctness of the Ewald formulas for various multipole orders by comparing its energy to brute force evaluation of eq 1. Analogous plots for FFP and PME with suitably dense FFT grids and B-spline order are indistinguishable from the Ewald results shown in Figure 1. The errors in the energy flatten to a constant in the range of 10–13–10–14, because the large number of calculations involved in the brute force evaluation suffer from a loss of numerical precision. As the multipole order is increased, the energy converges more quickly, because the Coulomb interactions decay by r–. When L ≥ 4, the range of the interactions become so small that the energy is well-approximated by the minimum-image interactions.

Errors and Timings as a Function of Multipole Order

Before comparing the Ewald, PME, and FFP RFEs, it is useful to briefly discuss the relative magnitude of the forces associated with increasing the multipole order. As shown in Figure 2, the charge–charge interactions represent 40%–50% of the force when higher-order multipoles are used in the ad hoc water model. Furthermore, higher-order multipoles progressively contribute less to the force, since their interactions become shorter-ranged. All interactions involving the octupoles within the Lmax = 3 water model, for example, contribute <5% to the force, whereas all interactions involving L = 5 contribute <0.5%.

The wall-clock time required to evaluate the real- and reciprocal-space contributions to the multipolar potentials, as a function of multipole expansion order, are displayed in Figures 3 and 4, respectively. The real-space corrections evaluated from eqs 29 and 30 and the algorithm in ref (56) are particularly efficient for Lmax ≤ 3. For higher-order expansions, techniques based on rotations into an internal-coordinate system can be used to further improve performance;[52,79] however, considering the short-range and relative insignificance of these interactions, it can be verily questioned if practical models for condensed-phase simulations will consider their additional cost worthwhile. Including quadrupoles on the oxygen causes the real-space evaluation timings to slow down by a factor of 1.8, relative to a charge-only model. Comparison of the reciprocal-space timings in Figure 4 suggests that, for this system of 1024 waters, FFP and Ewald are ∼2 and ∼100 times slower than PME, respectively, for comparable error levels. Furthermore, the evaluation of the real-space corrections at Lmax = 2 with a cutoff of 9 Å are approximately twice as slow as the sixth-order B-spline PME reciprocal-space evaluation times.

The optimized Gaussian exponents shown in Figure 4 are far more sensitive to the choice of real-space cutoff than they are to multipole order. As the real-space cutoff is increased, the Gaussian exponent decreases to become better resolved in the plane wave basis, and this consequently decreases the errors of the reciprocal-space forces. Similarly, increasing the PME B-spline order affords the opportunity to increase the Gaussian exponent to reduce the errors in the real-space correction without sacrificing accuracy in the reciprocal-space potential.

Comparisons within a Linear-Scaling Quantum Force Field Framework

Much of this manuscript uses the word “error” to describe a difference between an accurate electrostatics protocol and a more approximate treatment. However, this phrasing should not be misconstrued to incorrectly suggest that the PME gradients are inconsistent with the PME energy, for example. To emphasize this, Figure 5 demonstrates conservation of total energy in a NVE simulation of DMG for 1 million time steps evaluated with the mDC quantum force field, whose electrostatics are computed with PME. If there were an inconsistency between the energy and forces, then we would observe a significant drift in the energy as time is propagated; however, we do not observe a drift.

The mDC method is linear scaling when using either PME or FFP, as illustrated in Figure 6. When evaluated with PME, a box of 4096 waters can be SCF converged within <0.7 s. Furthermore, our timings indicate that the use of quadrupoles slows the mDC calculation by a factor of 1.5, relative to a charge-only model. This is slightly less than the ratio observed in Figure 3, because the mDC and charge-only mDC evaluations share common operations that are independent of the multipole order used in their intermolecular electrostatics (for example, their Fock matrix diagonalizations). The number of SCF cycles required to reach convergence is not particularly sensitive to the size of the system, which can be inferred from the linear-scaling shown in Figure 6. The SCF convergence rate is largely determined by the quality of the initial guess orbitals. In practice, the difference between MD steps is so small that convergence can typically be reached within 4–6 cycles.

Ewald, PME, and FFP are approximate, but their accuracy can be systematically improved by adjusting the number of plane waves, the Gaussian exponent, and the real-space correction cutoff. In the analysis of Figure 4, we examine how these parameters alter the electrostatic energy protocol by quantifying their effect on atomic forces (for a fixed charge density), −∂E/∂R|, while ignoring their effect on electrostatic potentials ∂E/∂q|. However, the parameters do effect the electrostatic potentials. Furthermore, the electrostatic potentials enter the Fock matrix via eq 80. Therefore, the parameters indirectly effect the charge density through the propagation of electrostatic potential “errors” within the SCF, which cause its convergence to a different density matrix. In other words, not only will ab initio methods (or any polarizable model) suffer from the “primary error” associated with the partial derivative −∂E/∂R|, they also incur a “secondary error” when the model and reference forces are evaluated about different charge densities. The dotted lines in Figure 7 are the primary electrostatic errors; that is, they are the RFEs associated with the electrostatic forces evaluated from the atomic multipoles generated from the reference calculation. The solid lines in Figure 7 contain both the primary and secondary errors produced when the reference and approximate electrostatic protocols are used throughout the SCF procedures to yield two slightly different densities.

The reader may first notice that the primary FFP errors decrease and converge, with respect to FFT grid density quickly; however, the secondary errors trail off at 10–11 to 10–10. This limit is ultimately caused by our SCF numerical convergence criteria. The errors in Figure 7 were generated with a 10–12 tolerance on the maximum value appearing in the “error matrix”, that is, the commutator between the Fock matrix and the density matrix.

As expected, one observes a decrease in error as the PME B-spline order is increased, the FFT grid density is increased, and as the real-space cutoff is increased. The mDC RFEs observed in water and DMG parallel the primary errors; however, the DMG secondary errors are significantly larger than those observed in water. The reason for this is because a DMG molecule contains 16 atoms, whereas water contains 3. By having more AOs, the DMG molecules are more likely to converge to a slightly different density matrix.

Before PME was popularized in MD simulations,[37,38] electrostatic cutoffs were often used, and many of the original water models were thus parametrized using a cutoff protocol. It was later found that application of those water models with PME caused significant changes in the thermodynamic properties of liquid water,[80−82] and it has been suggested that new water models not be parametrized using cutoffs,[81] because their application within simulations of biomolecules often necessitates the use of an Ewald treatment.[83,84] Nonetheless, some recently developed water models have continued to be parametrized using cutoffs.[85] Therefore, we are obliged to investigate by how much this may effect condensed-phase properties. First, in Figure 7, we note that a cutoff of 9 Å produces an RFE value of 0.1, in comparison to the RFE value observed with a sixth-order B-spline PME with a grid spacing of 1 point/Å FFT (2 × 10–5). Furthermore, even if all N2 minimum image interactions were computed without PME, the RFE is still 0.04. Finally, we performed MD simulations using a parametrized mDC water model to obtain an equilibrated density and heat of vaporization, as shown in Table 1, and we reperformed the simulations using switched-cutoff electrostatics to make comparison. In agreement with previous works,[80,81] we observe that cut-off electrostatics cause the density of water to increase, and because the interactions in this range are attractive overall, this also causes the heat of vaporization to increase.

Conclusion

This work presented extensions of the Ewald, Particle Mesh Ewald (PME), and Fast Fourier–Poisson (FFP) methodologies to systems composed of point multipole expansions to arbitrary order by making use of the spherical tensor gradient operator. The timings and errors inherent to these methods were compared using ad hoc water models and with a parametrized water based on the modified divide-and-conquer (mDC) linear-scaling quantum force field. These comparisons lead us to conclude that (i) the FFP method is approximately twice as slow as the PME method at comparable error levels; (ii) the inclusion of quadrupoles in the linear-scaling force field slow the calculations by 1.5, relative to a charge-only model; and (iii) with the exception of the Ewald method, the real-space corrections are more expensive than the reciprocal-space calculations for typical cutoff values. Furthermore, our results suggest that the evaluation of multipolar electrostatics involving orders greater than 3 could likely be computed to an acceptable error using an electrostatic cutoff, because of their overall short-range and relative insignificance, in comparison to lower-order interactions.

The relative force errors exhibited within mDC were decomposed into primary and secondary errors, where the primary errors directly result from the approximations within the PME or FFP algorithms for a given density, and the secondary errors are the propagation of the model’s electrostatic potential within the self-consistent field (SCF) procedure, resulting in a different converged density matrix. It is found that the force errors closely follow the primary errors, and the magnitude of the secondary errors is related to the number of atomic orbitals (AOs). Nevertheless, the presence of these “errors” does not imply that the mDC forces are inconsistent with its energy, which was demonstrated with an NVE simulation that was devoid of an energy drift. Instead, these differences merely reflect how similar an electrostatic protocol is to another reference protocol.

Finally, the importance of using an Ewald treatment in simulations, as opposed to using electrostatic cutoffs, was emphasized by comparing the density and heat of vaporization of water. The electrostatic cut-off method was found to artificially increase the density and heat of vaporization of water.