Observations concerning the treatment of long-range interactions in molecular dynamics simulations

Molecular dynamics simulations of pure water employing two different empirical water models have been used to study the effects of different methods for truncation of long-range interactions in molecular mechanics calculations. As has been observed previously in integral equation studies, “shifting” these interactions on an atom-by-atom basis was found to produce artificial structuring in the water and affect diffusion rates. In cases where some form of short-range truncation must be used, the ST2 switching function applied on a group-by-group basis was found to be the most realistic procedure. If atom-based shifting must be employed, a cutoff distance greater than or equal to 12.0 Å was found to be required to produce realistic results. © 1993 John Wiley & Sons, Inc.     

Conformational and energetic effects of truncating nonbonded interactions in an aqueous protein dynamics simulation

In a previous aqueous protein dynamics study, we compared the rms deviation relative to the crystal structure for distance-dependent and constant dielectric models with and without a nonbonded cutoff. The structures obtained from a constant dielectric simulation with a cutoff were substantially different from the structures obtained from a distance-dependent dielectric simulation, with and without cutoff, and a constant dielectric model without a cutoff. In fact, structures from the distance-dependent dielectric simulations were insensitive to the nonbonded cutoff and in good agreement with the structures generated from the constant dielectric simulation without a cutoff. In addition, the solute-solvent temperature differential and solvent evaporation artifacts, characteristic of the constant dielectric simulation with a cutoff, were not present for the distance-dependent dielectric simulations. In this current work, we explore whether this dielectric-dependent cutoff-sensitive behavior for a constant dielectric model arises from the discontinuities in the forces at the nonbonded cutoff or from neglecting the structure-stabilizing interactions beyond the nonbonded cutoff. We also examine the origin of the dielectric-dependent artifacts, and its potential influence on the structural disparity. Several protocols for protein dynamics simulations are compared using both constant and distance-dependent dielectric models, including implementation of a switching function and a nonbonded cutoff and two different temperature coupling algorithms. We show that the distance-dependent dielectric model conserves energy in the SPASMS molecular mechanics and dynamics software for the time steps and nonbonded cutoffs commonly used in macromolecule simulations. Although the switching function simulation also conserved energy over a range of commonly used cutoffs, the constant dielectric model with a switching function yielded conformational results more similar to a constant dielectric simulation without a switching function than to a constant dielectric model without a nonbonded cutoff. Therefore, the conformational disparity between the dielectric models arises from neglecting important structure-stabilizing interactions beyond the cutoff, rather than differences in energy conservation. © 1993 John Wiley & Sons, Inc.     

The effect of different cutoff schemes in molecular simulations of proteins

Molecular simulations of nanoscale systems invariably involve assumptions and approximations to describe the electrostatic interactions, which are long-ranged in nature. One approach is the use of cutoff schemes with a reaction-field contribution to account for the medium outside the cutoff scheme. Recent reports show that macroscopic properties may depend on the exact choice of cutoff schemes in modern day simulations. In this work, a systematic analysis of the effects of different cutoff schemes was performed using a set of 52 proteins. We find no statistically significant differences between using a twin-range or a single-range cutoff scheme. Applying the cutoff based on charge groups or based on atomic positions, does lead to significant differences, which is traced to the cutoff noise for energies and forces. While group-based cutoff schemes show increased cutoff noise in the potential energy, applying an atomistic cutoff leads to artificial structure in the solvent at the cutoff distance. Carefully setting the temperature control, or using an atomistic cutoff for the solute and a group-based cutoff for the solvent significantly reduces the effects of the cutoff noise, without introducing structure in the solvent. This study aims to deepen the understanding of the implications different cutoffs have on molecular dynamics simulations.     

Molecular dynamics simulations of the isolated domain 1 of annexin I

