Integral equation theory for correcting truncation errors in molecular simulations

Various strategies for correcting structural and energetic artefacts of molecular simulations with truncated potentials based on integral equation theory are described and applied to liquid water. The performance of the methods is examined for a range of cutoff distances and different shifted-force potentials. With the recently enhanced damped Coulomb potential (D. Zahn, B. Schilling, S.M. Kast, J. Phys. Chem. B, 106 (2002) 10725), parameterised and corrected by integral equation theory, radial distribution functions and excess internal energy very close to the Ewald simulation limit are obtained from a simulation with a cutoff distance of only 6 Å.     

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.     

A molecular-dynamics simulation study on the dependence of Lennard-Jones gas-liquid phase diagram on the long-range part of the interactions

The particle-transfer molecular-dynamics technique is adopted to construct the Lennard-Jones fluid gas-liquid phase diagram. Detailed study of the dependence of the simulation results on the system size and the cutoff distance is performed to test the validity of the simulation technique. Both the traditional cutoff plus long-range correction (CPC) and Ewald summation methods are used in the simulations to calculate the interactions. In the intermediate range of temperatures, the results with the Ewald summation method are almost the same as those with the CPC method. However, in the range close to the critical point, the results with the CPC method deviate from those with the Ewald summation. Compared with the results obtained via the Ewald summation in a smaller system, simply increasing the system size in the CPC scheme may not give better results.     

Multiple "time step" Monte Carlo simulations: application to charged systems with Ewald summation

Recently, we have proposed an efficient scheme for Monte Carlo simulations, the multiple "time step" Monte Carlo (MTS-MC) [J. Chem. Phys. 117, 8203 (2002)] based on the separation of the potential interactions into two additive parts. In this paper, the structural and thermodynamic properties of the simple point charge water model combined with the Ewald sum are compared for the MTS-MC real-/reciprocal-space split of the Ewald summation and the common Metropolis Monte Carlo method. We report a number of observables as a function of CPU time calculated using MC and MTS-MC. The correlation functions indicate that speedups on the order of 4.5-7.5 can be obtained for systems of 108-500 waters for n=10 splitting parameter.     

Errors in the Ewald summation method with electrostatic layer correction for Coulomb interactions in slab geometry

New theoretical expressions for cut-off errors in 2D reciprocal-space summation of the electrostatic layer correction (ELC) term in energy and forces are derived, and a procedure to determine optimal parameters of the method is proposed. The procedure is tested in numerical calculations for charges distributed uniformly in a cubic box and charges located in two layers near the box basis. The summation errors for conventional Ewald method can be used to find out optimal values of the convergence parameter, and real- and reciprocal-space cut-off radii, whereas the ELC errors give possibility to choose an optimal value of an empty space gap in the simulation box.     

Performance evaluation of the zero‐multipole summation method in modern molecular dynamics software

The zero‐multiple summation method (ZMM) is a cutoff‐based method for calculating electrostatic interactions in molecular dynamics simulations, utilizing an electrostatic neutralization principle as a physical basis. Since the accuracies of the ZMM have been revealed to be sufficient in previous studies, it is highly desirable to clarify its practical performance. In this paper, the performance of the ZMM is compared with that of the smooth particle mesh Ewald method (SPME), where the both methods are implemented in molecular dynamics software package GROMACS. Extensive performance comparisons against a highly optimized, parameter‐tuned SPME implementation are performed for various‐sized water systems and two protein–water systems. We analyze in detail the dependence of the performance on the potential parameters and the number of CPU cores. Even though the ZMM uses a larger cutoff distance than the SPME does, the performance of the ZMM is comparable to or better than that of the SPME. This is because the ZMM does not require a time‐consuming electrostatic convolution and because the ZMM gains short neighbor‐list distances due to the smooth damping feature of the pairwise potential function near the cutoff length. We found, in particular, that the ZMM with quadrupole or octupole cancellation and no damping factor is an excellent candidate for the fast calculation of electrostatic interactions. © 2018 Wiley Periodicals, Inc.     

Amorphous silica modeled with truncated and screened Coulomb interactions: a molecular dynamics simulation study

