Hybrid Hamiltonian replica exchange molecular dynamics simulation method employing the Poisson-Boltzmann model

A hybrid Hamiltonian replica exchange molecular dynamics simulation scheme based on explicit water model hybrided with Poisson-Boltzmann model is brought out. In this method the motions of atoms are governed by potential energy obtained from explicit water model. However, the exchanges between different replicas under different temperatures are controlled by the solvation energies of the solute calculated using the Poisson-Boltzmann model. In order to get the correct canonical ensembles, the van der Waals radii, which are used to define the dielectric boundary, have to be optimized. The conformational spaces of three distinct pentapeptides, Met-enkephalin, alanine 5, and glycine 5, are explored. We find that with the optimized radii the structural ensembles are nearly identical to those obtained by standard replica exchange simulations while the number of replica needed is reduced greatly.     

Fluctuating hydrodynamics for multiscale simulation of inhomogeneous fluids: mapping all-atom molecular dynamics to capillary waves

We introduce a multiscale framework to simulate inhomogeneous fluids by coarse-graining an all-atom molecular dynamics (MD) trajectory onto sequential snapshots of hydrodynamic fields. We show that the field representation of an atomistic trajectory is quantitatively described by a dynamic field-theoretic model that couples hydrodynamic fluctuations with a Ginzburg-Landau free energy. For liquid-vapor interfaces of argon and water, the parameters of the field model can be adjusted to reproduce the bulk compressibility and surface tension calculated from the positions and forces of atoms in an MD simulation. These optimized parameters also enable the field model to reproduce the static and dynamic capillary wave spectra calculated from atomistic coordinates at the liquid-vapor interface. In addition, we show that a density-dependent gradient coefficient in the Ginzburg-Landau free energy enables bulk and interfacial fluctuations to be controlled separately. For water, this additional degree of freedom is necessary to capture both the bulk compressibility and surface tension emergent from the atomistic trajectory. The proposed multiscale framework illustrates that bottom-up coarse-graining and top-down phenomenology can be integrated with quantitative consistency to simulate the interfacial fluctuations in nanoscale transport processes.     

Simple and effective application of the Wang-Landau method for multicanonical molecular dynamics simulation

We propose a novel application of the Wang-Landau method (WLM) for multicanonical molecular dynamics (McMD) simulations. Originally, WLM was developed for Monte Carlo (MC) simulations. Fundamentally, WLM remarkably reduces simulation efforts because it estimates the optimal multicanonical energy function automatically. When WLM is applied to McMD, not only the multicanonical energy but also energy gradient must be estimated adequately. However, because of the rugged multicanonical energy function at the early simulation stage, applications of WLM for MD simulations are difficult and require a smoothing procedure: simulation efforts such as cubic-spline extrapolation and gathering multiple preruns are utilized for smoothing. We propose a simple and effective smoothing method that requires only one additional equation and two time-dependent parameters. As a result, our method produced the correct multicanonical energy function and succeeded in the flat sampling of a small biomolecule with reduced simulation effort.     

A new parallel method for molecular dynamics simulation of macromolecular systems

Short-range molecular dynamics simulations of molecular systems are commonly parallelized by replicated-data methods, in which each processor stores a copy of all atom positions. This enables computation of bonded 2-, 3-, and 4-body forces within the molecular topology to be partitioned among processors straightforwardly. A drawback to such methods is that the interprocessor communication scales as N (the number of atoms) independent of P (the number of processors). Thus, their parallel efficiency falls off rapidly when large numbers of processors are used. In this article a new parallel method for simulating macromolecular or small-molecule systems is presented, called force-decomposition. Its memory and communication costs scale as N/√P, allowing larger problems to be run faster on greater numbers of processors. Like replicated-data techniques, and in contrast to spatial-decomposition approaches, the new method can be simply load balanced and performs well even for irregular simulation geometries. The implementation of the algorithm in a prototypical macromolecular simulation code ParBond is also discussed. On a 1024-processor Intel Paragon, ParBond runs a standard benchmark simulation of solvated myoglobin with a parallel efficiency of 61% and at 40 times the speed of a vectorized version of CHARMM running on a single Cray Y-MP processor. © 1996 by John Wiley & Sons, Inc.     

An efficient molecular dynamics simulation method for calculating the diffusion-influenced reaction rates