Annexin molecules consist of a symmetrical arrangement of four domains of identical folds but very different sequences. Nuclear magnetic resonance (NMR) experiments on the isolated domains of annexin I in aqueous solution have indicated that domain 1 retains its native structure whereas domain 2 unfolds. Therefore these two domains constitute interesting models for comparative simulations of structural stability using molecular dynamics. Here we present the preliminary results of molecular dynamics simulations of the isolated domain 1 in explicit water at 300 K, using two different simulation protocols. For the first, domain 1 was embedded in a 46 ? cubic box of water. A group-based non-bonded cut-off of 9 ? with a 5–9 ? non-bonded switching function was used and a 2 fs integration step. Bonds containing hydrogens were constrained with the SHAKE algorithm. These conditions led to unfolding of the domain within 400 ps at 300 K. In the second protocol, the domain was embedded in a 62 ? cubic box of water. An atom-based non-bonded cut-off of 8–12 ? using a force switching function for electrostatics and a shifting function for van der Waals interactions were used with a 1 fs integration step. This second protocol led to a native-like conformation of the domain in accord with the NMR data which was stable over the whole trajectory (∼2 ns). A small, but well-defined relaxation of the structure, from that observed for the same domain in the entire protein, was observed. This structural relaxation is described and methodological aspects are discussed. Received: 10 May 1998 / Accepted: 4 August 1998 / Published online: 2 November 1998     

Harmonic analysis of large systems. II. Comparison of different protein models

A series of normal mode analyses of bovine pancreatic trypsin inhibitor (BPTI) has been performed. The results of modifying the long-range truncation of electrostatics, reducing the conformational space of the system (reduced basis normal mode analysis), and using different parameter sets and models for the potential function are reported. Both explicit (904 atoms) and polar hydrogen (580 atoms) representations of BPTI were examined and produced nearly identical normal mode vectors but slightly modified vibrational frequencies. The truncation methods—no cutoff, shift, and switch—were examined, and the use of a short switching function was found to alter harmonic motion greatly. A table relating the different cutoff methods to several previously published frequencies for BPTI indicates that the diversity of published lowest frequencies is due to the use of different electrostatic models rather than to inherent differences in the models or energy parameters. Examining reduced basis results demonstrates that a dihedral basis yields similar normal mode vectors, though the vibrational frequencies are shifted to higher values. The analysis of BPTI harmonic dynamics using a spherical harmonic reduced basis set yields significantly altered dynamics, indicating that BPTI is not well represented as a homogeneous object at low temperatures. © 1995 John Wiley & Sons, Inc.

1 This article is a U.S. Government work and, as such, is in the public domain in the United States of America.

     

Effects of system net charge and electrostatic truncation on all‐atom constant pH molecular dynamics

Constant pH molecular dynamics offers a means to rigorously study the effects of solution pH on dynamical processes. Here, we address two critical questions arising from the most recent developments of the all‐atom continuous constant pH molecular dynamics (CpHMD) method: (1) What is the effect of spatial electrostatic truncation on the sampling of protonation states? (2) Is the enforcement of electrical neutrality necessary for constant pH simulations? We first examined how the generalized reaction field and force‐shifting schemes modify the electrostatic forces on the titration coordinates. Free energy simulations of model compounds were then carried out to delineate the errors in the deprotonation free energy and salt‐bridge stability due to electrostatic truncation and system net charge. Finally, CpHMD titration of a mini‐protein HP36 was used to understand the manifestation of the two types of errors in the calculated pK_a values. The major finding is that enforcing charge neutrality under all pH conditions and at all time via cotitrating ions significantly improves the accuracy of protonation‐state sampling. We suggest that such finding is also relevant for simulations with particle mesh Ewald, considering the known artifacts due to charge‐compensating background plasma. © 2014 Wiley Periodicals, Inc.     

Accurate and efficient corrections for missing dispersion interactions in molecular simulations

In simulations, molecular dispersion interactions are frequently neglected beyond a cutoff of around 1 nm. In some cases, analytical corrections appropriate for isotropic systems are applied to the pressure and/or the potential energy. Here, we show that in systems containing macromolecules, either of these approaches introduce statistically significant errors in some observed properties; for example, the choice of cutoff can affect computed free energies of ligand binding to proteins by 1 to 2 kcal/mol. We review current methods for eliminating this cutoff-dependent behavior of the dispersion energy and identify some situations where they fail. We introduce two new formalisms, appropriate for binding free energy calculations, which overcome these failings, requiring minimal computational effort beyond the time required to run the original simulation. When these cutoff approximations are applied, which can be done after all simulations are completed, results are consistent across simulations run with different cutoffs. In many situations, simulations can be run with even shorter cutoffs than typically used, resulting in increased computational efficiency.     

Simple and accurate scheme to compute electrostatic interaction: Zero-dipole summation technique for molecular system and application to bulk water

