Communication: Shifted forces in molecular dynamics

Simulations involving the Lennard-Jones potential usually employ a cutoff at r = 2.5σ. This communication investigates the possibility of reducing the cutoff. Two different cutoff implementations are compared, the standard shifted potential cutoff and the less commonly used shifted forces cutoff. The first has correct forces below the cutoff, whereas the shifted forces cutoff modifies Newton's equations at all distances. The latter is nevertheless superior; we find that for most purposes realistic simulations may be obtained using a shifted forces cutoff at r = 1.5σ, even though the pair force is here 30 times larger than at r = 2.5σ.     

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.     

A coarse-grained protein-protein potential derived from an all-atom force field

In order to study protein-protein nonbonded interactions, we present the development of a new reduced protein model that represents each amino acid residue with one to three coarse grains, whose physical properties are derived in a consistent bottom-up procedure from the higher-resolution all-atom AMBER force field. The resulting potential energy function is pairwise additive and includes distinct van-der-Waals and Coulombic terms. The van-der-Waals effective interactions are deduced from preliminary molecular dynamics simulations of all possible amino acid homodimers. They are best represented by a soft 1/r6 repulsion and a Gaussian attraction, with parameters obeying Lorentz-Berthelot mixing rules. For the Coulombic interaction, coarse grain charges are optimized for each separate protein in order to best represent the all-atom electrostatic potential outside the protein core. This approach leaves the possibility of using any implicit solvent model to describe solvation effects and electrostatic screening. The coarse-grained force field is tested carefully for a small homodimeric complex, the magainin. It is shown to reproduce satisfactorily the specificity of the all-atom underlying potential, in particular within a PB/SA solvation model. The coarse-grained potential is applied to the redocking prediction of three different protein-protein complexes: the magainin dimer, the barnase-barstar, and the trypsin-BPTI complexes. It is shown to provide per se an efficient and discriminating scoring energy function for the protein-protein docking problem that remains pertinent at both the global and refinement stage.     

Iso-energy cutoff for the calculation of interionic potential of mean force in water

An iso-energy cutoff scheme is introduced for the calculation of the potential of mean force between two ions in water. The cutoff criterion is based on the optimal interaction of the water dipole with the ion pair, for which analytical expressions are derived. Formulas are also derived to characterize the solvent reorganization contribution to the potential of mean force. Treatment of the contributions from waters outside the cutoff is also discussed. © 1994 John Wiley & Sons, Inc.     

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.     

New spherical-cutoff methods for long-range forces in macromolecular simulation

New atom- and group-based spherical-cutoff methods have been developed for the treatment of nonbonded interactions in molecular dynamics (MD) simulation. A new atom-based method, force switching, leaves short-range forces unaltered by adding a constant to the potential energy, switching forces smoothly to zero over a specified range. A simple improvement to group-based cutoffs is presented: Switched group-shifting shifts the group–group potential energy by a constant before being switched smoothly to zero. Also introduced are generalizations of atom-based force shifting, which adds a constant to the Coulomb force between two charges. These new approaches are compared to existing methods by evaluating the energy of a model hydrogen-bonding system consisting of two N-methyl acetamide molecules and by full MD simulation. Thirty-five 150 ps simulations of carboxymyoglobin (MbCO) hydrated by 350 water molecules indicate that the new methods and atom-based shifting are each able to approximate no-cutoff results when a cutoff at or beyond 12 Å is used. However, atom-based potential-energy switching and truncation unacceptably contaminate group–group electrostatic interactions. Group-based potential truncation should not be used in the presence of explicit water or other mobile electrostatic dipoles because energy is not a state function with this method, resulting in severe heating (about 4 K/ps in the simulations of hydrated MbCO). The distance-dependent dielectric (? ∝? r) is found to alter the temperature dependence of protein dynamics, suppressing anharmonic motion at high temperatures. Force switching and force shifting are the best atom-based spherical cutoffs, whereas switched group-shifting is the preferred group-based method. To achieve realistic simulations, increasing the cutoff distance from 7.5 to 12 Å or beyond is much more important than the differences among the three best cutoff methods. © 1994 by John Wiley & Sons, Inc.

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

     

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.     

Monte Carlo versus molecular dynamics simulations in heterogeneous systems: an application to the n-pentane liquid-vapor interface

The Monte Carlo (MC) and molecular dynamics (MD) methodologies are now well established for computing equilibrium properties in homogeneous fluids. This is not yet the case for the direct simulation of two-phase systems, which exhibit nonuniformity of the density distribution across the interface. We have performed direct MC and MD simulations of the liquid-gas interface of n-pentane using a standard force-field model. We obtained density and pressure components profiles along the direction normal to the interface that can be very different, depending on the truncation and long range correction strategies. We discuss the influence on predicted properties of different potential truncation schemes implemented in both MC and MD simulations. We show that the MD and MC profiles can be made in agreement by using a Lennard-Jones potential truncated via a polynomial function that makes the first and second derivatives of the potential continuous at the cutoff distance. In this case however, the predicted thermodynamic properties (phase envelope, surface tension) deviate from experiments, because of the changes made in the potential. A further readjustment of the potential parameters is needed if one wants to use this method. We conclude that a straightforward use of bulk phase force fields in MD simulations may lead to some physical inconsistencies when computing interfacial properties.     

