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.     

On mesh-based Ewald methods: optimal parameters for two differentiation schemes

The particle-particle particle-mesh Ewald method for the treatment of long-range electrostatics under periodic boundary conditions is reviewed. The optimal Green's function for exact (real-space differentiation), which differs from that for reciprocal-space differentiation, is given. Simple analytic formulas are given to determine the optimal Ewald screening parameter given a differentiation scheme, a real-space cutoff, a mesh spacing, and an assignment order. Simulations of liquid water are performed to examine the effect of the accuracy of the electrostatic forces on calculation of the static dielectric constant. A target dimensionless root-mean-square error of 10(-4) is sufficient to obtain a well-converged estimate of the dielectric constant.     

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.     

Prediction, fine tuning, and temperature extrapolation of a vapor liquid equilibrium using COSMOtherm

As our entry for the third industrial fluid property simulation challenge, the COSMO-RS method in its COSMOtherm implementation has been used to predict the vapor liquid equilibrium (VLE) of ethanol and 1,1,1,2,3,3,3-heptafluoropropane (HFC-227ea) at 343.13 K. The absolute prediction of the VLE without any system specific adjustments already yields a reasonable agreement with the experimental data provided for the binary mixture at 283.17 K. Because the special emphasis of the this challenge is state conditions transferability, we also considered two special ways of fine tuning to the experimental data provided for this VLE at 283.17 K. The first way of fine tuning was by fitting of correction charges, the second was by adjustment of a single van der Waals (vdW) interaction parameter. Since vdW parameters used in COSMOtherm are the weakest part of the COSMO-RS parameterization, the fine tuning of the vdW parameters is considered as physically most plausible. Therefore our final prediction of the VLE ethanol and HFC-227ea at 343.13 K is based on the vdW fine tuning.     

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.     

Direct numerical simulation of aerosol coagulation with van der Waals forces

The van der Waals interaction energy has been computed for a large number of silver and water particle pairs in the range of 1 to 15 nm radius by first calculating the exact, retarded Hamaker factor for semi-infinite bodies and then applying a geometric correction for sphericity. The enhancement of particle collision rates has been deduced from these calculations. These values for the collisional enhancement were used in a direct numerical simulation of the coagulation of silver and water aerosols. Despite the size dependence of the van der Waals enhancement factors, the aerosol size distributions quickly attained nearly self-preserving forms.     

交联SU-8光刻胶与Ni基底结合性的分子动力学模拟

应用分子动力学模拟软件Materials Studio构建SU-8光刻胶与Ni基底的界面结构,研究后烘温度对界面结合性的影响.结合工艺中所采用的后烘温度,模拟计算了338～368K时Ni基底上SU-8胶的交联反应,在经过反复的能量最小化和分子动力学模拟后,对最终得到的平衡结构进行了界面结合能的计算.计算结果表明界面结合能随着后烘温度的升高而增大,在368K时结合能达到最大值,说明此时界面结合最好.对分子体系进行了能量分析,结果表明界面分子间的范德华力作用能是影响界面结合的主要因素.对体系界面原子间进行了径向分布函数分析,发现范德华力作用范围内(0.31～0.60nm)出现两组Ni—O的强峰,也证实了上述结论。     

α-quartz型GeO2高压相变的分子动力学研究

采用在经典离子晶体作用势中附加Morse势,并进行必要的量子化修正,对α-quartz型GeO2结构随压力变化特性,进行分子动力学计算模拟,获得了在压力高于6.0 GPa, α-quartz型GeO2从晶相向非晶相相变的模拟结果,并利用其摩尔体积变化、键角、径向分布、配位数等重要信息对模拟结果作了深入的探讨;相变后的非晶相,由占体积66%的八面体结构和33%的四面体结构组成的非晶体,其中还有极少量的α-quartz型GeO2存在。     

MOIL-opt: Energy-Conserving Molecular Dynamics on a GPU/CPU system

We report an optimized version of the molecular dynamics program MOIL that runs on a shared memory system with OpenMP and exploits the power of a Graphics Processing Unit (GPU). The model is of heterogeneous computing system on a single node with several cores sharing the same memory and a GPU. This is a typical laboratory tool, which provides excellent performance at minimal cost. Besides performance, emphasis is made on accuracy and stability of the algorithm probed by energy conservation for explicit-solvent atomically-detailed-models. Especially for long simulations energy conservation is critical due to the phenomenon known as "energy drift" in which energy errors accumulate linearly as a function of simulation time. To achieve long time dynamics with acceptable accuracy the drift must be particularly small. We identify several means of controlling long-time numerical accuracy while maintaining excellent speedup. To maintain a high level of energy conservation SHAKE and the Ewald reciprocal summation are run in double precision. Double precision summation of real-space non-bonded interactions improves energy conservation. In our best option, the energy drift using 1fs for a time step while constraining the distances of all bonds, is undetectable in 10ns simulation of solvated DHFR (Dihydrofolate reductase). Faster options, shaking only bonds with hydrogen atoms, are also very well behaved and have drifts of less than 1kcal/mol per nanosecond of the same system. CPU/GPU implementations require changes in programming models. We consider the use of a list of neighbors and quadratic versus linear interpolation in lookup tables of different sizes. Quadratic interpolation with a smaller number of grid points is faster than linear lookup tables (with finer representation) without loss of accuracy. Atomic neighbor lists were found most efficient. Typical speedups are about a factor of 10 compared to a single-core single-precision code.     

