3D RISM theory with fast reciprocal-space electrostatics Jochen Heil and Stefan M. Kast Citation: The Journal of Chemical Physics 142, 114107 (2015); doi: 10.1063/1.4914321 View online: http://dx.doi.org/10.1063/1.4914321 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/142/11?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Communication: Accurate hydration free energies at a wide range of temperatures from 3D-RISM J. Chem. Phys. 142, 091105 (2015); 10.1063/1.4914315 Molecular acidity: A quantitative conceptual density functional theory description J. Chem. Phys. 131, 164107 (2009); 10.1063/1.3251124 The role of atomic quadrupoles in intermolecular electrostatic interactions of polar and nonpolar molecules J. Chem. Phys. 119, 2192 (2003); 10.1063/1.1585016 Erratum: “Theory of interaction between helical molecules” [J. Chem. Phys. 107, 3656 (1997)] J. Chem. Phys. 108, 7035 (1998); 10.1063/1.476117 Theory of interaction between helical molecules J. Chem. Phys. 107, 3656 (1997); 10.1063/1.475320

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.30.242.61 On: Sat, 30 May 2015 16:09:51

THE JOURNAL OF CHEMICAL PHYSICS 142, 114107 (2015)

3D RISM theory with fast reciprocal-space electrostatics Jochen Heil and Stefan M. Kasta) Physikalische Chemie III, Technische Universität Dortmund, Otto-Hahn-Str. 6, 44227 Dortmund, Germany

(Received 9 December 2014; accepted 25 February 2015; published online 17 March 2015) The calculation of electrostatic solute-solvent interactions in 3D RISM (“three-dimensional reference interaction site model”) integral equation theory is recast in a form that allows for a computational treatment analogous to the “particle-mesh Ewald” formalism as used for molecular simulations. In addition, relations that connect 3D RISM correlation functions and interaction potentials with thermodynamic quantities such as the chemical potential and average solute-solvent interaction energy are reformulated in a way that calculations of expensive real-space electrostatic terms on the 3D grid are completely avoided. These methodical enhancements allow for both, a significant speedup particularly for large solute systems and a smoother convergence of predicted thermodynamic quantities with respect to box size, as illustrated for several benchmark systems. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4914321] ργ hγ (r) =

I. INTRODUCTION



cγ′ ∗ χγγ′(r).

(1)

γ′

Solvation effects are recognized to be essential for understanding and controlling molecular-scale properties in condensed-phase systems. Here, the excess chemical potential, µex, is a key quantity for understanding solvation phenomena since it represents the component of the solvation free energy that is associated with the perturbation of the liquid’s structure due to the presence of the solute. µex can be calculated more or less exactly (limited by statistical noise) for a given Hamiltonian from molecular simulations using, e.g., thermodynamic integration (TI1) or free energy perturbation (FEP2) methods. Since these direct approaches are frequently rather involved and time consuming, various levels of approximations have been developed as alternatives, ranging from dielectric continuum methods based on the Poisson-Boltzmann (PB) equation, that model the solvent as polarizable but structureless, to statistical-mechanical liquid state theories3 that retain the solvent’s granularity in order to account, for instance, for Hbonding in an aqueous environment. State-of-the-art methods for quantitatively describing granular solute-solvent systems in terms of thermodynamic properties are classical density functional theory (cDFT4,5) and integral equation (IE) theory in various flavors. IE theories of liquids and solutions aim at analytical representations of liquid structure in terms of molecular or site correlation functions computed directly from interaction potentials and can be derived from a general cDFT perspective. The particular appeal of IE methods is related to the additional benefit that thermodynamic quantities can be derived analytically from these correlation functions. For studying solvation phenomena, the three-dimensional reference interaction site model (3D RISM) represents a promising IE approach.6,7 The 3D RISM equations connect the total correlation functions hγ (r) between solute molecule and solvent site γ and the direct correlation function cγ (r) on a 3D grid specified by r by the convolution product a)Electronic mail: [email protected]

0021-9606/2015/142(11)/114107/13/$30.00

The molecule-site pair distribution function is then given by gγ (r) = hγ (r) + 1. ργ is the solvent site bulk density, χγγ′ represents the solvent site susceptibility that is typically precomputed by a dielectrically corrected 1D RISM calculation (DRISM)8,9 for the pure solvent. Equation (1) is solved together with the “closure” relation hγ (r) = exp[t γR(r) + Bγ (r)] − 1

(2)

with the “fully renormalized” (for details of this terminology, see below) indirect correlation function t γR(r) = t γ (r) − βuγ (r) = hγ (r) − cγ (r) − βuγ (r),

(3)

where uγ (r) is typically (but not necessarily restricted to) the Lennard-Jones and Coulomb potential between the solute molecule and a solvent site, and β is the inverse temperature. B is the so-called bridge function which cannot be computed from first principles with high precision and is therefore often neglected, leading to the “hypernetted chain” (HNC) closure10 hγ (r) = exp[t γR(r)] − 1.

(4)

In practice, the HNC closure often causes numerical instability, for instance, due to a violation of the Lambert W condition connecting bridge and direct correlation function.11 Numerically more stable are the partial series expansions to n-th order (PSE-n),12 n   [t R(r)]i /i! − 1 ⇔ t γR(r) > 0  i=0 γ (5) hγ (r) =    exp[t γR(r)] − 1 ⇔ t γR(r) ≤ 0  that specialize to the Kovalenko-Hirata (KH) closure13 for n = 1 and to the HNC approximation for infinite order. In water or aqueous solutions, PSE-3 and PSE-4 typically yield an accuracy that is comparable to HNC,12,14 making these closures suitable for stable routine calculations. From the solution of the nonlinear 3D RISM/closure system of equations, the excess chemical potential can be

142, 114107-1

© 2015 AIP Publishing LLC

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.30.242.61 On: Sat, 30 May 2015 16:09:51

114107-2

J. Heil and S. M. Kast

J. Chem. Phys. 142, 114107 (2015)

obtained analytically for HNC and PSE-n, obeying the general form,12 µex = µex,(1) + µex,(2)     1 = − β −1 ργ dr hγ (r) cγ (r) + µex,(2), 2 γ with µex,(2) HNC

= β

−1



 ργ

γ

and µex,(2) PSE−n



1 dr hγ2 (r) − cγ (r) 2







ργ

(7)



1 2 h (r) − cγ (r) 2 γ γ  − Θ(hγ (r))[t γR(r)]n+1/(n + 1)! .

=β

−1

(6)

dr

(8)

For non-zero bridge functions, an integral over a functional derivative of B multiplied by the pair distribution function appears which has to be added to Eq. (6)15,16 such that the derivations below, which concern the efficient calculation of the interaction potential and µex,(1), are valid even for a hypothetical exact theory (note that the integration of the formally non-zero B within the PSE-n approximation is contained in the last term of Eq. (8) as an analytical result). With the availability of an analytical free energy expression, the scope of 3D RISM applications is broad, ranging from small molecules in solution to biomolecular complexes,17 ion distributions,18 and the computational characterization of potassium channel proteins19,20 with particular emphasis on ion selectivity.21 A particular strength arises from the combination with quantumchemical calculations. For instance, 3D RISM can replace dielectric continuum methods for the solvent description22,23 where solvent granularity is crucial,24,25 as applied to tautomer ratios,26 hydration free energies,27 or solvatochromic shifts.28 A typical 3D RISM calculation is composed of three consecutive steps: (a) potential energy functions are calculated on a 3D grid, (b) 3D RISM/closure equations are solved for this potential yielding the total and direct correlation functions, (c) the excess chemical potential is calculated from the correlation functions obtained after convergence of step (b). For small solute molecules, the most time-consuming part is typically step (b). Therefore, the development of improved 3D RISM solvers has in the past mostly focused on this part and several efficient numerical schemes have been applied.29–33 Probably the most widely used convergence acceleration method is the “modified direct inversion of iterative subspace” (MDIIS) algorithm,29,30,34 which improves upon the simple fixed-point (“Picard”) iteration by modifying Pulay’s selfconsistent field calculation accelerator34 often used in quantum chemistry. For large solute species on the nanoscopic scale and the consequential requirement of not only a large number of grid points but also high-resolution grids (as indicated by recent results for finite-difference PB solvers35), these improvements shift a large part of the remaining computational burden towards the efficient computation of the potential functions in step (a). The accuracy-maintaining acceleration of step (a) (and the subsequently resulting implications for the computation of step (c), as we will see shortly) has not

