Interplay of local hydrogen-bonding and long-ranged dipolar forces in simulations of confined water

Bulk Gaussian-Truncated Water

In LMF theory, we split the 1/r potential associated with each∼charge into rapidly varying (v₀) and slowly varying (v₁) components as

This potential separation isolates strong short-ranged and rapidly varying interactions in v₀(r), and the remaining slowly varying long-ranged forces are encompassed by v₁(r). This interaction v₁(r) is proportional to the electrostatic potential arising from a smooth normalized Gaussian charge distribution with width σ, and is defined via the convolution

Thus, v₁(r) is slowly varying in r-space over the smoothing length σ and is simultaneously localized to small wave vectors in reciprocal space, as can be seen from its Fourier transform v̂₁(k) = 4πk^-2exp[-(kσ)²/4].

Hence, v₀(r) ≡ 1/r - v₁(r) is the screened potential resulting from a point charge surrounded by a neutralizing Gaussian charge distribution whose width σ also sets the scale for the smooth truncation of v₀. For small r < σ, the force due to v₀(r) approaches that of the bare point charge. By increasing σ we increase the effective range of essentially unscreened Coulomb interactions included in v₀(r).

LMF theory introduces a general mapping from the original system with full Coulomb interactions and an applied electrostatic potential V(r) arising from fixed charges or an applied electric field to a well-chosen “mimic system” with Coulomb interactions replaced by the short-ranged v₀(r):

Here, V_R(r) is a restructured electrostatic potential in the mimic system that accounts for the averaged effects of the neglected long-ranged forces represented by v₁(r), in essence by solving a modified Poisson's equation as explained in a later section.

In this article, we simulate the SPC/E model of water as shown in Fig. 1A. All data are results from molecular dynamics simulations as described in Methods. The point charges in SPC/E water play a dual role, generating both strong short-ranged forces that lead to the local hydrogen-bond network and the longer-ranged dipole–dipole interactions. Neither the point charges nor the Lennard–Jones (LJ) potential describing the molecular cores are perfect representations of actual water interactions, but adjusting the constants makes the total potential quite reasonable for temperatures and densities of interest.

Because the charges generating the local hydrogen-bonding network in water are very strongly coupled, a reasonable first approximation is simply to replace all point–charge interactions by the truncated v₀(r) without employing a V_R. LMF theory suggests that Coulomb interactions are best truncated on a site–site basis, because it is the basic 1/r potential that we attenuate about each site. Such a site–site scheme has been implemented in other spherical truncations of 1/r interactions as well (10). When the v₀(r) truncations are made, we have Gaussian-truncated (GT) water.

Crucial to the success of the Gaussian-truncation scheme is the choice of a σ that will be large enough to accurately describe the nearest-neighbor hydrogen bonds. By using a σ of 0.4 nm or larger, GT water can give a remarkably accurate description of all O and H pair correlation functions in bulk water, as shown in Fig. 1B. The difference between g_OO from a full electrostatic treatment (indicated by the labels “Full” in the figures) and g_OO from GT water is shown in Fig. 1C to emphasize the good agreement. We expect that certain earlier truncation schemes (10) with very slowly varying v₁(r) could give comparably accurate results.

The local hydrogen-bond network is often viewed as the most important qualitative feature of water (24), and GT water captures that feature with a very good description of nearest-neighbor hydrogen-bond energetics and local tetrahedral order, while ignoring the longer-ranged electrostatic effects of the bound charges. As argued elsewhere (Z. Hu, J.M.R., and J.D.W., unpublished data), the accuracy of GT water in bulk arises from a strong cancellation of the long-ranged electrostatic forces exerted by surrounding molecules in the uniform bulk on any given molecular site. These site-based spherical truncations can also yield highly accurate angular correlation functions, a property that is missed with molecule-based truncations that do not neglect only slowly varying forces (8). LMF theory also suggests simple analytic corrections to the energy and pressure such that the thermodynamics of bulk GT water are accurate as well (J.M.R. and J.D.W., unpublished data).

LMF-Corrected GT Water in an Applied Field

Following previous work (25, 26), we confine SPC/E water between smoothed hydrophobic LJ walls and apply a constant electric field E₀ normal to the walls, as sketched in Fig 2A. The confined water should behave like a dielectric slab from introductory electrostatics, resulting in a polarization field E_pol opposing the applied electric field (27). The oxygen density profile and charge density profile using only GT water are completely incorrect when we apply an electric field of 10.0 V/nm as shown in Fig. 2B and C.