We show that finite-range alternatives to the standard long-range pair potential for silica by van Beest et al. [Phys. Rev. Lett. 64, 1955 (1990)] might be used in molecular dynamics simulations. We study two such models that can be efficiently simulated since no Ewald summation is required. We first consider the Wolf method, where the Coulomb interactions are truncated at a cutoff distance rc such that the requirement of charge neutrality holds. Various static and dynamic quantities are computed and compared to results from simulations using Ewald summations. We find very good agreement for rc approximately 10 A. For lower values of rc, the long-range structure is affected which is accompanied by a slight acceleration of dynamic properties. In a second approach, the Coulomb interaction is replaced by an effective Yukawa interaction with two new parameters determined by a force fitting procedure. The same trend as for the Wolf method is seen. However, slightly larger cutoffs have to be used in order to obtain the same accuracy with respect to static and dynamic quantities as for the Wolf method.     

Fragment-based quantum mechanical methods for periodic systems with Ewald summation and mean image charge convention for long-range electrostatic interactions

We describe an Ewald-summation method to incorporate long-range electrostatic interactions into fragment-based electronic structure methods for periodic systems. The present method is an extension of the particle-mesh Ewald technique for combined quantum mechanical and molecular mechanical (QM/MM) calculations, and it has been implemented into the explicit polarization (X-Pol) potential to illustrate the computational details. As in the QM/MM-Ewald method, the X-Pol-Ewald approach is a linear-scaling electrostatic method, in which the short-range electrostatic interactions are determined explicitly in real space and the long-range Ewald pair potential is incorporated into the Fock matrix as a correction. To avoid the time-consuming Fock matrix update during the self-consistent field procedure, a mean image charge (MIC) approximation is introduced, in which the running average with a user-chosen correlation time is used to represent the long-range electrostatic correction as an average effect. Test simulations on liquid water show that the present X-Pol-Ewald method takes about 25% more CPU time than the usual X-Pol method using spherical cutoff, whereas the use of the MIC approximation reduces the extra costs for long-range electrostatic interactions by 15%. The present X-Pol-Ewald method provides a general procedure for incorporating long-range electrostatic effects into fragment-based electronic structure methods for treating biomolecular and condensed-phase systems under periodic boundary conditions.     

A simple algorithm to accelerate the computation of non-bonded interactions in cell-based molecular dynamics simulations

Cell lists are ubiquitous in molecular dynamics simulations--be it for the direct computation of short-range inter-atomic potentials, the short-range direct part of a long-range interaction or for the periodic construction of Verlet lists. The conventional approach to computing pairwise interactions using cell lists leads to a large number of unnecessary interparticle distance calculations. In this paper, an algorithm is presented which reduces the number of spurious distance calculations by first sorting the particles along the cell pair axis and then only interacting two particles if their distance along the axis is smaller than the cutoff distance of the interaction. This approach is shown to be more efficient than the conventional approach and similar approaches using smaller cells.     

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.     

A new method applicable to study solid compounds with multiple polyhedral structures

A new direct summation method, named as polyhedron method, is proposed to calculate Madelung energy. This method calculates sums of electrostatic interactions over sets of neutral polyhedron unit pairs rather than conventional ion pairs; this gives Madelung constant in a matrix. With robustly rapid convergence, polyhedron method is generally applicable for complex compounds containing multiple polyhedral building‐blocks and numerical polyhedral connection modes. The matrical analysis suggests face‐sharing between octahedral pairs and edge‐sharing between tetrahedral pairs can be electrostatically stable, against Pauling's third rule. Further, the matrical calculation of Madelung energies offers a unique advantage to evaluate enormous configurations of cation distributions in a given lattice in a high‐throughput manner. That is applicable to study solid solution composites, polymorphism, and defect structures, including but not limited to intermediate phase of delithiated cathode compounds, charge order or antisite defects, and extensively magnetic order. © 2016 Wiley Periodicals, Inc.     

Some comments and corrections regarding the calculation of electrostatic potential derivatives using the Ewald summation technique

A review of the literature on the calculation of electrostatic potentials, fields, and field gradients in systems consisting of charges and dipoles using the Ewald summation technique is presented. Discrepancies between the previous formulas are highlighted, and an error in the derivation of the reciprocal contributions to the electrostatic field and field gradient is corrected. The new formulas for the field and field gradient are shown to exhibit a termwise identity with the ones for the electrostatic energy.     

