The flexible, polarizable, thole-type interaction potential for water (TTM2-F) revisited

We present a revision of the flexible, polarizable, Thole-type interaction potential for water [J. Chem. Phys.2002, 116, 5115], which allows for condensed-phase simulations. The revised version (TTM2.1-F) of the potential correctly describes the individual water molecular dipole moment and alleviates problems arising at short intermolecular separations that can be sampled in the course of molecular dynamics and Monte Carlo simulations of condensed environments. Furthermore, its parallel implementation under periodic boundary conditions enables the efficient calculation of the macroscopic structural and thermodynamic properties of liquid water, as its performance scales superlinearly with up to a number of 64 processors for a simulation box of 512 molecules. We report the radial distribution functions, average energy, internal geometry, and dipole moment in the liquid as well as the density, dielectric constant, and self-diffusion coefficient at T = 300 K from (NVT) and (NPT) classical molecular dynamics simulations by using the revised version of the potential.     

Stress tensor in model polymer systems with periodic boundaries

The calculation of the stress tensor from molecular simulations of atomistic model polymer systems employing periodic boundary conditions is discussed. Starting from the dynamical equations governing the motion of sites, correct double summation forms of the atomic and the molecular virial equations are derived, which are valid for flexible, infinitely stiff and rigid chain models even in the presence of interactions between different images of the same parent macromolecule. A new expression for the true instantaneous stress (flux of momentum through the faces of the simulation box) is derived and shown to exhibit large fluctuations when applied in molecular dynamics simulations. A new equation for the thermodynamic stress, cast exclusively in terms of intermolecular forces on interaction sites, is also derived. Application to Monte Carlo simulations shows that the molecular virial expression exhibits the smallest fluctuations among all stress expressions discussed, and thus allows computation of the thermodynamic stress with least uncertainty. A scheme is developed for the calculation of surface tension from intermolecular forces only.     

A theoretical study of aqueous solvation of K comparing ab initio, polarizable, and fixed-charge models

The hydration of K(+) is studied using a hierarchy of theoretical approaches, including ab initio Born-Oppenheimer molecular dynamics and Car-Parrinello molecular dynamics, a polarizable force field model based on classical Drude oscillators, and a nonpolarizable fixed-charge potential based on the TIP3P water model. While models based more directly on quantum mechanics offer the possibility to account for complex electronic effects, polarizable and fixed-charges force fields allow for simulations of large systems and the calculation of thermodynamic observables with relatively modest computational costs. A particular emphasis is placed on investigating the sensitivity of the polarizable model to reproduce key aspects of aqueous K(+), such as the coordination structure, the bulk hydration free energy, and the self diffusion of K(+). It is generally found that, while the simple functional form of the polarizable Drude model imposes some restrictions on the range of properties that can simultaneously be fitted, the resulting hydration structure for aqueous K(+) agrees well with experiment and with more sophisticated computational models. A counterintuitive result, seen in Car-Parrinello molecular dynamics and in simulations with the Drude polarizable force field, is that the average induced molecular dipole of the water molecules within the first hydration shell around K(+) is slightly smaller than the corresponding value in the bulk. In final analysis, the perspective of K(+) hydration emerging from the various computational models is broadly consistent with experimental data, though at a finer level there remain a number of issues that should be resolved to further our ability in modeling ion hydration accurately.     

Ab initio molecular dynamics using hybrid density functionals

Ab initio molecular dynamics simulations with hybrid density functionals have so far found little application due to their computational cost. In this work, an implementation of the Hartree-Fock exchange is presented that is specifically targeted at ab initio molecular dynamics simulations of medium sized systems. We demonstrate that our implementation, which is available as part of the CP2K/Quickstep program, is robust and efficient. Several prescreening techniques lead to a linear scaling cost for integral evaluation and storage. Integral compression techniques allow for in-core calculations on systems containing several thousand basis functions. The massively parallel implementation respects integral symmetry and scales up to hundreds of CPUs using a dynamic load balancing scheme. A time-reversible multiple time step scheme, exploiting the difference in computational efficiency between hybrid and local functionals, brings further time savings. With extensive simulations of liquid water, we demonstrate the ability to perform, for several tens of picoseconds, ab initio molecular dynamics based on hybrid functionals of systems in the condensed phase containing a few thousand Gaussian basis functions.     

