Local molecular field theory for the treatment of electrostatics
Abstract
We examine in detail the theoretical underpinnings of previous successful applications of local molecular field (LMF) theory to charged systems. LMF theory generally accounts for the averaged effects of longranged components of the intermolecular interactions by using an effective or restructured external field. The derivation starts from the exact YvonBorn Green hierarchy and shows that the approximation can be very accurate when the interactions averaged over are slowly varying at characteristic nearestneighbor distances. Application of LMF theory to Coulomb interactions alone allows for great simplifications of the governing equations. LMF theory then reduces to a single equation for a restructured electrostatic potential that satisfies Poisson’s equation defined with a smoothed charge density. Because of this charge smoothing by a Gaussian of width , this equation may be solved more simply than the detailed simulation geometry might suggest. Proper choice of the smoothing length plays a major role in ensuring the accuracy of this approximation. We examine the results of a basic confinement of water between corrugated wall and justify the simple LMF equation used in a previous publication. We further generalize these results to confinements that include fixed charges in order to demonstrate the broader impact of charge smoothing by . The slowlyvarying part of the restructured electrostatic potential will be more symmetric than the local details of confinements.
I Introduction
The treatment of longranged Coulomb interactions remains a substantial challenge in nearly all classical molecular simulations. The basic problem is that interactions from even very distant charges remain important, as indicated by the divergence of the integral of 1/r from any finite truncation distance to infinity. This causes problems in computer simulations of systems with charges, where contributions from distant periodic images of the simulation cell must be taken into account, typically by Ewald sum techniques. This also hampers the development of simple, intuitive local pictures and analytic theory when charges are present.
Here we discuss basic features of a new approach for charged systems, local molecular field (LMF) theory. LMF theory was originally applied to LennardJones (LJ) systems for both simulation and analytical work (1); (2); (3); (4); (5); (6), and the LMF approach has been used recently to analyze many charged systems, both uniform and nonuniform, with even greater success (7); (8); (9); (10); (11).
LMF theory can be generally characterized as a mapping that relates structural and thermodynamic properties of a nonuniform system with longranged and slowlyvarying intermolecular interactions in a given external potential energy field, representing, e.g., fixed solutes or confining walls, to those of a simpler “mimic system” with properly chosen shortranged intermolecular interactions in an effective or restructured field. The strong shortranged interactions in the mimic system are chosen to give an accurate description of the forces between typical nearestneighbor molecules in the full system, e.g., the repulsive molecular cores in a dense LJ fluid, local hydrogen bonds in water, or ion pairs in a strongly coupled ionic system. Because of the short range of the interactions in the mimic system, fast and efficient simulations scaling linearly with system size are possible. The restructured field is determined selfconsistently and is the sum of the bare external field in the original system and a densityweighted meanfield average of the remaining longranged and slowlyvarying parts of the intermolecular interactions. The restructured field corrects major errors that can arise, especially in nonuniform systems, from applying solely a spherical truncation of the intermolecular pair interactions, exemplified by shifted force truncations for nonuniform Coulomb systems (12); (13). More advanced spherical truncations of exist, such as sitesite reaction field techniques (14), but the difficulties in nonuniform situations remain the same.
The LMF approach offers a very general perspective and the averaging can be carried out for any slowlyvarying component of the intermolecular interactions. We show here that LMF theory takes an especially powerful and simple form, with strong analogies to classical electrostatics, when it is applied uniformly to the basic Coulomb interaction alone (and not, e.g., to LJ interactions as well). We can then take advantage of symmetries between charges of different signs and magnitudes, interacting with the same spherically symmetric potential, and determine the effective field using only the total charge density, and not separate number densities for each charged species or interaction site as would be required for a general mixture. We find that the contribution to the restructured field in LMF theory from the slowlyvarying parts of the Coulomb interactions can be exactly determined from a restructured electrostatic potential that satisfied Poisson’s equation, but with a Gaussiansmoothed charge distribution. A proper choice of the smoothing length in the Gaussian plays a key role in the accuracy of LMF theory, since it also determines the size of the local strongcoupling region within which strong forces from nearestneighbor sites must be accurately captured by the shortranged interactions in the mimic system.
Previous discussions of LMF theory for Coulomb or LJ systems relied on general arguments and did not exploit all the important simplifications that arise by applying LMF theory only to with the same potential splitting regardless of point charge details. Here we seek to justify and document LMF theory as completely as possible, accounting for Coulomb symmetries. We show how this leads to a simple, broadly applicable equation for a restructured electrostatic potential, and then we extend results from (11) to show that smooth LMF solutions can still arise from a general class of molecularlycorrugated confinement potentials.
Ii Derivation of the LMF Equation
The LMF equation may be derived for mixtures and standard sitesite molecular models. However the basic features of the derivation are identical regardless of whether we consider a mixture, a sitesite molecular fluid, or a singlecomponent system. Therefore, we present the derivation for a singlecomponent system in order to highlight the physical approximations in the derivation without the clutter of the multiple indices necessary in the derivations for mixtures or molecules. Furthermore, symmetries available to us when treating charge interactions will allow us to simplify the LMF equation for mixtures into a single electrostatic potential equation, as we will show later in section III.
ii.1 Motivation
We assume that the intermolecular interactions in the full system of interest are slowly varying at large separations. The first step in LMF theory is to divide into properly chosen short and longranged parts:
(1) 
such that contains all the the slowlyvarying longranged interactions in but remains slowly varying at characteristic nearestneighbor distances where there are strong forces between neighboring molecules. By construction, the shortranged component then captures those strong shortranged forces but vanishes quickly at larger separations.
LMF theory then defines a mapping from the full system with pair potential and an external singleparticle potential energy function to a mimic system defined by the shortranged pair interactions and an effective or restructured external potential energy function :
(2) 
We call the final combination a mimic system when and are properly chosen because we seek a system that can capture many relevant structural and thermodynamic properties of the full system. In addition, we seek this mapping to be closed, in the sense that the mapping should not require any information about the full system other than the original and . It seems plausible that some effective could account for averaged effects of the neglected longranged interactions in the mimic system, but we seek a more careful development of this intuition.
ii.2 Exact Starting Point for the LMF Derivation
The derivation of the LMF equation starts with an exact statistical mechanical equation involving structure and forces, the YvonBornGreen (YBG) hierarchy of equations (15); (16). In one form, the first equation in the hierarchy reads
(3) 
This equation relates the gradient of the (singlet) density profile of the system, , to the force from the external single particle potential, , and the averaged force between particles, weighted by the conditional singlet density, . This function is defined as the density at position given that a particle is at , and may be expressed mathematically as the ratio of the twoparticle density at and and the singleparticle density at :
(4) 
Thus, the YBG hierarchy expresses the rapidlyvarying singlet density in terms of the more complex and even more rapidlyvarying twoparticle density. There are two related difficulties in trying to extract practically useful information about the singlet density from (3). First, the density profile is related to more complicated rapidlyvarying functions. Second, there is no obvious manner to create a selfcontained loop to determine the singlet density profile, again due to its relation to the higherorder density profile.
Typical superposition closures of the YBG hierarchy attempt to resolve both difficulties by approximating , leading to an equation where the gradient over can be pulled out of (3) and is defined in terms of itself. However, even for a bulk fluid, this is a poor approximation for any moderately dense liquid (17); (18). The superposition approximation neglects any excluded volume or specificbinding interactions that alter the small nature of relative to . Replacing the conditional singlet density with the singlet density simplifies the equation but also effectively includes with large weight too many configurations where two particles are improbably close to each other, leading to great inaccuracies.
We assume that there exists a good choice of that captures those shortranged excludedvolume and specificbinding contributions. We will describe such a choice (7) for splitting using the smoothing length in section III. Given , we write the YBG equations for both the full sytem and the LMFmapped system, denoted by the subscript R, from (2) as:
(5) 
The density profiles and conditional singlet densities are written as functionals of and to emphasize the restructuring of the external potential energy in the LMF mapping.
As written with an arbitrary choice of , (5) simply displays formally exact but not particularly useful equations for two arbitrary and unconnected systems with different intermolecular interactions in different external fields. As in (2), we now connect these systems and achieve substantial simplifications and cancellations by requiring that be chosen so that
(6) 
It seems very plausible that such a exists, especially when captures the strong shortranged intermolecular forces. We then take the difference between the equations in (5). The terms involving the singlet densities on the lefthand sides cancel and the result can be exactly written in a form useful for further analysis as
(7a)  
(7b)  
(7c) 
Henceforth, for simplicity, we denote as .
Equation (7) formally holds for any choice of and , and because conditional densities still appear, it generally has most of the problems of the standard hierarchy equations noted after (4). Thus it may seem we have made little progress. But this is not the case for proper choice of and with the properties sketched above. This will allow considerable simplification of (7) and defines a mimic system that