A modified TIP3P water potential for simulation with Ewald summation

The charges and Lennard-Jones parameters of the TIP3P water potential have been modified to improve its performance under the common condition for molecular dynamics simulations of using Ewald summation in lieu of relatively short nonbonded truncation schemes. These parameters were optimized under the condition that the hydrogen atoms do not have Lennard-Jones parameters, thus making the model independent of the combining rules used for the calculation of nonbonded, heteroatomic interaction energies, and limiting the number of Lennard-Jones calculations required. Under these conditions, this model provides accurate density (rho = 0.997 g/ml) and heat of vaporization (DeltaH(vap) = 10.53 kcal/mol) at 25 degrees C and 1 atm, but also provides improved structure in the second peak of the O-O radial distribution function and improved values for the dielectric constant (epsilon(0) = 89) and the diffusion coefficient (D = 4.0 x 10(-5) cm(2)/s) relative to the original parametrization. Like the original parameterization, however, this model does not show a temperature density maximum. Several similar models are considered with the additional constraint of trying to match the performance of the optimized potentials for liquid simulation atom force field to that obtained when using the simulation conditions under which it was originally designed, but no model was entirely satisfactory in reproducing the relative difference in free energies of hydration between the model compounds, phenol and benzene. Finally, a model that incorporates a long-range correction for truncated Lennard-Jones interactions is presented, which provides a very accurate dielectric constant (epsilon(0) = 76), however, the improvement in this estimate is on the same order as the uncertainty in the calculation.     

Ewald summation and multiple time step methods for molecular dynamics simulation of biological molecules

Methods by which to determine conditions for a molecular dynamics (MD) simulation of biological molecules were investigated. Derivation of the optimal parameters of the Ewald summation was described so as to give same precision to the real space, the reciprocal space summations and the van der Waals interaction. Later, the procedure by which to determine the condition of the multiple time step method by RESPA (REference System Propagator Algorithm; Tuckerman et al., 1992, J. Chem. Phys., 97, 1990) was described as exemplified by MD simulations of a solvated β-sheet peptide. The conservation of the total energy in a microcanonical ensemble was measured to investigate the stability of the simulation conditions. The most feasible respective combinations of the time steps were: 0.25 fs for bond, angle and torsion interactions; 2 fs for van der Waals interaction and Ewald real-space summation; and 4 fs for Ewald reciprocal-space summation. Though it retained an acceptable accuracy, this condition accelerated the simulation ten-fold compared to that in which a simple velocity-Verlet method with a time step of 0.25 fs was used. The update of the correction term due to excluded neighbors was then investigated. Better results were obtained when the correction was updated with the real-space than when it was updated with the reciprocal-space summation. Finally, an MD simulation as long as 50 ps performed under the optimal Ewald and RESPA parameters was thus determined. The trajectory showed a good stability, indicating the feasibility of the parameters.     

Pairwise long-range compensation for strongly ionic systems

We propose a pairwise compensation method for long-range electrostatics, as an alternative to traditional infinite lattice sums. The approach represents the third generation in a series beginning with the shifted potential corresponding to counterions surrounding a cutoff sphere. That simple charge compensation scheme resulted in pairwise potentials that are continuous at the cutoff, but forces that are not. A second-generation approach modified both the potential and the force such that both are continuous at the cutoff. Here we introduce another layer of softening such that the derivative of the force is also continuous at the cutoff. In strongly ionic liquids, this extension removes structural artifacts associated with the earlier pairwise compensation schemes, and provides results that compare well with Ewald sums.     

Ab initio based polarizable force field generation and application to liquid silica and magnesia

We extend the program potfit, which generates effective atomic interaction potentials from ab initio data, to electrostatic interactions and induced dipoles. The potential parametrization algorithm uses the Wolf direct, pairwise summation method with spherical truncation. The polarizability of oxygen atoms is modeled with the Tangney-Scandolo interatomic force field approach. Due to the Wolf summation, the computational effort in simulation scales linearly in the number of particles, despite the presence of electrostatic interactions. Thus, this model allows to perform large-scale molecular dynamics simulations of metal oxides with realistic potentials. Details of the implementation are given, and the generation of potentials for SiO(2) and MgO is demonstrated. The approach is validated by simulations of microstructural, thermodynamic, and vibrational properties of liquid silica and magnesia.     

