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.     

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.     

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 Å.     

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.     

Cutoff radius effect of the isotropic periodic sum and Wolf method in liquid-vapor interfaces of water

As a more economical but similarly accurate computation method than the Ewald sum, the isotropic periodic sum (IPS) method for nonpolar molecules (IPSn) and polar molecules (IPSp), along with the Wolf method are of interest, but the cutoff radius dependence is an important issue. To evaluate the cutoff radius effect of the three methods, a water-vapor interfacial system has been studied by molecular dynamics. The Wolf method can produce adequate results for surface tension compared to that of the Ewald sum (within 2.9%) at a long enough cutoff radius, r(c). However, the estimation of the electrostatic potential profile and dipole orientational function is poor. The Wolf method cannot estimate electrostatic configuration at r(c) ≤ L(z)∕2 (L(z) is the longest lattice of the system). We have found that the convergence of the surface tension and the electrostatic configuration of the IPSn method is faster than that of the IPSp method. Moreover, the IPSn method is most accurate among the three methods for the same cutoff radius. Furthermore, the behavior of the surface tension against the cutoff radius shows a greater difference for the IPSn and IPSp method. The surface tension of the IPSp method fluctuates and presents a similar result to that of the Ewald sum, but the surface tension for the IPSn method greatly deviates near r(c) = L(z)∕3. The cause of this deviation is the difference between the interfacial configuration of the water surface and the cutoff treatment of the IPS method. The deviation becomes insignificant far from r(c) = L(z)∕3. In spite of this shortcoming, the IPSn method gives the most accurate result in estimating the surface tension at r(c) = L(z)∕2. From all the results in this work, the IPSn and IPSp method have been found to be more accurate than the Wolf method. In conclusion, the surface tension and structure of water-vapor interface can be calculated by the IPSn method when r(c) is greater than or equal to the longest lattice of the system. The IPSp method and the Wolf method require a longer cutoff radius than the longest lattice of the system to estimate interfacial properties.     

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.     

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.     

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.     

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.     

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.     

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.     

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.     

Application of MDGRAPE-3, a special purpose board for molecular dynamics simulations, to periodic biomolecular systems

We describe the application of a special purpose board for molecular dynamics simulations, named MDGRAPE-3, to the problem of simulating periodic bio-molecular systems. MDGRAPE-3 is the latest board in a series of hardware accelerators designed to calculate the nonbonding long-range interactions much more rapidly than normal processors. So far, MDGRAPEs were mainly applied to isolated systems, where very many nonbonded interactions were calculated without any distance cutoff. However, in order to regulate the density and pressure during simulations of membrane embedded protein systems, one has to evaluate interactions under periodic boundary conditions. For this purpose, we implemented the Particle-Mesh Ewald (PME) method, and its approximation with distance cutoffs and charge neutrality as proposed by Wolf et al., using MDGRAPE-3. When the two methods were applied to simulations of two periodic biomolecular systems, a single MDGRAPE-3 achieved 30-40 times faster computation times than a single conventional processor did in the both cases. Both methods are shown to have the same molecular structures and dynamics of the systems.     

A new approach for efficient simulation of Coulomb interactions in ionic fluids

We propose a simplified version of local molecular field (LMF) theory to treat Coulomb interactions in simulations of ionic fluids. LMF theory relies on splitting the Coulomb potential into a short-ranged part that combines with other short-ranged core interactions and is simulated explicitly. The averaged effects of the remaining long-ranged part are taken into account through a self-consistently determined effective external field. The theory contains an adjustable length parameter sigma that specifies the cutoff distance for the short-ranged interaction. This can be chosen to minimize the errors resulting from the mean-field treatment of the complementary long-ranged part. Here we suggest that in many cases an accurate approximation to the effective field can be obtained directly from the equilibrium charge density given by the Debye theory of screening, thus eliminating the need for a self-consistent treatment. In the limit sigma-->0, this assumption reduces to the classical Debye approximation. We examine the numerical performance of this approximation for a simple model of a symmetric ionic mixture. Our results for thermodynamic and structural properties of uniform ionic mixtures agree well with similar results of Ewald simulations of the full ionic system. In addition, we have used the simplified theory in a grand-canonical simulation of a nonuniform ionic mixture where an ion has been fixed at the origin. Simulations using short-ranged truncations of the Coulomb interactions alone do not satisfy the exact condition of complete screening of the fixed ion, but this condition is recovered when the effective field is taken into account. We argue that this simplified approach can also be used in the simulations of more complex nonuniform systems.     

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.     

