Direct summation of dipole-dipole interactions using the Wolf formalism , Björn Stenqvist , Martin Trulsson, Alexei I. Abrikosov, and Mikael Lund
Citation: The Journal of Chemical Physics 143, 014109 (2015); doi: 10.1063/1.4923001 View online: http://dx.doi.org/10.1063/1.4923001 View Table of Contents: http://aip.scitation.org/toc/jcp/143/1 Published by the American Institute of Physics
THE JOURNAL OF CHEMICAL PHYSICS 143, 014109 (2015)
Direct summation of dipole-dipole interactions using the Wolf formalism Björn Stenqvist,1,a) Martin Trulsson,2 Alexei I. Abrikosov,3 and Mikael Lund1
1
Department of Theoretical Chemistry, Lund University, P.O. Box 124, SE-22100 Lund, Sweden Universite Paris-Sud, CNRS, LPTMS, UMR 8626, Orsay F-91405, France 3 Department of Physical Chemistry, Lund University, P.O. Box 124, SE-22100 Lund, Sweden 2
(Received 16 April 2015; accepted 11 June 2015; published online 2 July 2015) We present an expanded Wolf formalism for direct summation of long-range dipole-dipole interactions and rule-of-thumbs how to choose optimal spherical cutoff (Rc ) and damping parameter (α). This is done by comparing liquid radial distribution functions, dipole-dipole orientation correlations, particle energies, and dielectric constants, with Ewald sums and the Reaction field method. The resulting rule states that ασ < 1 and αRc > 3 for reduced densities around ρ∗ = 1 where σ is the particle size. Being a pair potential, the presented approach scales linearly with system size and is applicable to simulations involving point dipoles such as the Stockmayer fluid and polarizable water models. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4923001] I. INTRODUCTION
The computation of long-range electrostatic interactions in molecular simulations is often the time-limiting step and subject to current research.1–3 Ewald summation techniques (Ewald)4 and its derivatives, e.g., Particle Mesh Ewald5 (PME), are popular and with optimized parameters,6 the scaling with number of particles, N, is O(N 3/2) while for large systems, PME scales as O(N log N). Both methods are vastly faster than direct summation, O(N 2), but require special attention when implemented on parallel architectures.7 Besides that the Ewald method introduces numerical errors (as the sum, in practice, needs to be truncated), it also relies on periodic boundary condition (PBC), introducing intrinsic long-ranged order while avoiding surface effects.8 These artefacts can be avoided with, e.g., hyperspherical geometry9–12 which on the other hand introduces curvature effects and O(N 2) complexity. An alternative scheme for handling electrostatics is the Wolf method14 and derivatives thereof,13,15–17 based on truncation of electrostatics giving linear scaling. Other O(N) methods are the reaction field method (RF)18,19 and the Fast Multipole Method (FMM).20 The above mentioned methods, except the hypersphere, require some form of parametrization. Ewald and Wolf require a damping parameter, α, and a spherical cutoff distance, Rc , (Ewald further requires a second reciprocal cutoff, k c ) and RF needs a dielectric constant, ε RF, of the surrounding medium as well as a cutoff. The aim of this study is to provide a parametrization for the Wolf method for dipolar systems and to give rules of how to choose Rc and α parameters as only a few studies have focused on this task.15 Recent studies have explored different approaches for dipole-dipole interactions,21,22 required in, e.g., simulations with polarizable force fields, though parametrization was not the prime focus. Optimal parameters for Ewald methods have already been thoroughly studied23–25 and shall not be discussed. For the purpose of parametrization, we chose to study Stockmayer fluids at a)Electronic mail: [email protected]
0021-9606/2015/143(1)/014109/9/$30.00
different reduced dipole moments. We concentrate on radial distribution functions (RDF’s), dipole-dipole correlation functions, energies, and the collective behavior of the dipoles giving rise to the dielectric constant. The results are compared with the standard Ewald and the RF method for dipolar systems.
II. METHODS
The total interaction energy of a PBC system is given by N
ETot =
N
′
1 u(ri , r j , n), 2 i=1 j=1 n
(1)
where ri denotes the position of particle i, n is a vector, with elements nγ Lγ , from the original box to a distant replica, and u the pair-potential between two particles. Lγ is the simulation box-lengths in the γ-direction, nγ ∈ Z, and the prime indicates ¯ The sum involves infinitely many terms that i , j when n = 0. which is unmanageable for any computer. Some truncation is thus required which, for long-ranged interactions, clearly is a delicate issue. The dipole-dipole interaction energy between two point dipoles [ri , µ i ] and [r j , µ j ] is ui j = −µ i T2(ri j )µ j ,
(2)
where T2 denotes a dipolar interaction tensor and ri j = ri − r j . Different approaches for handling long-ranged electrostatic interactions are captured by this tensor and for an infinite, nonPBC space, ( ) 1 1 1 1 ∇∇ = (3ˆrrˆ T − I), (3) T2(r) = 4πε 0 |r| 4πε 0 |r|3 where rˆ = r/|r| and I is the identity tensor. A. Reaction field
The reaction field method, developed in the 1930s,18,19 uses the concept of dipoles in a cavity immersed in a surrounding dielectric continuum medium. The dipoles in the cavity
143, 014109-1
© 2015 AIP Publishing LLC
014109-2
Stenqvist et al.
J. Chem. Phys. 143, 014109 (2015)
creates a total dipole moment that polarizes the surrounding continuum which in turn affects the dipoles. With the advent of computer simulations, the method was included in calculations by Barker and Watts.26 The implementation is straightforward as the interaction tensor is only slightly altered compared to the infinite case as 1 2(ε RF − 1) I if |ri j | ≤ Rc . (4) TRF 2 (ri j ) = T2(ri j ) − 4πε 0 (2ε RF + 1)Rc3 For other ri j , the tensor is equal to the zero tensor. B. Ewald
The Ewald method is well-known4,23 and only a summary of the theory will be presented here. The total system energy, Eq. (1), is only conditionally convergent and thus cannot be evaluated directly. To manage this problem, the sum is split into two parts, a short-ranged interaction part handled in real space and a long-ranged part summed in reciprocal space. The splitting is performed by a unity multiplication with the error-function and its complement, erf(α|ri j |) + erfc(α|ri j |) = 1, where after the two parts converge rapidly, giving welldefined energies.23 The physical interpretation of the errorfunction is to introduce Gaussian moment distributions around each particle, but in practice, one can chose any splitting function. The final energy is given by a sum of terms, Ewald Ewald Ewald Ewald Ewald ETot = EReal + EReciprocal + ESelf + ESurface
(5)
Ewald which are given in the Appendix B. The term EReal runs over all particles in the simulation box and its images, but in practice, one chooses a spherical cutoff, Rc , smaller or equal to half the box-length (L/2). A similar approach is used Ewald Ewald for the EReciprocal term in reciprocal space. The term ESelf includes the interaction between the point particles and their corresponding surrounding Gaussian distributions and is thus Ewald independent of particle positions. The last term ESurface is sometimes introduced in Ewald and corresponds to a solvation energy in a surrounding dielectric medium.
C. Wolf
The first version of the Wolf method14 was based on a shift of the interaction potential at the cutoff, Rc , such that the potential is zero at and, per definition, outside of this cutoff sphere. A simple truncation of the potential generally gives fairly inaccurate results, something that Wolf and coworkers acknowledged to stem from non-electroneutrality within the spherical system. By adding oppositely signed charges at the cutoff, effectively shifting the potential (SP), results improve significantly. By employing the same logic for dipoles, i.e., adding dipoles so as to cancel the dipole moment, we get ( ) N 1 SP µ i −T2(ri j ) + lim T2(ri j ) µ j . ETot = ri j → rˆ i j R c 2 i=1 i, j |ri j |< R c
(6) The shift in Eq. (6) corresponds to introducing a dipole at the cutoff, linearly extrapolated along the ri j -vector and having
opposite direction of the dipole-vector. This approach thus shows some analogies with the RF method which too can be viewed as a fully or partially compensating dipole on the cutoff sphere. Wolf et al. found that simply cancelling the monopole is insufficient to get fast converging energies. To remedy this, a damping-function akin the real-space part of the Ewald summation, f (|r|, α) = erfc(|r|α), was proposed (DSP). The introduced error is given by the complementary erf-part and can be tuned to be arbitrarily small with the damping factor α. In the case of dipoles, the damped (D) interaction tensor becomes ) ( f (|r|, α) D T2 (r, α) = ∇∇ |r| 1 (7) = 3 (3ˆrrˆ TC(|r|, α) − IB(|r|, α)), |r| where B(|r|) and C(|r|) are given in the Appendix B by Eqs. (B2) and (B3), respectively. To make the Wolf potential compatible with molecular dynamics simulations, Fennell and Gezelter13 introduced damped force shifting (DSF), similar to Zahn.15 Extending this approach to dipoles, we get D D T2 (ri j , α) TDSF lim 2 (ri j , α) = T2 (ri j , α) − ri j → rˆ i j R c
− (|ri j | − Rc )
lim
ri j → rˆ i j R c
∂T D (r , α) 2 i j . ∂ri j
(8)
Wolf et al. also found that adding the self-term, i.e., an electrostatic cancellation of the central particles moment, achieves faster convergence when calculating Madelung energies. In the Appendix A, we have derived the self-energy for the dipolar version. The result is presented in Eq. (9), 2 2
N
1 erfc(αRc ) 2α e−α R c 4 α 3 + 2 DSF +√ + √ µ. ESelf =− * 2, 3 π - i=1 i π Rc2 Rc3
(9)
The present approach, essentially a dipolar version of the Fennel-potential, needs both a parametrization of α and Rc . These are exactly what we aim to investigate. The presented potential has been independently derived with regard to previously presented work21 though the self-energy term differs slightly.
III. DIELECTRIC CONSTANT
For simulations using the reaction field, the general equation for the relative dielectric constant, ε r , is27 −1 εr − 1 εr − 1 ˜ 1 ⟨M 2⟩ 1− T(0) = , (10) εr + 2 εr + 2 3ε 0 3V k BT ˜ where T(0) = 2(ε RF − 1)/(2ε RF + 1), ⟨M 2⟩ = ⟨| i µ i |2⟩ is the mean squared dipole moment of the simulation box, k B is Boltzmann’s constant, T the temperature, ε 0 is the vacuum permittivity, ε RF is the dielectric constant of the surrounding medium, and V is the volume. Thus, for conducting, insulating, and vacuum boundaries,
014109-3
Stenqvist et al.
εr − 1 3 (2ε r + 1)(ε r − 1) 1 ⟨M 2⟩ = 9ε r 3ε 0 3V k BT εr − 1 εr + 2
J. Chem. Phys. 143, 014109 (2015)
if ε RF = ∞, if ε RF = ε r ,
(11)
if ε RF = 1.
For Ewald with conducting boundaries,28 the expression for the dielectric constant is identical to that of the reaction field with ε RF = ∞. For the present dipolar case, we use the approach of Neumann and Steinhauser29–31 but with the interaction tensor introduced by Eq. (8). We thereby arrive at a similar expression as Zahn et al.,15 Eq. (10) with 2[(αRc )4 + 2(αRc )2 + 3]αRc e−α ˜ T(0) = erf(αRc ) − √ 3 π
2 R2 c
.
(12)
˜ Note that T(0) ≈ 1 in the limit αRc ≫ 1 whereby the conducting Ewald and RF results are retrieved. IV. COMPUTER SIMULATIONS
Metropolis Monte Carlo (MC) simulations were used to generate configurations of the Stockmayer fluid32 in the canonical (NVT) ensemble and were performed using the Faunus software.33 In reduced units, the density was ρ∗ = N σ 3/V = 0.924, temperature T ∗ = k BT/ϵ LJ = 1.333, µ∗2 = µ2/(4πϵ 0σ 3ϵ LJ ) = {1.663, 3.470}, and N = 3000 particles. When we have studied different densities, ρ∗ = {0.231, 0.462, 0.693, 0.924, 1.155, 1.386}, the system size has been constant, L = 10.265σ. The different, spherical cutoffs, Rc , were chosen to be smaller than or equal to half the box length. Simulations were equilibrated in 104 sweeps and ten times longer productions runs using translational and rotational MC moves. In each sweep, N moves were attempted. We used Stockmayer potentials consisting of the electrostatic potential given by the different interaction tensors defined above, combined with a Lennard-Jones (LJ) potential. All Ewald simulations were performed with conducting boundary conditions, a real-space cutoff of half the box length, α = 0.404σ −1 and 686 wavevectors. V. RESULTS
Unless otherwise stated, the reduced, squared dipole moment is µ∗2 = 3.470 but the main conclusions reached are valid also for µ∗2 = 1.663. To simplify the analysis, we address results as acceptable if they are within approximately 1–2 standard deviations (std’s) from the corresponding Ewald result. The presented electrostatic energies are including the self-energy.
FIG. 1. RDF’s (a) and dipole-dipole correlations (b) using different α −1, see legend. The black lines in the main figures correspond to Ewald simulations, often overlapping with the Wolf results. The insets show the corresponding errors; grey area is confined by ±1 Ewald std.
1. Effect of the damping
We initially compare five simulations with different α-values for a constant cutoff, Rc ≃ 7.4σ (L/2), testing the role of the smooth damping. Figure 1(a) shows RDF’s as well as their errors (i.e., the difference, see inset figure) between Ewald and Wolf RDF’s. As can be seen, g(r) is seemingly insensitive to changes of the damping parameter and dominated chiefly by the Lennard-Jones part of the pair interaction. Still, according to the error plot, α −1 = 0.9σ do not satisfy a limit about ±1–2 Ewald std and show oscillating tendencies, even if the relative error is negligible. While oscillating, the remaining α values fall within the error range of the Ewald results. Fig. 1(b) shows that for the dipole-dipole correlation, acceptable α −1 values are also obtained between 1.2σ and 2.1σ. Table I shows the average energy per particle converging to that of Ewald for all α −1 but 0.9σ. The same trend is seen for the dielectric constants.
A. Structural properties
In this section, we evaluate RDF’s of the Stockmayer fluid, and for a more critical measure of the anisotropic electrostatics, ˆ · also dipole-dipole correlation functions, defined as ⟨ µ(0) ˆ ˆ µ(r)⟩ where µ(0) is a central, normalized dipole surrounded ˆ by dipoles µ(r) in a spherical volume element at distance r + dr. The angular brackets denote an ensemble average.
2. Effect of the cutoff
We now study the opposite case: constant α and varying cutoffs. Results for α −1 = 1.5σ and α −1 = 1.8σ, are presented in Figs. 2 and 3 and in Table II. For α −1 = 1.5σ, it is clear that Rc = 5σ is a saturated cutoff, i.e., the results will not improve by increasing it further. In the α −1 = 1.8σ case, half the
014109-4
Stenqvist et al.
J. Chem. Phys. 143, 014109 (2015)
TABLE I. Dipolar, Lennard-Jones, and total particle energies as well as dielectric constants. The std is 0.002 for βu µ µ and βu LJ where β = 1/k BT . α −1
βu µ µ
βu LJ
βu Tot
εr
0.9σ 1.2σ 1.5σ 1.8σ 2.1σ Ewald
−4.552 −4.523 −4.520 −4.518 −4.517 −4.523
−4.000 −3.996 −3.998 −3.996 −3.996 −3.995
−8.552 −8.520 −8.518 −8.514 −8.513 −8.518
100 ± 30 120 ± 20 120 ± 40 130 ± 50 140 ± 40 140 ± 40
box-length is needed to get accurate results. This cutoff may be a saturated cutoff—more on this Sec. V B. B. Energies
Fig. 4(a) shows that the LJ energy is nearly constant when α −1 ≥ σ for every cutoff, except Rc = 2σ. The electrostatic energy, see Fig. 4(b), has a stationary or quasi-stationary point (indicated by dots) in the interval σ ≤ α −1 ≤ 2.3σ for every cutoff but Rc = 2σ. By using the least-square method, the fitted lines corresponding to the quasi-stationary points are described by Rc ≈ 3.2α −1 + 0.1σ for µ∗2 = 3.470 and Rc ≈ 3.1α −1 − 0.3σ for µ∗2 = 1.663. Fig. 4(b) shows that the larger the cutoff, the wider the width of the “stationary”
FIG. 3. RDF’s (a) and dipole-dipole correlations (b) with α −1 = 1.8σ using different cutoff’s, cf. legend; black lines the Ewald results, often overlapping with the Wolf results. The insets show the corresponding errors; black indicates ±1 Ewald std.
segment around the quasi-stationary point. Furthermore, in this section, we see that for α −1 = 1.8σ, a mere cutoff of 6σ is needed to get the correct energy, thus in this aspect Rc = 6σ is a saturated cutoff. TABLE II. Dipolar, Lennard-Jones, and total energies per particle as well as dielectric constants for α −1 = 1.5σ and α −1 = 1.8σ. The std is ∼ 0.002 for βu µ µ and βu L J .
FIG. 2. RDF’s (a) and dipole-dipole correlations (b) with α −1 = 1.5σ using different cutoff’s, cf. legend; black lines are Ewald results, often overlapping with the Wolf results. The insets show the corresponding errors; grey area is confined by ±1 Ewald std.
Rc
βu µ µ
2σ 3σ 4σ 5σ 6σ L/2 Ewald
−3.021 −4.138 −4.457 −4.518 −4.520 −4.520 −4.523
2σ 3σ 4σ 5σ 6σ L/2 Ewald
−2.876 −3.964 −4.349 −4.485 −4.514 −4.518 −4.523
βu LJ α −1 = 1.5σ −4.039 −3.991 −3.993 −3.995 −3.995 −3.998 −3.995 α −1 = 1.8σ −4.043 −3.995 −3.987 −3.995 −3.998 −3.996 −3.995
βu Tot
εr
−7.060 −8.130 −8.450 −8.513 −8.515 −8.518 −8.513
10 ± 10 10 ± 10 40 ± 10 130 ± 40 120 ± 30 120 ± 30 140 ± 40
−6.919 −7.959 −8.337 −8.479 −8.512 −8.514 −8.513
10 ± 10 10 ± 10 10 ± 10 60 ± 10 110 ± 40 130 ± 50 140 ± 40
014109-5
Stenqvist et al.
J. Chem. Phys. 143, 014109 (2015)
FIG. 4. βu LJ (a) and βu µ µ (b) at different cutoff values. The insets show the corresponding energy as a function of the cutoff for α −1 = 2.1σ and the Ewald energy. The energy plateau in (a) starts at α −1 = σ for every cutoff and ends with the equally colored star.
C. Dielectric constant
In this section, we investigate how damped electrostatic interactions influence the dielectric constant which is a measure of the collective polarization of the liquid. Throughout the analysis, the general behavior as a function of Rc and α −1 is preserved, see Figs. 5(a) and 6(a). Some information can readily be extracted from Fig. 5(a): The dielectric constant is a macroscopic property and to capture this behavior, a limit on the damping and minimal cutoff are required, α −1 ≥ σ and Rc ≥ 3σ. At the point Rc = 3σ and α −1 = σ, a triangle starts to form, inside of which acceptable results emerge. There are two interesting lines: α −1 = σ and Rc ≈ 3α −1 + 0.7σ. The area between these are statistically equal (around ±1 std) to the Ewald result. A large interval easy to target is α −1 = σ − 2.2σ when Rc is half the box length. From the analysis of the RDF’s and dipole-correlations, this is too a satisfactory interval. An extended cutoff increases the triangle and widens the possible margin of error in the choice of α such that acceptable results are still obtained. This is in close agreement with our earlier conclusion coming from the analysis of the energies. From this (limited) parameter space, we deduce that with increasing α −1 must come an even faster increase in cutoff, satisfying Rc ≥ 3α −1 + 0.7σ. Also, the dots in Fig. 4(b) can be seen in corresponding color in Fig. 5(a).
∆ε (a), and ε std in % (b) for µ ∗2 = 3.470—values outside the FIG. 5. ε Ewald Ewald color-bars are shown as white. The solid, red lines represent α −1 = σ and −1 R c = 3α + 0.7σ and form a triangular area within which valid results (±1 std) are obtained. The colored dots correspond to (quasi-) stationary points from Fig. 4(b). The dielectric constant for Ewald is ε Ewald = 142.1 ± 36.1.
The arguments established so far hold also for the less coupled µ∗2 = 1.663 case, see Fig. 6(a), where the lines are α −1 = σ and Rc ≈ 2.7α −1 + 0.4σ. The main difference is an effectively larger triangle which is reasonable since µ∗2 → 0 should trivially result in the entire parameter space being valid. D. Densities
So far we have studied a single density, ρ∗ = 0.924 and here we explore more dilute as well as denser systems for constant L = 29.629 Å. A cutoff of half the box-length was used in all simulations in this section. The results, condensed in Fig. 7, show that the dipole-dipole correlations become less accurate for the lowest and highest densities while the RDF only suffers from low densities. This may be due to an appearance of non-isotropic or crystalline phases; networks, ring, and chain formations in the dilute regime, and ferro-electric liquids/solids in the dense limit.34 The long-range order in these phases is suppressed when electrostatics are truncated as in the Wolf formalism which has already been pointed out.14
014109-6
Stenqvist et al.
J. Chem. Phys. 143, 014109 (2015)
FIG. 7. RMS of error between Wolf and Ewald for RDF’s (a) and dipoledipole correlations (b) using different reduced densities, cf. the legends. ∆ε (a), and ε std in % (b) for µ ∗2 = 1.663—values outside the FIG. 6. ε Ewald Ewald color-bars are shown as white. The solid, red lines represent α −1 = σ and R c = 2.7α −1 + 0.4σ and form a triangular area within which valid results (±1 std) are obtained. The colored dots correspond to (quasi-) stationary points in the electrostatic and self-energy. The dielectric constant for Ewald is ε Ewald = 16.4 ± 2.3.
In between these extreme densities, correct results emerge. The triangular area—compare, e.g., Figs. 5(a) and 6(a)—in parameter space giving correct results tends to narrow as the dielectric constant increases. At constant dipolar coupling, increasing the density leads to a larger dielectric constant which may explain the smaller window (in α −1) of accuracy in Fig. 7(b) for high densities; vice versa for low densities. The low-density precision windows are however here most likely limited by the small cutoff (determined by the box-length). Although the original Wolf method for charges is expected to work best for liquids, reasonable Madelung energies are obtained for a number of crystal lattices14 which is encouraging for treating dipolar crystals even if it is not guarantied to provide accurate results due to the anisotropic nature of the potential. E. Reaction field
The reaction field energies are shown in Table III. For ε RF = 1, energies are generally lower than for Ewald. The
dipole-dipole energies under these conditions are non-monotonic, narrowing the gap to the Ewald dipole-dipole energy with increasing cutoff, as expected, until Rc = 4σ, where after and the gap is widened again. This behavior makes the TABLE III. Reaction field dipole-dipole and LJ energies in units of k BT for different cutoff and dielectric boundaries. The Ewald dipole-dipole and total energies are −4.523 k BT and −3.995 k BT , respectively. The std is ∼ 0.002 for both βu µ µ and βu L J . Rc
ε RF = 1
2σ 3σ 4σ 5σ 6σ 7σ L/2
−4.710 −4.658 −4.644 −4.666 −4.742 −4.816 −4.839
2σ 3σ 4σ 5σ 6σ 7σ L/2
ε RF = ∞
Dipole-dipole energy −4.703 −4.520 −4.555 −4.514 −4.512 −4.512 −4.513 LJ energy −3.910 −3.929 −3.997 −3.998 −3.982 −3.995 −4.003 −3.996 −4.010 −3.997 −4.020 −3.997 −4.024 −3.997
ε RF = ε r
−4.695 −4.513 −4.547 −4.509 −4.512 −4.512 −4.512 −3.928 −3.997 −3.994 −3.997 −3.995 −3.996 −3.995
014109-7
Stenqvist et al.
J. Chem. Phys. 143, 014109 (2015)
TABLE IV. Dielectric constants ± their std’s from reaction field simulations. The corresponding Ewald dielectric constant is 140 ± 40. Rc
ε RF = 1
ε RF = ∞
ε RF = ε r
2σ 3σ 4σ 5σ 6σ 7σ L/2
−50 ± 240 90 ± 110 70 ± 100 50 ± 30 100 ± 90 40 ± 30 50 ± 30
390 ± 160 150 ± 30 200 ± 80 150 ± 60 110 ± 30 110 ± 30 130 ± 40
420 ± 220 140 ± 30 220 ± 50 120 ± 40 120 ± 20 130 ± 30 120 ± 30
dipole-dipole energy (and thus the total energy) seem to diverge with increasing cutoff. Since the reaction field method is identical to both the Wolf and Ewald methods in the limit Rc → ∞, we deduce that our parameter space for the cutoff is too small to capture the convergence. The inaccuracy when ε RF = 1 is also seen in Table IV as proportionally large deviations from Ewald for any tested cutoff. When ε RF = ε r or ε RF = ∞, energies do converge with increasing cutoff, while the dielectric constants are indistinguishable within one std. This is true for any cutoff larger than 5σ after which the std’s remain unaffected by further increase in Rc . The contrast
to the vacuum case can be rationalized by the fact that the more polarizable surroundings tend to dampen electrostatic correlations, whereby convergence is reached even for small system sizes. Figure 8(a) shows that, as previously stated, ε RF = 1 gives deviating results compared to Ewald. Surprisingly, the smaller cutoff gives noticeable more accurate results than for half the box length. Another abnormal result is a bias on the L/2 ε RF = ∞ case. The remaining cases are all acceptable. When including the more sensitive dipole-dipole correlations in Fig. 8(b), we see that there is no evidence of any of the cases to be accurate. An increasing cutoff seems to damp the fluctuations, although the increasing trend is unaffected.
VI. CONCLUSIONS
We have presented an expanded Wolf formalism for long range dipolar interactions, evaluated on Stockmayer fluids where structural and dielectric properties are investigated. Comparisons are made with Ewald summation and a reaction field method. In short, the α-parameter in the Wolf method should respect α −1 ≥ σ and αRc > 3 to get values in agreement with Ewald results. The latter limit ensures sufficiently soft damping before the hard cutoff is introduced. This points to the fact that even if both the potential and force are shifted, a hard cutoff alters the system fluctuations, while soft damping does not suffer from similar effects. Combining the two rules leads to a cutoff criteria of Rc > 3σ. Increasing Rc beyond αRc ∼ 3 at fixed α gives little or no increase in accuracy—at least for densities around ρ∗ = 1. Finding the (quasi-) stationary electrostatic energy with regard to α seems to give parameters with accurate results. If no such stationary point can be found, or its width is small, it is an indication that the chosen cutoff is too short. Previous work22 on multipolar fluids (dipoles, and quadrupoles) in the micro-canonical ensemble is qualitatively consistent with the damping parameter interval deduced in this work. Finally, Wolf with optimal parameters presents more or equally accurate results as the equivalently scaling RF method used with either conducting, insulating, or vacuum boundaries. ACKNOWLEDGMENTS
We thank LUNARC for computational resources and for financial support: Södra Research foundation; the Swedish Research Council; and the Swedish Foundation for Strategic Research. Finally, we thank Per Linse for useful discussions and comments.
APPENDIX A: SELF-ENERGY
FIG. 8. Reaction field RDF’s (a) and dipole-dipole correlations (b) using different ε RF-values (labeled); the black line is from Ewald. The cutoffs are 5σ for the three top legends and half the box length (L/2) for the three bottom legends. Note that most RDF’s and correlations are overlapping. The insets show the corresponding errors; grey area is confined by ±1 Ewald std.
In this section, we will use a direct expansion of the point charge approach in Wolf14 to dipoles, in order to get the expanded DSF self-energy. Since we in this section are referring to similar equations in Wolf,14 we will, when possible, use similar notation. The damped version of Wolf uses a unity multiplication, erfc(α|ri j + n|) + erf(α|ri j + n|) = 1, of
014109-8
Stenqvist et al.
J. Chem. Phys. 143, 014109 (2015)
the total energy of a periodic boundary condition system (see Wolf14 Eq. (5.1)). For dipoles, we get
where (see Wolf14 Eq. (5.12))
) ( N ∞ 1 T T erfc(α|ri j + n|) µj ETot = − µ ∇∇ 2 i=1 j,i=1 i |ri j + n|
DSP ESelf =−
) ( N ∞ 1 T T erf(α|ri j + n|) − µ j. µ ∇∇ 2 i=1 j,i=1 i |ri j + n|
(A1)
The second part of Eq. (A1) is large for large values of α. By extracting the diverging term, we end up with (see Wolf14 Eqs. (5.2)–(5.5)) N ∞ 1 T D µ T (ri j + n, α)µ j 2 i=1 j,i=1 i 2 ( ) N N 1 T T erf(α|ri j + n|) µ i ∇∇ µj − 2 i=1 j=1 |ri j + n|
ETot = −
N 2α 3 2 − √ µi , 3 π i=1
(A2)
where TD 2 (r, α) is described by Eq. (7). By making this extraction, the remainder (i.e., the second part of Eq. (A2)) is now a small correction. Hence it is neglected, giving (see Wolf14 Eq. (5.6)) Approx
ETot
=−
1 2
N i=1
µTi T2D (ri j + n, α)µ j −
j,i |r i j |< R c
N 3
2α √ 3 π
µ2i .
i=1
(A3) By using the condition about cancellation of dipole moment inside some cutoff (smaller than half the box-length), the expression becomes (see Wolf14 Eq. (5.7)) DSP ETot =−
i=1
N
(A8)
All described so far is a strict expansion of the approach for charges in Wolf14 for dipoles. The only assumption is that the compensation for the centred dipole is in the form of a smeared out dipole on the cutoff sphere. However, one can interpret Wolf14 as if the compensating charge is in the form of a smeared out charge on the surface, i.e., just like we have presented but for dipoles. Hence, the presented approach is backward compatible with Wolf.14
Ewald EReal
N N ′ ( µi · µ j 1 = B(|ri j + n|) 2 i=1 j=1 |ri j + n|3 3
) ((ri j + n) · µ i )((ri j + n) · µ j ) − 3C(|ri j + n|) , |ri j + n|5 (B1) where (A5)
which can be written as (see Wolf14 Eq. (5.10)) 2
2 2
1 erfc(αRc ) 2α e−α R c 4 α 3 + 2 +√ + √ µ. =− * 2, 3 π - i=1 i π Rc2 Rc3
n∈Z
*. +/ N 1 lim ..− µTi T2D (ri j , α)µ j // |ri j |→ R c . 2 / j,i i=1 |ri j |< R c , N 1 µTi T2D (r, α) r= R c r µ i + , + lim *− |ri i |→ R c |r i i | i i , 2 i=1 -
DSP ETot =−
DSP DSF ESelf = ESelf
(A4)
where (see Wolf14 Eq. (5.9))
N 1
First of all, we acknowledge that the term rii is ill-defined, the direction seems arbitrary. To remedy this we, instead of a point dipole use a smeared out dipole on the cutoff sphere. Since it is the cancellation of the zeroth order electrostatic moment (i.e., charges) that is the principle of Wolf14 which we extrapolate to first order moment (i.e., dipoles), this is still within the framework of the expanded formalism. The dipole moment density is chosen such that ρ(θ) = ρ0 cos(θ) where ρ(θ) is the dipole density in spherical coordinates (θ ∈ [0, π]), using a ρ0 such that the integrated dipole moment density on the sphere is oppositely signed but equal in magnitude to the dipole moment of the centred dipole. The DSP self-energy one then gets is actually the same as the DSF selfenergy,
Here, we show the equations used in the Ewald summation method. The real-space (short ranged) part shows similarities with the Wolf method,
|ri j |< R c
Cancel ETot =
(A7)
APPENDIX B: EWALD EQUATIONS
N 1 µTi T2D (ri j , α)µ j 2 i=1 i, j
N 2α 3 2 Cancel µi − ETot , − √ 3 π i=1
N N 1 2α 3 2 µi . µ i T2D (r, α) r= R c r µ i + √ |ri i | i i 2 i=1 3 π i=1
µTi T2D (ri j , α)
i, j |ri j |< R c
) |r i j | i j
DSP , µ j − ESelf
(B2)
( ) 2αr 3 3 2 2 C(r, α) = erfc(αr) + √ 2α 2 + 2 e−α r . r 3 π
(B3)
and
(
− T2D (r, α) r= R c r
2αr 2 2 B(r, α) = erfc(αr) + √ e−α r π
(A6)
The reciprocal-space part takes care of the long-ranged interactions,
014109-9
Stenqvist et al.
Ewald EReciprocal
J. Chem. Phys. 143, 014109 (2015)
( 2 2) 1 4π π |k| = exp − 2 2 3 2 2L k,0 |k| α L
7D.
k∈Z3
) ( N N 2πik · rmn . × (µ m · k)(µ n · k) exp L m=1 n=1 (B4) Finally, Eqs. (B5) and (B6) give the self-energy and the surface energy, arising from an interaction between the enclosed system and a surrounding dielectric medium with dielectric constant ε Surface, N 2α 3 Ewald ESelf =− √ |µ i |2 3 π i=1
(B5)
and Ewald ESurface =
2π (2ε Surface + 1)L
N N (µ i · µ j ), 3
(B6)
i=1 j=1
where we used L = L x = L y = L z . 1B.
W. McCann and O. Acevedo, J. Chem. Theory Comput. 9, 944 (2013). Levitt, M. Hirshberg, R. Sharon, and V. Daggett, Comput. Phys. Commun. 91, 215 (1995). 3B. Linse and P. Linse, J. Chem. Phys. 141, 184114 (2014). 4P. Ewald, Ann. Phys. 369, 253287 (1921). 5T. Darden, D. York, and L. Pedersen, J. Chem. Phys. 98, 10089 (1993). 6J. W. Perram, H. G. Petersen, and S. W. De Leeuw, Mol. Phys. 65, 875 (1988). 2M.
E. Shaw, M. M. Deneroff, R. O. Dror, J. S. Kuskin, R. H. Larson, J. K. Salmon, C. Young, B. Batson, K. J. Bowers, and J. C. Chao, ACM SIGARCH Comput. Archit. News 35, 1 (2007). 8J. Stenhammar, P. Linse, and G. Karlström, J. Phys.: Condens. Matter 20, 494204 (2008). 9J. Hansen, D. Levesque, and J. Weis, Phys. Rev. Lett. 43, 979 (1979). 10J. Caillol and D. Levesque, J. Chem. Phys. 94, 597 (1991). 11J. Caillol, J. Chem. Phys. 96, 1455 (1992). 12M. Trulsson, J. Chem. Phys. 133, 174105 (2010). 13C. J. Fennell and J. D. Gezelter, J. Chem. Phys. 124, 234104 (2006). 14D. Wolf, P. Keblinski, S. Phillpot, and J. Eggebrecht, J. Chem. Phys. 110, 8254 (1999). 15D. Zahn, B. Schilling, and S. M. Kast, J. Phys. Chem. B 106, 10725 (2002). 16Y. Yonezawa, J. Chem. Phys. 136, 244103 (2012). 17G. S. Fanourgakis, J. Phys. Chem. B 119, 1974 (2015). 18A. Martin, London, Edinburgh, Dublin Philos. Mag. J. Sci. 8, 547 (1929). 19R. P. Bell, Trans. Faraday Soc. 27, 797 (1931). 20L. Greengard and V. Rokhlin, J. Comput. Phys. 73, 325 (1987). 21M. Lamichhane, J. D. Gezelter, and K. E. Newman, J. Chem. Phys. 141, 134109 (2014). 22M. Lamichhane, K. E. Newman, and J. D. Gezelter, J. Chem. Phys. 141, 134110 (2014). 23S. de Leeuw, J. Perram, and E. Smith, Proc. R. Soc. A 373, 27 (1980). 24Z. Wang and C. Holm, J. Chem. Phys. 115, 6351 (2001). 25J. Kolafa and J. W. Perram, Mol. Simul. 9, 351 (1992). 26J. Barker and R. Watts, Mol. Phys. 26, 789 (1973). 27M. Neumann, Mol. Phys. 50, 841 (1983). 28E. Pollock and B. Alder, Phys. A 102, 1 (1980). 29M. Neumann and O. Steinhauser, Chem. Phys. Lett. 95, 417 (1983). 30M. Neumann, O. Steinhauser, and G. S. Pawley, Mol. Phys. 52, 97 (1984). 31M. Neumann, Mol. Phys. 57, 97 (1986). 32W. H. Stockmayer, J. Chem. Phys. 9, 398 (1941). 33B. Stenqvist, A. Thuresson, A. Kurut, R. Vácha, and M. Lund, Mol. Simul. 39, 1233 (2013). 34C. Holm and J.-J. Weis, Curr. Opin. Colloid Interface Sci. 10, 133 (2005).
Citation: The Journal of Chemical Physics 143, 014109 (2015); doi: 10.1063/1.4923001 View online: http://dx.doi.org/10.1063/1.4923001 View Table of Contents: http://aip.scitation.org/toc/jcp/143/1 Published by the American Institute of Physics
THE JOURNAL OF CHEMICAL PHYSICS 143, 014109 (2015)
Direct summation of dipole-dipole interactions using the Wolf formalism Björn Stenqvist,1,a) Martin Trulsson,2 Alexei I. Abrikosov,3 and Mikael Lund1
1
Department of Theoretical Chemistry, Lund University, P.O. Box 124, SE-22100 Lund, Sweden Universite Paris-Sud, CNRS, LPTMS, UMR 8626, Orsay F-91405, France 3 Department of Physical Chemistry, Lund University, P.O. Box 124, SE-22100 Lund, Sweden 2
(Received 16 April 2015; accepted 11 June 2015; published online 2 July 2015) We present an expanded Wolf formalism for direct summation of long-range dipole-dipole interactions and rule-of-thumbs how to choose optimal spherical cutoff (Rc ) and damping parameter (α). This is done by comparing liquid radial distribution functions, dipole-dipole orientation correlations, particle energies, and dielectric constants, with Ewald sums and the Reaction field method. The resulting rule states that ασ < 1 and αRc > 3 for reduced densities around ρ∗ = 1 where σ is the particle size. Being a pair potential, the presented approach scales linearly with system size and is applicable to simulations involving point dipoles such as the Stockmayer fluid and polarizable water models. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4923001] I. INTRODUCTION
The computation of long-range electrostatic interactions in molecular simulations is often the time-limiting step and subject to current research.1–3 Ewald summation techniques (Ewald)4 and its derivatives, e.g., Particle Mesh Ewald5 (PME), are popular and with optimized parameters,6 the scaling with number of particles, N, is O(N 3/2) while for large systems, PME scales as O(N log N). Both methods are vastly faster than direct summation, O(N 2), but require special attention when implemented on parallel architectures.7 Besides that the Ewald method introduces numerical errors (as the sum, in practice, needs to be truncated), it also relies on periodic boundary condition (PBC), introducing intrinsic long-ranged order while avoiding surface effects.8 These artefacts can be avoided with, e.g., hyperspherical geometry9–12 which on the other hand introduces curvature effects and O(N 2) complexity. An alternative scheme for handling electrostatics is the Wolf method14 and derivatives thereof,13,15–17 based on truncation of electrostatics giving linear scaling. Other O(N) methods are the reaction field method (RF)18,19 and the Fast Multipole Method (FMM).20 The above mentioned methods, except the hypersphere, require some form of parametrization. Ewald and Wolf require a damping parameter, α, and a spherical cutoff distance, Rc , (Ewald further requires a second reciprocal cutoff, k c ) and RF needs a dielectric constant, ε RF, of the surrounding medium as well as a cutoff. The aim of this study is to provide a parametrization for the Wolf method for dipolar systems and to give rules of how to choose Rc and α parameters as only a few studies have focused on this task.15 Recent studies have explored different approaches for dipole-dipole interactions,21,22 required in, e.g., simulations with polarizable force fields, though parametrization was not the prime focus. Optimal parameters for Ewald methods have already been thoroughly studied23–25 and shall not be discussed. For the purpose of parametrization, we chose to study Stockmayer fluids at a)Electronic mail: [email protected]
0021-9606/2015/143(1)/014109/9/$30.00
different reduced dipole moments. We concentrate on radial distribution functions (RDF’s), dipole-dipole correlation functions, energies, and the collective behavior of the dipoles giving rise to the dielectric constant. The results are compared with the standard Ewald and the RF method for dipolar systems.
II. METHODS
The total interaction energy of a PBC system is given by N
ETot =
N
′
1 u(ri , r j , n), 2 i=1 j=1 n
(1)
where ri denotes the position of particle i, n is a vector, with elements nγ Lγ , from the original box to a distant replica, and u the pair-potential between two particles. Lγ is the simulation box-lengths in the γ-direction, nγ ∈ Z, and the prime indicates ¯ The sum involves infinitely many terms that i , j when n = 0. which is unmanageable for any computer. Some truncation is thus required which, for long-ranged interactions, clearly is a delicate issue. The dipole-dipole interaction energy between two point dipoles [ri , µ i ] and [r j , µ j ] is ui j = −µ i T2(ri j )µ j ,
(2)
where T2 denotes a dipolar interaction tensor and ri j = ri − r j . Different approaches for handling long-ranged electrostatic interactions are captured by this tensor and for an infinite, nonPBC space, ( ) 1 1 1 1 ∇∇ = (3ˆrrˆ T − I), (3) T2(r) = 4πε 0 |r| 4πε 0 |r|3 where rˆ = r/|r| and I is the identity tensor. A. Reaction field
The reaction field method, developed in the 1930s,18,19 uses the concept of dipoles in a cavity immersed in a surrounding dielectric continuum medium. The dipoles in the cavity
143, 014109-1
© 2015 AIP Publishing LLC
014109-2
Stenqvist et al.
J. Chem. Phys. 143, 014109 (2015)
creates a total dipole moment that polarizes the surrounding continuum which in turn affects the dipoles. With the advent of computer simulations, the method was included in calculations by Barker and Watts.26 The implementation is straightforward as the interaction tensor is only slightly altered compared to the infinite case as 1 2(ε RF − 1) I if |ri j | ≤ Rc . (4) TRF 2 (ri j ) = T2(ri j ) − 4πε 0 (2ε RF + 1)Rc3 For other ri j , the tensor is equal to the zero tensor. B. Ewald
The Ewald method is well-known4,23 and only a summary of the theory will be presented here. The total system energy, Eq. (1), is only conditionally convergent and thus cannot be evaluated directly. To manage this problem, the sum is split into two parts, a short-ranged interaction part handled in real space and a long-ranged part summed in reciprocal space. The splitting is performed by a unity multiplication with the error-function and its complement, erf(α|ri j |) + erfc(α|ri j |) = 1, where after the two parts converge rapidly, giving welldefined energies.23 The physical interpretation of the errorfunction is to introduce Gaussian moment distributions around each particle, but in practice, one can chose any splitting function. The final energy is given by a sum of terms, Ewald Ewald Ewald Ewald Ewald ETot = EReal + EReciprocal + ESelf + ESurface
(5)
Ewald which are given in the Appendix B. The term EReal runs over all particles in the simulation box and its images, but in practice, one chooses a spherical cutoff, Rc , smaller or equal to half the box-length (L/2). A similar approach is used Ewald Ewald for the EReciprocal term in reciprocal space. The term ESelf includes the interaction between the point particles and their corresponding surrounding Gaussian distributions and is thus Ewald independent of particle positions. The last term ESurface is sometimes introduced in Ewald and corresponds to a solvation energy in a surrounding dielectric medium.
C. Wolf
The first version of the Wolf method14 was based on a shift of the interaction potential at the cutoff, Rc , such that the potential is zero at and, per definition, outside of this cutoff sphere. A simple truncation of the potential generally gives fairly inaccurate results, something that Wolf and coworkers acknowledged to stem from non-electroneutrality within the spherical system. By adding oppositely signed charges at the cutoff, effectively shifting the potential (SP), results improve significantly. By employing the same logic for dipoles, i.e., adding dipoles so as to cancel the dipole moment, we get ( ) N 1 SP µ i −T2(ri j ) + lim T2(ri j ) µ j . ETot = ri j → rˆ i j R c 2 i=1 i, j |ri j |< R c
(6) The shift in Eq. (6) corresponds to introducing a dipole at the cutoff, linearly extrapolated along the ri j -vector and having
opposite direction of the dipole-vector. This approach thus shows some analogies with the RF method which too can be viewed as a fully or partially compensating dipole on the cutoff sphere. Wolf et al. found that simply cancelling the monopole is insufficient to get fast converging energies. To remedy this, a damping-function akin the real-space part of the Ewald summation, f (|r|, α) = erfc(|r|α), was proposed (DSP). The introduced error is given by the complementary erf-part and can be tuned to be arbitrarily small with the damping factor α. In the case of dipoles, the damped (D) interaction tensor becomes ) ( f (|r|, α) D T2 (r, α) = ∇∇ |r| 1 (7) = 3 (3ˆrrˆ TC(|r|, α) − IB(|r|, α)), |r| where B(|r|) and C(|r|) are given in the Appendix B by Eqs. (B2) and (B3), respectively. To make the Wolf potential compatible with molecular dynamics simulations, Fennell and Gezelter13 introduced damped force shifting (DSF), similar to Zahn.15 Extending this approach to dipoles, we get D D T2 (ri j , α) TDSF lim 2 (ri j , α) = T2 (ri j , α) − ri j → rˆ i j R c
− (|ri j | − Rc )
lim
ri j → rˆ i j R c
∂T D (r , α) 2 i j . ∂ri j
(8)
Wolf et al. also found that adding the self-term, i.e., an electrostatic cancellation of the central particles moment, achieves faster convergence when calculating Madelung energies. In the Appendix A, we have derived the self-energy for the dipolar version. The result is presented in Eq. (9), 2 2
N
1 erfc(αRc ) 2α e−α R c 4 α 3 + 2 DSF +√ + √ µ. ESelf =− * 2, 3 π - i=1 i π Rc2 Rc3
(9)
The present approach, essentially a dipolar version of the Fennel-potential, needs both a parametrization of α and Rc . These are exactly what we aim to investigate. The presented potential has been independently derived with regard to previously presented work21 though the self-energy term differs slightly.
III. DIELECTRIC CONSTANT
For simulations using the reaction field, the general equation for the relative dielectric constant, ε r , is27 −1 εr − 1 εr − 1 ˜ 1 ⟨M 2⟩ 1− T(0) = , (10) εr + 2 εr + 2 3ε 0 3V k BT ˜ where T(0) = 2(ε RF − 1)/(2ε RF + 1), ⟨M 2⟩ = ⟨| i µ i |2⟩ is the mean squared dipole moment of the simulation box, k B is Boltzmann’s constant, T the temperature, ε 0 is the vacuum permittivity, ε RF is the dielectric constant of the surrounding medium, and V is the volume. Thus, for conducting, insulating, and vacuum boundaries,
014109-3
Stenqvist et al.
εr − 1 3 (2ε r + 1)(ε r − 1) 1 ⟨M 2⟩ = 9ε r 3ε 0 3V k BT εr − 1 εr + 2
J. Chem. Phys. 143, 014109 (2015)
if ε RF = ∞, if ε RF = ε r ,
(11)
if ε RF = 1.
For Ewald with conducting boundaries,28 the expression for the dielectric constant is identical to that of the reaction field with ε RF = ∞. For the present dipolar case, we use the approach of Neumann and Steinhauser29–31 but with the interaction tensor introduced by Eq. (8). We thereby arrive at a similar expression as Zahn et al.,15 Eq. (10) with 2[(αRc )4 + 2(αRc )2 + 3]αRc e−α ˜ T(0) = erf(αRc ) − √ 3 π
2 R2 c
.
(12)
˜ Note that T(0) ≈ 1 in the limit αRc ≫ 1 whereby the conducting Ewald and RF results are retrieved. IV. COMPUTER SIMULATIONS
Metropolis Monte Carlo (MC) simulations were used to generate configurations of the Stockmayer fluid32 in the canonical (NVT) ensemble and were performed using the Faunus software.33 In reduced units, the density was ρ∗ = N σ 3/V = 0.924, temperature T ∗ = k BT/ϵ LJ = 1.333, µ∗2 = µ2/(4πϵ 0σ 3ϵ LJ ) = {1.663, 3.470}, and N = 3000 particles. When we have studied different densities, ρ∗ = {0.231, 0.462, 0.693, 0.924, 1.155, 1.386}, the system size has been constant, L = 10.265σ. The different, spherical cutoffs, Rc , were chosen to be smaller than or equal to half the box length. Simulations were equilibrated in 104 sweeps and ten times longer productions runs using translational and rotational MC moves. In each sweep, N moves were attempted. We used Stockmayer potentials consisting of the electrostatic potential given by the different interaction tensors defined above, combined with a Lennard-Jones (LJ) potential. All Ewald simulations were performed with conducting boundary conditions, a real-space cutoff of half the box length, α = 0.404σ −1 and 686 wavevectors. V. RESULTS
Unless otherwise stated, the reduced, squared dipole moment is µ∗2 = 3.470 but the main conclusions reached are valid also for µ∗2 = 1.663. To simplify the analysis, we address results as acceptable if they are within approximately 1–2 standard deviations (std’s) from the corresponding Ewald result. The presented electrostatic energies are including the self-energy.
FIG. 1. RDF’s (a) and dipole-dipole correlations (b) using different α −1, see legend. The black lines in the main figures correspond to Ewald simulations, often overlapping with the Wolf results. The insets show the corresponding errors; grey area is confined by ±1 Ewald std.
1. Effect of the damping
We initially compare five simulations with different α-values for a constant cutoff, Rc ≃ 7.4σ (L/2), testing the role of the smooth damping. Figure 1(a) shows RDF’s as well as their errors (i.e., the difference, see inset figure) between Ewald and Wolf RDF’s. As can be seen, g(r) is seemingly insensitive to changes of the damping parameter and dominated chiefly by the Lennard-Jones part of the pair interaction. Still, according to the error plot, α −1 = 0.9σ do not satisfy a limit about ±1–2 Ewald std and show oscillating tendencies, even if the relative error is negligible. While oscillating, the remaining α values fall within the error range of the Ewald results. Fig. 1(b) shows that for the dipole-dipole correlation, acceptable α −1 values are also obtained between 1.2σ and 2.1σ. Table I shows the average energy per particle converging to that of Ewald for all α −1 but 0.9σ. The same trend is seen for the dielectric constants.
A. Structural properties
In this section, we evaluate RDF’s of the Stockmayer fluid, and for a more critical measure of the anisotropic electrostatics, ˆ · also dipole-dipole correlation functions, defined as ⟨ µ(0) ˆ ˆ µ(r)⟩ where µ(0) is a central, normalized dipole surrounded ˆ by dipoles µ(r) in a spherical volume element at distance r + dr. The angular brackets denote an ensemble average.
2. Effect of the cutoff
We now study the opposite case: constant α and varying cutoffs. Results for α −1 = 1.5σ and α −1 = 1.8σ, are presented in Figs. 2 and 3 and in Table II. For α −1 = 1.5σ, it is clear that Rc = 5σ is a saturated cutoff, i.e., the results will not improve by increasing it further. In the α −1 = 1.8σ case, half the
014109-4
Stenqvist et al.
J. Chem. Phys. 143, 014109 (2015)
TABLE I. Dipolar, Lennard-Jones, and total particle energies as well as dielectric constants. The std is 0.002 for βu µ µ and βu LJ where β = 1/k BT . α −1
βu µ µ
βu LJ
βu Tot
εr
0.9σ 1.2σ 1.5σ 1.8σ 2.1σ Ewald
−4.552 −4.523 −4.520 −4.518 −4.517 −4.523
−4.000 −3.996 −3.998 −3.996 −3.996 −3.995
−8.552 −8.520 −8.518 −8.514 −8.513 −8.518
100 ± 30 120 ± 20 120 ± 40 130 ± 50 140 ± 40 140 ± 40
box-length is needed to get accurate results. This cutoff may be a saturated cutoff—more on this Sec. V B. B. Energies
Fig. 4(a) shows that the LJ energy is nearly constant when α −1 ≥ σ for every cutoff, except Rc = 2σ. The electrostatic energy, see Fig. 4(b), has a stationary or quasi-stationary point (indicated by dots) in the interval σ ≤ α −1 ≤ 2.3σ for every cutoff but Rc = 2σ. By using the least-square method, the fitted lines corresponding to the quasi-stationary points are described by Rc ≈ 3.2α −1 + 0.1σ for µ∗2 = 3.470 and Rc ≈ 3.1α −1 − 0.3σ for µ∗2 = 1.663. Fig. 4(b) shows that the larger the cutoff, the wider the width of the “stationary”
FIG. 3. RDF’s (a) and dipole-dipole correlations (b) with α −1 = 1.8σ using different cutoff’s, cf. legend; black lines the Ewald results, often overlapping with the Wolf results. The insets show the corresponding errors; black indicates ±1 Ewald std.
segment around the quasi-stationary point. Furthermore, in this section, we see that for α −1 = 1.8σ, a mere cutoff of 6σ is needed to get the correct energy, thus in this aspect Rc = 6σ is a saturated cutoff. TABLE II. Dipolar, Lennard-Jones, and total energies per particle as well as dielectric constants for α −1 = 1.5σ and α −1 = 1.8σ. The std is ∼ 0.002 for βu µ µ and βu L J .
FIG. 2. RDF’s (a) and dipole-dipole correlations (b) with α −1 = 1.5σ using different cutoff’s, cf. legend; black lines are Ewald results, often overlapping with the Wolf results. The insets show the corresponding errors; grey area is confined by ±1 Ewald std.
Rc
βu µ µ
2σ 3σ 4σ 5σ 6σ L/2 Ewald
−3.021 −4.138 −4.457 −4.518 −4.520 −4.520 −4.523
2σ 3σ 4σ 5σ 6σ L/2 Ewald
−2.876 −3.964 −4.349 −4.485 −4.514 −4.518 −4.523
βu LJ α −1 = 1.5σ −4.039 −3.991 −3.993 −3.995 −3.995 −3.998 −3.995 α −1 = 1.8σ −4.043 −3.995 −3.987 −3.995 −3.998 −3.996 −3.995
βu Tot
εr
−7.060 −8.130 −8.450 −8.513 −8.515 −8.518 −8.513
10 ± 10 10 ± 10 40 ± 10 130 ± 40 120 ± 30 120 ± 30 140 ± 40
−6.919 −7.959 −8.337 −8.479 −8.512 −8.514 −8.513
10 ± 10 10 ± 10 10 ± 10 60 ± 10 110 ± 40 130 ± 50 140 ± 40
014109-5
Stenqvist et al.
J. Chem. Phys. 143, 014109 (2015)
FIG. 4. βu LJ (a) and βu µ µ (b) at different cutoff values. The insets show the corresponding energy as a function of the cutoff for α −1 = 2.1σ and the Ewald energy. The energy plateau in (a) starts at α −1 = σ for every cutoff and ends with the equally colored star.
C. Dielectric constant
In this section, we investigate how damped electrostatic interactions influence the dielectric constant which is a measure of the collective polarization of the liquid. Throughout the analysis, the general behavior as a function of Rc and α −1 is preserved, see Figs. 5(a) and 6(a). Some information can readily be extracted from Fig. 5(a): The dielectric constant is a macroscopic property and to capture this behavior, a limit on the damping and minimal cutoff are required, α −1 ≥ σ and Rc ≥ 3σ. At the point Rc = 3σ and α −1 = σ, a triangle starts to form, inside of which acceptable results emerge. There are two interesting lines: α −1 = σ and Rc ≈ 3α −1 + 0.7σ. The area between these are statistically equal (around ±1 std) to the Ewald result. A large interval easy to target is α −1 = σ − 2.2σ when Rc is half the box length. From the analysis of the RDF’s and dipole-correlations, this is too a satisfactory interval. An extended cutoff increases the triangle and widens the possible margin of error in the choice of α such that acceptable results are still obtained. This is in close agreement with our earlier conclusion coming from the analysis of the energies. From this (limited) parameter space, we deduce that with increasing α −1 must come an even faster increase in cutoff, satisfying Rc ≥ 3α −1 + 0.7σ. Also, the dots in Fig. 4(b) can be seen in corresponding color in Fig. 5(a).
∆ε (a), and ε std in % (b) for µ ∗2 = 3.470—values outside the FIG. 5. ε Ewald Ewald color-bars are shown as white. The solid, red lines represent α −1 = σ and −1 R c = 3α + 0.7σ and form a triangular area within which valid results (±1 std) are obtained. The colored dots correspond to (quasi-) stationary points from Fig. 4(b). The dielectric constant for Ewald is ε Ewald = 142.1 ± 36.1.
The arguments established so far hold also for the less coupled µ∗2 = 1.663 case, see Fig. 6(a), where the lines are α −1 = σ and Rc ≈ 2.7α −1 + 0.4σ. The main difference is an effectively larger triangle which is reasonable since µ∗2 → 0 should trivially result in the entire parameter space being valid. D. Densities
So far we have studied a single density, ρ∗ = 0.924 and here we explore more dilute as well as denser systems for constant L = 29.629 Å. A cutoff of half the box-length was used in all simulations in this section. The results, condensed in Fig. 7, show that the dipole-dipole correlations become less accurate for the lowest and highest densities while the RDF only suffers from low densities. This may be due to an appearance of non-isotropic or crystalline phases; networks, ring, and chain formations in the dilute regime, and ferro-electric liquids/solids in the dense limit.34 The long-range order in these phases is suppressed when electrostatics are truncated as in the Wolf formalism which has already been pointed out.14
014109-6
Stenqvist et al.
J. Chem. Phys. 143, 014109 (2015)
FIG. 7. RMS of error between Wolf and Ewald for RDF’s (a) and dipoledipole correlations (b) using different reduced densities, cf. the legends. ∆ε (a), and ε std in % (b) for µ ∗2 = 1.663—values outside the FIG. 6. ε Ewald Ewald color-bars are shown as white. The solid, red lines represent α −1 = σ and R c = 2.7α −1 + 0.4σ and form a triangular area within which valid results (±1 std) are obtained. The colored dots correspond to (quasi-) stationary points in the electrostatic and self-energy. The dielectric constant for Ewald is ε Ewald = 16.4 ± 2.3.
In between these extreme densities, correct results emerge. The triangular area—compare, e.g., Figs. 5(a) and 6(a)—in parameter space giving correct results tends to narrow as the dielectric constant increases. At constant dipolar coupling, increasing the density leads to a larger dielectric constant which may explain the smaller window (in α −1) of accuracy in Fig. 7(b) for high densities; vice versa for low densities. The low-density precision windows are however here most likely limited by the small cutoff (determined by the box-length). Although the original Wolf method for charges is expected to work best for liquids, reasonable Madelung energies are obtained for a number of crystal lattices14 which is encouraging for treating dipolar crystals even if it is not guarantied to provide accurate results due to the anisotropic nature of the potential. E. Reaction field
The reaction field energies are shown in Table III. For ε RF = 1, energies are generally lower than for Ewald. The
dipole-dipole energies under these conditions are non-monotonic, narrowing the gap to the Ewald dipole-dipole energy with increasing cutoff, as expected, until Rc = 4σ, where after and the gap is widened again. This behavior makes the TABLE III. Reaction field dipole-dipole and LJ energies in units of k BT for different cutoff and dielectric boundaries. The Ewald dipole-dipole and total energies are −4.523 k BT and −3.995 k BT , respectively. The std is ∼ 0.002 for both βu µ µ and βu L J . Rc
ε RF = 1
2σ 3σ 4σ 5σ 6σ 7σ L/2
−4.710 −4.658 −4.644 −4.666 −4.742 −4.816 −4.839
2σ 3σ 4σ 5σ 6σ 7σ L/2
ε RF = ∞
Dipole-dipole energy −4.703 −4.520 −4.555 −4.514 −4.512 −4.512 −4.513 LJ energy −3.910 −3.929 −3.997 −3.998 −3.982 −3.995 −4.003 −3.996 −4.010 −3.997 −4.020 −3.997 −4.024 −3.997
ε RF = ε r
−4.695 −4.513 −4.547 −4.509 −4.512 −4.512 −4.512 −3.928 −3.997 −3.994 −3.997 −3.995 −3.996 −3.995
014109-7
Stenqvist et al.
J. Chem. Phys. 143, 014109 (2015)
TABLE IV. Dielectric constants ± their std’s from reaction field simulations. The corresponding Ewald dielectric constant is 140 ± 40. Rc
ε RF = 1
ε RF = ∞
ε RF = ε r
2σ 3σ 4σ 5σ 6σ 7σ L/2
−50 ± 240 90 ± 110 70 ± 100 50 ± 30 100 ± 90 40 ± 30 50 ± 30
390 ± 160 150 ± 30 200 ± 80 150 ± 60 110 ± 30 110 ± 30 130 ± 40
420 ± 220 140 ± 30 220 ± 50 120 ± 40 120 ± 20 130 ± 30 120 ± 30
dipole-dipole energy (and thus the total energy) seem to diverge with increasing cutoff. Since the reaction field method is identical to both the Wolf and Ewald methods in the limit Rc → ∞, we deduce that our parameter space for the cutoff is too small to capture the convergence. The inaccuracy when ε RF = 1 is also seen in Table IV as proportionally large deviations from Ewald for any tested cutoff. When ε RF = ε r or ε RF = ∞, energies do converge with increasing cutoff, while the dielectric constants are indistinguishable within one std. This is true for any cutoff larger than 5σ after which the std’s remain unaffected by further increase in Rc . The contrast
to the vacuum case can be rationalized by the fact that the more polarizable surroundings tend to dampen electrostatic correlations, whereby convergence is reached even for small system sizes. Figure 8(a) shows that, as previously stated, ε RF = 1 gives deviating results compared to Ewald. Surprisingly, the smaller cutoff gives noticeable more accurate results than for half the box length. Another abnormal result is a bias on the L/2 ε RF = ∞ case. The remaining cases are all acceptable. When including the more sensitive dipole-dipole correlations in Fig. 8(b), we see that there is no evidence of any of the cases to be accurate. An increasing cutoff seems to damp the fluctuations, although the increasing trend is unaffected.
VI. CONCLUSIONS
We have presented an expanded Wolf formalism for long range dipolar interactions, evaluated on Stockmayer fluids where structural and dielectric properties are investigated. Comparisons are made with Ewald summation and a reaction field method. In short, the α-parameter in the Wolf method should respect α −1 ≥ σ and αRc > 3 to get values in agreement with Ewald results. The latter limit ensures sufficiently soft damping before the hard cutoff is introduced. This points to the fact that even if both the potential and force are shifted, a hard cutoff alters the system fluctuations, while soft damping does not suffer from similar effects. Combining the two rules leads to a cutoff criteria of Rc > 3σ. Increasing Rc beyond αRc ∼ 3 at fixed α gives little or no increase in accuracy—at least for densities around ρ∗ = 1. Finding the (quasi-) stationary electrostatic energy with regard to α seems to give parameters with accurate results. If no such stationary point can be found, or its width is small, it is an indication that the chosen cutoff is too short. Previous work22 on multipolar fluids (dipoles, and quadrupoles) in the micro-canonical ensemble is qualitatively consistent with the damping parameter interval deduced in this work. Finally, Wolf with optimal parameters presents more or equally accurate results as the equivalently scaling RF method used with either conducting, insulating, or vacuum boundaries. ACKNOWLEDGMENTS
We thank LUNARC for computational resources and for financial support: Södra Research foundation; the Swedish Research Council; and the Swedish Foundation for Strategic Research. Finally, we thank Per Linse for useful discussions and comments.
APPENDIX A: SELF-ENERGY
FIG. 8. Reaction field RDF’s (a) and dipole-dipole correlations (b) using different ε RF-values (labeled); the black line is from Ewald. The cutoffs are 5σ for the three top legends and half the box length (L/2) for the three bottom legends. Note that most RDF’s and correlations are overlapping. The insets show the corresponding errors; grey area is confined by ±1 Ewald std.
In this section, we will use a direct expansion of the point charge approach in Wolf14 to dipoles, in order to get the expanded DSF self-energy. Since we in this section are referring to similar equations in Wolf,14 we will, when possible, use similar notation. The damped version of Wolf uses a unity multiplication, erfc(α|ri j + n|) + erf(α|ri j + n|) = 1, of
014109-8
Stenqvist et al.
J. Chem. Phys. 143, 014109 (2015)
the total energy of a periodic boundary condition system (see Wolf14 Eq. (5.1)). For dipoles, we get
where (see Wolf14 Eq. (5.12))
) ( N ∞ 1 T T erfc(α|ri j + n|) µj ETot = − µ ∇∇ 2 i=1 j,i=1 i |ri j + n|
DSP ESelf =−
) ( N ∞ 1 T T erf(α|ri j + n|) − µ j. µ ∇∇ 2 i=1 j,i=1 i |ri j + n|
(A1)
The second part of Eq. (A1) is large for large values of α. By extracting the diverging term, we end up with (see Wolf14 Eqs. (5.2)–(5.5)) N ∞ 1 T D µ T (ri j + n, α)µ j 2 i=1 j,i=1 i 2 ( ) N N 1 T T erf(α|ri j + n|) µ i ∇∇ µj − 2 i=1 j=1 |ri j + n|
ETot = −
N 2α 3 2 − √ µi , 3 π i=1
(A2)
where TD 2 (r, α) is described by Eq. (7). By making this extraction, the remainder (i.e., the second part of Eq. (A2)) is now a small correction. Hence it is neglected, giving (see Wolf14 Eq. (5.6)) Approx
ETot
=−
1 2
N i=1
µTi T2D (ri j + n, α)µ j −
j,i |r i j |< R c
N 3
2α √ 3 π
µ2i .
i=1
(A3) By using the condition about cancellation of dipole moment inside some cutoff (smaller than half the box-length), the expression becomes (see Wolf14 Eq. (5.7)) DSP ETot =−
i=1
N
(A8)
All described so far is a strict expansion of the approach for charges in Wolf14 for dipoles. The only assumption is that the compensation for the centred dipole is in the form of a smeared out dipole on the cutoff sphere. However, one can interpret Wolf14 as if the compensating charge is in the form of a smeared out charge on the surface, i.e., just like we have presented but for dipoles. Hence, the presented approach is backward compatible with Wolf.14
Ewald EReal
N N ′ ( µi · µ j 1 = B(|ri j + n|) 2 i=1 j=1 |ri j + n|3 3
) ((ri j + n) · µ i )((ri j + n) · µ j ) − 3C(|ri j + n|) , |ri j + n|5 (B1) where (A5)
which can be written as (see Wolf14 Eq. (5.10)) 2
2 2
1 erfc(αRc ) 2α e−α R c 4 α 3 + 2 +√ + √ µ. =− * 2, 3 π - i=1 i π Rc2 Rc3
n∈Z
*. +/ N 1 lim ..− µTi T2D (ri j , α)µ j // |ri j |→ R c . 2 / j,i i=1 |ri j |< R c , N 1 µTi T2D (r, α) r= R c r µ i + , + lim *− |ri i |→ R c |r i i | i i , 2 i=1 -
DSP ETot =−
DSP DSF ESelf = ESelf
(A4)
where (see Wolf14 Eq. (5.9))
N 1
First of all, we acknowledge that the term rii is ill-defined, the direction seems arbitrary. To remedy this we, instead of a point dipole use a smeared out dipole on the cutoff sphere. Since it is the cancellation of the zeroth order electrostatic moment (i.e., charges) that is the principle of Wolf14 which we extrapolate to first order moment (i.e., dipoles), this is still within the framework of the expanded formalism. The dipole moment density is chosen such that ρ(θ) = ρ0 cos(θ) where ρ(θ) is the dipole density in spherical coordinates (θ ∈ [0, π]), using a ρ0 such that the integrated dipole moment density on the sphere is oppositely signed but equal in magnitude to the dipole moment of the centred dipole. The DSP self-energy one then gets is actually the same as the DSF selfenergy,
Here, we show the equations used in the Ewald summation method. The real-space (short ranged) part shows similarities with the Wolf method,
|ri j |< R c
Cancel ETot =
(A7)
APPENDIX B: EWALD EQUATIONS
N 1 µTi T2D (ri j , α)µ j 2 i=1 i, j
N 2α 3 2 Cancel µi − ETot , − √ 3 π i=1
N N 1 2α 3 2 µi . µ i T2D (r, α) r= R c r µ i + √ |ri i | i i 2 i=1 3 π i=1
µTi T2D (ri j , α)
i, j |ri j |< R c
) |r i j | i j
DSP , µ j − ESelf
(B2)
( ) 2αr 3 3 2 2 C(r, α) = erfc(αr) + √ 2α 2 + 2 e−α r . r 3 π
(B3)
and
(
− T2D (r, α) r= R c r
2αr 2 2 B(r, α) = erfc(αr) + √ e−α r π
(A6)
The reciprocal-space part takes care of the long-ranged interactions,
014109-9
Stenqvist et al.
Ewald EReciprocal
J. Chem. Phys. 143, 014109 (2015)
( 2 2) 1 4π π |k| = exp − 2 2 3 2 2L k,0 |k| α L
7D.
k∈Z3
) ( N N 2πik · rmn . × (µ m · k)(µ n · k) exp L m=1 n=1 (B4) Finally, Eqs. (B5) and (B6) give the self-energy and the surface energy, arising from an interaction between the enclosed system and a surrounding dielectric medium with dielectric constant ε Surface, N 2α 3 Ewald ESelf =− √ |µ i |2 3 π i=1
(B5)
and Ewald ESurface =
2π (2ε Surface + 1)L
N N (µ i · µ j ), 3
(B6)
i=1 j=1
where we used L = L x = L y = L z . 1B.
W. McCann and O. Acevedo, J. Chem. Theory Comput. 9, 944 (2013). Levitt, M. Hirshberg, R. Sharon, and V. Daggett, Comput. Phys. Commun. 91, 215 (1995). 3B. Linse and P. Linse, J. Chem. Phys. 141, 184114 (2014). 4P. Ewald, Ann. Phys. 369, 253287 (1921). 5T. Darden, D. York, and L. Pedersen, J. Chem. Phys. 98, 10089 (1993). 6J. W. Perram, H. G. Petersen, and S. W. De Leeuw, Mol. Phys. 65, 875 (1988). 2M.
E. Shaw, M. M. Deneroff, R. O. Dror, J. S. Kuskin, R. H. Larson, J. K. Salmon, C. Young, B. Batson, K. J. Bowers, and J. C. Chao, ACM SIGARCH Comput. Archit. News 35, 1 (2007). 8J. Stenhammar, P. Linse, and G. Karlström, J. Phys.: Condens. Matter 20, 494204 (2008). 9J. Hansen, D. Levesque, and J. Weis, Phys. Rev. Lett. 43, 979 (1979). 10J. Caillol and D. Levesque, J. Chem. Phys. 94, 597 (1991). 11J. Caillol, J. Chem. Phys. 96, 1455 (1992). 12M. Trulsson, J. Chem. Phys. 133, 174105 (2010). 13C. J. Fennell and J. D. Gezelter, J. Chem. Phys. 124, 234104 (2006). 14D. Wolf, P. Keblinski, S. Phillpot, and J. Eggebrecht, J. Chem. Phys. 110, 8254 (1999). 15D. Zahn, B. Schilling, and S. M. Kast, J. Phys. Chem. B 106, 10725 (2002). 16Y. Yonezawa, J. Chem. Phys. 136, 244103 (2012). 17G. S. Fanourgakis, J. Phys. Chem. B 119, 1974 (2015). 18A. Martin, London, Edinburgh, Dublin Philos. Mag. J. Sci. 8, 547 (1929). 19R. P. Bell, Trans. Faraday Soc. 27, 797 (1931). 20L. Greengard and V. Rokhlin, J. Comput. Phys. 73, 325 (1987). 21M. Lamichhane, J. D. Gezelter, and K. E. Newman, J. Chem. Phys. 141, 134109 (2014). 22M. Lamichhane, K. E. Newman, and J. D. Gezelter, J. Chem. Phys. 141, 134110 (2014). 23S. de Leeuw, J. Perram, and E. Smith, Proc. R. Soc. A 373, 27 (1980). 24Z. Wang and C. Holm, J. Chem. Phys. 115, 6351 (2001). 25J. Kolafa and J. W. Perram, Mol. Simul. 9, 351 (1992). 26J. Barker and R. Watts, Mol. Phys. 26, 789 (1973). 27M. Neumann, Mol. Phys. 50, 841 (1983). 28E. Pollock and B. Alder, Phys. A 102, 1 (1980). 29M. Neumann and O. Steinhauser, Chem. Phys. Lett. 95, 417 (1983). 30M. Neumann, O. Steinhauser, and G. S. Pawley, Mol. Phys. 52, 97 (1984). 31M. Neumann, Mol. Phys. 57, 97 (1986). 32W. H. Stockmayer, J. Chem. Phys. 9, 398 (1941). 33B. Stenqvist, A. Thuresson, A. Kurut, R. Vácha, and M. Lund, Mol. Simul. 39, 1233 (2013). 34C. Holm and J.-J. Weis, Curr. Opin. Colloid Interface Sci. 10, 133 (2005).