The zero-dipole summation method was extended to general molecular systems, and then applied to molecular dynamics simulations of an isotropic water system. In our previous paper [I. Fukuda, Y. Yonezawa, and H. Nakamura, J. Chem. Phys. 134, 164107 (2011)], for evaluating the electrostatic energy of a classical particle system, we proposed the zero-dipole summation method, which conceptually prevents the nonzero-charge and nonzero-dipole states artificially generated by a simple cutoff truncation. Here, we consider the application of this scheme to molecular systems, as well as some fundamental aspects of general cutoff truncation protocols. Introducing an idea to harmonize the bonding interactions and the electrostatic interactions in the scheme, we develop a specific algorithm. As in the previous study, the resulting energy formula is represented by a simple pairwise function sum, enabling facile applications to high-performance computation. The accuracy of the electrostatic energies calculated by the zero-dipole summation method with the atom-based cutoff was numerically investigated, by comparison with those generated by the Ewald method. We obtained an electrostatic energy error of less than 0.01% at a cutoff length longer than 13 A for a TIP3P isotropic water system, and the errors were quite small, as compared to those obtained by conventional truncation methods. The static property and the stability in an MD simulation were also satisfactory. In addition, the dielectric constants and the distance-dependent Kirkwood factors were measured, and their coincidences with those calculated by the particle mesh Ewald method were confirmed, although such coincidences are not easily attained by truncation methods. We found that the zero damping-factor gave the best results in a practical cutoff distance region. In fact, in contrast to the zero-charge scheme, the damping effect was insensitive in the zero-charge and zero-dipole scheme, in the molecular system we treated. We discussed the origin of this difference between the two schemes and the dependence of this fact on the physical system. The use of the zero damping-factor will enhance the efficiency of practical computations, since the complementary error function is not employed. In addition, utilizing the zero damping-factor provides freedom from the parameter choice, which is not trivial in the zero-charge scheme, and eliminates the error function term, which corresponds to the time-consuming Fourier part under the periodic boundary conditions.     

Solution structure of papain as studied by molecular mechanics and molecular dynamics techniques

In this paper, we report the molecular mechanics and molecular dynamics studies of the hydration of papain using the AMBER and CHARMm programs. We studied papain in an environment with minimal hydration involving only the X-ray waters and also in the hydrated environment by adding an extra 9 Å layer of water around the residues. The effect of nonbond cutoff was studied by performing minimizations with 8 Å and 15 Å nonbond cutoffs using the program AMBER. Two different solvent models—a constant dielectric and a distance-dependent dielectric—were considered. The AMBER-minimized structure and the average structure obtained from the CHARMm simulations for papain solvated with X-ray waters are presented and compared with the X-ray crystal structure results. Results of a similar comparison of the hydrated structures were also presented. The calculated root mean square deviation between the minimized structure and the X-ray structure is smaller for the hydrated system than for the system hydrated with only the X-ray waters. Results of the molecular mechanics and molecular dynamics simulations were presented for the various regions of papain. The hydration of the active site of papain and the effect of hydration on the torsional motion of the active site residues are presented. © 1996 by John Wiley & Sons, Inc.     

Conserving energy during molecular dynamics simulations of water,proteins, and proteins in water

Molecular dynamics simulations have been carried out for a series of systems of increasing complexity including: pure water, a model polypeptide (α-helical decaglycine) in vacuo, a protein (Pancreatic Trypsin Inhibitor, PTI) in vacuo, and a fully solvated protein (PTI in water). The equations of motion were integrated using Andersen's velocity version of the Verlet algorithm with internal contraints (the RATTLE algorithm). The accuracy with which the equations of motion are integrated has been analyzed for several different simulation conditions. The effects of various nonbonded interaction truncation schemes on the conservation of energy have been examined, including the use of atomic cutoffs, and (neutral group) residue cutoffs. The use of a smoothing function to eliminate the discontinuities in the potential at the cutoff leads to a significant improvement in the accuracy of the integration for each of the systems studied. The accuracy with which the equations of motion are integrated using the RATTLE algorithm for pure water and for the solvated protein are found to be comparable when the nonbonded interactions are tapered with a smoothing function at the cutoff.     

Liquid water simulation: a critical examination of cutoff length