yet been investigated in a rigorous fashion. As to step (a), recently Luchko et al. introduced an interpolation scheme where, outside a given radius, the potential is calculated only for every eighth grid point and interpolated onto the remaining ones.36 While it was shown that stable implicit solvent molecular dynamics (MD) simulations could be performed with this method, the effect on the accuracy of single structure energetics is not known. In practice, long range terms on finite grids pose a computational and numerical challenge. In the 3D RISM context, such terms can originate from electrostatic contributions to the molecule-solvent site interactions that require proper treatment within both, the closure approximation (2) and the chemical potential expression (6). The standard procedure to solve the 3D RISM equations involves the computation of electrostatic interactions in reciprocal (k) space, implying 3D translational periodicity of the grid data (utilizing the minimum image convention for short range contributions). In this way, the convolution product (1) becomes algebraic, and there is in principle no need to compute a real-space representation of electrostatics for step (b) since long range interactions can be “absorbed” in the indirect or direct correlation functions that enter the closure (2), which can be written alternatively as hγ (r) = exp[− βuγ (r) + t γ (r) + Bγ (r)] − 1 = exp[− βuγS(r) + t γS(r) + Bγ (r)] − 1 = exp[− βuγS(r) − cγS(r) + hγ (r) + Bγ (r)] − 1,

(9)

where we define the process of “renormalization” to split the total potential into short range (superscript “S”) and long range part (superscript “L”), uγ (r) = uγS(r) + uγL(r),

(10)

giving rise to (partially) renormalized direct and indirect correlation functions cγS(r) = cγ (r) + βuγL(r),

(11)

t γS(r) = t γ (r) − βuγL(r),

(12)

from which the equivalent forms of Eq. (9) become immediately clear. Note that correlation functions must be expressed in reciprocal space by a Fourier transform in order to compute the convolution product in Eq. (1) efficiently. This also applies to long range potentials and automatically implies translational periodicity of all functions in 3D space (known as “KovalenkoHirata supercell method” in the literature) and therefore the minimum image convention for computing real-space short range potentials. The 3D box (where real- and reciprocalspace resolution are mutually dependent) should be chosen large enough for a fluid, i.e., intrinsically aperiodic system such that periodicity artifacts are negligible. In contrast to an MD framework, 3D RISM does not require access to the total potential energy of a solute-solvent system but only to solute-solvent potentials. Our derivations below agree with this general model perspective. For net-charged solutes, the finite 3D box size has consequences for describing the long range asymptotics correctly. To this end, we37 developed a

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.30.242.61 On: Sat, 30 May 2015 16:09:51

114107-3

J. Heil and S. M. Kast

J. Chem. Phys. 142, 114107 (2015)

concept to handle the monopolar net charge term on a larger radial grid, a methodology later extended by Perkyns et al.38 to an exact treatment of site-site electrostatics. Renormalization offers the principal possibility to accelerate step (a) by employing a strategy analogous to the “particle-mesh Ewald” (PME) approximation39 to calculate the Ewald reciprocal-space long range potential with a computational complexity of O(Nα + nxnynz log(n x n y nz )) instead of O(Nα n x n y nz ) where Nα is the number of solute sites, and nx, ny, and nz are the number of grid points in each direction. In what follows, we call this the “particle-mesh supercell” (PMS) method which implies a charge spreading procedure to transfer the solute’s site charge distribution to a regular 3D grid such that an approximate Ewald reciprocal-space potential results from a 3D fast Fourier transform (FFT). Formulation and implementation of this approach represent one major goal of the present article. However, the evaluation of the chemical potential in step (c) poses an additional barrier. Although µex,(2) in Eqs. (6)–(8) is unproblematic since cγ can be safely replaced by cγS due to the electroneutrality condition, µex,(1) still involves the full real-space potential which would reduce the significance of the speed advantage gained by the PMS approximation. To this end, we express µex,(1) (and the average solute-solvent energy as well) completely in reciprocal space, termed the “no real-space supercell” (NRS) method in the following, and discuss various peculiarities. This represents the second major goal of this paper, facilitating complete circumvention of real-space electrostatics computations (except for the unproblematic short-range component). Contrary to the approach by Luchko et al.,36 NRS by itself does not involve significant additional approximations. We first develop the theory in full detail starting with a derivation of the NRS approach followed by the PMS approximation, after which we will show for a number of numerical examples that the implementation leads to substantial time savings, is accurate and reduces 3D grid artifacts. Furthermore, the convergence of the excess chemical potential with respect to the box size is accelerated and memory consumption is reduced. Therefore, the combination of NRS and PMS significantly speeds up 3D RISM calculations for large solutes and high-resolution grids. II. THEORY A. Renormalized 3D RISM equations and excess energies

For numerical calculation, the convolution product in Eq. (1) is discretized using the 3D discrete Fourier transform that is evaluated by FFT. In practice, we solve the renormalized form as the fixed-point problem S S S cγ,i+1 (r) = exp[uγS(r) + t γ,i (r) + Bγ (r)] − 1 − t γ,i (r), (13)  S S ′ S t γ,i (r) = ℑ−1[{ℑ[cγ,i (r)] − βuˆγL′(k)} χˆ γγ ′(|k|)] − cγ,i (r), ′ γ

(14) where χˆ ′ = ρ−1χˆ (bold letters indicating matrices containing all pairs of site-site functions, the density matrix being diagonal) and i denotes the iteration counter. Throughout, we define the

Fourier transform of a real-space function by  fˆ(k) = ℑ[ f (r)] = dr f (r) exp(−ik · r)

(15)

and the inverse Fourier transform as  1 dk fˆ(k) exp(ik · r). (16) f (r) = ℑ−1[ fˆ(k)] = (2π)3 The long and short range potential functions uγL and uγS are constructed using a suitably chosen splitting function for the Coulomb part of the total potential. Typically, one employs Ewald summation-type splitting for the distance (r) dependent site-site potential,40 i.e., 1 erf(κr) erfc(κr) = + , (17) r r r where erf and erfc denote the error function and its complement, and κ is the Ewald charge smearing parameter. The total reciprocal-space potential arising from all solute site charges interacting with a specific solvent site then follows from ˆ sˆ(k) = qγ ϕ(k) uˆγL (k) = qγ a(k) ˆ b(k) ˆ sˆ(k),

(18)

where 4π 4π = 2, 2 |k| k ( ) 2 ˆb(k) = exp − k 4κ 2

a(k) ˆ =

(19) (20)

are the Fourier-transformed continuum Coulomb Green’s function and a “charge smearing function” that accounts for the Gaussian shape of the Ewald charge distributions as opposed to point charges, respectively, such that ϕ(k) ˆ is the smeared continuum Coulomb Green’s function in reciprocal space. The key quantity that describes the solute’s charge distribution is  Nα sˆ(k) = qα exp(−ik · Rα ), (21) α=1