We present a molecular dynamics (MD) simulation method for calculating the diffusion-influenced reaction rates in the limit of low reactant concentrations. To calculate the reaction rate coefficient, we use MD trajectories of a nonreactive equilibrium system that are initiated with a pair of reactant molecules in reactive configuration. Hence reaction systems involving complicated reactant molecules with geometrically restricted reactivities can be treated with comparable efficiency as the simple hard-sphere reaction system. Compared to the similar MD method proposed by Van Beijeren, Dong, and Bocquet [J. Chem. Phys. 114, 6265 (2001)], the present method has a couple of advantages. First, reactions involving more general sink functions can be treated. Second, more accurate results can be obtained when the reaction probability upon collision is less than unity. As an application, we investigate the effects of nondiffusive dynamics and hydrodynamic interaction of reactants on the reaction rate.     

Accelerated molecular dynamics: a promising and efficient simulation method for biomolecules

Many interesting dynamic properties of biological molecules cannot be simulated directly using molecular dynamics because of nanosecond time scale limitations. These systems are trapped in potential energy minima with high free energy barriers for large numbers of computational steps. The dynamic evolution of many molecular systems occurs through a series of rare events as the system moves from one potential energy basin to another. Therefore, we have proposed a robust bias potential function that can be used in an efficient accelerated molecular dynamics approach to simulate the transition of high energy barriers without any advance knowledge of the location of either the potential energy wells or saddle points. In this method, the potential energy landscape is altered by adding a bias potential to the true potential such that the escape rates from potential wells are enhanced, which accelerates and extends the time scale in molecular dynamics simulations. Our definition of the bias potential echoes the underlying shape of the potential energy landscape on the modified surface, thus allowing for the potential energy minima to be well defined, and hence properly sampled during the simulation. We have shown that our approach, which can be extended to biomolecules, samples the conformational space more efficiently than normal molecular dynamics simulations, and converges to the correct canonical distribution.     

A hybrid fluctuating hydrodynamics and kinetic Monte Carlo method for modeling chemically-powered nanoscale motion

We present a stochastic multiscale method for modeling heterogeneous catalysis at the nanoscale. The system is decomposed into the fluid domain and the catalyst-fluid interface. We implemented the fluctuating hydrodynamics framework to model the diffusion of the chemical species in the fluid domain, and the chemical master equation to describe the catalytic activity at the interface. The coupling between the domains occurs simultaneously. Using a simple one-dimensional (1D) linear model, we showed that the predictions of our scheme are in excellent agreement with deterministic simulations. The method was specifically developed to model the spatially asymmetric catalysis on the surface of self-propelled nanoswimmers. Numerical simulations showed that our approach can estimate the uncertainty in the swimming velocity resulting from inherent stochastic nature of the chemical reactions at the catalytic interface. Although the method has been applied to simple 1D and 2D models, it can be generalized to handle different geometries and more sophisticated chemical reactions. Therefore, it can serve as a practical mathematical tool to study how the efficiency of chemically powered nanomachines is affected by the interplay between structural complexity, nonlinear reactivity, and nonequilibrium fluctuations.     

Local elevation: A method for improving the searching properties of molecular dynamics simulation

Summary The concept of memory has been introduced into a molecular dynamics algorithm. This was done so as to persuade a molecular system to visit new areas of conformational space rather than be confined to a small number of low-energy regions. The method is demonstrated on a simple model system and the 11-residue cyclic peptide cyclosporin A. For comparison, calculations were also performed using simulated temperature annealing and a potential energy annealing scheme. Although the method can only be applied to systems with a small number of degrees of freedom, it offers the chance to generate a multitude of different low-energy structures, where other methods only give a single one or few. This is clearly important in problems such as drug design, where one is interested in the conformational spread of a system.     

Three-dimensional Brownian diffusion of rod-like macromolecules in the presence of randomly distributed spherical obstacles: molecular dynamics simulation

Brownian diffusion of rod-like polymers in the presence of randomly distributed spherical obstacles is studied using molecular dynamics simulations. It is observed that dependence of the reduced diffusion coefficient of these macromolecules on the available volume fraction can be described reasonably by a power law function. Despite the case of obstructed diffusion of flexible polymers in which reduced diffusion coefficient has a weak dependence on the polymer length, this dependence is noticeably strong in the case of rod-like polymers. Diffusion of these macromolecules in the presence of obstacles is observed that is anomalous at short time scales and normal at long times. Duration time of the anomalous diffusion regime is found that increases very rapidly with increasing both the polymer length and the obstructed volume fraction. Dynamics of diffusion of these polymers is observed that crosses over from Rouse to reptation type with increasing the density of obstacles.     