Use of spherical truncations alone neglects the long-ranged electrostatic consequences of the local ordering of the water molecules, which act to attenuate and weaken the applied electric field in the central region as suggested by Le Chatelier's principle. Thus, the GT water molecules spectacularly overorder with the unattenuated electric field because the externally applied field induces an incorrect, strongly polarized state throughout the GT water slab that happens to be compatible with an energetically favored arrangement of hydrogen bonding. Inclusion of a self-consistent V_R(z) to represent the long-ranged net effect of water molecule polarization makes the density profiles again agree very well with the results for the full system as shown in Fig. 2B and C. With the inclusion of net long-ranged effects, the highly ordered state does not occur because it is entropically disfavored and thus is highly unlikely in the true attenuated field.

The inclusion of a restructured V_R(z) also corrects many other properties of confined water when an electric field is applied such as the dielectric constants and details of molecular orientation profiles (J.M.R. and J.D.W., unpublished data). However, the density profile is the simplest and most spectacular example of the failure of GT water in this system, and LMF theory handily corrects that failure.

LMF Theory for Electrostatics

LMF theory has been applied to many different ionic systems, including charge mixtures (23, 28–30). General derivations of the approach are available in the literature (20–22). However, when the LMF truncation is applied to the general 1/r functional form and the same σ is applied to every charge–charge interaction in the system, regardless of identity, we may write a substantially simpler formulation of the LMF mixture equation than previously stated (28). An explanation of these symmetries and the validity of the LMF equation below is given in reference 31 for charge mixtures and can be extended to site-to-site molecular models as well (Z. Hu, J.M.R., and J.D.W., unpublished data).

Provided that the σ chosen is at least large enough for v₀(r) to incorporate the nearest-neighbor interactions such as hydrogen bonding in water and for v₁(r) to be slowly varying over those nearby correlations, LMF theory predicts that the restructured potential V_R is accurately given by

The notation ρ_R^q(r;[V_R]) simply indicates that we now use the charge density profile of the mimic system (indicated by R) in the presence of the restructured potential V_R. The charge density ρ_R^q encompasses the effect of v₀(r) and V_R as well as all other nonelectrostatic interactions like the hydrophobic walls and the LJ cores in SPC/E water. We use the constant C to set the zero of our potential. In the systems examined in this article, this formula may be rewritten as a 1-dimensional equation as detailed in ref. 28.

V_R can be understood as a mean-field electrostatic potential representing the long-ranged forces that may be reasonably averaged over. Alternatively, harnessing the convolution definition of v₁(r) in Eq. 2, we may rewrite the LMF equation as the convolution of a smoothed charge density with the full Coulomb potential. Thus, V_R exactly satisfies a modified Poisson's equation

defined with a Gaussian-smoothed charge density

This again emphasizes that all charge interactions in the system must be identically truncated for this modified electrostatics to be valid.

The form of the short-ranged potential v₀(r) might suggest an analogy to Ewald summation. However, LMF theory focuses on the short-ranged mimic system and how best to use and correct its properties based on results from equilibrium theory, whereas Ewald summation focuses on the full system with exact periodic boundary conditions. As a result, V_R is a static external potential based on equilibrium properties that is determined self-consistently as described in Methods, whereas Ewald summation evaluates fluctuating forces arising from both the short- and long-ranged components of the Coulomb interactions at each time step.

Classical Failure of GT Water Explained

Although we believe that the failure above of GT water alone in an electric field has not been noted before, problems with electrostatic properties of confined truncated water in zero field are well known (13). Shown in Fig. 3 is the classic example of the failure of GT water in nonuniform situations—the incorrect induced potential profile of water molecules. SPC/E water is again confined between the hydrophobic walls, now with no applied electric field. As well-established by Rossky and coworkers (25), water molecules next to the wall tend to point one hydroxyl group into the wall so that one hydrogen bond per molecule on average is broken rather than two. This induces a dipole layer near each hydrophobic surface with the dipole pointing toward the walls. This atomically detailed dipole layer leads to a semidiscrete jump at the surface in the electrostatic potential defined as

