Angular resolution and range of dipole–dipole correlations in water

How Thick is the Air-Water Interface?─A Direct Experimental Measurement of the Decay Length of the Interfacial Structural Anisotropy

The air-water interface is a highly prevalent phase boundary impacting many natural and artificial processes. The significance of this interface arises from the unique properties of water molecules within the interfacial region, with a crucial parameter being the thickness of its structural anisotropy, or "healing depth". This quantity has been extensively assessed by various simulations which have converged to a prediction of a remarkably short length of ∼6 Å. Despite the absence of any direct experimental measurement of this quantity, this predicted value has surprisingly become widely accepted as fact. Using an advancement in nonlinear vibrational spectroscopy, we provide the first measurement of this thickness and, indeed, find it to be ∼6-8 Å, finally confirming the prior predictions. Lastly, by combining the experimental results with depth-dependent second-order spectra calculated from ab initio parametrized molecular dynamics simulations, which are also in excellent agreement with this experimental result, we shed light on this surprisingly short correlation length of molecular orientations at the interface.

Bottom-up derived flexible water model with dipole and quadrupole moments for coarse-grained molecular simulations

A bottom-up CG water model is developed to capture the electrostatic multipoles, structural correlation and thermodynamics of water.

A fast and accurate computational method for the linear-combination-based isotropic periodic sum

An isotropic periodic sum (IPS) is a powerful technique to reasonably calculate intermolecular interactions for wide range of molecular systems under periodic boundary conditions. A linear-combination-based IPS (LIPS) has been developed to attain computational accuracy close to an exact lattice sum, such as the Ewald sum. The algorithm of the original LIPS method has a high computational cost because it needs long-range interaction calculations in real space. This becomes a performance bottleneck for long-time molecular simulations. In this work, the combination of an LIPS and fast Fourier transform (FFT) was developed, and evaluated on homogeneous and heterogeneous molecular systems. This combinational approach of LIPS/FFT attained computational efficiency close to that of a smooth particle mesh Ewald while maintaining the same high accuracy as the original LIPS. We concluded that LIPS/FFT has great potential to extend the capability of IPS techniques for the fast and accurate computation of many types of molecular systems.

Critical test of isotropic periodic sum techniques with group-based cut-off schemes

Truncation is still chosen for many long-range intermolecular interaction calculations to efficiently compute free-boundary systems, macromolecular systems and net-charge molecular systems, for example. Advanced truncation methods have been developed for long-range intermolecular interactions. Every truncation method can be implemented as one of two basic cut-off schemes, namely either an atom-based or a group-based cut-off scheme. The former computes interactions of “atoms” inside the cut-off radius, whereas the latter computes interactions of “molecules” inside the cut-off radius. In this work, the effect of group-based cut-off is investigated for isotropic periodic sum (IPS) techniques, which are promising cut-off treatments to attain advanced accuracy for many types of molecular system. The effect of group-based cut-off is clearly different from that of atom-based cut-off, and severe artefacts are observed in some cases. However, no severe discrepancy from the Ewald sum is observed with the extended IPS techniques.

The "Hot-Solvent/Cold-Solute" Problem Revisited

The temperature steers the equilibrium and nonequilibrium conformational dynamics of macromolecules in solution. Therefore, corresponding molecular dynamics simulations require a strategy for temperature control which should guarantee that the experimental statistical ensemble is also sampled in silico. Several algorithms for temperature control have been proposed. All these thermostats interfere with the macromolecule's "natural" dynamics as given by the Newtonian mechanics. Furthermore, using a single thermostat for an inhomogeneous solute-solvent system can lead to stationary temperature gradients. To avoid this "hot solvent/cold solute" problem, two separate thermostats are frequently applied, one to the solute and one to the solvent. However, such a separate temperature control will perturb the dynamics of the macromolecule much more strongly than a global one and, therefore, can introduce large artifacts into its conformational dynamics. Based on the concept that an explicit solvent environment represents an ideal thermostat concerning the magnitude and time correlation of temperature fluctuations of the solute, we propose a temperature control strategy that, on the one hand, provides a homogeneous temperature distribution throughout the system together with the correct statistical ensemble for the solute molecule while, on the other hand, minimally perturbing its dynamics.

Fluctuations of local electric field and dipole moments in water between metal walls