Cutoff treatment is the simplest approach for evaluating intermolecular interactions in molecular dynamics simulations. It has been believed that increasing cutoff length makes simulation results better. On the contrary, our results of the bulk water simulations studied within the range of cutoff lengths, 9-18 A, showed an opposite tendency: the artifact was enhanced by increasing the cutoff length. Especially, in terms of the distance dependent Kirkwood factor GK(r), it was clearly shown that the orientational behavior of water molecules becomes gradually worse as the cutoff length becomes longer. The artifact enhanced by the increased cutoff length led to a reported spurious artifact, i.e., phase transition [Y. Yonetani, Chem. Phys. Lett. 406, 49 (2005)]. Though the cutoff artifact was largely reduced by adopting a force switching technique, it did not completely remove the anomalous cutoff length dependence of the artifact. These results suggest that increasing the cutoff should not be attempted regardless of whether the switching force is adopted or not.     

Free energy derivatives: A new method for probing the convergence problem in free energy calculations

The convergence behavior of free energy calculations has been explored in more detail than in any previously reported work, using a model system of two neon atoms in a periodic box of water. We find that for thermodynamic integration-type free energy calculations as much as a nanosecond or more molecular dynamics sampling is required to obtain a fully converged value for a single λ point of the integrand. The concept of “free energy derivatives” with respect to the individual parameters of the force field is introduced. This formalism allows the total convergence of the simulation to be deconvoluted into components. A determination of the statistical “sampling ratio” from these simulations indicates that for window-type free energy calculations carried out in a periodic waterbox of typical size at least 0.6 ps of sampling should be performed at each window (0.7 ps if constraint contributions to the free energy are being determined). General methods to estimate and reduce the error in thermodynamic integration and free energy perturbation calculations are discussed. We show that the difficulty in applying such methods is determining a reliable estimate of the correlation length from a short series of data. © 1994 by John Wiley & Sons, Inc.     

Assessing group-based cutoffs and the Ewald method for electrostatic interactions in clusters and in saturated,superheated, and supersaturated vapor phases of dipolar molecules

Monte Carlo simulations in the canonical, isobaric-isothermal, grand canonical, and Gibbs ensembles were used to assess whether the computationally expensive Ewald summation method for the computation of the first-order electrostatic energy can be replaced with a simpler truncation approach for accurate simulations of the saturated, superheated, and supersaturated vapor phases of dipolar and hydrogen-bonding molecules. Rotationally averaged hydrogen fluoride dimer and trimer energies, thermophysical properties and aggregation in the superheated vapor phase of hydrogen fluoride, nucleation free energy barriers for water, and the vapor–liquid coexistence properties of hydrogen fluoride and water were investigated over a wide range of state points. We find that for densities not too close to the critical density, results obtained from simulations using a spherical potential truncation based on neutral groups (molecules or fragments) for the Coulomb interactions are statistically identical to those obtained using the Ewald summation method. Use of the simpler spherical truncation results in a significant reduction of the computational effort for simulations employing molecular mechanics force fields and also allows for straightforward implementation of many-body expansion methods to compute the potential energy from electronic structure calculations of subsystems of the entire vapor-phase system.     

Hydrated electron yields in the heavy ion radiolysis of water

Experimental measurements coupled with Monte Carlo track simulations have been used to examine the yields of hydrated electrons in the radiolysis of water with protons, helium ions, and carbon ions. Glycylglycine, in concentrations ranging from 10(-4) to 1 M, was employed as a scavenger and the production of the ammonium cation used as a probe of hydrated electron yields from about 2 ns to 20 mus. Monte Carlo track simulations employing diffusion-kinetic calculations of product yields are found to reproduce experimental observations satisfactorily. Model details are used to elucidate the heavy ion track physics and chemistry. Comparison of the heavy ion results with those found in gamma radiolysis shows intratrack reactions are significant on the nanosecond to microsecond time scale as the ion track relaxes, and that a constant (escape) yield is never attained on this time scale. Numerical interpolation techniques are used to obtain both track average and track segment yields for use in practical applications or comparison with other models. The model results give the first hints that initial ( approximately 5 ps) hydrated electron yields, and possibly other water decomposition products, are dependent on the type and energy of the incident radiation.     

Methodological problems in pressure profile calculations for lipid bilayers

