Structural Anisotropy in Polar Fluids Subjected to Periodic Boundary Conditions.

A heuristic model based on dielectric continuum theory for the long-range solvation free energy of a dipolar system possessing periodic boundary conditions (PBCs) is presented. The predictions of the model are compared to simulation results for Stockmayer fluids simulated using three different cell geometries. The boundary effects induced by the PBCs are shown to lead to anisotropies in the apparent dielectric constant and the long-range solvation free energy of as much as 50%. However, the sum of all of the anisotropic energy contributions yields a value that is very close to the isotropic one derived from dielectric continuum theory, leading to a total system energy close to the dielectric value. It is finally shown that the leading-order contribution to the energetic and structural anisotropy is significantly smaller in the noncubic simulation cell geometries compared to when using a cubic simulation cell.

Introduction

Molecular dynamics (MD) and Monte Carlo (MC) simulations have grown to become a central tool in physics, chemistry, and biology over the past three decades.[1,2] However, in spite of the huge advancement of both algorithms and hardware, there are still some unresolved methodological issues. Arguably, the most persistent of these is the question of how to handle long-range electrostatic (Coulomb and dipole–dipole) interactions in a simulation.[3−5] The basic problem is that the integraldiverges for all finite values of the cutoff radius rcut as long as the intermolecular potential v(r) does not decay faster than r–3. Thus, applying a simple (spherical or cubic) cutoff to the electrostatic potentials may often lead to serious artifacts in the structure and thermodynamics of the system under study.

Although several solutions to the infinite-range interaction problem have been proposed, the most common way to circumvent this problem is the use of lattice-based summation techniques, or periodic boundary conditions (PBCs). These methods compose a plethora of different algorithms that all rest on the same basic assumption, namely that the (finitely sized) simulation cell is duplicated in all directions to create an infinite lattice. The original implementation of this idea was developed by Ewald[6] and is built upon a separation of the interaction into short-range and long-range parts, where the former is summed up in real space and the latter in reciprocal space. The original Ewald method has since been developed in many ways, and today different mesh-based methods[7−10] are numerically faster alternatives to the classical Ewald summation.

When simulating a fluid phase, the assumption of periodicity is clearly not a correct description of the real system. This criticism has been put forward several times in the literature but was originally noted by Valleau and Whittington,[11] who gave a qualitative argument about the inability of lattice summation methods to correctly reproduce long-range fluctuations in fluid systems. Furthermore, several studies have addressed the issue of periodicity effects on the properties of Lennard-Jones fluids,[12,13] ionic solutions,[14−17] and biomolecules.[18−21] In the context of dipolar systems, Boresch and Steinhauser[22] conducted a careful study of dipole fluctuations and correlations in SPC water simulated using the Ewald summation technique. In particular, they addressed the importance of the so-called surface term,[23] which describes the solvation from the dielectric surroundings of the infinite lattice on structural properties such as the dielectric permittivity, dipole time correlation functions, and the Kirkwood g factor. However, the total dipole moment of the simulation box is a special property, in the sense that its total interaction with all its periodic images is identically zero, as long as the contributions are summed in spherical shells.[24−27] Therefore, the periodicity effects on the fluctuating dipole moment of the whole simulation box (and related properties) are expected to be small. In a recent contribution,[28] we showed, however, that the fluctuations of higher order electric multipole moments of the whole simulation box are greatly influenced by the interaction between each instantaneous multipole and all of its periodic images. This effect is manifested through a difference of as much as 50% between the dielectric permittivities calculated from different multipole components, depending on whether the multipole component has an attractive or a repulsive (or, in some cases, zero) interaction with its neighbors. A schematic picture of the coupling of different multipole components in a system under PBCs is given in Figure 1.

Figure 1

Schematic picture of the coupling of the total dipole (left) and higher multipoles (right) of a simulation cell subjected to PBCs. The dipole does not “see” its neighbors since its self-interaction energy is zero but is solvated by the dielectric response from the surrounding medium through the surface term. In contrast, higher multipoles (Q, l > 1) couple to its neighbors through their nonzero self-interaction but are not affected by the surface term. In addition, the dipole as well as the higher multipoles interact with the set of “unconstrained” multipoles Q, l′ ≠ l and/or m′ ≠ m, (not depicted) in the surrounding cells, giving an (approximately) isotropic contribution to the solvation.