This potential Φ_pol should plateau in the central bulk region; however, the use of GT water fails dramatically in this respect. Because of the initial overorientation of water molecules at the surface, driven by the maintenance of hydrogen bonds in GT water, the polarization potential never plateaus. This failure is well documented (13, 14) and is attributable to the neglect of net-additive long-ranged forces due to the formation of dipole layers at the surface. As also shown in Fig. 3, LMF theory correctly includes those neglected forces with tremendous success.

The form of V_R(z) shown in Fig. 4 highlights its physical effect on GT water molecules. Without V_R, the truncated molecules overorient near the wall. The restructured electrostatic potential V_R from LMF theory reflects the net long-ranged forces due to the ordering of water near the surfaces by applying a reorienting torque on surface water molecules. Because V_R is defined via a self-consistent equation, it represents a statistical mechanical balance between the favorability of maintaining local hydrogen bonds and the electrostatic penalty for creating an overly severe dipole layer. The potential V_R varies smoothly with z because it is determined by averaging over the slowly varying v₁(r) in Eq. 4. This suggests that a simple theory yielding this general form may be possible.

Although the physics of this failure was convincingly established by Spohr (14), its structural roots in the charge density profile remained mysterious. Brooks and coworkers (13) observed an overorientation of molecular dipoles near the walls in spherically truncated water. But as shown in Fig. 5A, the charge density profile for GT water tracks quite closely with the results from full Ewald summation techniques for SPC/E water. A dipole layer near the wall seems quite clear, and there are only slight deviations in the peaks next to the walls. However, the polarization potential Φ_pol belies the apparent accuracy of the atomic-level charge density profiles for GT water. The polarization potential may be fit accurately by a quadratic function in the central region between the walls. This fit suggests the existence of a constant residual charge density of approximately —0.04 e₀/nm³ throughout the “bulk” region neutralizing the overaccumulation of positive charge next to each wall.

LMF theory allows us to highlight exactly that relevant structural problem. The writing of the LMF equation as a modified Poisson's equation in Eq. 5 suggests that the charge density that best reflects the long-ranged electrostatics is not ρ^q, but rather the Gaussian-smoothed charge density ρ^qσ. The smoothed charge density ρ^qσ in Fig. 5B clearly shows the dipolar layer at each wall for the water confined with no electric field. Beyond that, GT water alone very obviously spreads the dipole layer into the bulk. The residual negative charge density is approximately the value suggested by the form of Φ_pol for GT water. Because the smoothed charge density is 2 orders of magnitude smaller than the atomic level charge density, this underlying electrostatics is buried alongside the slight simulation noise in ρ^q. On Gaussian smoothing, the random simulation noise in ρ^q cancels, uncovering the smaller but coherent electrostatic features. In a sense, the strong charge layering over 0.2 nm displayed in Fig. 5A is not the long-ranged dipole near the wall. The true dipole layer is due to average behavior over a few molecular layers. Classical electrostatic textbooks often posit some sort of smoothing over such local order (27).

This smoothed charge density allows for a more direct analogy between molecular simulations and classical electrostatics. The LMF approach naturally distinguishes the long-ranged charge–charge interactions representing overall dipolar ordering from the short-ranged charge–charge interactions in simulations that capture hydrogen-bonding and short-ranged polar attraction. Given that any σ greater than some σ_min will describe relevant short-ranged interactions like hydrogen bonding and can be used in simulations, too much significance should not be assigned to the particular chosen σ value or the exact magnitudes of the resulting ρ^qσ. However, we strengthen a previous suggestion (32) that charge density profiles reveal electrostatic effects in a less-biased fashion than molecule-based profiles: perhaps a smoothed profile is even less biased.

LMF Theory Is More Broadly Applicable

We now illustrate the application of LMF theory to moderately more realistic molecular surfaces by replacing the hydrophobic walls with empirical Pt(111) surfaces that order the water molecules at the surface, attracting the oxygen atoms to localized binding sites as shown by Berkowitz and coworkers (33). This surface is meant to represent specific interactions and detailed ordering of water molecules at an atomically corrugated surface, thus mimicking a hydrophilic interaction. The particular wall potential we use does not have explicitly fixed wall charges, but LMF theory can readily handle these as well (31).