From molecular dynamics simulations of a dipalmitoyl-phosphatidyl-choline (DPPC) lipid bilayer in the liquid crystalline phase, pressure profiles through the bilayer are calculated by different methods. These profiles allow us to address two central and unresolved problems in pressure profile calculations: The first problem is that the pressure profile is not uniquely defined since the expression for the local pressure involves an arbitrary choice of an integration contour. We have investigated two different choices leading to the Irving-Kirkwood (IK) and Harasima (H) expressions for the local pressure tensor. For these choices we find that the pressure profile is almost independent of the contour used, which indicates that the local pressure is well defined for a DPPC bilayer in the liquid crystalline phase. This may not be the case for other systems and we therefore suggest that both the IK and H profiles are calculated in order to test the uniqueness of the profile. The second problem is how to include electrostatic interactions in pressure profile calculations when the simulations are conducted without truncating the electrostatic potential, i.e., using the Ewald summation technique. Based on the H expression for the local pressure, we present a method for calculating the contribution to the lateral components of the local pressure tensor from electrostatic interactions evaluated by the Ewald summation technique. Pressure profiles calculated with an electrostatic potential truncation (cutoff) from simulations conducted with Ewald summation are shown to depend on the cutoff in a subtle manner which is attributed to the existence of long-ranged charge ordering in the system. However, the pressure profiles calculated with relatively long cutoffs are qualitatively similar to the Ewald profile for the DPPC bilayer studied here.     

Mechanism and kinetics of hydrated electron diffusion

Molecular dynamics simulations are used to study the mechanism and kinetics of hydrated electron diffusion. The electron center of mass is found to exhibit Brownian-type behavior with a diffusion coefficient considerably greater than that of the solvent. As previously postulated by both experimental and theoretical works, the instantaneous response of the electron to the librational motions of surrounding water molecules constitutes the principal mode of motion. The diffusive mechanism can be understood within the traditional framework of transfer diffusion processes, where the diffusive step is akin to the exchange of an extramolecular electron between neighboring water molecules. This is a second-order process with a computed rate constant of 5.0 ps(-1) at 298 K. In agreement with experiment the electron diffusion exhibits Arrhenius behavior over the temperature range of 298-400 K. We compute an activation energy of 8.9 kJ mol(-1). Through analysis of Arrhenius plots and the application of a simple random walk model it is demonstrated that the computed rate constant for exchange of an excess electron is indeed the phenomenological rate constant associated with the diffusive process.     

On "the complete basis set limit" and plane-wave methods in first-principles simulations of water

Water structure, measured by the height of the first peak in oxygen-oxygen radial distributions, is converged with respect to plane-wave basis energy cutoffs for ab initio molecular dynamics simulations, confirming the reliability of plane-wave methods.     

Hybrid integral equation/simulation model for enhancing free energy computations

Integral equation theory is used for extrapolating free energy data from molecular simulations of a reference state with respect to a modification of the interaction potential. The methodology is applied to the correction of artefacts arising from potential shifting and truncation. Corrective contributions for the hydration free energy with respect to the full potential are analysed for the case that both the solute-solvent as well as the solvent-solvent potentials are truncated and modified by a shifted-force term, reaching beyond the range of the dielectric continuum approximation and simple long-range correction expressions. The model systems argon in water and pure water are used as examples for apolar and polar solutes, revealing significant correction contributions even for the short-ranged dispersive interactions and the magnitude of solute-solvent and solvent-solvent components. In comparison with simulation-based extrapolation techniques the integral equation method is shown to be capable of quantitatively predicting truncation artefacts at negligible computational overhead.     

The influence of different treatments of electrostatic interactions on the thermodynamics of folding of peptides

Replica exchange molecular dynamics simulations were performed to investigate the effects of different electrostatic treatments on the structure and thermodynamics of a small beta-hairpin forming peptide. Three different electrostatic schemes were considered: regular cutoffs, generalized reaction field (GRF), and particle mesh Ewald (PME), with the peptide modeled using OPLS/AA all-atom force field with explicit TIP3P water. Both the GRF and PME methods yielded results consistent with experiment, with free energy surfaces displaying a single minimum corresponding to the native beta-hairpin structure. In contrast, use of straight cutoffs led to the population of an additional local minimum corresponding to nonhairpin conformations that compete with the formation of the native beta-hairpin at low temperatures. This extra minimum would not be apparent in conventional constant-temperature molecular dynamics simulations run for a few nanoseconds. This result points to the critical need of careful sampling of conformational space to assess the quality of different numerical treatments of long-range forces. While differences emerged in the nature of the unfolded states populated using PME and GRF approaches, simulations on the beta-hairpin forming peptide and on two additional control peptides indicate that the GRF treatment of electrostatics offers a satisfactory compromise between accuracy and computational speed for the identification of low-energy conformations. A GRF-based approach emerges as a viable means for treating larger biological systems that would be prohibitively costly to simulate using PME methods.