In addition to the cubic simulation cell used in the majority of computer simulations, some alternative simulation cell geometries have been suggested and implemented,[29−33] most notably the rhombic dodecahedron (RD) and the truncated octahedron (TO). These two bodies have the appealing property of more closely resembling the geometry of solvated spherical solutes, in the sense that they have larger inscribed spheres than a cubic simulation cell of the same volume. Even though these alternative geometries are implemented in major simulation packages, there are only a few studies[29,34,35] probing the effect of changing the cell geometry on the thermodynamic properties of the system under study. Because these cells pack in lattice structures different from that of the cube, it seems reasonable to expect their periodicity effects to differ qualitatively from those of a cubic cell.

In the present contribution, we will extend our previous analysis[28] of the periodicity effects on a dipolar model system from the qualitative to the quantitative level, as well as from cubic to noncubic simulation cells. We will develop a heuristic model describing the solvation of and electrostatic fluctuations in a spherical subvolume of a dielectric medium exhibited to PBCs. This model will be compared to values of the dielectric constant calculated from simulations of a simple dipolar model system.

Theory

General ansatz

In the following, the electrostatic fluctuations in a spherical subvolume of a dipolar model system treated using the Ewald summation technique will be described by dividing the long-range solvation energy of this subvolume into two contributions:

An approximately isotropic part, coming from the interaction between the instantaneous multipole Q of a spherical volume in the central simulation cell and its noncorrelated neighbors, i.e., Q with l′ ≠ l and/or m′ ≠ m, in the other cells. This interaction is, at least partly, Boltzmann-weighted in a simulation, and we will thus attempt to describe it using formulas valid for an isotropic dielectric medium.

A strongly anisotropic part, coming from the “self-interaction” between Q in the central cell and its fully correlated replicas (Q with l′ = l and m′ = m) in the rest of the lattice. This part of the interaction is not Boltzmann-weighted, because of the perfect periodicity imposed by the PBCs. We will thus describe this interaction using the reduced lattice-interaction tensors introduced previously.[28]

On the basis of this description, we will present a heuristic derivation of the long-range solvation free energy of the spherical subvolume. This will be compared to the behavior expected from a spherical subvolume inside an infinite isotropic dielectric medium and the analysis will thus enable us to directly probe the magnitude of the periodicity effect introduced by the PBCs.

Periodic Boundary Conditions

In the present study, the term “periodic boundary conditions” refers to a system with a potential energy Upot of the formwhere n = (n, n, n) is a vector that runs over all lattice points in the particular (unit length) lattice and a denotes the side length of the unit cell. Furthermore, the primed sum indicates that the term with i = j for n = 0 should be excluded, and v(r, ω, ω) denotes the intermolecular potential between particles i and j, depending in general on their separation r and orientations ω and ω. In practice, v(r, ω, ω) is usually long-range in the sense that it decays no faster than r–3, the two most important examples being the Coulomb and dipole–dipole potentials.

Since the sum in eq 2 is slowly (and conditionally) convergent, more elaborate methods to evaluate the potential energy in a PBC system need to be used in practice. The by far most popular technique to achieve a fast convergence of the potential energy is the technique originally due to Ewald[6] and different mesh-based variants[7−10] thereof. Within the Ewald-based methods, the short-range (n = 0) part of Upot is screened through the addition of a Gaussian charge (dipole) cloud and is thereafter summed within a, usually spherical, cutoff after considering the nearest image convention.[1] The long-range (n ≠ 0) part of the potential energy is summed up in Fourier space, leading to a quickly (and absolutely) convergent sum.

Simulation Cells with Noncubic Geometries