GROMACS: fast, flexible, and free

This article describes the software suite GROMACS (Groningen MAchine for Chemical Simulation) that was developed at the University of Groningen, The Netherlands, in the early 1990s. The software, written in ANSI C, originates from a parallel hardware project, and is well suited for parallelization on processor clusters. By careful optimization of neighbor searching and of inner loop performance, GROMACS is a very fast program for molecular dynamics simulation. It does not have a force field of its own, but is compatible with GROMOS, OPLS, AMBER, and ENCAD force fields. In addition, it can handle polarizable shell models and flexible constraints. The program is versatile, as force routines can be added by the user, tabulated functions can be specified, and analyses can be easily customized. Nonequilibrium dynamics and free energy determinations are incorporated. Interfaces with popular quantum-chemical packages (MOPAC, GAMES-UK, GAUSSIAN) are provided to perform mixed MM/QM simulations. The package includes about 100 utility and analysis programs. GROMACS is in the public domain and distributed (with source code and documentation) under the GNU General Public License. It is maintained by a group of developers from the Universities of Groningen, Uppsala, and Stockholm, and the Max Planck Institute for Polymer Research in Mainz. Its Web site is http://www.gromacs.org.     

The ice-vapor interface and the melting point of ice I(h) for the polarizable POL3 water model

We use molecular dynamics simulations to determine the melting point of ice I(h) for the polarizable POL3 water force field (Dang, L. X. J. Chem. Phys.1992, 97, 2659). Simulations are performed on a slab of ice I(h) with two free surfaces at several different temperatures. The analysis of the time evolution of the total energy in the course of the simulations at the set of temperatures yields the melting point of the POL3 model to be T(m) = 180 ± 10 K. Moreover, the results of the simulations show that the degree of hydrogen-bond disorder occurring in the bulk of POL3 ice is larger (at the corresponding degree of undercooling) than in ice modeled by nonpolarizable water models. These results demonstrate that the POL3 water force field is rather a poor model for studying ice and ice-liquid or ice-vapor interfaces. While a number of polarizable water models have been developed over the past years, little is known about their performance in simulations of supercooled water and ice. This study thus highlights the need for testing of the existing polarizable water models over a broad range of temperatures, pressures, and phases, and developing a new polarizable water force field, reliable over larger areas of the phase diagram.     

Elastic bag model for molecular dynamics simulations of solvated systems: application to liquid water and solvated peptides

The fluctuating elastic boundary (FEB) model for molecular dynamics has recently been developed and validated through simulations of liquid argon. In the FEB model, a flexible boundary which consists of particles connected by springs is used to confine the solvated system, thereby eliminating the need for periodic boundary conditions. In this study, we extend this model to the simulation of bulk water and solvated alanine dipeptide. Both the confining potential and boundary particle interaction functions are modified to preserve the structural integrity of the boundary and prevent the leakage of the solute-solvent system through the boundary. A broad spectrum of structural and dynamic properties of liquid water are computed and compared with those obtained from conventional periodic boundary condition simulations. The applicability of the model to biomolecular simulations is investigated through the analysis of conformational population distribution of solvated alanine dipeptide. In most cases we find remarkable agreement between the two simulation approaches.     

Hybrid particle‐field molecular dynamics simulations: Parallelization and benchmarks