Because this surface induces a charge density that is explicitly a function of x, y, and z, one might reasonably conclude that we now must solve the full 3-dimensional LMF equation given in Eq. 4 to obtain a V_R that is a function of r. However, the charge smoothing inherent in the LMF equation justifies a much simpler treatment of V_R. The lateral variations in charge density are washed out when smoothed with the Gaussian of width σ, leaving only a slowly varying ρ^qσ(z) in the Poisson-like formulation of the LMF equation given in Eq. 5. This ρ^qσ(z) again reflects the formation of dipole layers near each wall. Thus, solving a 1-dimensional LMF equation yields remarkably accurate results for the polarization potential profile Φ_pol(z), as shown in Fig. 6. Such findings give us confidence that the LMF treatment of interesting and more complex simulations, such as water–membrane interfaces and liquid–liquid interfaces, will prove to be just as simple.

Conclusions and Outlook

Here, we have demonstrated the utility of the LMF approach both in simulating a molecular system in a slab geometry and also in analyzing the structural properties of this system. The inclusion of a restructured electrostatic potential V_R allows for remarkably accurate results in this slab geometry. The spectacular failure of the spherical truncation of 1/r interactions to describe even the density profile of water in an applied field is corrected through the use of LMF theory. Earlier known failures of truncated water models to describe the electrostatic potential for confined water in zero field are also corrected by LMF theory regardless of the specific details of the confining surfaces. Furthermore, the splitting of 1/r into short-ranged and long-ranged components provides a natural way to disentangle the impact of local hydrogen bonding and of long-ranged dipolar interactions in various aqueous systems. In light of LMF theory's success in this slab geometry, we certainly believe that it could be of use in other physically relevant slab-geometry simulations such as biological membranes and liquid–liquid interfaces. Perhaps it also could help resolve a recent dispute over whether or not depletion of water between capacitor plates occurs in the simulation of an open system (34–36). Furthermore, given previous work treating primitive models of ionic solutions with LMF theory (23, 28, 30), the study of aqueous ionic solutions is another natural step forward. Truncating all 1/r interactions via v₀(r) with the same σ would lead to a natural mimic system for such ionic solutions.

Examination of either V_R or the smoothed charge density ρ^qσ prescribed by LMF theory also clarifies the underlying electrostatics and the failure of spherical truncations. The Gaussian-smoothed perspective for long-ranged electrostatics greatly elucidates the perceived contradiction between seemingly good charge density profiles but very poor polarization potential profiles for water confined between walls in zero field. Analysis of long-ranged electrostatic properties by using the smoothed charge density seems more plausible because the point charges in classical site–site models approximate the real quantum-mechanical electron cloud, and longer-ranged electrostatic effects should be relatively insensitive to the exact placement and magnitude of point charges in models with the same dipole moment. Thus, in GT water, the short-ranged cores are altered to reflect both hydrogen-bonding and polar attractions represented by v₀(r) as well as the LJ cores. Simulation by using GT water without a self-consistent V_R will only be reasonable if the resulting ρ^qσ reflects the underlying electrostatics. Because V_R is the solution to a modified Poisson's equation, we should not expect this to be generally true.

Overall, the use of LMF theory for electrostatics in simulations seems very promising, but certainly this approach needs more investigation. Future work should extend the derivation of the simple LMF electrostatic potential given in Eq. 4 to site-to-site molecular models (Z. Hu, J.M.R. and J.D.W., unpublished data) and describe simple analytic corrections for thermodynamics of these mimic systems (J.M.R. and J.D.W., unpublished data). Further investigation into optimizing the solution of the LMF equation is called for, and the slowly varying nature of the smoothed charge density suggests that simple approximations could prove quite useful (30). At present, the LMF mapping is for equilibrium properties only. However, current work has shown that LMF theory used in tandem with with a Born–Oppenheimer-like approximation can describe the slow dynamics and the associated free-energy landscape of charged polymers folding in the presence of mobile salt ions (N.A. Denesyuk and J.D.W., unpublished data). In general, however, a full extension of these ideas to dynamics remains a challenging open question. Its resolution may clarify the relationship between the fluctuations of these short-ranged systems and those in experimental systems.

Although all nonuniform systems examined herein have a slab geometry, LMF theory is valid for general nonuniformity. As such, we view the work presented in this article as a demonstration of the utility of LMF theory to treat and to understand electrostatic interactions in molecular systems far more complex than those studied here.