Simulation of crack propagation in alumina with ab initio based polarizable force field

We present an effective atomic interaction potential for crystalline α-Al(2)O(3) generated by the program potfit. The Wolf direct, pairwise summation method with spherical truncation is used for electrostatic interactions. The polarizability of oxygen atoms is included by use of the Tangney-Scandolo interatomic force field approach. The potential is optimized to reproduce the forces, energies, and stresses in relaxed and strained configurations as well as {0001}, {1010}, and {1120} surfaces of Al(2)O(3). Details of the force field generation are given, and its validation is demonstrated. We apply the developed potential to investigate crack propagation in α-Al(2)O(3) single crystals.     

Many-body force field models based solely on pairwise Coulomb screening do not simultaneously reproduce correct gas-phase and condensed-phase polarizability limits

It is demonstrated that many-body force field models based solely on pairwise Coulomb screening cannot simultaneously reproduce both gas-phase and condensed-phase polarizability limits. Several many-body force field model forms are tested and compared with basis set-corrected ab initio results for a series of bifurcated water chains. Models are parameterized to reproduce the ab initio polarizability of an isolated water molecule, and pairwise damping functions are set to reproduce the polarizability of a water dimer as a function of dimer separation. When these models are applied to extended water chains, the polarization is over-predicted, and this over-polarization increased as a function of the overlap of molecular orbitals as the chains are compressed. This suggests that polarizable models based solely on pairwise Coulomb screening have some limitations, and that coupling with non-classical many-body effects, in particular exchange terms, may be important.     

Trimer based polarization as a multibody molecular model. Application to hydrogen fluoride

A molecular modeling approach is introduced as a way to treat multibody (more than two molecules) contributions to the intermolecular potential. There are two key features to the method. First, it employs polarizable electrostatics on the molecules, but converges the charges and fields for only three molecules at a time, taken separately for all trimers (three molecules falling within a cutoff distance) in the system. This feature introduces significant computational savings when applied in Monte Carlo simulation (in comparison to a full N-body polarization treatment), as movement of a single molecule does not require re-converging of the polarization of all molecules, and it achieves this without approximations that cause the value of the energy to depend on the history of the simulation. Second, the approach defines the polarization energy in excess of the pairwise contribution, meaning that the trimer energy has subtracted from it the sum of the energies obtained by converging the polarization of each molecule pair in the trimer. This feature is advantageous because it removes the need (often found in polarizable models) to stiffen inappropriately the repulsive part of the pair potential. The polarization contribution is thus a purely three-body potential. The approach is applied to model hydrogen fluoride, which in experiments exhibits unusual properties that have proven difficult to capture well by molecular models. The new HF model is shown to be much more successful than previous modeling efforts in obtaining agreement with a broad range of experimental data (volumetric properties, heat effects, molecular structure, and vapor-liquid equilibria).     

Melting line of the Lennard-Jones system, infinite size, and full potential

Literature estimates of the melting curve of the Lennard-Jones system vary by as much as 10%. The origin of such discrepancies remains unclear. We present precise values for the Lennard-Jones melting temperature, and we examine possible sources of systematic errors in the prediction of melting points, including finite-size and interaction-cutoff effects. A hypothetical thermodynamic integration path is used to find the relative free energies of the solid and liquid phases, for various system sizes, at constant cutoff radius. The solid-liquid relative free energy and melting temperature scale linearly as the inverse of the number of particles, and it is shown that finite-size effects can account for deviations in the melting temperature (from the infinite-size limit) of up to 5%. An extended-ensemble density-of-states method is used to determine free energy changes for each phase as a continuous function of the cutoff radius. The resulting melting temperature predictions exhibit an oscillatory behavior as the cutoff radius is increased. Deviations in the melting temperature (from the full potential limit) arising from a finite cutoff radius are shown to be of comparable magnitude as those resulting from finite-size effects. This method is used to identify melting temperatures at five different pressures, for the infinite-size and full potential Lennard-Jones system. We use our simulation results as references to connect the Lennard-Jones solid equation of state of van der Hoef with the Lennard-Jones fluid equation of state of Johnson. Once the references are applied the two equations of state are used to identify a melting curve. An empirical equation that fits this melting curve is provided. We also report a reduced triple point temperature T(tr)=0.694.     

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.     

