A combination of the tree-code and IPS method to simulate large scale systems by molecular dynamics

An IPS/Tree method which is a combination of the isotropic periodic sum (IPS) method and tree-based method was developed for large-scale molecular dynamics simulations, such as biological and polymer systems, that need hundreds of thousands of molecules. The tree-based method uses a hierarchical tree structure to reduce the calculation cost of long-range interactions. IPS/Tree is an efficient method like IPS/DFFT, which is a combination of the IPS method and FFT in calculating large-scale systems that require massively parallel computers. The IPS method has two different versions: IPSn and IPSp. The basic idea is the same expect for the fact that the IPSn method is applied to calculations for point charges, while the IPSp method is used to calculate polar molecules. The concept of the IPS/Tree method is available for both IPSn and IPSp as IPSn/Tree and IPSp/Tree. Even though the accuracy of the Coulomb forces with tree-based method is well known, the accuracy for the combination of the IPS and tree-based methods is unclear. Therefore, in order to evaluate the accuracy of the IPS/Tree method, we performed molecular dynamics simulations for 32,000 bulk water molecules, which contains around 10(5) point charges. IPSn/Tree and IPSp/Tree were both applied to study the interaction calculations of Coulombic forces. The accuracy of the Coulombic forces and other physical properties of bulk water systems were evaluated. The IPSp/Tree method not only has reasonably small error in estimating Coulombic forces but the error was almost the same as the theoretical error of the ordinary tree-based method. These facts show that the algorithm of the tree-based method can be successfully applied to the IPSp method. On the other hand, the IPSn/Tree has a relatively large error, which seems to have been derived from the interaction treatment of the original IPSn method. The self-diffusion and radial distribution functions of water were calculated each by both the IPSn/Tree and IPSp/Tree methods, where both methods showed reasonable agreement with the Ewald method. In conclusion, the IPSp/Tree method is a potentially fast and sufficiently accurate technique for predicting transport coefficients and liquid structures of water in a homogeneous system.     

The Wolf method applied to the liquid-vapor interface of water

The Wolf method for the calculation of electrostatic interactions is applied in a liquid phase and at the liquid-vapor interface of water and its results are compared with those from the Ewald sums method. Molecular dynamics simulations are performed to calculate the radial distribution functions at room temperature. The interface simulations are used to obtain the coexisting densities and surface tension along the coexistence curve. The water model is a flexible version of the extended simple point charge model. The Wolf method gives good structural results, fair coexistence densities, and poor surface tensions as compared with those obtained using the Ewald sums method.     

A long-range electrostatic potential based on the Wolf method charge-neutral condition

Molecular simulations rely heavily on a long range electrostatic Coulomb interaction. The Coulomb potential decays inversely with distance, indicating infinite effective range. In practice, molecular simulations do not directly take into account such an infinite interaction. Therefore, the Ewald, fast multipole, and cutoff methods are frequently used. Although cutoff methods are implemented easily and the calculations are fast, it has been pointed out that they produce serious artifacts. Wolf and coworkers recently discovered one source of the artifacts. They found that when the total charge in a cutoff sphere disappeared, the cutoff error is dramatically suppressed. The Wolf method uses the charge-neutral principle combined with a potential damping that is realized using a complementary error function. To date, many molecular simulation studies have demonstrated the accuracy and reliability of the Wolf method. We propose a novel long-range potential that is constructed only from the charge-neutral condition of the Wolf method without potential damping. We also show that three simulation systems, in which involve liquid sodium-chloride, TIP3P water, and a charged protein in explicit waters with neutralized ions using the new potential, provide accurate statistical and dielectric properties when compared with the particle mesh Ewald method.     

Is the Ewald summation still necessary? Pairwise alternatives to the accepted standard for long-range electrostatics

We investigate pairwise electrostatic interaction methods and show that there are viable computationally efficient (O(N)) alternatives to the Ewald summation for typical modern molecular simulations. These methods are extended from the damped and cutoff-neutralized Coulombic sum originally proposed by Wolf et al. [J. Chem. Phys. 110, 8255 (1999)]. One of these, the damped shifted force method, shows a remarkable ability to reproduce the energetic and dynamic characteristics exhibited by simulations employing lattice summation techniques. Comparisons were performed with this and other pairwise methods against the smooth particle-mesh Ewald summation to see how well they reproduce the energetics and dynamics of a variety of molecular simulations.     

Electrostatic interactions in dissipative particle dynamics using the Ewald sums