Van der Waals surface graphs and molecular shape

A van der Waals surface graph is the graph defined on a van der Waals surface by the intersections of the atomic van der Waals spheres. A van der Waals shape graph has a vertex for each atom with a visible face on the van der Waals surface, and edges between vertices representing atoms with adjacent faces on the van der Waals surface. These are discrete invariants of three‐dimensional molecular shape. Some basic properties of van der Waals surface graphs are studied, including their relationship with the Voronoi diagram of the atom centres, and a class of molecular embeddings is identified for which the dual of the van der Waals surface graph coincides with the van der Waals shape graph. This revised version was published online in July 2006 with corrections to the Cover Date.     

Computationally efficient canonical molecular dynamics simulations by using a multiple time‐step integrator algorithm combined with the particle mesh Ewald method and with the fast multipole method

An efficient implementation of the canonical molecular dynamics simulation using the reversible reference system propagator algorithm (r‐RESPA) combined with the particle mesh Ewald method (PMEM) and with the macroscopic expansion of the fast multipole method (MEFMM) was examined. The performance of the calculations was evaluated for systems with 3000, 9999, 30,000, 60,000, and 99,840 particles. For a given accuracy, the optimal conditions for minimizing the CPU time for the implementation of the Ewald method, the PMEM, and the MEFMM were first analyzed. Using the optimal conditions, we evaluated the performance and the reliability of the integrated methods. For all the systems examined, the r‐RESPA with the PMEM was about twice as fast as the r‐RESPA with the MEFMM. The difference arose from the difference in the numerical complexities of the fast Fourier transform in the PMEM and from the transformation of the multipole moments into the coefficients of the local field expansion in the MEFMM. Compared with conventional methods, such as the velocity‐verlet algorithm with the Ewald method, significant speedups were obtained by the integrated methods; the speedup of the calculation was a function of system size, and was a factor of 100 for a system with 3000 particles and increased to a factor of 700 for a system with 99,840 particles. These integrated calculations are, therefore, promising for realizing large‐scale molecular dynamics simulations for complex systems. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 201–217, 2000     

Energy ranking of molecular crystals using density functional theory calculations and an empirical van der waals correction

By combination of high level density functional theory (DFT) calculations with an empirical van der Waals correction, a hybrid method has been designed and parametrized that provides unprecedented accuracy for the structure optimization and the energy ranking of molecular crystals. All DFT calculations are carried out using the VASP program. The van der Waals correction is expressed as the sum over atom-atom pair potentials with each pair potential for two atoms A and B being the product of an asymptotic C(6,A,B)/r(6) term and a damping function d(A,B)(r). Empirical parameters are provided for the elements H, C, N, O, F, Cl, and S. Following Wu and Yang, the C(6) coefficients have been determined by least-squares fitting to molecular C(6) coefficients derived by Meath and co-workers from dipole oscillator strength distributions. The damping functions d(A,B)(r) guarantee the crossover from the asymptotic C(6,A,B)/r(6) behavior at large interatomic distances to a constant interaction energy at short distances. The careful parametrization of the damping functions is of crucial importance to obtain the correct balance between the DFT part of the lattice energy and the contribution from the empirical van der Waals correction. The damping functions have been adjusted to yield the best possible agreement between the unit cells of a set of experimental low temperature crystal structures and their counterparts obtained by lattice energy optimization using the hybrid method. On average, the experimental and the calculated unit cell lengths deviate by 1%. To assess the performance of the hybrid method with respect to the lattice energy ranking of molecular crystals, various crystal packings of ethane, ethylene, acetylene, methanol, acetic acid, and urea have been generated with Accelrys' Polymorph Predictor in a first step and optimized with the hybrid method in a second step. In five out of six cases, the experimentally observed low-temperature crystal structure corresponds to the most stable calculated structure.     

乙酰胆碱酯酶AChE与1,7-二氮杂咔唑抑制剂的作用机理的分子动力学模拟

通过分子对接、分子动力学(MD)模拟以及成键自由能分析方法,从原子水平上模拟研究了3种1,7-二氮杂咔唑衍生物(分别记为M1、M2和M3)与ACh E的结合模式及相互作用机理,分析和讨论了研究体系的静电相互作用和范德华相互作用(vd W)。用MM-PBSA方法计算的3种抑制剂与ACh E之间的结合自由能与抑制剂的实验生物活性数据(IC50值)相对应。分析结果表明,残基S286与抑制剂之间形成的氢键作用有利于抑制剂与ACh E之间的结合。范德华相互作用,尤其是抑制剂与关键残基W279和Y334的作用,对抑制剂与ACh E之间的结合自由能有较大的贡献,在区分抑制剂M1(或M2)和M3的生物活性上发挥着重要的作用。     

乙酰胆碱酯酶AChE与1,7-二氮杂咔唑抑制剂的作用机理的分子动力学模拟