Although cubic simulation cells are used for the majority of simulation studies, the use of alternative simulation cell geometries has started to become increasingly popular. In total, five classes of geometrical bodies are translationally space-filling and can thus be used for simulating a periodic system;[31] however, due to their relatively “sphere-like” symmetry, the two most useful alternatives to the cube, at least for the simulation of bulk systems, are the rhombic dodecahedron (RD) and truncated octahedron (TO). While the cube, of course, packs in a simple cubic (SC) lattice structure, the natural choice for the lattice structures of the RD and TO are face-centered cubic (FCC) and body-centered cubic (BCC), respectively.[31] However, Smith and Fincham[36] showed that the use of a body-centered tetragonal (BCT) lattice structure for the RD, with one side of the unit cell elongated by a factor √2 compared to the other two, greatly facilitates the implementation of Ewald summation for this geometry, by simply excluding k-space terms of certain parity. Thus, we will use the BCC and BCT lattice structures as the basis of our analysis. In Figure 2, the RD and the TO are shown, inscribed in their respective unit cells.

Figure 2

The rhombic dodecahedron (left) and truncated octahedron (right) inscribed in their BCT and BCC unit cells, respectively.

Solvation in Dielectric Media

In the following subsections, we will treat relevant parts of the theory of solvation and fluctuations in dielectric media. First, we will review the theory for the solvation of a polarizable dipole in a dielectric medium (section ). Subsequently, in section , we will treat electrostatic fluctuations and solvation in isotropic dielectric media. Finally, in section , we will use tools from the two preceding parts to develop a heuristic model for the solvation of a dielectric subvolume in a PBC system.

Generally, the collective electrostatic fluctuations will be quantified through the spherical multipole moments Q, defined throughwhere ρ(r) denotes the charge density in a point r = (r, Ω) = (r, φ, θ) ∈ V and C(Ω) represents Racah’s unnormalized spherical harmonics. The index l denotes the order of the multipole, whereas m describes its orientation in an external coordinate frame. Just as for the spherical harmonics, m takes on all integer values between –l and +l. However, the –m and +m components are related according towhere * denotes complex conjugation; thus, Q and Q are not independent degrees of freedom. Instead, we will adopt the approach taken previously[37,38] and treat separately the real and imaginary parts of Q, denoted respectively by superscripts R and I, for m ≥ 0. Since Q is real, these l + 1 real and l imaginary multipole components form 2l + 1 linearly independent fluctuation modes. Furthermore, we will use the bracketed superscripts (R) and (I) to denote quantities that are somehow related to the real and imaginary parts of Q, although not themselves complex quantities.

A Polarizable Dipole in a Dielectric Medium

The solvation energy Usolv of a polarizable point dipole of magnitude μ, radius R, and polarizability α embedded in a dielectric medium of dielectric permittivity ε is given by[39]wherequantifies the reaction field, parallel to the dipole, coming from the surrounding dielectric medium. A physical interpretation of the expression for Usolv is facilitated by expanding eq 5 in a geometric series, i.e.From this expression, we can identify the prefactor – gμ2/2 as the solvation energy of a permanent dipole immersed in a dielectric medium, whereas the factor ∑(gα) takes into account the increase of the solvation energy due to the additional polarization of the particle by the reaction field. The infinite sum is due to the incremental nature of this process; the reaction field increases the total dipole moment of the particle, which in turn polarizes the dielectric to yield a larger reaction field, etc. In section , we will show how this partitioning of the solvation energy can be mapped onto the solvation of a dielectric subvolume exhibited to PBCs.

Electrostatic Fluctuations in a Dielectric Medium

The (unnormalized) probability distribution Pvac(Q) of the 2-pole moment of a spherical dielectric volume with radius R and dielectric permittivity ε in a vacuum is given by the Gaussian function[37,40]where Uvac†, given bydenotes the (free) energy cost for creating an instantaneous multipole moment Q in the dielectric volume, β = (kBT)−1 is the inverse thermal energy, and X ∈ {R, I}. If the dielectric volume is immersed in an infinite dielectric medium with the same value of ε as the sphere itself, U† is decreased due to the depolarizing reaction field from the surrounding medium, changing the probability distribution to[37,40]whereand the second equality is accurate for not too small values of ε. The energy expression in eq 9 is roughly independent of ε for high-dielectric media and thus not numerically useful for determining ε from a computer simulation. In contrast, the right-hand-side of eq 11 shows that Udiel† ∼ ε–1 for large and intermediate values of ε, and thus determining the width of Pdiel in a computer simulation can be used to determine the dielectric permittivity of the system under study. More specifically, eq 10 can be transformed into a formula for the mean-square quantity ⟨(Q)2⟩ by noting that the Gaussian form implies that ⟨(Q)2⟩ = (2βc)−1. After some rearrangements, still using the simplified form for high and intermediate ε, we get