Methods

For this research, we have used the DLPOLY 2.16 molecular dynamics package (37) modified to include the Lennard–Jones wall potential (25), the Pt(111) wall potential (33), the LMF theory v₀(r) potential, and the V_R(z) mean field potential, as well as the 3-dimensional corrected Ewald sum for slab systems (EW3DC) (38). SPC/E water (3) was simulated at constant volume and temperature by using the Berendsen thermostat (39) with a time step of 1 fs. For the bulk water pair correlation functions, we simulated 1,728 SPC/E molecules at bulk density for 1.5 ns after 500 ps of equilibration at a temperature of 300 K. For each nonuniform system and electrostatic technique we also have carried out at least 500 ps of equilibration and, subsequently, 1.5 ns of simulation at 298 K. The lateral simulation box size for the smooth Lennard–Jones walls was 2.772 nm × 2.772 nm with 1,024 SPC/E molecules, and the walls were placed at ±2.25 nm to obtain the bulk density of water in the central region. When the electric field of 10.0 V/nm was applied, the walls were moved to 2.153 nm to maintain the bulk density of water. (Therefore, effects of electrostriction were not accounted for.) The Pt(111) wall area was 3.05 nm × 2.88 nm with 1,054 SPC/E molecules, and the z_c for the walls (as defined in ref. 33) were placed at ±1.56 nm. In both cases the simulation cell spanned 14.0 nm in the z direction, though this was only crucial for the Ewald summations to maintain sufficient vacuum space between periodically replicated slabs.

The three electrostatic techniques used for the slab geometry were Ewald summation (15), minimum image simulations using a spherically truncated v₀(r), and minimum image simulations using a spherically truncated v₀(r) and V_R(z). The Ewald technique used was the slab-corrected 3-dimensional Ewald scheme of Yeh and Berkowitz (38). We used an α of 3.40 nm^-1 and k-vectors spanning (12,12,60) for LJ walls and (13,12,60) for Pt(111) walls. For v₀(r) we used a σ of 0.60 nm for slab simulations and 0.40 nm for the bulk simulations. For bulk, the use of a larger σ of 0.60 nm produced equivalent results, though the δg will be slightly different.

The LMF equation must be solved self-consistently so that the ρ^q determined via simulation by using V_R yields the same LMF equation prediction of V_R. Self-consistently solving the LMF equation in previous simulation work (29) proved relatively straightforward with a linear mixing parameter λ as V_Rⁱ⁺¹ = λ · V_R^LMF + (1-λ) · V_Rⁱ, where V_R^LMF is the restructured potential determined by inserting the ρ^q due to V_Rⁱ into Eq. 4. However, for these systems, the solution proved more challenging because of the need to distinguish between inherent equilibrium fluctuations in molecular orientations in the short-ranged system and the induced response of the short-ranged system to a variation in the small and slowly varying V_R.

We may still solve the LMF equation fairly quickly, given a reasonable initial V_R⁰, by carrying out N simultaneous short-time simulations with distinct starting points during the evolution to a self-consistent solution. The averaged ρ^q of the N systems simulated in V_Rⁱ at the i^th solution step better revealed the response of the system to the change in V_R from V_R^i-1 to V_Rⁱ. These short-time simulations each ran for 25 ps of equilibration and 50 ps of accumulation with a larger Δt of 2.5 fs. Then the ρ^q(z) resulting from each of the N short-time simulations are averaged together and inserted in the LMF equation to yield V_R^LMF. Just as for the simpler solution schemes, we then used linear mixing with λ ranging from 0.05 to 0.10 to predict the subsequent V_Rⁱ⁺¹. For σ = 0.60 nm, V_R converged within 3 steps of 10 parallel simulations with no applied E₀ and within 10 steps of 10 parallel simulations with an applied field of 10.0 V/nm. However, this process was employed simply to demonstrate the utility of a self-consistent LMF treatment of electrostatics in molecular simulations. An optimized approach will likely appear quite different. Currently, by using the interpretation of the LMF equation given in Eq. 5, we are exploring options such as fast Poisson solvers and Car–Parinello techniques in solving the equation adaptively during simulation. Such solutions for V_R would still be static in nature once self-consistency was achieved.