We examine the thermal fluctuations of the local electric field Ek (loc) and the dipole moment μk in liquid water at T = 298 K between metal walls in electric field applied in the perpendicular direction. We use analytic theory and molecular dynamics simulation. In this situation, there is a global electrostatic coupling between the surface charges on the walls and the polarization in the bulk. Then, the correlation function of the polarization density pz(r) along the applied field contains a homogeneous part inversely proportional to the cell volume V. Accounting for the long-range dipolar interaction, we derive the Kirkwood-Fröhlich formula for the polarization fluctuations when the specimen volume v is much smaller than V. However, for not small v/V, the homogeneous part comes into play in dielectric relations. We also calculate the distribution of Ek (loc) in applied field. As a unique feature of water, its magnitude |Ek (loc)| obeys a Gaussian distribution with a large mean value E0 ≅ 17 V/nm, which arises mainly from the surrounding hydrogen-bonded molecules. Since |μk|E0 ∼ 30kBT, μk becomes mostly parallel to Ek (loc). As a result, the orientation distributions of these two vectors nearly coincide, assuming the classical exponential form. In dynamics, the component of μk(t) parallel to Ek (loc)(t) changes on the time scale of the hydrogen bonds ∼5 ps, while its smaller perpendicular component undergoes librational motions on time scales of 0.01 ps.

Toward the correction of effective electrostatic forces in explicit-solvent molecular dynamics simulations: restraints on solvent-generated electrostatic potential and solvent polarization

Despite considerable advances in computing power, atomistic simulations under nonperiodic boundary conditions, with Coulombic electrostatic interactions and in systems large enough to reduce finite-size associated errors in thermodynamic quantities to within the thermal energy, are still not affordable. As a result, periodic boundary conditions, systems of microscopic size and effective electrostatic interaction functions are frequently resorted to. Ensuing artifacts in thermodynamic quantities are nowadays routinely corrected a posteriori, but the underlying configurational sampling still descends from spurious forces. The present study addresses this problem through the introduction of on-the-fly corrections to the physical forces during an atomistic molecular dynamics simulation. Two different approaches are suggested, where the force corrections are derived from special potential energy terms. In the first approach, the solvent-generated electrostatic potential sampled at a given atom site is restrained to a target value involving corrections for electrostatic artifacts. In the second approach, the long-range regime of the solvent polarization around a given atom site is restrained to the Born polarization, i.e., the solvent polarization corresponding to the ideal situation of a macroscopic system under nonperiodic boundary conditions and governed by Coulombic electrostatic interactions. The restraints are applied to the explicit-water simulation of a hydrated sodium ion, and the effect of the restraints on the structural and energetic properties of the solvent is illustrated. Furthermore, by means of the calculation of the charging free energy of a hydrated sodium ion, it is shown how the electrostatic potential restraint translates into the on-the-fly consideration of the corresponding free-energy correction terms. It is discussed how the restraints can be generalized to situations involving several solute particles. Although the present study considers a very simple system only, it is an important step toward the on-the-fly elimination of finite-size and approximate-electrostatic artifacts during atomistic molecular dynamics simulations.

Design of a reaction field using a linear-combination-based isotropic periodic sum method

In our previous study (Takahashi et al., J. Chem. Theory Comput. 2012, 8, 4503), we developed the linear-combination-based isotropic periodic sum (LIPS) method. The LIPS method is based on the extended isotropic periodic sum theory that produces a ubiquitous interaction potential function to estimate homogeneous and heterogeneous systems. The LIPS theory also provides the procedure to design a periodic reaction field. To demonstrate this, in the present work, a novel reaction field of the LIPS method was developed. The novel reaction field was labeled LIPS-SW, because it provides an interaction potential function with a shape that resembles that of the switch function method. To evaluate the ability of the LIPS-SW method to describe in homogeneous and heterogeneous systems, we carried out molecular dynamics (MD) simulations of bulk water and water-vapor interfacial systems using the LIPS-SW method. The results of these simulations show that the LIPS-SW method gives higher accuracy than the conventional interaction potential function of the LIPS method. The accuracy of simulating water-vapor interfacial systems was greatly improved, while that of bulk water systems was maintained using the LIPS-SW method. We conclude that the LIPS-SW method shows great potential for high-accuracy, high-performance computing to allow large scale MD simulations. © 2014 Wiley Periodicals, Inc.

Hyper-Rayleigh scattering from correlated molecules

The polarization dependence of hyper-Rayleigh scattering has been calculated for spherical domains of orientation correlated molecules. Distributions with radial or azimuthal mean polar orientation of the molecules are found that give results consistent with experimental observations, and expressions for the polarization ratios in terms of the product of correlation strength and correlated domain size are derived for these distributions. Assuming a plausible correlation strength, it is estimated that the correlated domain size in typical polar liquids is of order 100 molecular diameters.