Hybrid quantum mechanics/molecular mechanics-based molecular dynamics simulation of acid-catalyzed dehydration of polyols in liquid water

We use the conversion of protonated glycerol to acrolein for a case study of the mechanism of acid-catalyzed dehydration of polyols in aqueous environments. We employ hybrid Quamtum Mechanics/Molecular Mechanics Molecular Dynamics (QM/MM MD) simulations with biased sampling and perform free energy calculations for the elementary steps of the reaction. We investigate the effects of solvent dynamics and in particular the role of quantum mechanical water in the dehydration mechanism. We present results supporting a mechanism that proceeds via water-mediated proton transfers and thus through an enol intermediate. We find that the first dehydration may take place by two, low-energy pathways requiring, respectively, 20.9 and 18.8 kcal/mol of activation free energy. The second dehydration requires 19.9 kcal/mol of activation free energy while for the overall reaction we compute a free energy change of -8 kcal/mol.     

Normal mode analysis using the driven molecular dynamics method. II. An application to biological macromolecules

The driven molecular-dynamics (DMD) method, recently proposed by Bowman, Zhang, and Brown [J. Chem. Phys. 119, 646 (2003)], has been implemented into the TINKER molecular modeling program package. The DMD method yields frequencies and normal modes without evaluation of the Hessian matrix. It employs an external harmonic driving term that can be used to scan the spectrum and determine resonant absorptions. The molecular motions, induced by driving at resonant frequencies, correspond to the normal-mode vibrations. In the current work we apply the method to a 20-residue protein, Trp-cage, that has been reported by Neidigh, Fesinmeyer, and Andersen [Nature Struct. Biol. 9, 425 (2002)]. The structural and dynamical properties of this molecule, such as B-factors, root-mean square fluctuations, anisotropies, vibrational entropy, and cross-correlations coefficients, are calculated using the DMD method. The results are in very good agreement with ones calculated using standard normal-mode analysis method. Thus, the DMD method provides a viable alternative to the standard Hessian-based method and has considerable potential for the study of large systems, where the Hessian-based method is not feasible.     

Fluctuating hydrodynamics for multiscale modeling and simulation: Energy and heat transfer in molecular fluids

This work illustrates that fluctuating hydrodynamics (FHD) simulations can be used to capture the thermodynamic and hydrodynamic responses of molecular fluids at the nanoscale, including those associated with energy and heat transfer. Using all-atom molecular dynamics (MD) trajectories as the reference data, the atomistic coordinates of each snapshot are mapped onto mass, momentum, and energy density fields on Eulerian grids to generate a corresponding field trajectory. The molecular length-scale associated with finite molecule size is explicitly imposed during this coarse-graining by requiring that the variances of density fields scale inversely with the grid volume. From the fluctuations of field variables, the response functions and transport coefficients encoded in the all-atom MD trajectory are computed. By using the extracted fluid properties in FHD simulations, we show that the fluctuations and relaxation of hydrodynamic fields quantitatively match with those observed in the reference all-atom MD trajectory, hence establishing compatibility between the atomistic and field representations. We also show that inclusion of energy transfer in the FHD equations can more accurately capture the thermodynamic and hydrodynamic responses of molecular fluids. The results indicate that the proposed MD-to-FHD mapping with explicit consideration of finite molecule size provides a robust framework for coarse-graining the solution phase of complex molecular systems.     

Hybrid molecular dynamics simulations of living filaments

We propose a hybrid molecular dynamics/multi-particle collision dynamics model to simulate a set of self-assembled semiflexible filaments and free monomers. Further, we introduce a Monte Carlo scheme to deal with single monomer addition (polymerization) or removal (depolymerization), satisfying the detailed balance condition within a proper statistical mechanical framework. This model of filaments, based on the wormlike chain, aims to represent equilibrium polymers with distinct reaction rates at both ends, such as self-assembled adenosine diphosphate-actin filaments in the absence of adenosine triphosphate (ATP) hydrolysis and other proteins. We report the distribution of filament lengths and the corresponding dynamical fluctuations on an equilibrium trajectory. Potential generalizations of this method to include irreversible steps like ATP-actin hydrolysis are discussed.     