The parallel implementation of a recently developed hybrid scheme for molecular dynamics (MD) simulations (Milano and Kawakatsu, J Chem Phys 2009, 130, 214106) where self‐consistent field theory (SCF) and particle models are combined is described. Because of the peculiar formulation of the hybrid method, considering single particles interacting with density fields, the most computationally expensive part of the hybrid particle‐field MD simulation can be efficiently parallelized using a straightforward particle decomposition algorithm. Benchmarks of simulations, including comparisons of serial MD and MD‐SCF program profiles, serial MD‐SCF and parallel MD‐SCF program profiles, and parallel benchmarks compared with efficient MD program GROMACS 4.5.4 are tested and reported. The results of benchmarks indicate that the proposed parallelization scheme is very efficient and opens the way to molecular simulations of large scale systems with reasonable computational costs. © 2012 Wiley Periodicals, Inc.     

Parallel replica dynamics with a heterogeneous distribution of barriers: application to n-hexadecane pyrolysis

Parallel replica dynamics simulation methods appropriate for the simulation of chemical reactions in molecular systems with many conformational degrees of freedom have been developed and applied to study the microsecond-scale pyrolysis of n-hexadecane in the temperature range of 2100-2500 K. The algorithm uses a transition detection scheme that is based on molecular topology, rather than energetic basins. This algorithm allows efficient parallelization of small systems even when using more processors than particles (in contrast to more traditional parallelization algorithms), and even when there are frequent conformational transitions (in contrast to previous implementations of the parallel replica algorithm). The parallel efficiency for pyrolysis initiation reactions was over 90% on 61 processors for this 50-atom system. The parallel replica dynamics technique results in reaction probabilities that are statistically indistinguishable from those obtained from direct molecular dynamics, under conditions where both are feasible, but allows simulations at temperatures as much as 1000 K lower than direct molecular dynamics simulations. The rate of initiation displayed Arrhenius behavior over the entire temperature range, with an activation energy and frequency factor of E(a) = 79.7 kcal/mol and log A/s(-1) = 14.8, respectively, in reasonable agreement with experiment and empirical kinetic models. Several interesting unimolecular reaction mechanisms were observed in simulations of the chain propagation reactions above 2000 K, which are not included in most coarse-grained kinetic models. More studies are needed in order to determine whether these mechanisms are experimentally relevant, or specific to the potential energy surface used.     

Comparison of selected polarizable and nonpolarizable water models in molecular dynamics simulations of ice I(h)

We present a molecular dynamics simulation study in which we determined the melting point of ice I(h) for the polarizable SWM4-NDP water model (Lamoureux et al., Chem. Phys. Lett., 2006, 418, 245-249) and compared the performance of several popular water force fields, both polarizable and nonpolarizable, in terms of melting temperature, stability and orientational structuring of ice. The simulations yield the melting temperature of SWM4-NDP ice as low as T(m) = 185 ± 10 K, despite the quadrupole moment of a SWM4-NDP water molecule being close to the experimental gas phase value. The results thus show that the dependence of T(m) on the molecular quadrupole, observed for the three- and four-site water models, is generally lost if polarization is explicitly included. The study also shows that adding polarizability to a planar three-charge water model increases orientational disorder in hexagonal ice. In addition, analysis of the tetrahedral order in bulk ice reveals a correlation between the pre-existing degree of orientational disorder in ice simulated using different polarizable and nonpolarizable models and the melting temperature of the models. Our findings thus suggest some new considerations regarding the role of polarization forces in a crystalline solid that may guide future development of reliable polarizable water models for ice.     

Midpoint cell method for hybrid (MPI+OpenMP) parallelization of molecular dynamics simulations