通过分子对接、分子动力学（MD）模拟以及成键自由能分析方法,从原子水平上模拟研究了3种1,7-二氮杂咔唑衍生物（分别记为M1、M2和M3）与AChE的结合模式及相互作用机理,分析和讨论了研究体系的静电相互作用和范德华相互作用（vdW）。用MM-PBSA方法计算的3种抑制剂与AChE之间的结合自由能与抑制剂的实验生物活性数据（IC₅₀值）相对应。分析结果表明,残基S286与抑制剂之间形成的氢键作用有利于抑制剂与AChE之间的结合。范德华相互作用,尤其是抑制剂与关键残基W279和Y334的作用,对抑制剂与AChE之间的结合自由能有较大的贡献,在区分抑制剂M1（或M2）和M3的生物活性上发挥着重要的作用。     

Implementations of Nosé-Hoover and Nosé-Poincaré thermostats in mesoscopic dynamic simulations with the united-residue model of a polypeptide chain

Molecular dynamics (MD) simulations generate a canonical ensemble only when integration of the equations of motion is coupled to a thermostat. Three extended phase space thermostats, one version of Nose-Hoover and two versions of Nose-Poincare, are compared with each other and with the Berendsen thermostat and Langevin stochastic dynamics. Implementation of extended phase space thermostats was first tested on a model Lennard-Jones fluid system; subsequently, they were implemented with our physics-based protein united-residue (UNRES) force field MD. The thermostats were also implemented and tested for the multiple-time-step reversible reference system propagator (RESPA). The velocity and temperature distributions were analyzed to confirm that the proper canonical distribution is generated by each simulation. The value of the artificial mass constant, Q, of the thermostat has a large influence on the distribution of the temperatures sampled during UNRES simulations (the velocity distributions were affected only slightly). The numerical stabilities of all three algorithms were compared with each other and with that of microcanonical MD. Both Nose-Poincare thermostats, which are symplectic, were not very stable for both the Lennard-Jones fluid and UNRES MD simulations started from nonequilibrated structures which implies major changes of the potential energy throughout a trajectory. Even though the Nose-Hoover thermostat does not have a canonical symplectic structure, it is the most stable algorithm for UNRES MD simulations. For UNRES with RESPA, the "extended system inside-reference system propagator algorithm" of the RESPA implementation of the Nose-Hoover thermostat was the only stable algorithm, and enabled us to increase the integration time step.     

Characterization of V5+/TiO2-Rutile Using Molecular Dynamics Simulations

Molecular dynamics (MD) simulations with a quantum correction were used to study the 0.5 mol% V5+/TiO2 rutile under conditions of 300 K and 101 kPa. The interatomic potential function in MD simulations is composed of Coulomb, short-range repulsion, van der Waals, and Morse interactions. The topology of rutile was found to undergo a large deformation when Ti4+ is replaced by V5+, which can be related to the difference of valence and ionic radius between V5+ and Ti4+. The graphed distributions of Ti-O and O-O bond lengths, and O-Ti-O angles are all broaden. The simulations showed that when Ti4+ is replaced by V5+, the V5+ moves out of the O6 polyhedron and toward the interstice between TiO6 octahedra, which results in serious distortion of the octahedra near the interstice, but the phase with 0.5 mol% V5+ doping still remains in rutile. The structural features, such as the bond lengths, migration of the dopant, and crystal phase in the MD simulation are consistent with the experimental observations by FTIR, Raman, XRD and ESR.     

Exact Non-Linear Relationship Between Exponential-6 and Lennard-Jones (12-6) Potential Functions

The relationship between two of the most frequently adopted van der Waals potential functions – Exponential-6 and Lennard-Jones (12-6) – is shown to be faulty upon comparison of the repulsive terms. By using the Maclaurin's series expansions, an exact relationship between both the potential functions is obtained by means of a non-linear correction factor. A purely Exponential-6 form is recovered when the correction factor is taken at zeroth order, while a purely Lennard-Jones (12-6) form is attained as the order of the correction factor tends to infinity. For real van der Waals systems that are bounded by the Exponential-6 and Lennard-Jones (12-6) potential functions, a positive integer order of the correction factor may possibly provide good curve-fitting to experimental data.     

Ionic Liquids: Dissecting the Enthalpies of Vaporization

We calculate the heats of vaporisation for imidazolium‐based ionic liquids [C_nmim][NTf₂] with n=1, 2, 4, 6, 8 by means of molecular dynamics (MD) simulations and discuss their behavior with respect to temperature and the alkyl chain length. We use a force field developed recently. The different cohesive energies contributing to the overall heats of vaporisations are discussed in detail. With increasing alkyl chain length, the Coulomb contribution to the heat of vaporisation remains constant at around 80 kJ mol^?1, whereas the van der Waals interaction increases continuously. The calculated increase of about 4.7 kJ mol^?1 per CH₂‐group of the van der Waals contribution in the ionic liquid exactly coincides with the increase in the heats of vaporisation for n‐alcohols and n‐alkanes, respectively. The results support the importance of van der Waals interactions even in systems completely composed of ions.