Understanding the Origins of Dipolar Couplings and Correlated Motion in the Vibrational Spectrum of Water

The combination of vibrational spectroscopy and molecular dynamics simulations provides a powerful tool to obtain insights into the molecular details of water structure and dynamics in the bulk and in aqueous solutions. Applying newly developed approaches to analyze correlations of charge currents, molecular dipole fluctuations, and vibrational motion in real and k-space, we compare results from nonpolarizable water models, widely used in biomolecular modeling, to ab initio molecular dynamics. For the first time, we unfold the infrared response of bulk water into contributions from correlated fluctuations in the three-dimensional, anisotropic environment of an average water molecule, from the OH-stretching region down to the THz regime. Our findings show that the absence of electronic polarizability in the force field model not only results in differences in dipolar couplings and infrared absorption but also induces artifacts into the correlated vibrational motion between hydrogen-bonded water molecules, specifically at the intramolecular bending frequency. Consequently, vibrational motion is partially ill-described with implications for the accuracy of non-self-consistent, a posteriori methods to add polarizability.

Long-distance correlations in molecular orientations of liquid water and shape-dependent hydrophobic force

Liquid water, at ambient conditions, has short-range density correlations which are well known in literature. Surprisingly, large scale molecular-dynamics simulations reveal an unusually long-distance correlation in "longitudinal" part of dipole-dipole orientational correlations. It is nonvanishing even at 75 A and falls off exponentially with a correlation length of about 24 A beyond solvation region. Numerical evidence suggests that the long-range nature of dipole-dipole correlation is due to underlying fluctuating network of hydrogen bonds in the liquid phase. This correlation is shown to give a shape dependant attraction between two hydrophobic surfaces at large distances of separation and the range of this attractive force is in agreement with experiments. In addition it is seen that quadrupolar fluctuations vanish within the first solvation peak (3 A) .

Nondielectric long-range solvation of polar liquids in cubic symmetry

Long-range solvation properties of strongly coupled dipolar systems simulated using the Ewald and reaction field methods are assessed by using electric fluctuation formulas for a dielectric medium. Some components of the fluctuating electric multipole moments are suppressed, whereas other components are favored as the boundary of the simulation box is approached. An analysis of electrostatic interactions in a periodic cubic system suggests that these structural effects are due to the periodicity embedded in the Ewald method. Furthermore, the results obtained using the reaction field method are very similar to those obtained using the Ewald method, an effect which we attribute to the use of toroidal boundary conditions in the former case. Thus, the long-range solvation properties of polar liquids simulated using either of the two methods are nondielectric in their character.

Molecular dynamics simulations of a reversibly folding beta-heptapeptide in methanol: influence of the treatment of long-range electrostatic interactions

Eight 100-ns molecular dynamics simulations of a beta-heptapeptide in methanol at 340 K (within cubic periodic computational boxes of about 6-nm edge) are reported and compared. These simulations were performed with three different charge-state combinations at the peptide termini, one of them with or without a neutralizing chloride counterion, and using either the lattice-sum (LS) or reaction-field (RF) scheme to handle electrostatic interactions. The choice of the electrostatic scheme has essentially no influence on the folding-unfolding equilibrium when the peptide termini are uncharged and only a small influence when the peptide is positively charged at its N-terminus (with or without inclusion of a neutralizing chloride counterion). However, when the peptide is zwitterionic, the LS scheme leads to preferential sampling of the high-dipole folded helical state, whereas the RF scheme leads to preferential sampling of a low-dipole unfolded salt-bridged state. A continuum electrostatics analysis based on the sampled configurations (zwitterionic case) suggests that the LS scheme stabilizes the helical state through artificial periodicity, but that the magnitude of this perturbation is essentially negligible (compared to the thermal energy) for the large box size and relatively polar solvent considered. The results thus provide clear evidence (continuum electrostatics analysis) for the absence of LS artifacts and some indications (still not definitive because of the limited sampling of the folding-unfolding transition) for the presence of RF artifacts in this specific system.

Molecular mechanism of water bridge buildup: field-induced formation of nanoscale menisci

We perform molecular dynamics calculations to describe, at the molecular level, the formation of a water bridge induced by an electric field. Restriction of orientational degrees of freedom (confinement) of water dipoles at the interfaces leads to a polarizability that depends on the shape of the water system, that is, droplet versus pillar. Above a threshold field of 1.2 V nm(-1), the competition between orientational confinement and electric field leads to the sudden formation of a water pillar. The formation of a water bridge is marked by a first order discontinuity in the total energy of the system. The simulations offer a molecular explanation for the threshold voltage and hysteresis behavior observed in the formation of nanoscale liquid bridges with a force microscope.