We have developed a new hybrid (MPI+OpenMP) parallelization scheme for molecular dynamics (MD) simulations by combining a cell‐wise version of the midpoint method with pair‐wise Verlet lists. In this scheme, which we call the midpoint cell method, simulation space is divided into subdomains, each of which is assigned to a MPI processor. Each subdomain is further divided into small cells. The interaction between two particles existing in different cells is computed in the subdomain containing the midpoint cell of the two cells where the particles reside. In each MPI processor, cell pairs are distributed over OpenMP threads for shared memory parallelization. The midpoint cell method keeps the advantages of the original midpoint method, while filtering out unnecessary calculations of midpoint checking for all the particle pairs by single midpoint cell determination prior to MD simulations. Distributing cell pairs over OpenMP threads allows for more efficient shared memory parallelization compared with distributing atom indices over threads. Furthermore, cell grouping of particle data makes better memory access, reducing the number of cache misses. The parallel performance of the midpoint cell method on the K computer showed scalability up to 512 and 32,768 cores for systems of 20,000 and 1 million atoms, respectively. One MD time step for long‐range interactions could be calculated within 4.5 ms even for a 1 million atoms system with particle‐mesh Ewald electrostatics. © 2014 Wiley Periodicals, Inc.     

Structural and dynamical properties of methane clathrate hydrates

Equilibrium molecular dynamics (MD) simulations have been performed in both the NVT and NPT ensembles to study the structural and dynamical properties of fully occupied methane clathrate hydrates at 50, 125, and 200 K. Five atomistic potential models were used for water, ranging from fully flexible to rigid polarizable and nonpolarizable. A flexible and a rigid model were utilized for methane. The phonon densities of states were evaluated and the localized rattling modes for the methane molecules were found to couple to the acoustic phonons of the host lattice. The calculated methane density of states was found to be in reasonable agreement with available experimental data.     

Analytic energy gradient in combined second-order Møller-Plesset perturbation theory and polarizable force field calculation

Second-order M?ller-Plesset perturbation theory (MP2) is used to describe electronic correlation on the basis of Hartree-Fock (HF) variational calculations that incorporate induced dipole polarizable force fields (i.e., QM/MMpol style HF and MP2). The Z-vector equations for regular closed shell and open shell MP2 methods (RMP2, ZAPT2, and UMP2) are extended to include induced dipole contributions to determine the MP2 response density so that nuclear gradient and other properties can be efficiently evaluated. A better estimation of the induced dipole polarization energy can be obtained using the MP2 relaxed density. QM/MMpol style MP2 molecular dynamics simulations are performed for the ground state and first triplet state of acetone solvated by 1024 polarizable water molecules. A switching function is used to ensure energy conservation in QM/MM simulation under periodic boundary condition.     

Ab initio based polarizable force field parametrization

Experimental and simulation studies of anion-water systems have pointed out the importance of molecular polarization for many phenomena ranging from hydrogen-bond dynamics to water interfaces structure. The study of such systems at molecular level is usually made with classical molecular dynamics simulations. Structural and dynamical features are deeply influenced by molecular and ionic polarizability, which parametrization in classical force field has been an object of long-standing efforts. Although when classical models are compared to ab initio calculations at condensed phase, it is found that the water dipole moments are underestimated by approximately 30%, while the anion shows an overpolarization at short distances. A model for chloride-water polarizable interaction is parametrized here, making use of Car-Parrinello simulations at condensed phase. The results hint to an innovative approach in polarizable force fields development, based on ab initio simulations, which do not suffer for the mentioned drawbacks. The method is general and can be applied to the modeling of different systems ranging from biomolecular to solid state simulations.     

On the pressure calculation for polarizable models in computer simulation

We present a short overview of pressure calculation in molecular dynamics or Monte Carlo simulations. The emphasis is given to polarizable models in order to resolve the controversy caused by the paper of M. J. Louwerse and E. J. Baerends [Chem. Phys. Lett. 421, 138 (2006)] about pressure calculation in systems with periodic boundaries. We systematically derive expressions for the pressure and show that despite the lack of explicit pairwise additivity, the pressure formula for polarizable models is identical with that of nonpolarizable ones. However, a strict condition for using this formula is that the induced dipole should be in perfect mechanical equilibrium prior to pressure calculation. The perfect convergence of induced dipoles ensures conservation of energy as well. We demonstrate using more cumbersome but exact methods that the derived expressions for the polarizable model of water provide correct numerical results. We also show that the inaccuracy caused by imperfect convergence of the induced dipoles correlates with the inaccuracy of the calculated pressure.     