The ENUF method—Ewald summation based on nonuniform fast Fourier transform: Implementation,parallelization, and application

Computer simulations of model systems are widely used to explore striking phenomena in promising applications spanning from physics, chemistry, biology, to materials science and engineering. The long range electrostatic interactions between charged particles constitute a prominent factor in determining structures and states of model systems. How to efficiently calculate electrostatic interactions in simulation systems subjected to partial or full periodic boundary conditions has been a grand challenging task. In the past decades, a large variety of computational schemes has been proposed, among which the Ewald summation method is the most reliable route to accurately deal with electrostatic interactions between charged particles in simulation systems. In addition, extensive efforts have been done to improve computational efficiencies of the Ewald summation based methods. Representative examples are approaches based on cutoffs, reaction fields, multi-poles, multi-grids, and particle-mesh schemes. We sketched an ENUF method, an abbreviation for the Ewald summation method based on the nonuniform fast Fourier transform technique, and have implemented this method in particle-based simulation packages to calculate electrostatic energies and forces at micro- and mesoscopic levels. Extensive computational studies of conformational properties of polyelectrolytes, dendrimer-membrane complexes, and ionic fluids demonstrated that the ENUF method and its derivatives conserve both energy and momentum to floating point accuracy, and exhibit a computational complexity of with optimal physical parameters. These ENUF based methods are attractive alternatives in molecular simulations where high accuracy and efficiency of simulation methods are needed to accelerate calculations of electrostatic interactions at extended spatiotemporal scales.     

Molecular dynamics simulation of trp-repressor/operator complex: analysis of hydrogen bond patterns of protein–DNA interaction

Molecular dynamics simulation of Escherichia coli trp-repressor/operator complex was performed to elucidate protein–DNA interactions in solution for 800 ps on special-purpose computer MD-GRAPE. The Ewald summation method was employed to treat the electrostatic interaction without cutoff. DNA kept stable conformation in comparison with the result of the conventional cutoff method. Thus, the trajectories obtained were used to analyze the protein–DNA interaction and to understand the role of dynamics of water molecules forming sequence specific recognition interface. The dynamical cross-correlation map showed a significant positive correlation between the helix-turn-helix DNA-binding motifs and the major grooves of operator DNA. The extensive contact surface was stable during the simulation. Most of the contacts consisted of direct interactions between phosphates of DNA and the protein, but several water-mediated polar contacts were also observed. These water-mediated interactions, which were also seen in the crystal structure (Z. Otwinowski, et al., Nature, 335 (1998) 321) emerged spontaneously from the randomized initial configuration of the solvent. This result suggests the importance of the water-mediated interaction in specific recognition of DNA by the trp-repressor, consistent with X-ray structural information.     

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.     

Effect of different treatments of long-range interactions and sampling conditions in molecular dynamic simulations of rhodopsin embedded in a dipalmitoyl phosphatidylcholine bilayer

The present study analyzes the effect of the simulation conditions on the results of molecular dynamics simulations of G-protein coupled receptors (GPCRs) performed with an explicit lipid bilayer. Accordingly, the present work reports the analysis of different simulations of bovine rhodopsin embedded in a dipalmitoyl phosphatidylcholine (DPPC) lipid bilayer using two different sampling conditions and two different approaches for the treatment of long-range electrostatic interactions. Specifically, sampling was carried out either by using the statistical ensembles NVT or NPT (constant number of atoms, a pressure of 1 atm in all directions and fixed temperature), and the electrostatic interactions were treated either by using a twin-cutoff, or the particle mesh Ewald summation method (PME). The results of the present study suggest that the use of the NPT ensemble in combination with the PME method provide more realistic simulations. The use of NPT during the equilibration avoids the need of an a priori estimation of the box dimensions, giving the correct area per lipid. However, once the system is equilibrated, the simulations are irrespective of the sampling conditions used. The use of an electrostatic cutoff induces artifacts on both lipid thickness and the ion distribution, but has no direct effect on the protein and water molecules.