Molecular dynamics scheme for precise estimation of electrostatic interaction via zero-dipole summation principle

We propose a novel idea, zero-dipole summation, for evaluating the electrostatic energy of a classical particle system, and have composed an algorithm for effectively utilizing the idea for molecular dynamics. It conceptually prevents the nonzero-charge and nonzero-dipole states artificially generated by a simple cutoff truncation. The resulting energy formula is nevertheless represented by a simple pairwise function sum, which enables facile application to high-performance computation. By following a heuristic approach to derive the current electrostatic energy formula, we developed an axiomatic approach to construct the method consistently. Explorations of the theoretical details of our method revealed the structure of the generated error, and we analyzed it by comparisons with other methods. A numerical simulation using liquid sodium chloride confirmed that the current method with a small damping factor yielded sufficient accuracy with a practical cutoff distance region. The current energy function also conducts stable numerical integration in a liquid MD simulation. Our method is an extension of the charge neutralized summation developed by Wolf et al. [J. Chem. Phys. 110, 8254 (1999)]. Furthermore, we found that the current method becomes a generalization of the preaveraged potential method proposed by Yakub and Ronchi [J. Chem. Phys. 119, 11556 (2003)], which is based on a viewpoint different from the neutrality. The current study presents these relationships and suggests possibilities for their further applications.     

液态RbCl的态密度

给出了在分子动力学模拟基础上Fumi-Tosi势离子液体的正则模式分析方法,用Fumi-Tosi势(包括长程势)代替Lennard-Jones势,并且用等效Coulomb势处理长程Coulomb作用.讨论了Hessian矩阵元的计算方法和Hessian矩阵特征值的计算方法.计算实践表明，取用余误差函数形式的等效库仑势,可以合理地得到Hessian矩阵和态密度.液态RbCl中构型平均态密度的数值结果表明，液态RbCl的态密度表现出与Lennard-Jones液体的态密度相仿的特点.     

Gaussian split Ewald: A fast Ewald mesh method for molecular simulation

Gaussian split Ewald (GSE) is a versatile Ewald mesh method that is fast and accurate when used with both real-space and k-space Poisson solvers. While real-space methods are known to be asymptotically superior to k-space methods in terms of both computational cost and parallelization efficiency, k-space methods such as smooth particle-mesh Ewald (SPME) have thus far remained dominant because they have been more efficient than existing real-space methods for simulations of typical systems in the size range of current practical interest. Real-space GSE, however, is approximately a factor of 2 faster than previously described real-space Ewald methods for the level of force accuracy typically required in biomolecular simulations, and is competitive with leading k-space methods even for systems of moderate size. Alternatively, GSE may be combined with a k-space Poisson solver, providing a conveniently tunable k-space method that performs comparably to SPME. The GSE method follows naturally from a uniform framework that we introduce to concisely describe the differences between existing Ewald mesh methods.     

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.     

Efficient lookup table using a linear function of inverse distance squared

The major bottleneck in molecular dynamics (MD) simulations of biomolecules exist in the calculation of pairwise nonbonded interactions like Lennard‐Jones and long‐range electrostatic interactions. Particle‐mesh Ewald (PME) method is able to evaluate long‐range electrostatic interactions accurately and quickly during MD simulation. However, the evaluation of energy and gradient includes time‐consuming inverse square roots and complementary error functions. To avoid such time‐consuming operations while keeping accuracy, we propose a new lookup table for short‐range interaction in PME by defining energy and gradient as a linear function of inverse distance squared. In our lookup table approach, densities of table points are inversely proportional to squared pair distances, enabling accurate evaluation of energy and gradient at small pair distances. Regardless of the inverse operation here, the new lookup table scheme allows fast pairwise nonbonded calculations owing to efficient usage of cache memory. © 2013 Wiley Periodicals, Inc.