Ewald summation of multipolar interactions at an arbitrary order on a two-dimensional lattice

There exist problems in condensed-matter theory that require evaluating infinite Bloch sums of multipolar potential r^?l?1Y_l,m(θ,φ) on a periodic lattice. For an arbitrary multipolar order l, tractable formulas are given for summing such interactions on a two-dimensional Bravais lattice and evaluating their Bloch sums at a point outside as well as inside the plane of the lattice. The approach used is the Ewald method, which consists of separating the original series in rapidly converging sums in reciprocal and real spaces. Computational aspects of the present formulation are briefly reviewed. © 1993 John Wiley & Sons, Inc.     

Aqueous solutions containing amino acids and peptides. Part 19: The enthalpic coefficients for the interactions of N-acetylsarcosinamide with 2-(N-acetylamino)acyl amides at 25°C

Enthalpies of dilution both of solutions of N-acetylsarcosinamide and of ternary solutions equimolal in N-acetylsarcosinamide and N-acetylglycinamide, N-acetyl-L-alaninamide, N-acetyl-L-valinamide or N-acetyl-L-leucinamide have been determined by a microcalorimetric method. The results were employed to calculate the pairwise enthalpic coefficients for both homotactic (like-like) and heterotactic (like-unlike) solute interactions. These pairwise interaction coefficients have been analyzed by means of a group additivity approach and some comments on the utility of this, when applied to such systems, are made.     

Molecular mechanics and the deformation of macromolecules: The use of a very short cutoff combined with a quadratic approximation

An approximation to the molecular mechanical treatment of structural deformations of macromolecules is presented. The method is based on a partitioning of the conformational energy into three parts. The first part is covered by the condensed potential functions which absorb many short-range nonbonding interactions. The second part consists of a few nonbonded interactions below a very short cutoff radius of 4 Å. The third part, consisting of the vast majority of pairwise interactions, is approximated by a quadratic expression confined to a subspace of the conformational space. A detailed computational example on LH-RH, including an analysis of the errors resulting from other conventional approximation methods, is given. A comparison to the conventional cutoff approximation used in x-ray refinement delivers a speedup factor of at least two orders of magnitude.     

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.     

Empirical force fields for complex hydrated calcio-silicate layered materials

The use of empirical force fields is now a standard approach in predicting the properties of hydrated oxides which are omnipresent in both natural and engineering applications. Transferability of force fields to analogous hydrated oxides without rigorous investigations may result in misleading property predictions. Herein, we focus on two common empirical force fields, the simple point charge ClayFF potential and the core-shell potential to study tobermorite minerals, the most prominent family of Calcium-Silicate-Hydrates that are complex hydrated oxides. We benchmark the predictive capabilities of these force fields against first principles results. While the structural information seem to be in close agreement with DFT results, we find that for higher order properties such as elastic constants, the core-shell potential quantitatively improves upon the simple point charge model, and shows a larger degree of transferability to complex materials. In return, to remedy the deficiencies of the simple point charge potential for hydrated calcio-silicates, we suggest using both structural data and elasticity data for potential calibration, a new force field potential, CSH-FF. This re-parameterized version of ClayFF is then applied to simulating an atomistic model of cement (Pellenq et al., PNAS, 2009). We demonstrate that this force field improves the predictive capabilities of ClayFF, being considerably less computational intensive than the core-shell model.     

Dynamics of benzene molecules situated in metal-organic frameworks

In this paper, we investigate the gyroscopic motion of a benzene molecule C₆H₆, which comprises an inner carbon ring and an outer hydrogen ring, and is suspended rigidly inside a metal-organic framework. The metal-organic framework provides a sterically unhindered environment and an electronic barrier for the benzene molecule. We model such gyroscopic motion from the inter-molecular interactions between the benzene ring and the metal-organic framework by both the Columbic force and the van der Waals force. We also capture additional molecular interactions, for example due to sterical compensations arising from the carboxylate ligands between the benzene molecule and the framework, by incorporating an extra empirical energy into the total molecular energy. To obtain a continuous approximation to the total energy of such a complicated atomic system, we assume that the atoms of the metal-organic framework can be smeared over the surface of a cylinder, while those for the benzene molecule are smeared over the contour line of the molecule. We then approximate the pairwise molecular energy between the molecules by performing line and surface integrals. We firstly investigate the freely suspended benzene molecule inside the framework and find that our theoretical results admit a two-fold flipping, with the possible maximum rotational frequency reaching the terahertz regime, and gigahertz frequencies at room temperature. We also show that the electrostatic interaction and the thermal energy dominate the gyroscopic motion of the benzene molecule, and we deduce that the extra energy term could possibly reduce the rotational frequency of the rigidly suspended benzene molecule from gigahertz to megahertz frequencies at room temperature, and even lower frequencies might be obtained when the strength of these interactions increases.