Electrostatics of proteins in dielectric solvent continua. II. First applications in molecular dynamics simulations

In the preceding paper by Stork and Tavan, [J. Chem. Phys. 126, 165105 (2007)], the authors have reformulated an electrostatic theory which treats proteins surrounded by dielectric solvent continua and approximately solves the associated Poisson equation [B. Egwolf and P. Tavan, J. Chem. Phys. 118, 2039 (2003)]. The resulting solution comprises analytical expressions for the electrostatic reaction field (RF) and potential, which are generated within the protein by the polarization of the surrounding continuum. Here the field and potential are represented in terms of Gaussian RF dipole densities localized at the protein atoms. Quite like in a polarizable force field, also the RF dipole at a given protein atom is induced by the partial charges and RF dipoles at the other atoms. Based on the reformulated theory, the authors have suggested expressions for the RF forces, which obey Newton's third law. Previous continuum approaches, which were also built on solutions of the Poisson equation, used to violate the reactio principle required by this law, and thus were inapplicable to molecular dynamics (MD) simulations. In this paper, the authors suggest a set of techniques by which one can surmount the few remaining hurdles still hampering the application of the theory to MD simulations of soluble proteins and peptides. These techniques comprise the treatment of the RF dipoles within an extended Lagrangian approach and the optimization of the atomic RF polarizabilities. Using the well-studied conformational dynamics of alanine dipeptide as the simplest example, the authors demonstrate the remarkable accuracy and efficiency of the resulting RF-MD approach.

Electrostatics of proteins in dielectric solvent continua. I. Newton’s third law marries qE forces

The authors reformulate and revise an electrostatic theory treating proteins surrounded by dielectric solvent continua [B. Egwolf and P. Tavan, J. Chem. Phys. 118, 2039 (2003)] to make the resulting reaction field (RF) forces compatible with Newton's third law. Such a compatibility is required for their use in molecular dynamics (MD) simulations, in which the proteins are modeled by all-atom molecular mechanics force fields. According to the original theory the RF forces, which are due to the electric field generated by the solvent polarization and act on the partial charges of a protein, i.e., the so-called qE forces, can be quite accurately computed from Gaussian RF dipoles localized at the protein atoms. Using a slightly different approximation scheme also the RF energies of given protein configurations are obtained. However, because the qE forces do not account for the dielectric boundary pressure exerted by the solvent continuum on the protein, they do not obey the principle that actio equals reactio as required by Newton's third law. Therefore, their use in MD simulations is severely hampered. An analysis of the original theory has led the authors now to a reformulation removing the main difficulties. By considering the RF energy, which represents the dominant electrostatic contribution to the free energy of solvation for a given protein configuration, they show that its negative configurational gradient yields mean RF forces obeying the reactio principle. Because the evaluation of these mean forces is computationally much more demanding than that of the qE forces, they derive a suggestion how the qE forces can be modified to obey Newton's third law. Various properties of the thus established theory, particularly issues of accuracy and of computational efficiency, are discussed. A sample application to a MD simulation of a peptide in solution is described in the following paper [M. Stork and P. Tavan, J. Chem. Phys., 126, 165106 (2007).

Molecular dynamics study of orientational cooperativity in water

Recent experiments on liquid water show collective dipole orientation fluctuations dramatically slower than expected (with relaxation time > 50 ns) [D.P. Shelton, Phys. Rev. B 72, 020201(R) (2005)]. Molecular dynamics simulations of extended simple point charge (SPC/E) water show a large vortexlike structure of the dipole field at ambient conditions surviving over [J. Higo, Proc. Natl. Acad. Sci. U.S.A. 98, 5961 (2001)]. Both results disagree with previous results on water dipoles in similar conditions, for which autocorrelation times are a few picoseconds. Motivated by these recent results, we study the water dipole reorientation using molecular dynamics simulations of the SPC/E model in bulk water for temperatures ranging from ambient 300 K down to the deep supercooled region of the phase diagram at 210 K. First, we calculate the dipole autocorrelation function and find that our simulations are well described by a stretched exponential decay, from which we calculate the orientational autocorrelation time t(alpha). Second, we define a second characteristic time, namely, the time required for the randomization of molecular dipole orientation, the self-dipole randomization time t(r), which is an upper limit on t(alpha); we find that t(r) is approximately equal to 5t(alpha). Third, to check if there are correlated domains of dipoles in water which have large relaxation times compared to the individual dipoles, we calculate the randomization time t(box) of the site-dipole field, the net dipole moment formed by a set of molecules belonging to a box of edge L(box). We find that the site-dipole randomization time t(box) is approximately equal to 2.5t(alpha) for L(box) approximately equal to 3 A, i.e., it is shorter than the same quantity calculated for the self-dipole. Finally, we find that the orientational correlation length is short even at low T.