The electrostatic interactions in dissipative particle dynamics (DPD) simulations are calculated using the standard Ewald [Ann. Phys. 64, 253 (1921)] sum method. Charge distributions on DPD particles are included to prevent artificial ionic pair formation. This proposal is an alternative method to that introduced recently by Groot [J. Chem. Phys. 118, 11265 (2003)] where the electrostatic field was solved locally on a lattice. The Ewald method is applied to study a bulk electrolyte and polyelectrolyte-surfactant solutions. The structure of the fluid is analyzed through the radial distribution function between charged particles. The results are in good agreement with those reported by Groot for the same systems. We also calculated the radius of gyration of a polyelectrolyte in salt solution as a function of solution pH and degree of ionization of the chain. The radius of gyration increases with the net charge of the polymer in agreement with the trend found in static light scattering experiments of polystyrene sulfonate solutions.     

Comparative study of the effect of tail corrections on surface tension determined by molecular simulation

We report results from a comparative study of the influence of tail corrections on the surface tension of the Lennard-Jones fluid. We find that cutoff-independent surface tensions can be obtained by applying a set of tail corrections recently introduced by Janecek at each step of an interfacial Monte Carlo (MC) or molecular dynamics (MD) simulation. The effect of tail corrections on an alternative methodology for calculating surface tension, the combination of finite-size scaling and grand-canonical transition-matrix Monte Carlo (FSS/GC-TMMC), was also investigated. Using this indirect method, surface tensions were calculated with standard (bulk-fluid) tail corrections and lattice sums, the latter usually considered more accurate but computationally more intensive than the former. With standard tail corrections, we find that the surface tension decreases with increasing cutoff distance, reaching a limiting value corresponding to the maximum cutoff possible, namely half the simulation box length. In contrast, surface tension values obtained with the lattice summation were cutoff-independent. More importantly, these values were equivalent to those surface tension values obtained using standard tail corrections and a cutoff distance of half the box length. We also find that the surface tension values obtained here are in agreement with those found in the literature. Last, we find that surface tension values obtained by MD and FSS/GC-TMMC are in decent agreement so long as the appropriate tail correction schemes are used, and that the relative uncertainties in the surface tensions calculated by MD are generally an order of magnitude greater than those calculated by FSS/GC-TMMC. However, the time required by MD on a single central processing unit is less than that required by FSS/GC-TMMC.     

Computer simulation of sedimentation of ionic systems using the Wolf method

We present computer simulation results for 1:1 and 2:1 electrolyte solutions in the presence of a gravitational field, using the Monte Carlo method in the NVT ensemble for the restrictive primitive model. Coulombic interactions were taken into account comparing the Ewald and Wolf methods. Three variations of Ewald summations were considered: the exact method for slab geometries (EW2D), and the three-dimensional (3D) versions with and without a dipolar correction (EW3DC and EW3D, respectively). The equivalent 3D Wolf protocols were applied under the same conditions (WF3DC and WF3D, respectively). The Wolf and Ewald methods agree accurately in the prediction of several thermodynamic and structural properties for these inhomogeneous systems: excess internal energies, isochoric heath capacities, and density and electrostatic potential profiles. The main advantage using the Wolf method is the significant saving in computing time, which is approximately six times faster than EW3D and EW3DC, and sixty times faster than EW2D.     

A comparison of Coulombic interaction methods in non-equilibrium studies of heat transfer in water

We investigate the impact of the treatment of electrostatic interactions on the heat conduction of liquid water. With this purpose, we report a series of non-equilibrium molecular dynamics computer simulations of the Modified Central Force Model of water. We consider both the Ewald summation approach, which includes the full range of the electrostatic interactions, and the Wolf method, which uses a cutoff to truncate the long range contributions. It is shown that the relaxation of the temperature profiles towards the stationary state solution and the equation of state of the liquid are not affected by the treatment of the electrostatic interactions. However, the truncation of the interactions results in lower internal energy fluxes as well as lower thermal conductivities. We also find that the anomalous increase of the thermal conductivity of water with temperature is reproduced by the different methods considered in this work, showing that this physical behavior is independent of the treatment of the long range electrostatic interactions.     

Ewald summation versus direct summation of shifted-force potentials for the calculation of electrostatic interactions in solids: A quantitative study