The l = 1 case of eq 12 applied to the total dipole moment has been widely used to determine ε from computer simulations of fluids, although care needs to be taken to use a form of the formula proper for the particular boundary conditions being used.[41,42] In an infinite, isotropic dielectric medium, ε is by definition independent of l, m, and R, but for a finite and/or molecular system, this does not necessarily hold. In particular, as we will show below, ε is not independent of l and m for a system exposed to PBCs.

In addition to the dielectric permittivity, the above formulas can be used to obtain the free energy change ΔAvac→diel of bringing the dielectric sphere from vacuum into its own medium. To this end, we will employ the standard relationship[43]where Zdiel and Zvac denote the configuration integrals in the solvated and nonsolvated states, respectively. Inserting eqs 8–11 and carrying out the integrations givesThus, ΔAvac→diel is (i) always negative, (ii) independent of m and R, and (iii) only weakly dependent on l, a dependence that disappears quickly in the limit l → ∞ . Finally, we note that ΔAvac→diel diverges logarithmically as ε → ∞ for all l.

Electrostatic Fluctuations in a System Subjected to PBCs

We will now propose a mapping of the energy expression in eq 7 for a polarizable dipole in a dielectric medium onto the solvation of a dielectric subvolume in a system exposed to PBCs. As a first assumption, we will describe the energy of creating an instantaneous multipole moment Q in the spherical volume, excluding the anisotropic part of the solvation, by the same expression as in an infinite dielectric medium. Using eq 11 and eq 7, we thus make the assignmentObviously, Udiel† is qualitatively different from the prefactor – gμ2/2 of eq 7; most importantly, it has a positive rather than a negative sign, since it also includes the energetic cost of creating the multipole moment in the dielectric medium, whereas the energy in eq 7 is valid for a permanent dipole, i.e., excluding the self-energy of the charge distribution.