Structural instability of the prion protein upon M205S/R mutations revealed by molecular dynamics simulations

The point mutations M205S and M205R have been demonstrated to severely disturb the folding and maturation process of the cellular prion protein (PrP(C)). These disturbances have been interpreted as consequences of mutation-induced structural changes in PrP, which are suggested to involve helix 1 and its attachment to helix 3, because the mutated residue M205 of helix 3 is located at the interface of these two helices. Furthermore, current models of the prion protein scrapie (PrP(Sc)), which is the pathogenic isoform of PrP(C) in prion diseases, imply that helix 1 disappears during refolding of PrP(C) into PrP(Sc). Based on molecular-dynamics simulations of wild-type and mutant PrP(C) in aqueous solution, we show here that the native PrP(C) structure becomes strongly distorted within a few nanoseconds, once the point mutations M205S and M205R have been applied. In the case of M205R, this distortion is characterized by a motion of helix 1 away from the hydrophobic core into the aqueous environment and a subsequent structural decay. Together with experimental evidence on model peptides, this decay suggests that the hydrophobic attachment of helix 1 to helix 3 at M205 is required for its correct folding into its stable native structure.

Combining the lattice-sum and reaction-field approaches for evaluating long-range electrostatic interactions in molecular simulations

A new scheme, the lattice-sum-emulated reaction-field (LSERF) method, is presented that combines the lattice-sum (LS) and reaction-field (RF) approaches for evaluating electrostatic interactions in molecular simulations. More precisely, the LSERF scheme emulates a RF calculation (based on an atomic cutoff) via the LS machinery. This is achieved by changing the form of the electrostatic interactions in a standard LS calculation (Coulombic) to the form corresponding to RF electrostatics (Coulombic plus quadratic reaction-field correction term, truncated at the cutoff distance). It is shown (both analytically and numerically) that in the limit of infinite reciprocal-space accuracy, (i) the LSERF scheme with a finite reaction-field cutoff and a given reaction-field permittivity is identical to the RF scheme with the same parameters (and an atomic cutoff), and (ii) the LSERF scheme is identical to the LS scheme in the limit of an infinite reaction-field cutoff, irrespective of the reaction-field permittivity. This new scheme offers two key advantages: (i) from a conceptual point of view, it shows that there is a continuity between the RF and LS schemes and unifies them into a common framework; (ii) from a practical point of view, it allows us to perform RF calculations with arbitrarily large reaction-field cutoff distances for the same computational costs as a corresponding LS calculation. The optimal choice for the cutoff will be the one that achieves the best compromise between artifacts arising from the dielectric heterogeneity of the system (short cutoff) and its artificial periodicity (long cutoff). The implementation of the LSERF method is extremely easy, requiring only very limited modifications of any standard LS code. For practical applications to biomolecular systems, the use of the LSERF scheme with large reaction-field cutoff distances is expected to represent a significant improvement over the current RF simulations involving comparatively much shorter cutoffs.

Molecular dynamics simulations indicate a possible role of parallel beta-helices in seeded aggregation of poly-Gln

The molecular structures of amyloid fibers characterizing neurodegenerative diseases such as Huntington's or transmissible spongiform encephalopathies are unknown. Recently, x-ray diffraction patterns of poly-Gln fibers and electron microscopy images of two-dimensional crystals formed from building blocks of prion rods have suggested that the corresponding amyloid fibers are generated by the aggregation of parallel beta-helices. To explore this intriguing concept, we study the stability of small beta-helices in aqueous solution by molecular dynamics simulations. In particular, for the Huntington aggregation nucleus, which is thought to be formed of poly-Gln polymers, we show that three-coiled beta-helices are unstable at the suggested circular geometries and stable at a triangular shape with 18 residues per coil. Moreover, we demonstrate that individually unstable two-coiled triangular poly-Gln beta-helices become stabilized upon dimerization, suggesting that seeded aggregation of Huntington amyloids requires dimers of at least 36 Gln repeats (or monomers of approximately 54 Gln) for the formation of sufficiently stable aggregation nuclei. An analysis of our results and of sequences occurring in native beta-helices leads us to the proposal of a revised model for the PrP(Sc) aggregation nucleus.