The calculated Madelung energies and Madelung forces of the electrostatic interaction for nine crystal structures are reported. The method of direct summation with two different shifted-force potentials is compared to the Ewald summation. There is a considerable difference in the convergence of the energy and the force for the two shifted-force potentials regarding the cutoff radius. The convergence depends not only on the potential itself, but also on the crystal structure. One of the shifted-force potentials used is implemented in the CHARMM force field. The energy calculated with this potential shows a good convergence for small cutoff radii. With the other shifted-force potential, the force shows a better convergence for small cutoff radii. The number of pair interactions for obtaining the Madelung limit using the Ewald summation and the direct summation of a shifted-force potential is also reported. For complex structures like zeolites, the number of relevant pair interactions is smaller using the direct summation of a shifted-force potential. For simple structures such as cesium chloride, the number of significant pair interactions is smaller using the Ewald summation. © 1997 by John Wiley & Sons, Inc.     

Electrostatic potentials in systems periodic in one, two, and three dimensions

We consider the electrostatic potential in a unit cell containing N point charges Q(j) with positions r(j) inside the cell. The cell is replicated periodically in one, two, and three dimensions. The purpose is to give representations for the potential which contain only lattice sums which are absolutely convergent and uniformly convergent in the sampling position r. These representations are derived using variants of the Ewald method and are primarily intended for use in evaluating the accuracy of any algorithm to evaluate electrostatic energies and forces in simulations of dense matter, rather than necessarily for use of themselves in simulations. In reduced dimensionality the Ewald representations can be numerically inefficient and other representations are also provided with careful specification which allows two forms to be used for the potential functions in order to improve numerical performance. These mixed representations may be satisfactory in simulations.     

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.     

Fast multipole method for three-dimensional systems with periodic boundary condition in two directions

We derived a new expression for the electrostatic interaction of three-dimensional charge-neutral systems with two-dimensional periodic boundary conditions (slab geometry) using a fast multipole method (FMM). Contributions from all the image cells are expressed as a sum of real and reciprocal space terms, and a self-interaction term. The reciprocal space contribution consists of two parts: zero and nonzero terms of the absolute value of the reciprocal lattice vector. To test the new expressions, electrostatic interactions were calculated for a randomly placed charge distribution in a cubic box and liquid water produced by molecular dynamics calculation. The accuracy could be controlled by the degree of expansion of the FMM. In the present expression, the computational complexity of the electrostatic interaction of N-particle systems is order N, which is superior to that of the conventional two-dimensional periodic Ewald method for a slab geometry and the particle mesh Ewald method with a large empty space at an interface of the unit cell. © 2020 Wiley Periodicals, Inc.     

Application of Ewald summations to long-range dispersion forces

Results illustrating the effects of using explicit summation terms for the r(-6) dispersion term on the interfacial properties of a Lennard-Jones fluid and SPC/E water are presented. For the Lennard-Jones fluid, we find that the use of long-range summations, even with a short "crossover radius," yields results that are consistent with simulations using large cutoff radii. Simulations of SPC/E water demonstrate that the long-range dispersion forces are of secondary importance to the Coulombic forces. In both cases, we find that the ratio of the box size L( parallel) to the crossover radius r(c) (k) plays an important role in determining the magnitude of the long-range dispersion correction, although its effect is secondary when Coulombic interactions are also present.     

Long-range Lennard-Jones and electrostatic interactions in interfaces: application of the isotropic periodic sum method

Molecular dynamics (MD) simulations of heptane/vapor, hexadecane/vapor, water/vapor, hexadecane/water, and dipalmitoylphosphatidylcholine (DPPC) bilayers and monolayers are analyzed to determine the accuracy of treating long-range interactions in interfaces with the isotropic periodic sum (IPS) method. The method and cutoff (rc) dependences of surface tensions, density profiles, water dipole orientation, and electrostatic potential profiles are used as metrics. The water/vapor, heptane/vapor, and hexadecane/vapor interfaces are accurately and efficiently calculated with 2D IPS (rc=10 A). It is demonstrated that 3D IPS is not practical for any of the interfacial systems studied. However, the hybrid method PME/IPS [Particle Mesh Ewald for electrostatics and 3D IPS for Lennard-Jones (LJ) interactions] provides an efficient way to include both types of long-range forces in simulations of large liquid/vacuum and all liquid/liquid interfaces, including lipid monolayers and bilayers. A previously published pressure-based long-range LJ correction yields results similar to those of PME/IPS for liquid/liquid interfaces. The contributions to surface tension of LJ terms arising from interactions beyond 10 A range from 13 dyn/cm for the hexadecane/vapor interface to approximately 3 dyn/cm for hexadecane/water and DPPC bilayers and monolayers. Surface tensions of alkane/vapor, hexadecane/water, and DPPC monolayers based on the CHARMM lipid force fields agree very well with experiment, whereas surface tensions of the TIP3P and TIP4P-Ew water models underestimate experiment by 16 and 11 dyn/cm, respectively. Dipole potential drops (DeltaPsi) are less sensitive to long-range LJ interactions than surface tensions. However, DeltaPsi for the DPPC bilayer (845+/-3 mV proceeding from water to lipid) and water (547+/-2 mV for TIP4P-Ew and 521+/-3 mV for TIP3P) overestimate experiment by factors of 3 and 5, respectively, and represent expected deficiencies in nonpolarizable force fields.     

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.     