has a wellchosen that mimics higher order correlations in the full system as well as the singlet density profile,

can be treated via statistical mechanics without any a priori knowledge of the full system aside from and , and

has an associated that depends only on singlet density profiles of the mimic system and not on the more complex conditional singlet densities.
ii.3 Approximations to Yield LMF Equation
Physicallymotivated approximations to the exact (7) will result in a final LMF equation below in (8) that has a simple meanfield form. The following discussion elucidates that the LMF equation is not a blind assertion of meanfield behavior but rather a controlled and accurate approximation, provided that we choose our mimic system carefully.
Recalling from (1) that , the LMF derivation focuses on making two connected and physicallyreasonable approximations based on choosing a shortranged that will induce the correct nearestneighbor structure, where the conditional singlet and singlet densities differ the most, and a corresponding that is slowlyvarying on that length scale. As we will see in section III.1, the form of the Coulomb potential grants us the freedom to fulfill both conditions by choosing a single smoothing length . When and are chosen correctly, we shall have a “mimic system.”
By making approximations to (7), the equality of densities expressed in (6) will only hold approximately, and we will find an equation based on the exact (7) that will also no longer be exact. But all will be accurately satisfied with proper choices of and . The general arguments below have been briefly sketched before (2), but the form of (7) allows us to make a more precise statement of the necessary approximations.
We note that (7a) depends only on the density response of the mimic system, and no information about the density response of the full system or any conditional density is needed. Thus if terms (7b) and (7c) could be neglected, then (7a) alone would define a closed equation for the mimic system. The right side of (7a) still transforms as a gradient and should display all the associated properties. In that sense, neglecting terms (7b) and (7c) together may well be more accurate than taking either one of the following approximations individually.
With proper choice of and , we may reasonably neglect the final two terms in (7) using the following arguments:

Term (7b) probes the difference between the conditional singlet densities for the full and mimic systems via convolution with . The integrand will be quickly forced to zero at larger by the vanishing gradient of the shortranged . Since both the full and mimic systems have the same strong shortranged core forces with an appropriately chosen , and these core forces should mainly determine the shortranged part of the conditional densities, it seems plausible that with proper choice of , term (7b) can be neglected. This is significantly better than any superposition approximation for the conditional densities and does not require us to numerically estimate either conditional density.

Term (7c) probes the difference between the conditional singlet density and the singlet density of the full system. As explained previously with relation to standard closures of the YBG equation, assuming their equality can be highly problematic at short distances. However, we are saved by the fact that this difference is paired with . Since has been chosen to encompass core interactions, will be simultaneously slowlyvarying over those nearestneighbor distances, so the associated force is essentially zero for exactly the range of distances where the conditional singlet density and singlet density differ significantly. Thus we expect for many liquids that the integrand in term (7c) may also be accurately approximated as zero.
These two approximations allow us to truncate the YBG hierarchy. No explicit knowledge of the conditional singlet density in either the full or the mimic system will be required, but we expect the conditional singlet density of the mimic system to closely resemble that of the full system at short range.
Notably, the second approximation will fail when there are longranged pair correlations arising from capillary waves at the liquidvapor interface or the divergence of the correlation length near the critical point. In those instances, the conditional singlet density will not be approximately the same as the singlet density well beyond nearestneighbor distances. However, we do expect this approximation to hold in stronglycoupled charged systems away from the critical point where there exist pair correlations that decay exponentially over those nearestneighbor distances.
With these two coupled approximations, only (7a) remains, and we can take the gradient with respect to outside of the integral. Then, integrating over yields the general LMF equation
(8) 
where the integration constant will be set by boundary conditions. The form of this equation highlights the connection with meanfield techniques. The function may be viewed as an averaging of the slowlyvarying, longranged portions of the pair potential over the singleparticle density . The hierarchy has been successfully truncated since is defined in terms of only . Provided sufficient computer time to simulate the shortranged mimic system and therefore close the selfconsistent loop between and , the LMF equation is exactly soluble.
This equation is not a naïve meanfield ansatz for since it arises from reasonable approximations to exact statistical mechanical equations. The crux of LMF theory, however, is choosing reasonable and such that the necessary approximations are accurate. The next section discusses this choice for .
Iii Electrostatics via LMF Theory
Here, we explain our choice of and for charge interactions. Then, we seek to highlight the symmetry and simplifications possible when applying the LMF equation to purely chargecharge interactions by developing the LMF equation for a single restructured electrostatic potential, rather than for separate potential energy functions for each atom type.
iii.1 Coulomb Potential Split
In CGS units, the pair interaction between two charges, and , is
(9) 
The factor is simply the dielectric constant of a uniform medium in which the charges are immersed. The details of charge magnitudes, signs, and choice of CGS or SI units will change from system to system, but the underlying functional form will remain constant. Thus, following (10), we introduce and to split this function as
(10) 
into shortranged and longranged components. Both and will yield either a potential energy or an electrostatic potential, with appropriate prefactors.
Work of previous researchers (7) revealed that a beneficial choice of is the functional form associated with the electrostatic potential due to a unit Gaussian charge distribution of width , defined as
(11) 
leading to a given by
(12) 
This is slowlyvarying in space over the smoothing length . As discussed in (7), this choice of also simultaneously isolates only the small wavevector components of the Coulomb potential. The Fourier transform of is
(13) 
and the Gaussian function attenuates any large components that would be rapidlyvarying in space. Thus the is particularly wellsuited for meanfield averaging.
The corresponding to our choice of will then be given as
(14) 
Therefore, is proportional to the electrostatic potential due to a point charge screened by a surrounding, neutralizing Gaussian charge distribution whose width also sets the scale for the smooth truncation of . This split is shown in figure 1 along with a demonstration of the impact of increasing the smoothing length on and . The form of is slowly varying over , and approximately represents a Coulomb core potential of range . Tuning to larger values increases the range of and simultaneously smooths the variation of over characteristic nearestneighbor interactions.
Since has no inherent length scale, the split between short and longranged by the smoothing length may be tailored to each situation. This length should not be viewed as a fitting parameter, but rather a consistency parameter. We expect that choosing too small will result in poor and rapidlyvarying results from full LMF theory as is changed since the meanfield average will be over rapidlyvarying interactions. However once reaches a that is related to characteristic nearestneighbor distances, the assumptions used in the previous section to derive the LMF equation from the YBG hierarchy become accurate, and essentially the same results should be found for larger values as well (19). Thus we expect that selfconsistently applied LMF theory with an appropriately chosen will yield a true mimic system. This requirement for reasonable core inclusion has been explored in previous papers (7); (8); (9).
For small , the force due to is quite similar to that of a bare point charge. By increasing , we increase the effective range of these essentially unscreened Coulomb interactions. Thus, can be thought of as a Gaussiantruncated (GT) or “short” Coulomb interaction. We refer to a model with genericallytruncated Coulomb interactions as a short model and with specifically used as a GT model. We showed in (11) that short GT water, where all the point charges in SPC/E water are replaced by with no account taken of the effective field, can give a surprisingly accurate description of all sitesite correlations for bulk water. But this same approach can fail spectacularly for electrostatic correlations in nonuniform systems unless corrected by the inclusion of the longranged forces via the LMFrestructured potential (11).
The functional form of suggests both Ewald sums (20) and Wolf sums (21). Ewald sums, however, seek to explicitly evaluate the forces from at each timestep using an exact sum over the fluctuating periodic images of the simulation cell, and thus have a very different philosophical approach. The analogy with Wolf sums is much more appropriate in that each results in a spherical truncation applied at each time step. However, the historical development of the Wolf summation technique relies on the observation that charges in ionic lattices are quickly neutralized by slightly shifted “mirror” lattices of neutralizing charges. In fluids, this is altered to assuming that charges within a sharp cutoff radius are exactly neutralized by a shell of charge at the cutoff radius, and the use of to damp the forces from is applied more as a numerical afterthought. In the context of LMF theory, we view such damping as much more crucial to the success of Wolf sums in reproducing the structure of uniform fluids than the imposition of strict neutrality within a sharp cutoff radius.
The GT also bears great similarity to sitesite reaction field (RF) approaches (14); (22); (23); (24), when the fluid is assumed to be conducting. The usual RF method generates a truncated from the electrostatic potential of a unit point charge surrounded by a uniform neutralizing spherical charge distribution. We discuss in the Appendix more details of the relation between the GT , and the usual RF as well as a smoother chargedclouds RF truncation (14). The choice of spheres as the neutralizing charge distribution in these RF methods leads to a potential that is exactly zero beyond a certain cutoff distance, but this induces a discontinuity in a higher order derivative at the sharp cutoff. The LMF choice of Gaussians leads to an infinitely smooth and soft truncation both in space and space. Moreover a GT model logically separates the choice of the smoothing length from the possible use for computational convenience of a sharp cutoff at a larger distance along with various local smoothing schemes that could be employed near the cutoff.
Despite the general similarities, the theoretical constructs of both reaction field techniques and Wolf sums allow for no obvious correction in nonuniform systems, such as the slab confinement to which we applied LMF theory in (11). As such, LMF theory is an important advance in the use of spherical truncations for .
iii.2 LMF Electrostatic Potential
Now we will derive a further simplified electrostatic form of the LMF equation. We need only consider the LMF mixture equation here because, with appropriate further approximations, the LMF equation for standard molecular models collapses onto the mixture LMF equation as will be shown in a paper dealing more broadly with sitesite molecular models (25). Using the YBG equation for mixtures, exactly the same approximations as in section II lead to the following LMF equation for a species in a mixture of different species indexed by and including , as originally given in (8):
(15) 
If the pair interactions being treated were of general form this would be the furthest simplification possible, with an LMF equation to be solved for each different species.
The observation in (8) that a single for all charged interactions proves most useful leads straightforwardly to a much more compact form. Using the LMF approach only for charges, each may be written as
(16) 
Thus (15) may be exactly written as
(17) 
Using the natural definition of charge density as
(18) 
we may reexpress the previous equation for as
(19) 
By noting as done in (8) that may be due to both electrostatic and nonelectrostatic components, we separate as
(20) 
where encompases nonelectrostatic confinements such as the smooth LJ walls used in (11) and represents the general external electrostatic potential.
This observation naturally leads to the following rewriting of each LMF equation for a species as
(21) 
with a single restructured electrostatic potential , defined as
(22) 
Thus all selfconsistent requirements on the diverse set of can be written in terms of a single restructured electrostatic potential defined by one LMF equation involving the equilibrium mobile charge density .
Given the convolution definition of as the potential due to a Gaussian charge density in (12), we may schematically represent the integral term in (22) as where indicates a convolution. Thus (22) can be exactly rewritten as
(23) 
where represents the smoothing of the charge density by in (11) via the convolution
(24) 
In almost all cases of interest we can picture as arising from a fixed, external charge distribution , and we should apply the same separation of the Coulomb potential to the fixed as well as mobile charges. This will allow us to express (23) in a suggestive and compact form using the sum of the fixed charge density and the equilibrated mobile charge density profile, . Thus we write
(25) 
Just as is the electrostatic potential due to a point charge convolved with a Gaussian of width , so is the electrostatic potential due to the convolution of the fixed external charge density with that same Gaussian as defined in (11) and used in (12) and (24). Given the definition of and , we write analogously
(26) 
thus defining , the slowlyvarying component of the restructured electrostatic potential.
Using (23), this slowlyvarying component is exactly given by
(27) 
We see that is the electrostatic potential arising from the Gaussian smoothing of the total charge density, .
Equation (27) appears as the standard definition of the electrostatic potential due to charge density from firstyear textbooks like (26). Thus the LMF equation for the slowlyvarying component of the restructured electrostatic potential may be viewed as simply doing classical electrostatics on an equilibrium, smoothed total charge density and defining a single external restructured electrostatic potential. In fact we may in general represent the slowlyvarying LMF equation (27) as a modified Poisson’s equation,
(28) 
based on the smoothed total charge density. This point was first made for the uniformlychargedwall model system (27), and it is immediately generalizable to the LMF equation dealt with here. The existence of this Poissonlike equation dictating a modified electrostatics emphasizes that, in order for LMF theory to be valid for electrostatics, all chargecharge interactions, whether fixed, bound, or mobile, should be mapped via LMF theory to shortranged interactions in the restructured electrostatic potential . In (11), we emphasized the value of this modifiedelectrostatics perspective in physically interpreting molecular simulations and in understanding why and when spherical Gaussian truncations, using alone, fail in nonuniform environments.
Iv Further Simplifications Due to Charge Smoothing
Understanding the LMF equation as standard electrostatics for a Gaussiansmoothed charge density also allows for simpler solutions of the LMF equation than might be initially expected. As an example, in (11), we were able to use a onedimensional LMF equation to treat a system that has a charge density profile with an explicitly threedimensional character. In this section, we give the reasoning behind such simplifications and argue that they should be expected generally in many physically relevant cases.
iv.1 Simple Corrugated Surface
In (11), we treated SPC/E water confined to a slab by two empirical Pt(111) surfaces (28) that order the water molecules at the surface, attracting the oxygen atoms to localized binding sites. Details of the simulation methods are available in (11). In brief, we simulated 1054 SPC/E molecules confined between the two walls with 500 ps of equilibration and 1.5 ns of subsequent data collection. Simulations using slabcorrected threedimensional Ewald sums, LMF theory with Å, and LMF theory with Å were carried out. During those runs, the timestep chosen was 1.0 fs. However, for selfconsistent solution of the LMF equation, we first simulated sets of 10 shorter runs with 25 ps of equilibration and 50 ps of accumulation with a timestep of 2.5 fs. These 10 shorter runs began with distinct initial conditions; we observed that averaging the resulting from such a set of simulations better isolated the system’s equilibrium response to changes in from the natural fluctuations in the charge density. For the larger , convergence was achieved within 3 iterations, and for the smaller , 10 iterations were required for selfconsistency. This approach certainly does not represent a numerically optimized solution of the LMF equation and is used here only to demonstrate the accuracy possible in principle when the full LMF theory is used. We will address the need for a more comprehensive and numerically efficient approach to the solution of the LMF equation and general timing issues further in section V.
In contrast to the smooth hydrophobic surface also dealt with in that paper, the empirical (111) surface was intended to represent detailed and specific ordering of water molecules at an atomicallycorrugated surface. A schematic of the rectangular unit cell of a (111) surface is shown in figure 2, demonstrating the spacing of the three topmost layers of atoms. The empirical Pt(111) potential (28) that we use does not have explicit atoms or charges. Rather the potential simply encompasses the presence of binding sites through exponential functions and various cosine functions to represent the periodicity of the surface.
As expected from the form of the potential, the charge density depends on and in addition to . In figure 2, the charge density in the 1.5 Å layer nearest the (111)wall is projected onto a unit cell. As previously demonstrated (28), the surface induces a charge density profile very nonuniform in the  and directions. The distinct periodicity is due to the fact that the spacing between binding sites (2.77 Å) is nearly commensurate with the first peak in the bulk water ( 2.8 Å). However despite a charge density that was explicitly a function of , we obtained accurate results for system properties using (11). We posited that this was the case due to the charge smoothing described in the previous section.
Here we seek to demonstrate this more quantitatively by probing the charge density profile away from the surface within hexagonal prisms defined by three distinct hexagonal binding regions – A, B, and C – as labeled in figure 2. Figure 3(a) shows the measured along each hexagonal prism and averaged laterally over the hexagonal area. These density profiles are quite distinct both from each other and from the full laterallysmoothed .
The Gaussian smoothing of charge density inherent in the LMF treatment allows us to distinguish the coherent and incoherent variations in charge density seen in figure 3(a). As shown in figure 3(b), convoluting with Gaussians in the  and directions and subsequently analyzing sites A, B, and C leads to a density, , that is indistinguishable from the calculated initially as only a function of . This equivalence is due to that fact that the corrugations in effect induce local random noise in the density profile that is therefore averaged out by the Gaussians of width greater than the corrugation. In contrast, the presence of the confining surface leads to a density profile that does not average out in the direction as seen in figure 3c. In essence, the LMF smoothing approach shows when variations in charge density profiles cancel, and when they do not.
iv.2 Simplifications Still Hold For Surfaces with Fixed Charges
However, a reasonable criticism of the previous surface is that it was not atomically realistic, even on a classical level. More realistic metallic surfaces would require a treatment of image charges, which we will not address here. However, molecular surfaces modeled with standard force fields would at least have charges at atomic sites near the surface. If the charges are free to move throughout the simulation as would generally be the case for biological membranes or liquidliquid interfaces, absolutely no change in tactic would be necessary. The charges associated with those molecules would also be included in the mobile .
If, instead, the charges are held fixed, as is more likely for solid surfaces in classical simulations of solidliquid interfaces, then . In such cases, we may no longer approximation as , but a similar simplification for the smooth will hold, as might be expected. Such a potential for silica surfaces (29) is currently being employed in our group to study the silicaacetonitrile interface (30).
For surfaces with fixed charges, the imposed external electrostatic potential would be
(29) 
where the prime indicates that these surface charges will need to be accounted for in the surface extending to infinity in the  and directions. The pair interactions representing the excluded volume of these atomic sites would be included in . The external electrostatic potential may be split into shortranged component and longranged component as
(30)  
(31) 
Now is a simple minimum image sum over shortranged pair interactions. The longranged component still must include the effect of surface sites extending out to infinity in the  and directions, but since may be interpreted as the electrostatic potential due to a Gaussian smoothing of the fixed charge sites, we expect that
(32) 
The sole requirement would be that should be larger than relevant lateral spacing of sites on the surface, such that the charge density will be sufficiently smoothed in the lateral directions. In such a case, is a function of , but all dependence on and is contained in , and we would expect the selfconsistent to hold to a very good approximation.
As a quick demonstration of the validity of these ideas, we construct a (111)surface with lattice constant Å as for the model (111)surface discussed in section IV.1. However, we arbitrarily assign charges of e for the first layer of atoms, e for the second layer of atoms, and e for the third layer of atoms.
The corresponding would be determined by smoothing out the charges in each layer laterally, such that the charged (111)surface is now approximated by three uniformly charged planes. With this representation, we may exactly express as the sum of three analytical functions,
(33) 
where is the surface charge density of each layer, and is the onedimensional Gaussiansmoothed Green’s function for a planar charge distribution (8); (9) analogous to in Eq. (12) for a point charge.
is plotted in figure 4 and is an important contributor to the forces near the surface. However, since the surface is net neutral, it decays essentially to zero within a distance of approximately from the surface.
We may calculate both and easily by constructing various configurations of one unit test charge above this (111)surface in the dlpoly2.16 molecular dynamics package. We exactly calculate to within the simulation energy precision, by carrying out single point energy calculations for a e test charge above various sites interacting with each particle in the surface via . This particle is a test charge in the sense that it measures the generated by the charges in the surface without disturbing the organization of charges within the surface. Thus the LMF estimation of the total is
(34) 
with given analytically by (33).
In order to determine independently, we employ slabcorrected threedimensional Ewald sums as a way to estimate that function (31). Since Ewald sums require net neutrality, we construct two surfaces with test charges – the surface we are interested in with a test charge above it and a second mirror surface with opposite charges with a test charge above it. The two (111)surfaces have lateral size 30.49 Å 28.806 Å with Å, and we use Å and for the corrected threedimensional Ewald sums. We separate the two surfaces within the simulation cell as much as possible in order to decouple each surface and testcharge pair from the other pair.
The energies yielded from these singlepointenergy calculations allow us to estimate , provided that we subtract the contribution from test chargetest charge interactions and the constant contribution from intrawall particle interactions () as follows,
(35) 
This is approximate since the Ewald sums will propagate the test charges laterally into all the periodic image cells, but given that the lateral area spans repeats of the minimum rectangular (111)surface unit cell of 2.77 Å 4.80 Å, interactions between each surface and the lateral periodic images of the nearest test charge should be minimal.
Shown in figure 5 is for fixed at Å as treated by (a) Ewald sums, (b) the LMF approximation with Å, and (c) the LMF approximation with Å. The similarity of the electrostatic potential surface for each to the full Ewald calculation demonstrates that 4.5 Å is greater than for the fixed surface at least. This stands to reason since, within a plane, particles are within 2.77 Å of each other and the planes of particles are separated by about 2.26 Å. Therefore Å is greater than relevant particle spacings. The agreement between and is expected but numerically nontrivial since at this distance above the surface. Therefore, achieving this degree of accuracy is a result of the accurate addition of two largemagnitude terms, and , to yield a much smaller magnitude .
In order to demonstrate this agreement more broadly, we again study sites A, B, and C as shown in figure 2. Figure 6 shows the difference
(36) 
laterally averaged above each hexagonal site. This variation is relatively small for a of either 4.5 Å or 6.0 Å since either is greater than any relevant spacings in the surface. Compared to the huge discrepancies that would arise from using solely and neglecting the several Volt contribution of , these differences are miniscule. At larger , there are slightly larger variations, which we believe are numerical artifacts arising from our use of a sharp cutoff radius (10.5 Å for Å and 13.5 Å for Å) in . Given the small scale on these graphs, even very tiny errors will visibly accumulate.
We fully expect that we should be able to treat as
(37) 
In fact, we expect this approximation for to be even better than for since the mobile charge density is included and will act to neutralize the fixed charged sites to an extent. The validity of approximating the Gaussiansmoothed equilibrium charge density profile as solely a function of for nonelectrostatic corrugated surface confinement discussed in section IV.1 should also hold for the mobile charge density profile in this model system with fixed point charges. Thus we conclude that this system should be treatable by LMF theory with the same onedimensional equation as used in (11) with the sole difference that is a function of .
V Conclusions and Outlook
Based on the success of LMF theory in treating a variety of charged and molecular systems (7); (8); (9); (10); (11); (25), we believe that LMF theory is a promising new approach that permits simple minimumimage simulations of interactions while still accurately accounting for the net additive longranged forces using the restructured external electrostatic potential . While the systems treated thus far have been simple ionic and molecular systems, the LMF framework is much more broadly applicable to simulation systems of biomolecular and experimental interest. Current research using the LMF approach in the Weeks research group includes the collapse of model charged polymers (32) and the behavior of acetonitrile near silica surfaces (30).
Here we have sought to present the statistical mechanical foundations of the LMF treatment of charged systems. We have highlighted the necessary but physically reasonable approximations leading to the LMF equation, and we have also derived the simple LMF equation for the electrostatic potential valid when treating pointcharge models, both ionic and molecular.
Furthermore, we have examined the useful consequences of understanding the smoothlyvarying LMF electrostatic potential as resulting from a Gaussiansmoothed charge density. In a previous paper (11), we showed that this charge density can be physically enlightening. Here we instead show the practical utility of this smoothing. The convolution of the charge density with a Gaussian of width allows us to solve much simpler LMF equations than might originally be believed. For confinement into a slab geometry, regardless of whether the confinement has nonelectrostatic or electrostatic corrugations, the relevant LMF equation will depend only on to a very good approximation.
In future work, we will examine the application of LMF theory to bulk sitesite fluids in much more detail (25) and also demonstrate very simple corrections that allow much more accurate thermodynamics from GT models (33). Another area which we hope to explore in greater detail is methods for optimally solving the LMF equation during simulation. As alluded to in (11), the writing of the LMF equation as a Poisson equation with a Gaussian smoothed charge density suggests that connection with work on efficient Poisson solvers would be numerically helpful, and that simple approximations may prove surprisingly accurate (10). Furthermore, since Poisson’s equation may be written as the minimum of a functional (34), using a modified CarParinello scheme (35) to slowly evolve to the averaged, equilibrium solution of the LMF equation seems like another fruitful path to pursue. In our (unoptimized) simulations of LMF systems, we observed a speedup of at least a factor of four relative to Ewald summation. However, such timing questions need to be studied in much greater detail once an optimized path to LMF solution is developed.
Regardless of such possibilities for solution of the LMF equation, LMF theory is still an inherently equilibrium theory. Given the success of a scaling principle of dynamics that states that equilibrium, dynamic properties will be wellcaptured in a reference system that reproduced static, structural properties (36), many dynamical properties for equilibrium systems could possibly be modeled quite well in our LMFmapped systems. But such a conjecture needs to be tested, and, regardless, much further theoretical work is required to develop the analog of LMF theory for dynamicallyevolving systems.
Acknowledgments
This work was supported by NSF grant # CHE0517818. JMR acknowledges the support of the University of Maryland Chemical Physics fellowship. The computations were supported in part by the University of Maryland.
Appendix: Relation between the GT and Reaction Field Truncations
The standard shortranged RF pair interaction is proportional to the potential energy between a unit point charge at and a unit point charge at surrounded by a uniform neutralizing spherical charge distribution terminating at the sharp cutoff radius of the RF model. Anticipating that we may want a more symmetric description of the charge distributions at and , we can also relate the RF to the difference between the potential energy of two unit point charges at and and the potential energy of a unit uniform spherical charge distribution at and a unit point charge at . But the subtracted term is not symmetric and the smoother chargedclouds interaction (14) can be found by subtracting the potential energy of two unit uniform spherical charge distributions at and (with cutoff radius ) from the potential energy of two unit point charges at and . The functional form was previously arrived at via a less symmetric argument (14). The functional form is different than that of .
In contrast, the chosen for LMF theory can be related to either the energy between a unit point charge at and a unit point charge at surrounded by a neutralizing Gaussian charge distribution of width or equivalently to the difference between the potential energy of two unit point charges at and and the potential energy of two unit Gaussian charge distributions of width at and (7). The same functional form results in either case.
References
 Weeks J D, Vollmayr K and Katsov K, 1997. Physica A 244 461
 Weeks J D, Katsov K and Vollmayr K, 1998. Phys. Rev. Lett. 81 4400
 VollmayrLee K, Katsov K and Weeks J D, 2001. J. Chem. Phys. 114 416
 Katsov K and Weeks J D, 2001. Phys. Rev. Lett. 86 440
 Katsov K and Weeks J D, 2001. J. Phys. Chem. B 105 6738
 Weeks J D, 2002. Annu. Rev. Phys. Chem. 53 533
 Chen Y G, Kaur C and Weeks J D, 2004. J. Phys. Chem. B 108 19874
 Chen Y G and Weeks J D, 2006. Proc. Nat. Acad. Sci. USA 103 7560
 Rodgers J M, Kaur C, Chen Y G and Weeks J D, 2006. Phys. Rev. Lett. 97 097801
 Denesyuk N A and Weeks J D, 2008. J. Chem. Phys. 128 124109
 Rodgers J M and Weeks J D, submitted. Proc. Nat. Acad. Sci. USA
 Feller S E, Pastor R W, Rojnuckarin A, Bogusz S and Brooks B R, 1996. J. Phys. Chem. 100 17011
 Spohr E, 1997. J. Chem. Phys. 107 6342
 Hummer G, Soumpasis D and Neumann M, 1994. J. Phys. – Condens. Matt. 6 A141
 McQuarrie D A, 2000. Statistical Mechanics (Sausalito, California: University Science Books)
 Hansen J P and McDonald I R, 2006. Theory of Simple Liquids (New York: Academic Press), 3rd ed.
 Krumhansl J and Wang S, 1972. J. Chem. Phys. 56 2179
 Wang S and Krumhansl J, 1972. J. Chem. Phys. 56 4287
 Since we have not bounded errors for the approximations in developing LMF theory, this is an approximate lower bound on , i.e. .
 Ewald P P, 1921. Ann. Phys. (Leipzig) 64 253
 Wolf D, Keblinski P, Phillpot S and Eggebrecht J, 1999. J. Chem. Phys. 110 82541
 Neumann M, 1983. Mol. Phys. 50 841
 Hummer G, Soumpasis D and Neumann M, 1992. Mol. Phys. 77 769
 Hummer G and Soumpasis D M, 1994. Phys. Rev. E 49 591
 Hu Z, Rodgers J M and Weeks J D, in preparation
 Purcell E M, 1985. Electricity and Magnetism, vol. 2 of Berkeley Physics Course (New York: McGrawHill), 2nd ed.
 Chen Y G, 2004. General Theory of Nonuniform Fluids: From Hard Spheres to Ionic Fluids. Ph.D. thesis, University of Maryland, College Park
 Raghavan K, Foster K, Motakabbir K and Berkowitz M, 1991. J. Chem. Phys. 94 2110
 Giovambattista N, Rossky P and Debenedetti P, 2006. Phys. Rev. E 73 041604
 Hu Z and Weeks J D, unpublished
 Yeh I and Berkowitz M, 1999. J. Chem. Phys. 111 3155
 Denesyuk N A and Weeks J D, unpublished
 Rodgers J M and Weeks J D, in preparation
 Jackson J D, 1999. Classical Electrodynamics (Hoboken, New Jersey: John Wiley and Sons, Inc.), 3rd ed.
 Frenkel D and Smit B, 2002. Understanding Molecular Simulation: From Algorithms to Applications (New York: Academic Press), 2nd ed.
 Young T and Andersen HC, 2005. J. Phys. Chem. B 109 2985