the electrostatic structure factor arising from solute sites α at positions Rα = (Rα, x , Rα, y , Rα, z )T. qα and qγ are the partial charges of solute site α and solvent site γ, respectively. Note that we exclusively use Gaussian units for electrostatics expressions in this work, i.e., all charges are given in units of √ e/ 4πε 0 (e: elementary charge, ε 0: vacuum permittivity). The short range potential is calculated in real space as, e.g., ) 12 ( ) Nα  ( σ  σαγ 6 αγ uγS(r) = 4 ε αγ  − |r − Rα |   |r − Rα | α=1  Nα  erfc(κ|r − Rα |) + qγ qα , (22) |r − Rα | α=1 where ε αγ and σαγ are the Lennard-Jones parameters and Nα is the number of solute sites (here and in Eq. (26), all Euclidean distances imply the minimum image convention). This part of the total potential can be safely truncated at a certain cutoff distance, which is not possible for real-space electrostatics without substantial quality loss. While the simple Ewald splitting is sufficient for uncharged systems, charged systems require another layer of renormalization. To this end, different methods have been devised. Kovalenko and Hirata proposed a 3D RISM adaptation

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.30.242.61 On: Sat, 30 May 2015 16:09:51

114107-4

J. Heil and S. M. Kast

J. Chem. Phys. 142, 114107 (2015)

of the well-known concept of an oppositely charged uniform background field for compensation41 and recently an extension for electrolyte solutions.30 In the first case, Eq. (14) becomes  ′ S S S (r)] − βuˆγL′(k)} χˆ γγ ℑ−1[{ℑ[cγ,i t γ, ′(|k|)] − cγ,i (r) i (r) = ′ γ

−

Na  qγ′ χˆ γγ′(k) 4π β * , qα + lim V ,α=1 - k→ 0 γ′ k2

(23)

where V is the cell volume. The monopole renormalization technique developed earlier by us has been shown to be applicable also to electrolytes.37 In this case, Eq. (14) is modified according to  S S ′ S t γ,i (r) = ℑ−1[{ℑ[cγ,i (r)] − βuˆγL′(k)} χˆ γγ ′(|k|)] − cγ,i (r) γ′  L(0) ˆ γγ′(k)], (24) ℑ−1 −β 1D[uˆγ ′ (k) χ γ′

where ℑ−1 1D represents an inverse 1D Fourier-Bessel transform that can be computed using the fast Fourier-Bessel transform algorithm devised by Talman,42 as frequently used for 1DRISM calculations; uˆγL′(0)(k) is the long range reciprocal-space supercell potential for a single monopole point-charge at an arbitrary position and of the same magnitude as the total charge of the solute. This kind of renormalization implies that the radial monopole term is subtracted from the total long range supercell potential, u˜ˆγL (k) = uˆγL (k) − uˆγL(0)(k).

(25)

The convolution of the monopole term with the solvent susceptibility is performed separately and the resulting contribution is added back.37 Note that, for charged solutes, use of either method (23), or (24) and (25) are necessary during iterations, meaning that the PMS approach described below applies to either the full long range supercell potential for background renormalization or to the neutralized form (25) for monopole renormalization, whereas the NRS energy expression always contains the full supercell potential. In the following derivation of the PMS technique, we therefore do not distinguish between uˆ˜γL (k) and uˆγL (k), implying a priori application of one of the renormalization methods. Even though we now see that the long range real-space potential, given formally by uγL(r) = qγ

Nα  qα erf(κ|r − Rα |) , |r − Rα | α=1

(26)

need not be known in order to compute converged distribution functions, the problem shows up by considering the first term of the chemical potential expression (6) combined with (11),     1 µex,(1) = − β −1 ργ dr hγ (r)cγ (r) 2 γ     1 −1 S L = −β ργ dr hγ (r){cγ (r) − βuγ (r)} (27) 2 γ giving rise to substantial computational effort for large solutes to calculate the real-space potential over all grid points. This is different from the situation encountered, for instance, in MD simulations where the evaluation of the reciprocal-space long

range Ewald potential uˆγL(k) by fast summation techniques completely replaces the time-consuming computation of its real-space equivalent. A similar situation occurs with the expression for the average excess solute-solvent interaction energy,     ex,L Uuexv = ργ dr gγ (r) {uγS(r) + uγL(r)} = Uuex,S v + Uu v , γ

(28) where we are again faced with the difficulty of computing a real-space long range potential for the second term. B. Reciprocal-space evaluation of energies: NRS formalism

The key to evaluate µex,(1) in reciprocal space is to make use of Parseval’s theorem43 (the star denotes the conjugate of a complex number),   1 dk [ y(k) ˆ zˆ∗(k)] . (29) dr [ y(r) z(r)] = (2π)3 Applied to the definition of µex,(1) in Eq. (6) and substituting the 3D RISM equation (1) for the total correlation function yields16,44  ργ    1 ∗ ˆ γ (k) µex,(1) = − β −1 dk c ˆ (k) h γ 2 (2π)3 γ     1  1  cˆ∗ (k) ′ ′ c ˆ (k) χ ˆ (k) dk = − β −1 γ γγ 3   γ 2 (2π) γ γ′   (30) which holds, strictly speaking, only in the limit of an infinitely large integration domain. Upon replacing the direct correlation function by the renormalization relation (11), we obtain a sum of three terms, µex,(1) = µex,(1a) + µex,(1b) + µex,(1c)     1 −1  1  S∗ S ′ =− β dk cˆγ (k) cˆγ′(k) χˆ γγ (k) 3 2 (2π)   γ γ′      1   S∗ L ′(k)  dk c ˆ (k) u ˆ ˆ + ′(k) χ γγ γ γ   (2π)3 γ γ′        1 1  L∗ L − β dk uˆγ (k) uˆγ′(k) χˆ γγ′(k) , 2 γ (2π)3   ′ γ (31) with uˆγL(k) = qγ ϕ(k) ˆ

Nα 

qα exp(ik · Rα )

(32)

α=1

being the Ewald long range potential function in reciprocal space. Since cγS is short-ranged, µex,(1a) and µex,(1b) can be directly evaluated in reciprocal space using the Ewald reciprocal-space potential already computed for the iteration procedure. The problematic asymptotic behavior of µex is dominated by µex,(1c). This last term can be recast in the form of a radially

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.30.242.61 On: Sat, 30 May 2015 16:09:51

114107-5

J. Heil and S. M. Kast

J. Chem. Phys. 142, 114107 (2015)

symmetric effective pairwise interaction between solute sites. To this end, we use the Fourier shift-theorem  dr[ f (r − r0) exp(−i k · r)] = exp(−i k · r0) ℑ [ f (r)] (33) to write µex,(1c) as 

Nα  Nα  dk  qα qα′qγ qγ′ α=1 α′=1  × exp{i k · (Rα − Rα′)} ϕˆ2(k) χˆ γγ′(k)   Nα  Nα 1  1 dk qα qα′ =− β 2 α=1 α′=1 (2π)3     2 qγ qγ′ χˆ γγ′(k) × exp{i k · (Rα − Rα′)} ϕˆ (k)   ′ γ,γ Nα  Nα 1  =− β qα qα′θ (|Rα − Rα′|). (34) 2 α=1 α′=1

1  1 µex,(1c) = − β 2 γ,γ′ (2π)3

θ(r) can be computed by a Fourier-Bessel transform via ∞  1 dk [ν(k,r)] qγ qγ′ χˆ γγ′(k)] = 2 θ(r) = ℑ [ϕˆ2(k) 2π γ,γ ′ 0

(35) with ν(k,r) = k 2 j0(k,r) ϕˆ2(k)



qγ qγ′ χˆ γγ′(k),

(36)

γ,γ ′

where j0(k,r) = sin(kr)/(kr) is the zero-order Bessel function of the first kind. Notice that the integrals occurring in Eqs. (34)–(36) are real which is easily seen if one applies Euler’s formula to the complex exponential and integrates each term in the resulting sum of trigonometric functions separately. Integrals over the odd sine functions vanish due to symmetry while the cosines are even functions and therefore their Fourier transforms are real. Numerical evaluation of Eq. (35) on logarithmic 1D grids, that are frequently employed and cannot start at k = 0, requires special attention. Focusing first on the off-diagonal elements with α , α ′, the integral must be split into two parts, the first one running from k = 0 to the first gridpoint of the logarithmic k grid, k0, and the second from k 0 to infinity k0 1 θ(r) = θ 1(r) + θ 2(r) = 2 dk[ν(k,r)] 2π 0

+

1 2π 2

∞ dk[ν(k,r)].

(37)

k0

For small enough grid spacing, θ 1(r) can be safely approximated by simple rectangular rule integration of ν(k) according to k0 ν(k0,r) . (38) θ 1(r) ≈ 2π 2 θ 2(r) can in principle be obtained in one step by computing the fast Fourier-Bessel transform according to Talman.42 However,

FIG. 1. θ(r ), defined in Eqs. (35)–(38), for modified SPC/E water (see Sec. III B) as a function of the site-site distance r from Talman’s fast FourierBessel transform and from analytical evaluation of the Fourier integral of a spline interpolant.

in practice, we find highly oscillatory behavior for small distances, as shown in Fig. 1. The numerically safe alternative is to construct a natural cubic B-spline interpolation of ν(k,r) and to compute the Fourier integral of the interpolant analytically, which yields θ 2(r) on the logarithmic real-space grid. This function is analogously interpolated and the resulting coefficients are stored in memory. This allows for very fast evaluation of the offdiagonal contribution to µex,(1c), since instead of using large 3D grids containing floating-point numbers, we need only to handle a few spline coefficients for each unique γγ ′ pair which can be precomputed for a given solvent model. For the diagonal elements with α = α ′, the argument in θ(r) approaches zero for small distances. In this case, we simply take the first 5 data-points of θ 2(r) (computed as described above) and use 4th order polynomial extrapolation to zero. Although evaluation of Eq. (34) implies quadratic scaling with the number of solute sites, we find in practice that the combination of precomputed data and spline evaluation generates negligible overhead compared to other steps of a complete 3D RISM calculation even for large systems. Further speedup can in principle be gained by truncating the function θ(r) and/or representing it on an adaptive irregular grid with less nodes. Analogous application of this formalism to the average interaction energy is also straightforward. The long range part of U ex is rewritten as     U ex,L = ργ dr (hγ (r) + 1) uγL(r) γ

 ργ    ˆ γ∗ (k) + 1) ˆ uˆγL(k) = dk ( h (2π)3 γ   1     L∗ S L = dk uˆγ (k) { cˆγ′(k)− βuˆγ′(k)} χˆ γγ′(k) 3 (2π)   γ γ′         1  S∗ L  ′   = dk c ˆ (k) u ˆ (k) χ ˆ (k) ′ γγ γ  γ  γ (2π)3   ′ γ           −β dk uˆγL∗(k) uˆγL′(k) χˆ γγ′(k)      γ γ′  = µex,(1b) + 2µex,(1c) (39)

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.30.242.61 On: Sat, 30 May 2015 16:09:51

114107-6

J. Heil and S. M. Kast

J. Chem. Phys. 142, 114107 (2015)

which allows the reuse of terms already computed for the chemical potential. The special appeal of the NRS method is that the number of memory accesses is typically reduced by orders of magnitude compared to evaluation of µex,(1) in Eq. (6) by direct evaluation/summation on a 3D grid due to the properties of the logarithmic grid which can also easily be made large enough that the assumption of an infinitely large box, as made for the Parseval relation, is warranted. The suitability of such grids for fluid integral equations has recently been highlighted by Löwen and coworkers.45 C. Fast evaluation of the reciprocal-space potential: PMS formalism

Since the calculation of the supercell long range realspace potential can be avoided as detailed above and the short range real-space Ewald potential can be truncated without further complications, the long range reciprocal-space supercell potential, Eq. (18), is the remaining time-consuming task in step (a) outlined in the Introduction. We now proceed to adapt the FFT-based PME algorithm39 in the context of the 3D RISM method, exploiting computational advantages similar to the massive efficiency gain in the context of MD simulations. For background information on particle-mesh methods, the reader is referred to the instructive treatise by Deserno and Holm46 or the book by Griebel et al.47 Our choice of PME over related schemes such as the particle-particleparticle mesh (P3M)48 or the “smooth” PME (SPME49) was guided by the reference formulation of the reciprocal-space potential we use for benchmarking the PMS approximation, where the (smeared) continuum Coulomb Green’s function specified in Eqs. (18)–(20) is expanded on the grid. In this case, as shown in Ref. 46, the Lagrange interpolation scheme associated with PME is the best approximation. However, calculations of energies and forces in MD simulations benefit from alternative interpolation variants (P3M or SPME) which in turn require an appropriate modification of the Coulomb Green’s function.46 Since we want to compare PMS results with data originating from a discretized continuum Coulomb Green’s function that also directly enters the NRS formalism, PME represents the proper framework. We start by decomposing the structure factor sˆ(k) according to sˆ(k, {Rα}) =

Nα 

qα exp(−ik x Rα, x )

α=1

× exp(−ik y Rα, y ) exp(−ik z Rα, z ),

(40)

where the complex exponentials can be approximated, separately for each direction, by a linear combination of Lagrange p basis polynomials L m of order p multiplied by exponentials evaluated at p equidistantly spaced grid points x m in the form of

exp(−ik x) ≈

p  m=0

with analogous expressions for other dimensions. The Lagrange basis polynomials of order p ∈ N0 are defined by p Lm (x) =

p  j=0, j,m

x − xj , xm − x j

(42)

with m = [0, p] ∈ N0 referencing grid points embedding a specific charge site, implicitly assuming that the empty product evaluates to one, and with ( p)  ′  x + d m − ⇔ p even  x   2 ( ) (43) xm =  p−1  ′   ⇔ p odd .  x + dx m − 2  Here, x ′ denotes the next nearest grid point to Rα, x if the polynomial order p is even, or the grid point corresponding to the next largest previous index if p is odd, and d x is the spacing between the grid points. That means that {x m } span a small subgrid around the charge sites Rα, x . Equation (41) effectively constitutes an “inverse” interpolation problem where we do not seek an interpolated function value between nodes but rather the function values at the nodes that, upon interpolation, give rise to some value at a given position between the nodes. In other words, the interpolation function serves to spread the irregularly placed solute site charges over the regular grid in order to allow for an FFT-based calculation of the structure factor. The approximate structure factor finally reads sˆ(k, {Rα}) ≈

Nα 



qα

α=1

p (Rα, y ) L lp (Rα, x )L m

l, m, n p L n (Rα, z ) exp(−i k

=

× 

· rl, m, n )

ψl, m, n ({Rα}) exp(−i k · rl, m, n ),

(44)

l, m, n

where ψl, m, n ({Rα}) is the grid-based “smeared” charge density of all solute sites defined as ψl, m, n ({Rα}) =

Nα  α=1

p qα L lp (Rα, x )L m (Rα, y )L np (Rα, z ).

(45)

ψl, m, n ({Rα}) can be quickly evaluated with a complexity of O(Nα ) due to the compact support of the Lagrange basis polynomials, i.e., we need only be concerned with those indices l, m, n that correspond to grid points near one of the solute sites, values for all others are zero. The basis polynomials themselves are not calculated directly from Eq. (42). Rather, the product is expanded as shown by Petersen50 who also tabulated orders 2 through 6. We additionally implemented orders up to 10, as explicitly worked out in the Appendix (Table III). In practice, only even orders are applicable due to a discontinuity for odd orders, as illustrated in Fig. 1 of Ref. 46. Note that the evaluation is subject to periodic boundary conditions. Upon insertion of Eqs. (44) and (45) into Eq. (18), we finally arrive at the expression  uˆγL (k) ≈ qγ ϕ(k) ˆ ψl, m, n ({Rα}) exp(−i k · rl, m, n ) (46) l, m, n

p Lm (x) exp(−ik x m ),

(41)

for the reciprocal-space supercell potential energy function. We choose the wave vectors k as centered around k = (0, 0, 0)T

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.30.242.61 On: Sat, 30 May 2015 16:09:51

114107-7

J. Heil and S. M. Kast

J. Chem. Phys. 142, 114107 (2015)

which results in ( )T (o − n x /2) (p − n y /2) (q − nz /2) ko, p,q = 2π , , , nx d x nyd y nz d z

(47)

with n x , n y , nz and d x , d y , d z being the number of grid points and grid spacing in x, y, and z direction, respectively, and o, p, q = [0, n x − 1] ∈ N0. When, typically, the solute is centered at the origin, we need to shift all solute sites by half the box length in all three directions such that the simulation cell resides completely inside the first quadrant. With rl, m, n = (l d x , m d y , n d z, )T, we finally recover the usual form of the discrete Fourier transform as a prerequisite for an application of the FFT  uˆγL(ko, p,q ) ≈ qγ ϕ(|k ˆ o, p,q |) ψl, m, n ({Rα}) l, m, n

× exp(−i ko, p,q · rl, m, n ) = qγ ϕ(|k ˆ o, p,q |)

y −1 n n z −1 x −1 n

in this case, we use 0.1 M aqueous KCl solution in addition to water as solvent in order to examine the consistency of the approach for both pure dielectric and electrolyte solvent environments. To assess the impact of the PMS approximation on µex as well as possible side-effects when applying PMS together with NRS for a larger system, we study statistically the conformational µex differences of 10 snapshots taken from a short simulated annealing MD simulation of the small protein thioredoxin in water. Here, the excess chemical potential was computed for each snapshot using PMS with even order Lagrange polynomials between 2 and 10, with and without NRS. In each case, mean absolute error (MAE) and (relative) “mean absolute percentage” error (MAPE) for the method studied with respect to reference data (real-space electrostatics with monopole renormalization,37 RS-M) were calculated by MAE(µex) =

ψl, m, n ({Rα})

) (l=0 m=0 n=0 ol pm qn + + × exp −2πi nx ny nz 

= qγ ϕ(|k ˆ o, p,q |)FFT[ψ(rl, m, n )].

2 Nsn (Nsn + 1) N sn  ex ex ex × | µex i − µ j − (µ i,ref − µ j,ref )|, (49) j >i

(48)

Thus, the interpolation and consequential suitability for an FFT treatment allows us to compute the long range potential with O (n x n y nz log(n x n y nz )) instead of O(Nα n x n y nz ) computational complexity required for direct evaluation of Eq. (18). In practice, for molecules with 104-105 atoms, the computational effort is barely significant in comparison with other potential terms, as we will show below.

III. NUMERICAL EXAMPLES A. Model systems

We now illustrate the performance characteristics of the NRS and PMS implementations with respect to speed and accuracy for several benchmark systems. We first conducted a test of the NRS method to show that, in the limit of infinitely large boxes and small grid spacing, the new formulation, Eqs. (31)–(38), is equivalent to direct evaluation of Eq. (6) on a 3D grid. To this end, the aqueous excess chemical potential of a spherical particle was calculated using both 3D and 1D RISM theory as a reference. Spherical symmetry implies that 1D and 3D RISM results for µex differ only due to different errors associated with finite grid resolution. These errors are typically very small for 1D RISM as long as suitable logarithmically spaced grids are used. Note, however, that using a single particle means that there are no off-diagonal elements in Eq. (34). Using this benchmark system, we investigated the behavior of NRS with respect to box size and grid spacing in comparison with real-space electrostatics combined with monopole renormalization37 (RS-M). As a more application-oriented example, we performed 3D RISM calculations on the acetate anion and the methylammonium as a function of box size for fixed grid spacing of 0.2 Å. Since molecular ions are anisotropic, the results are also influenced by off-diagonal elements in Eq. (34). Furthermore,

MAPE(µ ) = ex

2 Nsn (Nsn − 1) ex N s n µex − µex − (µex  j i i,ref − µ j,ref ) × | | (50) ex µex i,ref − µ j,ref j >i

for the excess chemical potential and analogously for the excess internal energy. Nsn is the number of simulation snapshots and µex i,ref is the excess chemical potential for snapshot i of the corresponding reference calculation without PMS, i.e., with RS-M. This was done for various grid spacings between 0.3 and 0.6 Å, which is the range typically employed in 3D RISM calculations, and using Lagrange polynomials of various orders for the PMS charge spreading procedure. Finally, timing experiments were carried out on a complex of the large poly-ADP-ribose-polymerase (PARP) protein with a small-molecule ligand in aqueous solution. This system was also used to visually inspect distribution functions obtained with and without the PMS approximation. B. Computational details

All benchmarks were run on a dual socket 6-core Intel Xeon X5690 system using the 1D/3D RISM software developed in our group. Large parts of the code are OpenMPparallel, notably the potential calculation and convergence accelerator routines, the 3D convolution products, and the closure evaluation. Fast Fourier transforms were computed using the shared-memory-parallel Intel MKL FFT routines using the FFTW wrapper functionality. The solvent site-site susceptibility needed as input for all RISM calculations was computed using the DRISM formalism8,9 with the HNC closure on a logarithmically spaced grid of 512 points reaching from 0.0059 to 164.02 Å. The equations were converged to a threshold of 10−8 for the maximum residual of the direct correlation functions. The temperature was set to T = 298.15 K and the dielectric

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.30.242.61 On: Sat, 30 May 2015 16:09:51

114107-8

J. Heil and S. M. Kast

constant of the solvents to ε r = 78. For pure water, the density was ρ = 0.033 34 Å−3. We used a modified SPC/E water model where the hydrogen Lennard-Jones parameter σ was set to 1.0 Å.51–53 For the 0.1 M KCl solution, we additionally used ion parameters from earlier studies19,20,54 and site densities of water and ions of 0.033 24 and 0.000 06 Å−3, respectively. The spherical test ion for studying NRS in comparison with 1D RISM data had a charge of q = ±1e and LennardJones parameters of ε = 0.301 11 kJ mol−1 and σ = 2 Å. For investigating the effect of the grid spacing, a fixed box size of 20 Å in each direction was used. Lennard-Jones parameters for the acetate anion and the methylammonium cation were taken from the general AMBER force field (GAFF)59 while the charges were derived from an (embedded cluster) EC-RISM25 calculation utilizing the MP2/6-31G/PSE-3 level of theory and with the same pure water solvent model as detailed above. All 1D/3D RISM calculations on the spherical ion, the acetate anion, and the methylammonium cation were performed with the PSE-3 closure. The MD snapshots of the thioredoxin protein (PDB code 1XOA) were generated with AMBER 12 PMEMD55 running a total of 1.2 ns with a 2 fs time step at 300 K using the Onufriev-Bashford-Case (OBC) variant of the generalized Born method56,57 with a Langevin-type thermostat and a collision frequency of 2 ps−1. After running for 100 ps, the system was cooled to 0 K in 20 ps and a sample was extracted. This was repeated 9 times. For subsequent 3D RISM calculations, the PSE-2 closure was used, and the total box size was kept constant at 72 Å in each direction while varying the number of grid points and the grid spacing. The AMBER 99 SB force field58 was used for the MD simulation as well as

J. Chem. Phys. 142, 114107 (2015)

for the subsequent 3D RISM calculations, and also for the timing experiments on the PARP complex (PDB code 4PAX, 5604 atoms, charge q = 2e), together with a GAFF/AM1BCC-parameterization for the ligand.59–61 Here, we used a grid of 240 × 202 × 218 points, a grid spacing of 0.4 Å (corresponding to a total of 30 Å solvent buffer in addition to the size of the complex) and the PSE-1(KH) closure. The reduction of the 3D distribution functions to 1D radial distribution functions was accomplished using an angular Lebedev-Laikov grid with 770 points62 on the logarithmic radial grid. For all 3D RISM calculations performed in this work, a tight convergence threshold of 10−6 in the maximum norm of the residual of the direct correlation functions was employed. Except for the protein systems, no real-space cutoff was applied to interaction potentials. All necessary data for reproducing the benchmarks are provided as supplementary material.63 C. Results and discussion

The results for the spherical ions in water (Fig. 2, Table I) and acetate and methylammonium ion cases in aqueous and electrolyte solution (Fig. 3) indicate that the NRS method offers an accuracy that is superior to the pure monopole renormalization37 for smaller boxes and is competitive with respect to the uniform charged background method.41 We furthermore see that, in particular, the real-space monopole method but, to lesser degree, also the real-space uniform charged background method exhibit strong oscillations, whereas their NRS counterparts show a smoother convergence behavior.

FIG. 2. µ ex ((a) and (b)) and deviation from reference values ((c) and (d)) as a function of box size (grid spacing: 0.05 Å) for the different renormalization schemes for spherical ions in water ((A)/(C): anion, (B)/(D): cation). “monop. renorm.” refers to the monopole renormalization method,37 “bg renorm” stands for the background charge method of Hirata and Kovalenko.41

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.30.242.61 On: Sat, 30 May 2015 16:09:51

114107-9

J. Heil and S. M. Kast

J. Chem. Phys. 142, 114107 (2015)

TABLE I. Deviations of µ ex from 1D RISM reference calculations as a function of grid spacing (fixed cubic box length of 28 Å) for the different renormalization schemes in the spherical ion case. Spacing/Å 0.05 0.1 0.2 0.4 0.8

ex ∆µ RS-M /(kJ mol−1)

ex /(kJ mol−1) ∆µ NRS

ex ∆URS−M /(kJ mol−1)

ex /(kJ mol−1) ∆UNRS

−0.13 −0.05 1.10 0.96 −20.58

0.05 0.13 1.29 1.15 −11.36

−0.28 −0.17 1.92 1.66 −40.09

0.09 0.20 2.29 2.04 −21.65

Whether the NRS method is coupled to the uniform charged background method or the monopole renormalization is obviously unimportant. However, coupling NRS to the monopole renormalization method gives us the opportunity to calculate correctly renormalized distribution functions and excess chemical potentials for electrolytes from a single calculation, which is not the case for the uniform charged background method. Therefore, it is important to test whether NRS also performs well for electrolytes, which is evident from Fig. 3(b)/3(d). From Table I, we see that the monopole renormalization/NRS method also systematically converges towards the 1D RISM reference value of the excess chemical potential when the grid spacing becomes smaller. A similar trend is observed for the solute-solvent interaction energy. The statistical assessment of the excess chemical potentials of the thioredoxin protein conformations shown in Table II allows us to investigate the effect of the PMS approximation. One can see that particularly for polynomial orders 8 and 10, the errors in µex are highly dependent on the grid spacing, being negligible for 0.4 Å and less and still small up to 0.5 Å. The error in the excess internal energy is typically twice as large as for µex, but still mostly acceptable.

Furthermore, we note the effect of error amplification between our two methods. The particle-mesh approximation enters into the NRS method by means of the 3D reciprocal space potential in the µex,(1b) term in Eq. (31) and strongly affects the accuracy of the combined NRS/PMS method on larger grids. Of course the question arises how large the impact of the particle-mesh approximation on the distribution functions is. Fig. 4 for 3D distributions and Fig. 5 for selected radially integrated site pair distribution functions show that, already for 2nd order Lagrange polynomials, the difference is very small. We also conclude that 10th order Lagrange polynomials offer good accuracy and that we could probably use even higher order polynomials. However, in 1901, Runge showed that the combination of high-order interpolating polynomials and equispaced interpolation points, as used for 3D RISM by necessity as a prerequisite for the Fast Fourier transform, can lead to serious interpolation errors under unfavorable circumstances.64 While we have not yet encountered a case where this specific problem arises in the context of the PMS method, we cannot rule out a priori that such a case exists, and the probability to encounter a situation where

FIG. 3. µ ex as a function of grid size for the different renormalization schemes for the acetate anion ((A): water, (B): 0.1 M KCl) and the methylammonium cation ((C): water, (D): 0.1 M KCl). For abbreviations, see Fig. 2.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.30.242.61 On: Sat, 30 May 2015 16:09:51

114107-10

J. Heil and S. M. Kast

J. Chem. Phys. 142, 114107 (2015)

TABLE II. MAE (relative MAPE in parentheses) for excess chemical potential and interaction energy for thioredoxin conformers (12.0 Å real-space cutoff for short-range interactions in all variants). Polynomial order 10 10 10 10 8 8 8 8 6 6 6 6 4 4 4 4 2 2 2 2

Spacing/Å 0.3 0.4 0.5 0.6 0.3 0.4 0.5 0.6 0.3 0.4 0.5 0.6 0.3 0.4 0.5 0.6 0.3 0.4 0.5 0.6

ex MAE(µ PMS/RS−M )/(kJ mol−1)

ex MAE(µ PMS/NRS )/(kJ mol−1)

ex MAE(UPMS/RS−M )/(kJ mol−1)

ex MAE(UPMS/NRS )/(kJ mol−1)

0.00 (0.01) 0.01 (0.02) 0.02 (0.04) 0.10 (0.14) 0.01 (0.01) 0.01 (0.01) 0.04 (0.08) 0.18 (0.26) 0.01 (0.05) 0.03 (0.04) 0.08 (0.14) 0.30 (0.40) 0.04 (0.33) 0.15 (0.30) 0.24 (0.22) 0.68 (0.58) 0.60 (3.48) 1.13 (1.88) 2.09 (2.27) 4.06 (6.28)

0.01 (0.01) 0.10 (0.27) 1.38 (2.86) 6.53 (15.03) 0.06 (0.08) 0.45 (1.36) 4.28 (7.97) 13.76 (32.17) 0.47 (0.63) 2.40 (7.34) 14.89 (24.96) 31.60 (73.94) 4.31 (5.93) 15.10 (52.80) 60.01 (94.66) 80.67 (182.16) 52.37 (78.40) 120.04 (559.76) 297.19 (461.17) 259.29 (497.60)

0.01 (0.01) 0.01 (0.02) 0.04 (0.05) 0.21 (0.10) 0.01 (0.02) 0.01 (0.01) 0.08 (0.10) 0.33 (0.19) 0.02 (0.02) 0.06 (0.08) 0.15 (0.13) 0.55 (0.36) 0.08 (0.07) 0.33 (0.55) 0.52 (0.34) 1.48 (0.74) 1.23 (0.91) 2.46 (4.22) 4.44 (3.11) 8.10 (4.58)

0.02 (0.06) 0.21 (0.10) 2.77 (5.31) 13.07 (10.16) 0.12 (0.29) 0.90 (0.66) 8.56 (14.14) 27.55 (21.97) 0.95 (1.37) 4.80 (5.23) 29.77 (42.37) 63.19 (50.75) 8.61 (5.12) 30.18 (41.14) 119.96 (157.26) 161.31 (126.32) 104.69 (144.25) 239.99 (375.03) 594.18 (771.38) 518.38 (366.07)

our PMS approach is inappropriate would certainly increase with rising polynomial order. More specifically, we should in principle be able to construct such a case by implementing and using polynomials of ever higher order which, of course, quickly becomes a tedious and error-prone task in practice. Additionally, one should bear in mind that at some point, the computational advantage gained through the local support of the polynomials is lost because the number of support grid points used for the charge spreading step grows cubically with the polynomial order.

FIG. 4. 3D water oxygen distribution function around the PARP complex (PDB code 4PAX). Purple: Ewald summation, no real-space cutoff, orange: PMS with 2nd order Lagrange polynomials, and a 12.0 Å real-space cutoff distance.

The timing experiments on the PARP complex, as illustrated in Fig. 6, provide a means to assess real-world speed improvements. On 6 CPU cores, the time needed

FIG. 5. Radial distribution functions of water oxygen (red) and hydrogen (blue) around (a) the ligand’s hydroxyl oxygen and (b) the N-terminal nitrogen atom in the PARP complex. Insets show the difference between reference data without PMS/cutoff and results obtained with PMS (2nd order) and a 12 Å short-range cutoff distance.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.30.242.61 On: Sat, 30 May 2015 16:09:51

114107-11

J. Heil and S. M. Kast

FIG. 6. Wall-time distribution between 3D RISM iteration effort and other contributions (dominated by steps (a) and (c) as outlined in the Introduction) for the PARP complex. A: Full real space potential/Ewald summation. B: NRS/PMS. C: NRS/PMS (10th order) with an additional 12 Å cutoff for Lennard-Jones potential and the short range real space-part of the Coulomb interaction.

for steps (a) and (c) outlined in the Introduction drops from ca. 18.5 min down to 6.7 min upon using NRS/PMS. If we additionally truncate the short-range part of the potential with a 12 Å cutoff, only ∼1.6 min are needed (though at the prize of some numerical deviations of thermodynamic quantities). The time required for computing electrostatic potentials drops by two orders of magnitude by using PMS. This represents a major reduction of total wall time, since the time needed for the 3D RISM iteration is about 21 min when using our very strict convergence threshold, which we could probably lower quite a bit for several purposes, e.g., qualitative visual analysis of distribution functions. In this case, the relative speedup achieved by accelerating the potential calculation would be even more important.

IV. CONCLUDING REMARKS

We have developed the theory for an efficient treatment of electrostatics in 3D RISM theory, consisting of a particlemesh technique (PMS) for reciprocal space potentials to be combined with a formalism (NRS) to avoid any reference

J. Chem. Phys. 142, 114107 (2015)

to costly real-space potential calculations for energy and chemical potential. The methodology was benchmarked on a number of test systems with respect to accuracy and computational efficiency. Given that potential calculation gives rise to a large fraction of the computational burden of 3D RISM calculations for large molecules and dense grids, this represents a substantial reduction of wall-time with acceptable loss of accuracy. The PMS method makes the computational effort for the calculation of the reciprocalspace supercell potential energy function insignificant, which benefits all kinds of large-scale 3D RISM calculations. 3D RISM solvation free energy calculations additionally benefit from the NRS technique, which removes nearly all overhead associated with computing the excess chemical potential. In combination with our monopole renormalization method developed earlier, the new methodology makes it possible to compute accurate excess chemical potentials together with the correct correlation functions in ionic solutions very efficiently, while also removing much of the spurious oscillation of µex with respect to the box size. It remains to be seen if more advanced fast grid summation schemes such as the P2NFFT65 method or SPME,49 which is basically equivalent to P3M,66 can be used in the context of 3D RISM calculations to gain an even more favorable accuracy to speed ratio and a potentially less pronounced grid spacing dependence. While it is in principle straightforward to implement alternative interpolation schemes, we will need to test the resulting modified reciprocal space potential with respect to both, liquid structure and solution thermodynamics. These questions will be addressed in forthcoming work. ACKNOWLEDGMENTS

We thank the Deutsche Forschungsgemeinschaft for financial and the IT and Media Center (ITMC) of the TU Dortmund as well as the Regionales Rechenzentrum Erlangen (RRZE) for computational support. Karen Kuhn and Thomas Kloss are greatly acknowledged for useful discussions, Christian Holm for important comments on the manuscript ahead of submission.

APPENDIX: LAGRANGE CHARGE ASSIGNMENT EXPRESSIONS

Table III shows expressions for the Lagrange charge assignment schemes for higher orders. Orders 2-6 can be found in Ref. 50. TABLE III. Lagrange charge assignment schemes for orders 7 through 10. Note that we define u to be the relative position of a unit charge with respect to the nearest grid point for even p or the previous grid point to the left for odd p, denoted by x ′, respectively. It follows, therefore, that u = x − x ′ and, for a given p p order p, L m (x) = Wm−p/2(u). For illustration, see also Fig. 8 in Ref. 46. Order 7

Expression 7 (u) = (−1575d 7 + 450d 6 u + 7252d 5 u 2 − 2072d 4 u 3 − 3920d 3 u 4 + 1120d 2 u 5 + 448d u 6 − 128u 7)/(645120d 7 ) W−7/2 x x x x x x x x 7 (u) = (2205d 7 − 882d 6 u − 9980d 5 u 2 + 3992d 4 u 3 + 4720d 3 u 4 − 1888d 2 u 5 − 320d u 6 + 128u 7)/(92160d 7 ) W−5/2 x x x x x x x x 7 (u) = (−3675d 7 + 2450d 6 u + 15588d 5 u 2 − 10392d 4 u 3 − 3600d 3 u 4 + 2400d 2 u 5 + 192d u 6 − 128u 7)/(30720d 7 ) W−3/2 x x x x x x x x 7 (u) = (11025d 7 − 22050d 6 u − 7564d 5 u 2 + 15128d 4 u 3 + 1328d 3 u 4 − 2656d 2 u 5 − 64d u 6 + 128u 7)/(18432d 7 ) W−1/2 x x x x x x x x 7 (u) = (11025d 7 + 22050d 6 u − 7564d 5 u 2 − 15128d 4 u 3 + 1328d 3 u 4 + 2656d 2 u 5 − 64d u 6 − 128u 7)/(18432d 7 ) W1/2 x x x x x x x x

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.30.242.61 On: Sat, 30 May 2015 16:09:51

114107-12

J. Heil and S. M. Kast

J. Chem. Phys. 142, 114107 (2015)

TABLE III. (Continued.) Order

Expression 7 (u) = (−3675d 7 − 2450d 6 u + 15588d 5 u 2 + 10392d 4 u 3 − 3600d 3 u 4 − 2400d 2 u 5 + 192d u 6 + 128u 7)/(30720d 7 ) W3/2 x x x x x x x x 7 (u) = (2205d 7 + 882d 6 u − 9980d 5 u 2 − 3992d 4 u 3 + 4720d 3 u 4 + 1888d 2 u 5 − 320d u 6 − 128u 7)/(92160d 7 ) W5/2 x x x x x x x x

8

7 (u) = (−1575d 7 − 450d 6 u + 7252d 5 u 2 + 2072d 4 u 3 − 3920d 3 u 4 − 1120d 2 u 5 + 448d 6 u 6 + 128u 7)/(645120d 7 ) W7/2 x x x x x x x x 8 (u) = (144d 7 u − 36d 6 u 2 − 196d 5 u 3 + 49d 4 u 4 + 56d 3 u 5 − 14d 2 u 6 − 4d u 7 + u 8)/(40320d 8 ) W−4 x x x x x x x x 8 (u) = (−192d 7 u + 64d 6 u 2 + 252d 5 u 3 − 84d 4 u 4 − 63d 3 u 5 + 21d 2 u 6 + 3d u 7 − u 8)/(5040d 8 ) W−3 x x x x x x x x 8 (u) = (288d 7 u − 144d 6 u 2 − 338d 5 u 3 + 169d 4 u 4 + 52d 3 u 5 − 26d 2 u 6 − 2d u 7 + u 8)/(1440d 8 ) W−2 x x x x x x x x 8 (u) = (−576d 7 u + 576d 6 u 2 + 244d 5 u 3 − 244d 4 u 4 − 29d 3 u 5 + 29d 2 u 6 + d u 7 − u 8)/(720d 8 ) W−1 x x x x x x x x

W08(u) = (576d 8x u − 820d 6x u 2 + 273d 4x u 4 − 30d 2x u 6 + u 8)/(576d 8x ) 8 W1 (u) = (576d 7x u + 576d 6x u 2 − 244d 5x u 3 − 244d 4x u 4 + 29d 3x u 5 + 29d 2x u 6 − d x u 7 − u 8)/(720d 8x ) W28(u) = (−288d 7x u − 144d 6x u 2 + 338d 5x u 3 + 169d 4x u 4 − 52d 3x u 5 − 26d 2x u 6 + 2d x u 7 + u 8)/(1440d 8x ) W38(u) = (192d 7x u + 64d 6x u 2 − 252d 5x u 3 − 84d 4x u 4 + 63d 3x u 5 + 21d 2x u 6 − 3d x u 7 − u 8)/(5040d 8x ) W48(u) = (−144d 7x u − 36d 6x u 2 + 196d 5x u 3 + 49d 4x u 4 − 56d 3x u 5 − 14d 2x u 6 + 4d x u 7 + u 8)/(40320d 8x ) 9

9 (u) = (99225d 9 − 22050d 8 u − 464976d 7 u 2 + 103328d 6 u 3 + 284256d 5 u 4 − 63168d 4 u 5 W−9/2 x x x x x x

− 48384d 3x u 6 + 10752d 2x u 7 + 2304d x u 8 − 512u 9)/(185794560d 9x ) 9 (u) = (−127575d 9 + 36450d 8 u + 593712d 7 u 2 − 169632d 6 u 3 − 346528d 5 u 4 + 99008d 4 u 5 + 51968d 3 u 6 W−7/2 x x x x x x x

− 14848d 2x u 7 − 1792d x u 8 + 512u 9)/(20643840d 9x ) 9 (u) = (178605d 9 − 71442d 8 u − 817200d 7 u 2 + 326880d 6 u 3 + 422240d 5 u 4 − 168896d 4 u 5 W−5/2 x x x x x x − 44800d 3x u 6 + 17920d 2x u 7 + 1280d x u 8 − 512u 9)/(5160960d 9x ) 9 (u) = (−297675d 9 + 198450d 8 u + 1277328d 7 u 2 − 851552d 6 u 3 − 353952d 5 u 4 + 235968d 4 u 5 W−3/2 x x x x x x

+ 29952d 3x u 6 − 19968d 2x u 7 − 768d x u 8 + 512u 9)/(2211840d 9x ) 9 (u) = (893025d 9 − 1786050d 8 u − 656784d 7 u 2 + 1313568d 6 u 3 + 137824d 5 u 4 − 275648d 4 u 5 W−1/2 x x x x x x

− 10496d 3x u 6 + 20992d 2x u 7 + 256d x u 8 − 512u 9)/(1474560d 9x ) 9 (u) = (893025d 9 + 1786050d 8 u − 656784d 7 u 2 − 1313568d 6 u 3 + 137824d 5 u 4 + 275648d 4 u 5 W1/2 x x x x x x − 10496d 3x u 6 − 20992d 2x u 7 + 256d x u 8 + 512u 9)/(1474560d 9x ) 9 (u) = (−297675d 9 − 198450d 8 u + 1277328d 7 u 2 + 851552d 6 u 3 − 353952d 5 u 4 − 235968d 4 u 5 W3/2 x x x x x x

+ 29952d 3x u 6 + 19968d 2x u 7 − 768d x u 8 − 512u 9)/(2211840d 9x ) 9 (u) = (178605d 9 + 71442d 8 u − 817200d 7 u 2 − 326880d 6 u 3 + 422240d 5 u 4 + 168896d 4 u 5 W5/2 x x x x x x

− 44800d 3x u 6 − 17920d 2x u 7 + 1280d x u 8 + 512u 9)/(5160960d 9x ) 9 (u) = (−127575d 9 − 36450d 8 u + 593712d 7 u 2 + 169632d 6 u 3 − 346528d 5 u 4 − 99008d 4 u 5 W7/2 x x x x x x + 51968d 3x u 6 + 14848d 2x u 7 − 1792d x u 8 − 512u 9)/(20643840d 9x ) 9 (u) = (99225d 9 + 22050d 8 u − 464976d 7 u 2 − 103328d 6 u 3 + 284256d 5 u 4 + 63168d 4 u 5 W9/2 x x x x x x

− 48384d 3x u 6 − 10752d 2x u 7 + 2304d x u 8 + 512u 9)/(185794560d 9x ) 10

10(u) = (−2880d 9 u + 576d 8 u 2 + 4100d 7 u 3 − 820d 6 u 4 − 1365d 5 u 5 + 273d 4 u 6 + 150d 3 u 7 − 30d 2 u 8 − 5d u 9 + u 10)/(3628800d 10) W−5 x x x x x x x x x x 10(u) = (3600d 9 u − 900d 8 u 2 − 5044d 7 u 3 + 1261d 6 u 4 + 1596d 5 u 5 − 399d 4 u 6 − 156d 3 u 7 + 39d 2 u 8 + 4d u 9 − u 10)/(362880d 10) W−4 x x x x x x x x x x 10(u) = (−4800d 9 u + 1600d 8 u 2 + 6492d 7 u 3 − 2164d 6 u 4 − 1827d 5 u 5 + 609d 4 u 6 + 138d 3 u 7 − 46d 2 u 8 − 3d u 9 + u 10)/(80640d 10) W−3 x x x x x x x x x x 10(u) = (7200d 9 u − 3600d 8 u 2 − 8738d 7 u 3 + 4369d 6 u 4 + 1638d 5 u 5 − 819d 4 u 6 − 102d 3 u 7 + 51d 2 u 8 + 2d u 9 − u 10)/(30240d 10) W−2 x x x x x x x x x x 10(u) = (−14400d 9 u + 14400d 8 u 2 + 6676d 7 u 3 − 6676d 6 u 4 − 969d 5 u 5 + 969d 4 u 6 + 54d 3 u 7 − 54d 2 u 8 − d u 9 + u 10)/(17280d 10) W−1 x x x x x x x x x x 8 2 6 4 4 6 2 8 10 10 W010(u) = (14400d 10 x − 21076d x u + 7645d x u − 1023d x u + 55d x u − u )/(14400d x ) 10 9 8 2 7 3 6 4 5 5 4 6 3 7 2 8 W1 (u) = (14400d x u + 14400d x u − 6676d x u − 6676d x u + 969d x u + 969d x u − 54d x u − 54d x u + d x u 9 + u 10)/(17280d 10 x) W210(u) = (−7200d 9x u − 3600d 8x u 2 + 8738d 7x u 3 + 4369d 6x u 4 − 1638d 5x u 5 − 819d 4x u 6 + 102d 3x u 7 + 51d 2x u 8 − 2d x u 9 − u 10)/(30240d 10 x)

W310(u) = (4800d 9x u + 1600d 8x u 2 − 6492d 7x u 3 − 2164d 6x u 4 + 1827d 5x u 5 + 609d 4x u 6 − 138d 3x u 7 − 46d 2x u 8 + 3d x u 9 + u 10)/(80640d 10 x) W410(u) = (−3600d 9x u − 900d 8x u 2 + 5044d 7x u 3 + 1261d 6x u 4 − 1596d 5x u 5 − 399d 4x u 6 + 156d 3x u 7 + 39d 2x u 8 − 4d x u 9 − u 10)/(362880d 10 x) W510(u) = (2880d 9x u + 576d 8x u 2 − 4100d 7x u 3 − 820d 6x u 4 + 1365d 5x u 5 + 273d 4x u 6 − 150d 3x u 7 − 30d 2x u 8 + 5d x u 9 + u 10)/(3628800d 10 x)

1J.

G. Kirkwood, J. Chem. Phys. 3, 300 (1935). W. Zwanzig, J. Chem. Phys. 22, 1420 (1954). 3J. P. Hansen and I. R. Mcdonald, Theory of Simple Liquids (Academic Press, London, 2006). 2R.

4R.

Ramirez, R. Gebauer, M. Mareschal, and D. Borgis, Phys. Rev. E 66, 031206 (2002). 5G. Jeanmairet, M. Levesque, R. Vuilleumier, and D. Borgis, J. Phys. Chem. Lett. 4, 619 (2013).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.30.242.61 On: Sat, 30 May 2015 16:09:51

114107-13 6D.

J. Heil and S. M. Kast

Beglov and B. Roux, J. Phys. Chem. B 101, 7821 (1997). Kovalenko and F. Hirata, Chem. Phys. Lett. 290, 237 (1998). 8J. S. Perkyns and B. M. Pettitt, Chem. Phys. Lett. 190, 626 (1992). 9J. S. Perkyns and B. M. Pettitt, J. Chem. Phys. 97, 7656 (1992). 10T. Morita and K. Hiroike, Prog. Theor. Phys. 23, 1003 (1960). 11S. M. Kast and D. Tomazic, J. Chem. Phys. 137, 171102 (2012). 12S. M. Kast and T. Kloss, J. Chem. Phys. 129, 236101 (2008). 13A. Kovalenko and F. Hirata, J. Chem. Phys. 110, 10095 (1999). 14S. Joung, T. Luchko, and D. A. Case, J. Chem. Phys. 138, 044103 (2013). 15R. Kjellander and S. Sarman, J. Chem. Phys. 90, 2768 (1989). 16S. M. Kast, Phys. Rev. E 67, 041203 (2003). 17S. Genheden, T. Luchko, S. Gusarov, A. Kovalenko, and U. Ryde, J. Phys. Chem. B 114, 8505 (2010). 18D. J. Sindhikara, N. Yoshida, and F. Hirata, J. Comput. Chem. 33, 1536 (2012). 19S. Tayefeh, T. Kloss, G. Thiel, B. Hertel, A. Moroni, and S. M. Kast, Biochemistry 46, 4826 (2007). 20S. Tayefeh, T. Kloss, M. Kreim, M. Gebhardt, D. Baumeister, B. Hertel, C. Richter, H. Schwalbe, A. Moroni, G. Thiel, and S. M. Kast, Biophys. J. 96, 485 (2009). 21S. M. Kast, T. Kloss, S. Tayefeh, and G. Thiel, J. Gen. Physiol. 138, 371 (2011). 22J. Tomasi, B. Mennucci, and R. Cammi, Chem. Rev. 105, 2999 (2005). 23A. Klamt and G. Schüürmann, J. Chem. Soc., Perkin Trans. 2, 799 (1993). 24N. Yoshida, Y. Kiyota, and F. Hirata, J. Mol. Liq. 159, 83 (2011). 25T. Kloss, J. Heil, and S. M. Kast, J. Phys. Chem. B 112, 4337 (2008). 26S. M. Kast, J. Heil, S. Güssregen, and K. F. Schmidt, J. Comput.-Aided Mol. Des. 24, 343 (2010). 27E. L. Ratkova and M. V. Fedorov, J. Chem. Theory Comput. 7, 1450 (2011). 28J. W. Kaminski, S. Gusarov, T. A. Wesolowski, and A. Kovalenko, J. Phys. Chem. A 114, 6082 (2010). 29A. Kovalenko, S. Ten-no, and F. Hirata, J. Comput. Chem. 20, 928 (1999). 30S. Gusarov, B. S. Pujari, and A. Kovalenko, J. Comput. Chem. 33, 1478 (2012). 31J. J. Howard, J. S. Perkyns, N. Choudhury, and B. M. Pettitt, J. Chem. Theory Comput. 4, 1928 (2008). 32Y. Maruyama and F. Hirata, J. Chem. Theory Comput. 8, 3015 (2012). 33V. P. Sergiievskyi and M. V. Fedorov, J. Chem. Theory Comput. 8, 2062 (2012). 34P. Pulay, Chem. Phys. Lett. 73, 393 (1980). 35R. C. Harris, A. H. Boschitsch, and M. O. Fenley, J. Chem. Theory Comput. 9, 3677 (2013). 36T. Luchko, S. Gusarov, D. R. Roe, C. Simmerling, D. A. Case, J. Tuszynski, and A. Kovalenko, J. Chem. Theory Comput. 6, 607 (2010). 37T. Kloss and S. M. Kast, J. Chem. Phys. 128, 134505 (2008). 38J. S. Perkyns, G. C. Lynch, J. J. Howard, and B. M. Pettitt, J. Chem. Phys. 132, 064106 (2010). 39T. Darden, D. York, and L. Pedersen, J. Chem. Phys. 98, 10089 (1993). 7A.

J. Chem. Phys. 142, 114107 (2015) 40P.

P. Ewald, Ann. Phys. 369, 253 (1921). Kovalenko and F. Hirata, J. Chem. Phys. 112, 10391 (2000). 42J. D. Talman, J. Comput. Phys. 29, 35 (1978). 43G. B. Arfken and H. J. Weber, Mathematical Methods for Physicists, 6th ed. (Academic Press, London, 2005). 44S. J. Singer and D. Chandler, Mol. Phys. 55, 621 (1985). 45M. Heinen, E. Allahyarov, and H. Löwen, J. Comput. Chem. 35, 275 (2014). 46M. Deserno and C. Holm, J. Chem. Phys. 109, 7678 (1998). 47M. Griebel, S. Knapek, and G. Zumbusch, Numerical Simulation in Molecular Dynamics (Springer, Berlin, Heidelberg, 2007). 48R. W. Hockney and J. W. Eastwood, Computer Simulation Using Particles (IOP, Bristol, 1988). 49U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Pedersen, J. Chem. Phys. 103, 8577 (1995). 50H. G. Petersen, J. Chem. Phys. 103, 3668 (1995). 51H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys. Chem. 91, 6269 (1987). 52S. Maw, H. Sato, S. Ten-no, and F. Hirata, Chem. Phys. Lett. 276, 20 (1997). 53H. Sato and F. Hirata, J. Chem. Phys. 111, 8545 (1999). 54D. Beglov and B. Roux, J. Chem. Phys. 100, 9050 (1994). 55D. A. Case, T. A. Darden, T. E. Cheatham III, C. L. Simmerling, J. Wang, R. E. Duke, R. Luo, R. C. Walker, W. Zhang, K. M. Merz, B. Roberts, S. Hayik, A. Roitberg, G. Seabra, J. Swails, A. W. Goetz, I. Kolossváry, K. F. Wong, F. Paesani, J. Vanicek, R. M. Wolf, J. Liu, X. Wu, S. R. Brozell, T. Steinbrecher, H. Gohlke, Q. Cai, X. Ye, J. Wang, M.-J. Hsieh, G. Cui, D. R. Roe, D. H. Mathews, M. G. Seetin, R. Salomon-Ferrer, C. Sagui, V. Babin, T. Luchko, S. Gusarov, A. Kovalenko, and P. A. Kollman, AMBER 12, University of California, San Francisco, 2012. 56A. Onufriev, D. Bashford, and D. A. Case, J. Phys. Chem. B 104, 3712 (2000). 57A. Onufriev, D. Bashford, and D. A. Case, Proteins: Struct., Funct., Bioinf. 55, 383 (2004). 58V. Hornak, R. Abel, A. Okur, B. Strockbine, A. Roitberg, and C. Simmerling, Proteins: Struct., Funct., Bioinf. 65, 712 (2006). 59J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman, and D. A. Case, J. Comput. Chem. 25, 1157 (2004). 60A. Jakalian, B. L. Bush, D. B. Jack, and C. I. Bayly, J. Comput. Chem. 21, 132 (2000). 61A. Jakalian, D. B. Jack, and C. I. Bayly, J. Comput. Chem. 23, 1623 (2002). 62V. I. Lebedev and D. N. Laikov, Dokl. Math. 59, 477 (1999). 63See supplementary material at http://dx.doi.org/10.1063/1.4914321 for structures and parameters. 64C. Runge, Z. Math. Phys. 46, 224 (1901). 65M. Pippig and D. Potts, SIAM J. Sci. Comput. 35, C411 (2013). 66V. Ballenegger, J. J. Cerdà, and C. Holm, J. Chem. Theory Comput. 8, 936 (2012). 41A.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.30.242.61 On: Sat, 30 May 2015 16:09:51