Extension of adaptive tree code and fast multipole methods to high angular momentum particle charge densities

The development and implementation of a tree code (TC) and fast multipole method (FMM) for the efficient, linear-scaling calculation of long-range electrostatic interactions of particle distributions with variable shape and multipole character are described. The target application of these methods are stochastic boundary molecular simulations with polarizable force fields and/or combined quantum mechanical/molecular mechanical potentials. Linear-scaling is accomplished through the adaptive decomposition of the system into a hierarchy of interacting particle sets. Two methods for effecting this decomposition are evaluated: fluc-splitting and box-splitting, for which the latter is demonstrated to be generally more accurate. In addition, a generalized termination criterion is developed that delivers optimal performance at fixed error tolerance that, in the case of quadrupole-represented Drude water, effects a speed-up by a factor of 2-3 relative to a multipole-independent termination criteria. The FMM is shown to be approximately 2-3 times faster than the TC, independent of the system size and multipole order of the particles. The TC and FMM are tested for a variety of static and polarizable water systems, and for the the 70S ribosome functional complex containing an assembly of transfer and messenger RNAs.     

Hydrogen-bond dynamics in the air-water interface

Hydrogen-bond (H-bond) dynamics in the air-water interface is studied by molecular dynamics simulations. The analysis reveals that the dynamics of breaking and forming hydrogen bonds in the air-water interface is faster than that in bulk water for the polarizable water models. This is in contrast to the results found on a protein surface. We show that the difference stems from more rapid translational diffusion in the interface. When the effect of pair diffusion is eliminated, the hydrogen-bond dynamics in the interface is observed to be slower than that in the bulk. This occurs because the number of water molecules adjacent to a hydrogen-bonded pair and available to accept or donate a hydrogen bond is smaller in the interface than in the bulk. The comparison between polarizable water models and fixed-charge models highlights the potential importance of the polarization effect in the water-vapor interface.     

High-performance scalable molecular dynamics simulations of a polarizable force field based on classical Drude oscillators in NAMD

Incorporating the influence of induced polarization in large-scale atomistic molecular dynamics (MD) simulations is a critical challenge in the progress toward computations of increased accuracy. One computationally efficient treatment is based on the classical Drude oscillator in which an auxiliary charged particle is attached by a spring to each nucleus. Here, we report the first implementation of this model in the program NAMD. An extended Lagrangian dynamics with a dual-Langevin thermostat scheme applied to the Drude-nucleus pairs is employed to efficiently generate classical dynamic propagation near the self-consistent field limit. Large-scale MD simulations based on the Drude polarizable force field scale very well on massively distributed supercomputing platforms, the computational demand being only about 50-100% higher than for nonpolarizable models. As an illustration, a large-scale 150 mM NaCl aqueous salt solution is simulated, and the calculated ionic conductivity is shown to be in excellent agreement with experiment.     

The QM‐MM interface for CHARMM‐deMon

We present a new QM/MM interface for fast and efficient simulations of organic and biological molecules. The CHARMM/deMon interface has been developed and tested to perform minimization and atomistic simulations for multi‐particle systems. The current features of this QM/MM interface include readability for molecular dynamics, tested compatibility with Free Energy Perturbation simulations (FEP) using the dual topology/single coordinate method. The current coupling scheme uses link atoms, but further extensions of the code to incorporate other available schemes are planned. We report the performance of different levels of theory for the treatment of the QM region, while the MM region was represented by a classical force‐field (CHARMM27) or a polarizable force‐field based on a simple Drude model. The current QM/MM implementation can be coupled to the dual‐thermostat method and the VV2 integrator to run molecular dynamics simulations. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010