Thermal analysis via molecular dynamics simulation

Thermal analysis by classical molecular dynamics simulations is discussed on hand of heat capacity of crystals of 9600 atoms. The differences between quantum mechanical and classical mechanical calculations are shown. Anharmonicity is proven to be an important factor. Finally, it is found that defects contribute to an increase in heat capacity before melting. The energy of conformational gauche defects within the crystal is only about 10% due to internal rotation. The other energy must be generated by cooperative strain. The conclusion is that the next generation of faster computers may permit wider use of molecular dynamics simulations in support of the interpretation of thermal analysis.Dedicated to Professor Bernhard Wunderlich on the occasion of his 65th birthday     

From molecular dynamics to hydrodynamics: a novel Galilean invariant thermostat

This paper proposes a novel thermostat applicable to any particle-based dynamic simulation. Each pair of particles is thermostated either (with probability P) with a pairwise Lowe-Andersen thermostat [C. P. Lowe, Europhys. Lett. 47, 145 (1999)] or (with probability 1-P) with a thermostat that is introduced here, which is based on a pairwise interaction similar to the Nosé-Hoover thermostat. When the pairwise Nosé-Hoover thermostat dominates (low P), the liquid has a high diffusion coefficient and low viscosity, but when the Lowe-Andersen thermostat dominates, the diffusion coefficient is low and viscosity is high. This novel Nosé-Hoover-Lowe-Andersen thermostat is Galilean invariant and preserves both total linear and angular momentum of the system, due to the fact that the thermostatic forces between each pair of the particles are pairwise additive and central. We show by simulation that this thermostat also preserves hydrodynamics. For the (noninteracting) ideal gas at P = 0, the diffusion coefficient diverges and viscosity is zero, while for P > 0 it has a finite value. By adjusting probability P, the Schmidt number can be varied by orders of magnitude. The temperature deviation from the required value is at least an order of magnitude smaller than in dissipative particle dynamics (DPD), while the equilibrium properties of the system are very well reproduced. The thermostat is easy to implement and offers a computational efficiency better than (DPD), with better temperature control and greater flexibility in terms of adjusting the diffusion coefficient and viscosity of the simulated system. Applications of this thermostat include all standard molecular dynamic simulations of dense liquids and solids with any type of force field, as well as hydrodynamic simulation of multiphase systems with largely different bulk viscosities, including surface viscosity, and of dilute gases and plasmas.     

Effective force field for liquid hydrogen fluoride from ab initio molecular dynamics simulation using the force-matching method

A recently developed force-matching method for obtaining effective force fields for condensed matter systems from ab initio molecular dynamics (MD) simulations has been applied to fit a simple nonpolarizable two-site pairwise force field for liquid hydrogen fluoride. The ab initio MD in this case was a Car-Parrinello (CP) MD simulation of 64 HF molecules at nearly ambient conditions within the Becke-Lee-Yang-Parr approximation to the electronic density functional theory. The force-matching procedure included a fit of short-ranged nonbonded forces, bonded forces, and atomic partial charges. The performance of the force-match potential was examined for the gas-phase dimer and for the liquid phase at various temperatures. The model was able to reproduce correctly the bent structure and energetics of the gas-phase dimer, while the results for the structural properties, self-diffusion, vibrational spectra, density, and thermodynamic properties of liquid HF were compared to both experiment and the CP MD simulation. The force-matching model performs well in reproducing nearly all of the liquid properties as well as the aggregation behavior at different temperatures. The model is computationally cheap and compares favorably to many more computationally expensive potential energy functions for liquid HF.     

Calculating the electric potential of macromolecules: A simple method for molecular surface triangulation

An improved version of the “marching-cube” method¹ is proposed for molecular surface triangulation. This new algorithm involves fewer and simpler basic building blocks and avoids the artificial gaps of the original one. Moreover, to make it applicable to the boundary element method, the procedures for the protein cavity identification and triangle reduction are also presented. The triangulation procedure was tested by incorporating it into the boundary element method (BEM) to estimate the pK_a values of subtilisin BPN′ and bovine trypsin inhibitor (BPTI). © 1995 by John Wiley & Sons, Inc.