In addition to the isotropic solvation, the instantaneous multipole moment induces a generalized reaction field coming from its own replicas in all of the surrounding boxes, which in turn polarizes the dielectric volume. This behavior is fully analogous to the polarization of a polarizable dipole by its own reaction field; however, in the case of PBCs the reaction field is not proportional to the factor g of eq 6 but rather to the lattice interaction tensor S( quantifying the interaction between the multipole component Q and all its replicas in the lattice. The use of eqs 44 and 47 of ref (38) leads us to the following definition of S(:where a is the side length of the unit cell, δ is the Kronecker delta, and the function f is defined bywith (···) representing the Wigner 3j symbol.[44] The two terms for m > 0 come from the interaction between Q in the central unit cell and (i) Q and (ii) Q in the surrounding cells. We furthermore note that (see Appendix A)where the sum runs over all multipole components with a given l. Thus, the (unweighted) mean value of the lattice interaction for any l ≥ 1 is zero for all multipoles and lattices. It should finally be clarified that, whereas g always yields an attractive coupling between the polarizable dipole and the reaction field, S( can represent attractive as well as repulsive couplings, depending on the symmetry properties of each multipole. Finally, we make the assumption that the polarizability α in eq 7 can be mapped according to α ∼ kselfR2, where kself is a positive constant related to the magnitude of the anisotropic solvation.

On the basis of the above discussion, we suggest that the (free) energy for creating an instantaneous multipole moment in a spherical subvolume of a dielectric exhibited to PBCs is given by

In accordance with eq 10, we also form the corresponding probability distribution PPBC(Q):

We note that, just like in the case of isotropic dielectric solvation, PPBC is Gaussian (as long as kselfR2S( < 1), but with the important difference that its exponent is now m-dependent through the dependence on S(. The Gaussian form with respect to Q implies that the mean-square multipole moment ⟨(Q)2⟩PBC can be expressed as

In analogy with eq 12, we now define the apparent dielectric permittivity εPBC,( aswhich from the above reasoning now becomes dependent on l and m due to the anisotropic polarization induced by the PBCs. Finally, inserting eq 21 into eq 22 gives the relationwhere we have added the subscript “diel” to ε, to stress that it represents the true (m-independent) dielectric permittivity that the fluid would have if it behaved as an isotropic dielectric medium. For a molecular system, εPBC,( can be obtained by sampling ⟨(Q)2⟩PBC in a computer simulation. By plotting εPBC,( as a function of S(, the two constants εdiel and kself appearing in eq 23 can be determined from the intercept and slope of a linear fit to the data points.

Just as in the case of dielectric solvation, we may use the analogy of eq 13 to define the free energy change ΔAvac→PBC of bringing a dielectric sphere from a vacuum into a system under PBCs. Using eqs 8, 9, and 19−20 gives, after performing the integrations,

Since ΔAvac→PBC describes the long-range part of the electrostatic free energy of the simulated system, it should ideally not differ too much from ΔAvac→diel, and therefore a comparison between these two quantities may be a good way of assessing the accuracy of the particular boundary conditions being used. In particular, for S( = 0 or kself = 0, ΔAvac→PBC reduces to eq 14 for an isotropic dielectric medium.

Computational Details

The molecular model system is composed by particles possessing a pairwise additive interparticle potential v(r, ωi, ωj), composed of a dipolar and a Lennard-Jones (LJ) part according towhereandIn the above equations, μ represents the dipole of particle i, r is the vector pointing from particle i to particle j, r = |r|, and εLJ and σLJ are the LJ parameters. Two different values of the molecular dipole moment μ = |μ| were employed: μ = 0.45 atomic units (0.23813e Å, μ* ≡ μ/(4πε0εLJσLJ3)1/2 = 1.290) and μ = 0.65 atomic units (0.34397e Å, μ* = 1.863). The LJ parameters were set to σLJ = 2.8863 Å and εLJ = 1.97023 kJ mol–1.

The thermodynamic properties of the model system were determined by performing MD simulations in the canonical (constant N,V,T) ensemble, using N = 1000 particles in a cell of volume V = 2.601 × 104 Å3 for all three simulation cell geometries. The temperature was kept constant at T = 315.78 K (T* ≡ kBT/εLJ = 1.333). Toroidal boundaries for the noncubic simulation cells were applied according to the procedures devised by Smith,[32] whereas Ewald summation with tinfoil boundaries were implemented using the formulas due to Smith and Fincham.[36] A spherical cutoff in real space of rcut = 14 Å was used in conjunction with the Ewald screening parameter α = 3.2/rcut. The cutoff ncut in reciprocal space was set to 7, 10, and 9 for the cube, RD, and TO geometries, respectively, to yield a constant relative error in the k-space energy of ∼10–5. For all simulations, the integrated MC/MD/Brownian dynamics simulation package Molsim[45] was used. For further details about the simulation parameters, the reader is referred to our previous study.[28]

The multipole moments Q, 1 ≤ l ≤ 4, of a sphere with radius R were evaluated after every 100th time step. The contribution Q from a molecular dipole Q1 located at r = (r, Ω) to the total multipole moment Q = ∑Q was calculated according to[28]where all terms containing C with |m| > l should be excluded. For each sampled configuration, Q was calculated with each particle used as the origin, giving in total N sampled values of Q per configuration. In addition, reference values of εdiel for various R values were calculated using eq 12 from a simulation in a cubic simulation cell and N = 105 particles, i.e., 100 times as large as the primary systems. This large (compared to R) system size was used in order to ensure that the values of εdiel thus obtained are unaffected by the boundary. This system is further described in conjunction with our previous study.[28]

Results and Discussion

In Figure 3, numerically calculated values of the reduced interaction tensor R2S( for l = 2 and 3 are given for the SC, BCT, and BCC lattices. We note that R2S( is highly dependent on m, taking on both positive, corresponding to repulsive net interaction energies, and negative, corresponding to attractive net interactions, values. We also note that there is no obvious correlation between either the sign or the magnitude of R2S( obtained from the three different lattices. For l = 2, the SC lattice yields significantly (≈ 100%) higher absolute values of R2S( than the other two lattices; however, for l = 3, the situation is the opposite. This means that the magnitude (and sign) of the coupling of a certain multipole depends strongly on its symmetry in relation to the symmetry of the particular lattice where it resides. We furthermore note (results not shown) that S( = 0 for all m and all lattices, meaning that the total dipole–dipole interaction is zero. This is a well-known fact for all cubic lattices, as long as the lattice sum is carried out in spherical shells.[24−27]

Figure 3

Values of the reduced interaction tensor R2S( for (a) l = 2 and (b) l = 3 relevant for the three different simulation cell geometries. The values were obtained from eq 16 using a spherical cutoff of nmax = 50.

In Figure 4, the probability distribution PPBC(Q) for two different octupole components, obtained from a simulation of the μ = 0.45 system in RD geometry, is shown. Clearly, there is a significant difference between the widths of the two probability distributions, although both distributions follow the predicted Gaussian form very well. The width of PPBC(Q) should be compared to the corresponding values of S( given in Figure 3. Obviously, the narrower one of the two probability distributions corresponds to a repulsive value of S( (R2S32(R) = 0.88), whereas the wider distribution corresponds to an attractive net interaction (R2S32(I) = −0.66). Thus, there is at least qualitative reason in our assumption that the width of PPBC(Q), and thus the magnitude of εPBC, can be described by the lattice interaction tensor S(.

Figure 4

Probability distribution PPBC(Q) of two octupole components obtained from a simulation with μ = 0.45 in RD geometry. The values of the corresponding reduced interaction tensors are R2S32(R) = 0.88 and R2S32(I) = −0.66.

In order to quantitatively assess the dependence of PPBC(Q) on S(, we evaluated the former quantity in terms of the apparent (m-dependent) dielectric permittivity εPBC,(, defined through eq 22. In Figure 5, plots of εPBC,( versus R2S( obtained for both dipole strengths and all three simulation cell geometries are presented. Clearly, the proposed linear relationship between εPBC,( and R2S( is very well reproduced by the simulation data in all cases. Furthermore, the results obtained from the μ = 0.45 systems exhibit slopes that are essentially independent of geometry and l. The slopes of the μ = 0.65 data show a larger variation, although no systematic dependence on geometry and l is apparent. We also note the perhaps somewhat nonintuitive fact that the magnitude of the self-interaction (quantified through R2S() does not decay with increasing l, at least not for l ≤ 4. Furthermore, the self-interaction magnitudes do not show any clear trend between the different cell geometries. As an example, we note that the cubic geometry exhibits the largest quadrupole (l = 2) self-interactions, whereas the octupole (l = 3) self-interaction has its largest magnitude in the RD and TO geometries.

Figure 5

Apparent dielectric permittivity εPBC,( obtained from simulations of dipoles with (a) μ = 0.45 and (b) μ = 0.65 versus R2S( for the corresponding lattice types. The results for l = 3 (l = 4) have been shifted vertically by 5 (10) units for μ = 0.45 and 20 (40) units for μ = 0.65 to enhance readability.

Because of the good linearity of the data, eq 23 can be used to obtain values of εdiel and kself from the intercept and slope, respectively, of the data in Figure 5. In Table 1, fitted values of εdiel and kself are given for both dipole strengths and all three cell geometries, together with values of εdiel independently calculated from a simulation with N = 105 using eq 12 and the same values of the sampling radius R. From this data, we note thatObservation 1 shows that our assumption that the interaction between noncorrelated multipoles (i.e. excluding the self-interaction) can be described using formulas for an isotropic dielectric medium (eq 15) is reasonable, perhaps with the exception for the l = 2 and μ = 0.65 case. Observation 2 indicates that there may be a slight variation of kself with εdiel (or μ), although the source of this variation is not clear. The larger variation in kself for μ = 0.65 we attribute to the larger statistical noise present in the more strongly coupled system.

Table 1

Fitted Values (Figure 5 and eq 23) of kself and εdiel for Dipoles with μ = 0.45 and 0.65 and Values of εdiel Obtained (eq 12) from a Simulation in Cubic Geometry with N = 105 Particles

μ	geometry	R [Å]	l	kself	εdiel(fit)	εdiel(ref)
0.45	cube	14.80	2	0.67	13.1	13.4
			3	0.75	12.6	12.5
			4	0.83	11.5	11.7
	RD	16.60	2	0.66	13.2	13.7
			3	0.74	13.1	12.9
			4	0.83	11.8	12.1
	TO	16.16	2	0.67	13.2	13.6
			3	0.76	12.9	12.8
			4	0.89	11.8	12.0
0.65	cube	14.80	2	0.96	59	67
			3	0.85	52	51
			4	1.16	39	39
	RD	16.60	2	0.87	59	73
			3	0.92	59	56
			4	1.15	43	44
	TO	16.16	2	0.73	59	72
			3	1.01	56	55
			4	1.40	42	43

The fitted values of εdiel are close (within 5%) to the ones calculated from eq 12 for all fittings except those with μ = 0.65 and l = 2, where the fittings generally yield too low values of εdiel.

The fitted values of kself are slightly larger and exhibit larger variations for μ = 0.65 than for μ = 0.45.

We furthermore note that εdiel decreases with increasing l, in line with what we have observed before.[28,46] The apparent geometry dependence of εdiel is not due directly to the geometry but rather to the slightly different radii of the inscribed spheres in the three simulation cells; this behavior is also consistent with our previous observation[28,46] that εdiel increases with increasing sampling sphere radius R for a given l.

In Figure 6, the data corresponding to Figure 5a, but obtained using sampling radii half as large, are presented. In this case, it is obvious that the magnitude of the self-interaction quickly becomes less significant for increasing l, due to its R–(2 dependence. The linear fits are however still satisfactory, even though the very small variation in εPBC,( over the range of R2S( values leads to larger statistical errors in the fittings, especially for l = 3 and 4. The apparent shift in the y direction between curves obtained using different geometries is merely due to the values used for the sampling radius R being geometry dependent, leading to different values of the intrinsic dielectric permittivity εdiel. In fact, the same shift is present in Figure 5, although it is not visible due to the much wider range of the ordinate axis. According to our theoretical assumptions, the value of kself should be independent of R, meaning that the slope of the lines in Figures 5a and 6 should be identical. In Table 2, the values of εdiel and kself obtained from the smaller sampling radii are presented. Although the variation in kself is larger than in Table 1 due to the larger statistical noise, our assumption for kself is not obviously contradicted. Using the smaller sampling radii for the μ = 0.65 system (data not shown), however, seems to yield somewhat larger values of kself than in Table 1, although the statistical significance of these values can be questioned.

Figure 6

Data corresponding to Figure 5a but using sampling radii R half as large. Note the different scale on both axes and that the data have not been shifted in the y direction.

Table 2

Data As in Table 1 but Obtained Using Sampling Radii Half As Large and Only for μ = 0.45

geometry	R [Å]	l	kself	εdiel(fit)	εdiel(ref)
cube	7.40	2	0.94	10.8	10.8
		3	1.84	9.4	9.4
		4	1.34	8.4	8.4
RD	8.80	2	0.77	11.5	11.5
		3	0.79	10.3	10.3
		4	1.33	9.4	9.4
TO	8.08	2	0.86	11.2	11.2
		3	0.59	10.0	10.0
		4	0.06	9.1	9.1

Figure 7 gives the free energy ΔAvac→PBC for l = 3, μ = 0.65, and RD geometry calculated using eq 24. Clearly, the anisotropy in εPBC discussed in the previous paragraphs also corresponds to a large anisotropy in the solvation free energy of the subvolume. Quantitatively, βΔAvac→PBC varies between −1.6 and −0.5, compared to the dielectric value βΔAvac→diel ≈ −1.30 (eq 14, red solid line in Figure 7). However, the average of βΔAvac→PBC over all seven octupole components is βΔAvac→PBC(avg) ≈ −1.26 (black dashed line in Figure 7), i.e., very close to the dielectric value. Thus, the total solvation free energy (at least on the octupolar level) of the simulation box is very close to that of an isotropic system, even though it is distributed in a highly anisotropic way. A possible key to understanding this behavior is to be found in eq 18, namely, that the (unweighted) interaction tensors for any given l ≥ 1 cancel out when summed over all multipole components. Thus, the suppression of some fluctuation modes is exactly compensated by the enhancement of others, leading to a reasonable “mean value” of the energy. This behavior is reproduced for all multipole components and geometries (results not shown), in the sense that βΔAvac→PBC(avg) and βΔAvac→diel are always within 10% of each other.

Figure 7

Solvation free energy ΔAvac→PBC for l = 3 (black solid line) obtained from a μ = 0.65 system in RD geometry using eq 24 and values of kself and εdiel(fit) from Table 1. The red solid line gives ΔAvac→diel obtained from eq 14 using εdiel(sim) from Table 1, and the black dashed line gives the mean value of ΔAvac→PBC, averaged over all seven octupole components (black symbols, two doubly degenerate values).

Conclusions

In the present study, we have presented a quantitative analysis of the periodicity effects induced in a dipolar system by the use of PBCs. Using classical electrostatics and statistical thermodynamics, we developed a heuristic model relating the apparent, anisotropic dielectric permittivity εPBC,( to the reduced lattice interaction tensors S(. The theory exhibits excellent agreement with results from MD simulations of Stockmayer fluids with two different dipole strengths and three different simulation cell geometries. Although the anisotropy in the electrostatic fluctuations is independent of l on the length scale of the simulation box, it is shown that the “range” of the boundary effects (i.e., the minimum value of R needed to induce significant boundary effects) decreases strongly with increasing l. Furthermore, it was shown that the large (∼200%) anisotropy in the solvation free energy on the length scale of the simulation box disappears when averaged over all fluctuation modes, leading to the total solvation free energy being practically identical to the value predicted for an isotropic system.

Even though the simulation part of our study is based on a Stockmayer model system, we argue that our use of a dielectric continuum model as the theoretical basis means that the effects are fully transferable to any polar system which may be described as a dielectric medium, in particular the many popular water models used in molecular simulations. We also expect that any structural property, i.e., not only the dielectric permittivity, evaluated on the length-scale of the full simulation box is equally affected by the boundary effects.

One of the most important observations from this study is that the total solvation free energy is, in spite of the large anisotropy of the individual contributions, very close to the correct, isotropic value. This observation is indeed closely analogous to the corresponding averaging in the anisotropy of the radial distribution function for a Lennard-Jones fluid under PBCs observed by Pratt and Haan.[13] We argue that this property explains the success of Ewald summation and related techniques, at least when it comes to evaluating energies and relatively short-range structural properties. Nevertheless, as we have also shown previously,[28] one should use caution when evaluating structural properties on length scales larger than half the length of the simulation box.

Another relevant question is whether there is any rationale behind using a noncubic (RD or TO) simulation cell in order to reduce periodicity effects, as has been suggested previously.[29] First of all, it is clear that the periodicity effects when taking the full simulation cell into account are as strong for all three cell geometries (Figure 5), albeit not identical for a given l. We note, however, that the influence from the periodicity on the quadrupolar (l = 2) fluctuations is significantly lower in the RD and TO geometries than for the cube. Since the magnitude of the boundary effect for a given l decays as R–(2, this means that the “leading-order” contribution to the boundary effects is about 50% smaller (Figure 6) in the RD and TO geometries compared to when using a cubic simulation cell. On the other hand, which simulation cell to be used also depends on the specific system under study. For example, if one wants to simulate a macromolecule (e.g., a protein) with a particularly large molecular octupole moment, a cubic box may be the most appropriate one, due to its lower lattice coupling for the system octupole moment. Furthermore, it may also be advantageous, when using rotational constraints, to orient the axis of the largest electrostatic moment along the axis with the lowest value of S(; for example, orienting the octupole moment in the S33(R) or S33(I) direction gives a lattice interaction of less than 3% of the value obtained when the octupole is oriented along the S32(I) axis.

The present study provides a quantitative understanding of the isotropic and anisotropic parts of the solvation in a polar system under PBCs and puts them in relation to the behavior of an isotropic system. This understanding is essential for the possibility to remedy the periodicity effects, for example by imposing suitable bias functions in an MC simulation in order to remove the anisotropic self-interaction.