A reoptimization of the five-site water potential (TIP5P) for use with Ewald sums

The five-site transferable interaction potential (TIP5P) for water is most accurate at reproducing experimental data when used with a simple spherical cutoff for the long-ranged electrostatic interactions. When used with other methods for treating long-ranged interactions, the model is considerably less accurate. With small modifications, a new TIP5P-like potential can be made which is very accurate for liquid water when used with Ewald sums, a more physical and increasingly more commonly used method for treating long-ranged electrostatic interactions. The new model demonstrates a density maximum near 4 degrees C, like the TIP5P model, and otherwise is similar to the TIP5P model for thermodynamic, dielectric, and dynamical properties of liquid water over a range of temperatures and densities. An analysis of this and other commonly used water models reveals how the quadrupole moment of a model can influence the dielectric response of liquid water.     

Ewald mesh method for quantum mechanical calculations

The Fourier transform Coulomb (FTC) method has been shown to be effective for the fast and accurate calculation of long-range Coulomb interactions between diffuse (low-energy cutoff) densities in quantum mechanical (QM) systems. In this work, we split the potential of a compact (high-energy cutoff) density into short-range and long-range components, similarly to how point charges are handled in the Ewald mesh methods in molecular mechanics simulations. With this linear scaling QM Ewald mesh method, the long-range potential of compact densities can be represented on the same grid as the diffuse densities that are treated by the FTC method. The new method is accurate and significantly reduces the amount of computational time on short-range interactions, especially when it is compared to the continuous fast multipole method.     

In silico prediction of drug solubility: 1. Free energy of hydration

As a first step in the computational prediction of drug solubility the free energy of hydration, DeltaG*(vw) in TIP4P water has been computed for a data set of 48 drug molecules using the free energy of perturbation method and the optimized potential for liquid simulations all-atom force field. The simulations were performed in two steps, where first the Coulomb and then the Lennard-Jones interactions between the solute and the water molecules were scaled down from full to zero strength to provide physical understanding and simpler predictive models. The results have been interpreted using a theory assuming DeltaG*(vw) = A(MS)gamma + E(LJ) + E(C)/2 where A(MS) is the molecular surface area, gamma is the water-vapor surface tension, and E(LJ) and E(C) are the solute-water Lennard-Jones and Coulomb interaction energies, respectively. It was found that by a proper definition of the molecular surface area our results as well as several results from the literature were found to be in quantitative agreement using the macroscopic surface tension of TIP4P water. This is in contrast to the surface tension for water around a spherical cavity that previously has been shown to be dependent on the size of the cavity up to a radius of approximately 1 nm. The step of scaling down the electrostatic interaction can be represented by linear response theory.     

Isotropic periodic sum: a method for the calculation of long-range interactions

This work presents an accurate and efficient approach to the calculation of long-range interactions for molecular modeling and simulation. This method defines a local region for each particle and describes the remaining region as images of the local region statistically distributed in an isotropic and periodic way, which we call isotropic periodic images. Different from lattice sum methods that sum over discrete lattice images generated by periodic boundary conditions, this method sums over the isotropic periodic images to calculate long-range interactions, and is referred to as the isotropic periodic sum (IPS) method. The IPS method is not a lattice sum method and eliminates the need for a reciprocal space sum. Several analytic solutions of IPS for commonly used potentials are presented. It is demonstrated that the IPS method produces results very similar to that of Ewald summation, but with three major advantages, (1) it eliminates unwanted symmetry artifacts raised from periodic boundary conditions, (2) it can be applied to potentials of any functional form and to fully and partially homogenous systems as well as finite systems, and (3) it is more computationally efficient and can be easily parallelized for multiprocessor computers. Therefore, this method provides a general approach to an efficient calculation of long-range interactions for various kinds of molecular systems.