Rigid-body dynamics in the isothermal-isobaric ensemble: a test on the accuracy and computational efficiency

We have developed a time-reversible rigid-body (rRB) molecular dynamics algorithm in the isothermal-isobaric (NPT) ensemble. The algorithm is an extension of rigid-body dynamics [Matubayasi and Nakahara, J Chem Phys 1999, 110, 3291] to the NPT ensemble on the basis of non-Hamiltonian statistical mechanics [Martyna, G. J. et al., J Chem Phys 1994, 101, 4177]. A series of MD simulations of water as well as fully hydrated lipid bilayer systems have been undertaken to investigate the accuracy and efficiency of the algorithm. The rRB algorithm was shown to be superior to the state-of-the-art constraint-dynamics algorithm SHAKE/RATTLE/ROLL, with respect to computational efficiency. However, it was revealed that both algorithms produced accurate trajectories of molecules in the NPT as well as NVT ensembles, as long as a reasonably short time step was used. A couple of multiple time-step (MTS) integration schemes were also examined. The advantage of the rRB algorithm for computational efficiency increased when the MD simulation was carried out using MTS on parallel processing computer systems; total computer time for MTS-MD of a lipid bilayer using 64 processors was reduced by about 40% using rRB instead of SHAKE/RATTLE/ROLL.     

Multibaric-multithermal ensemble molecular dynamics simulations

We present new generalized-ensemble molecular dynamics simulation algorithms, which we refer to as the multibaric-multithermal molecular dynamics. We describe three algorithms based on (1) the Nosé thermostat and the Andersen barostat, (2) the Nosé-Poincaré thermostat and the Andersen barostat, and (3) the Gaussian thermostat and the Andersen barostat. The multibaric-multithermal simulations perform random walks widely both in the potential-energy space and in the volume space. Therefore, one can calculate isobaric-isothermal ensemble averages in wide ranges of temperature and pressure from only one simulation run. We test the effectiveness of the multibaric-multithermal algorithm by applying it to a Lennard-Jones 12-6 potential system.     

Balancing an accurate representation of the molecular surface in generalized born formalisms with integrator stability in molecular dynamics simulations

Different integrator time steps in NVT and NVE simulations of protein and nucleic acid systems are tested with the GBMV (Generalized Born using Molecular Volume) and GBSW (Generalized Born with simple SWitching) methods. The simulation stability and energy conservation is investigated in relation to the agreement with the Poisson theory. It is found that very close agreement between generalized Born methods and the Poisson theory based on the commonly used sharp molecular surface definition results in energy drift and simulation artifacts in molecular dynamics simulation protocols with standard 2-fs time steps. New parameters are proposed for the GBMV method, which maintains very good agreement with the Poisson theory while providing energy conservation and stable simulations at time steps of 1 to 1.5 fs.     

A box-counting-based algorithm for computing Shannon entropy in molecular dynamics simulations

A box-counting-based algorithm (SEBC) has been developed for the numerical computation of the Shannon entropy from samples of continuous functions. Its performance was tested by applying it to several samples of known continuous distribution functions. The results obtained with SEBC reproduced those obtained by analytical or numerical integration. SEBC was also employed for computing the Shannon entropies of the steric energy, Sh(E(S)), of several amino acids from their in vacuo NVE molecular dynamics simulations using the AMBER-4 force field. The results obtained correlate linearly with the experimental standard thermodynamic entropies of these compounds. This work points to the possibility of introducing straightforward and reliable calculations of thermodynamic entropies from empirical linear relationships with Sh(E(S)) obtained from MD simulations.     

The distributed diagonal force decomposition method for parallelizing molecular dynamics simulations

Parallelization is an effective way to reduce the computational time needed for molecular dynamics simulations. We describe a new parallelization method, the distributed-diagonal force decomposition method, with which we extend and improve the existing force decomposition methods. Our new method requires less data communication during molecular dynamics simulations than replicated data and current force decomposition methods, increasing the parallel efficiency. It also dynamically load-balances the processors' computational load throughout the simulation. The method is readily implemented in existing molecular dynamics codes and it has been incorporated into the CHARMM program, allowing its immediate use in conjunction with the many molecular dynamics simulation techniques that are already present in the program. We also present the design of the Force Decomposition Machine, a cluster of personal computers and networks that is tailored to running molecular dynamics simulations using the distributed diagonal force decomposition method. The design is expandable and provides various degrees of fault resilience. This approach is easily adaptable to computers with Graphics Processing Units because it is independent of the processor type being used.     

Investigation of chiral recognition by molecular micelles with molecular dynamics simulations

Molecular dynamics simulations were used to characterize the binding of the chiral drugs chlorthalidone and lorazepam to the molecular micelle poly-(sodium undecyl-(L)-leucine-valine). The project’s goal was to characterize the nature of chiral recognition in capillary electrophoresis separations that use molecular micelles as the chiral selector. The shapes and charge distributions of the chiral molecules investigated, their orientations within the molecular micelle chiral binding pockets, and the formation of stereoselective intermolecular hydrogen bonds with the molecular micelle were all found to play key roles in determining where and how lorazepam and chlorthalidone enantiomers interacted with the molecular micelle.     

Nonlinear shear of entangled polymers from nonequilibrium molecular dynamics

This study aims to use molecular dynamics (MD) simulations of Kremer–Grest (KG) chains to inform future developments of models of entangled polymer dynamics. We perform nonequilibrium MD simulations, under shear flow, for well‐entangled KG chains. We study chains of 512 and 1000 KG beads, corresponding to 8 and 15 entanglements, respectively. We compute the linear rheological properties from equilibrium simulations of the stress autocorrelation and obtain from these data the tube model parameters. Under nonlinear shear flow, we compute the shear viscosity, the first and second normal stress differences, and chain contour length. For chains of 512 monomers, we obtain agreement with the results of Cao and Likhtman (ACS Macro. Lett. 2015, 4, 1376). We also compare our nonlinear results with the Graham, Likhtman and Milner‐McLeish (GLaMM) model. We identify some systematic disagreement that becomes larger for the longer chains. We made a comparison of the transient shear stress maximum from our simulations, two nonlinear models and experiments on a wide range of melts and solutions, including polystyrene (PS), polybutadiene, and styrene–butadiene rubber. This comparison establishes that the PS melt data show markedly different behavior to all other melts and solutions and KG simulations reproduce the PS data more closely than either the GLaMM or Xie and Schweizer models. We discuss the performance of these models against the data and simulations. Finally, by imposing a rapid reversing flow, we produce a method to extract the recoverable strain from MD simulations, valid for sufficiently entangled monodisperse polymers. We explore how the resulting data can probe the melt state just before the reversing flow. © 2019 The Authors. Journal of Polymer Science Part B: Polymer Physics published by Wiley Periodicals, Inc. J. Polym. Sci., Part B: Polym. Phys. 2019 , 57, 1692–1704     

Replica sub-permutation method for molecular dynamics and monte carlo simulations

We propose an improvement of the replica-exchange and replica-permutation methods, which we call the replica sub-permutation method (RSPM). Instead of considering all permutations, this method uses a new algorithm referred to as sub-permutation to perform parameter transition. The RSPM succeeds in reducing the number of combinations between replicas and parameters without the loss of sampling efficiency. For comparison, we applied the replica sub-permutation, replica-permutation, and replica-exchange methods to a β-hairpin mini protein, chignolin, in explicit water. We calculated the transition ratio and number of tunneling events in the parameter space, the number of folding–unfolding events, the autocorrelation function, and the autocorrelation time as measures of sampling efficiency. The results indicate that among the three methods, the proposed RSPM is the most efficient in both parameter and conformational spaces. © 2019 Wiley Periodicals, Inc.     

Systematic benchmarking of large molecular dynamics simulations employing GROMACS on massive multiprocessing facilities

The influence of the total number of cores, the number of cores dedicated to Particle mesh Ewald (PME) calculation and the choice of single vs. double precision on the performance of molecular dynamic (MD) simulations in the size of 70,000 to 1.7 million of atoms was analyzed on three different high‐performance computing facilities employing GROMACS 4 by running about 6000 benchmark simulations. Small and medium sized systems scaled linear up to 64 and 128 cores, respectively. Systems with half a million to 1.2 million atoms scaled linear up to 256 cores. The best performance was achieved by dedicating 25% of the total number of cores to PME calculation. Double precision calculations lowered the performance by 30–50%. A database for collecting information about MD simulations and the achieved performance was created and is freely available online and allows the fast estimation of the performance that can be expected in similar environments. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2011     

Pyrite: A blender plugin for visualizing molecular dynamics simulations using industry‐standard rendering techniques

Molecular dynamics (MD) simulations provide critical insights into many biological mechanisms. Programs such as VMD, Chimera, and PyMOL can produce impressive simulation visualizations, but they lack many advanced rendering algorithms common in the film and video‐game industries. In contrast, the modeling program Blender includes such algorithms but cannot import MD‐simulation data. MD trajectories often require many gigabytes of memory/disk space, complicating Blender import. We present Pyrite, a Blender plugin that overcomes these limitations. Pyrite allows researchers to visualize MD simulations within Blender, with full access to Blender's cutting‐edge rendering techniques. We expect Pyrite‐generated images to appeal to students and non‐specialists alike. A copy of the plugin is available at http://durrantlab.com/pyrite/ , released under the terms of the GNU General Public License Version 3. © 2017 Wiley Periodicals, Inc.     

Fast construction of assembly trees for molecular graphs

A number of modeling and simulation algorithms using internal coordinates rely on hierarchical representations of molecular systems. Given the potentially complex topologies of molecular systems, though, automatically generating such hierarchical decompositions may be difficult. In this article, we present a fast general algorithm for the complete construction of a hierarchical representation of a molecular system. This two-step algorithm treats the input molecular system as a graph in which vertices represent atoms or pseudo-atoms, and edges represent covalent bonds. The first step contracts all cycles in the input graph. The second step builds an assembly tree from the reduced graph. We analyze the complexity of this algorithm and show that the first step is linear in the number of edges in the input graph, whereas the second one is linear in the number of edges in the graph without cycles, but dependent on the branching factor of the molecular graph. We demonstrate the performance of our algorithm on a set of specifically tailored difficult cases as well as on a large subset of molecular graphs extracted from the protein data bank. In particular, we experimentally show that both steps behave linearly in the number of edges in the input graph (the branching factor is fixed for the second step). Finally, we demonstrate an application of our hierarchy construction algorithm to adaptive torsion-angle molecular mechanics.     

Charge distributions for molecular dynamics simulations from self-consistent polarization method

Partial atomic charges are important force field parameters. They are usually computed by applying quantum-chemical calculations and the assumed population scheme. In this study polarization consistent scheme of deriving a charge distribution inside solute molecule is proposed. The environment effect is explicitly taken into account by distributing solvent molecules around the solute target. The performed analysis includes a few computational schemes (HF, MP2, B3LYP, and M026X), basis sets (cc-pvnz, n = 2, 3, …, 6), and electrostatically derived charge distributions (KS, CHELP, CHELPG, and HLY). It is demonstrated that the environment effect is very important and cannot be disregarded. The second solvation shell should be included to achieve the charge convergence. Huge corrections to charge distribution are due to induction and dispersion. The B3LYP/cc-pvqz level of theory is recommended for deriving the charges within self-consistent polarization scheme.     

The flying ice cube: Velocity rescaling in molecular dynamics leads to violation of energy equipartition

This article describes an unexpected phenomenon encountered during MD simulations: velocity rescaling using standard protocols can systematically change the proportion of total kinetic energy (KE) found in motions associated with the various degrees of freedom. Under these conditions, the simulation violates the principle of equipartition of energy, which requires a mean kinetic energy of RT/2 in each degree of freedom. A particularly pathological form of this problem occurs if one does not periodically remove the net translation of (and rotation about) the center of mass. In this case, almost all of the kinetic energy is converted into these two kinds of motion, producing a system with almost no kinetic energy associated with the internal degrees of freedom. We call this phenomenon “the flying ice cube.” We present a mathematical analysis of a simple diatomic system with two degrees of freedom, to document the origin of the problem. We then present examples from three kinds of MD simulations, one being an in vacuo simulation on a diatomic system, one involving a low resolution model of DNA in vacuo, and the third using a traditional all-atom DNA model with full solvation, periodic boundary conditions, and the particle mesh Ewald method for treating long-range electrostatics. Finally, we discuss methods for avoiding the problem. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 726–740, 1998     

Force-field parametrization and molecular dynamics simulations of Congo red

Congo red, a diazo dye widely used in medical diagnosis, is known to form supramolecular systems in solution. Such a supramolecular system may interact with various proteins. In order to examine the nature of such complexes empirical force field parameters for the Congo red molecule were developed. The parametrization of bonding terms closely followed the methodology used in the development of the charmm22 force field, except for the calculation of charges. Point charges were calculated from a fit to a quantum mechanically derived electrostatic potential using the CHELP-BOW method. Obtained parameters were tested in a series of molecular dynamics simulations of both a single molecule and a micelle composed of Congo red molecules. It is shown that newly developed parameters define a stable minimum on the hypersurface of the potential energy and crystal and ab initio geometries and rotational barriers are well reproduced. Furthermore, rotations around C-N bonds are similar to torsional vibrations observed in crystals of diphenyl-diazene, which confirms that the flexibility of the molecule is correct. Comparison of results obtained from micelles molecular dynamics simulations with experimental data shows that the thermal dependence of micelle creation is well reproduced.     

On the use of mass scaling for stable and efficient simulated tempering with molecular dynamics

Simulated tempering (ST) is a generalized‐ensemble algorithm that employs trajectories exploring a range of temperatures to effectively sample rugged energy landscapes. When implemented using the molecular dynamics method, ST can require the use of short time steps for ensuring the stability of trajectories at high temperatures. To address this shortcoming, a mass‐scaling ST (MSST) method is presented in which the particle mass is scaled in proportion to the temperature. Mass scaling in the MSST method leads to velocity distributions that are independent of temperature and eliminates the need for velocity scaling after the accepted temperature updates that are required in conventional ST simulations. The homogeneity in time scales with changing temperature improves the stability of simulations and allows for the use of longer time steps at high temperatures. As a result, the MSST is found to be more efficient than the standard ST method, particularly for cases in which a large temperature range is employed. © 2016 Wiley Periodicals, Inc.     

Application of molecular dynamics simulations in molecular property prediction II: diffusion coefficient

In this work, we have evaluated how well the general assisted model building with energy refinement (AMBER) force field performs in studying the dynamic properties of liquids. Diffusion coefficients (D) have been predicted for 17 solvents, five organic compounds in aqueous solutions, four proteins in aqueous solutions, and nine organic compounds in nonaqueous solutions. An efficient sampling strategy has been proposed and tested in the calculation of the diffusion coefficients of solutes in solutions. There are two major findings of this study. First of all, the diffusion coefficients of organic solutes in aqueous solution can be well predicted: the average unsigned errors and the root mean square errors are 0.137 and 0.171 × 10(-5) cm(-2) s(-1), respectively. Second, although the absolute values of D cannot be predicted, good correlations have been achieved for eight organic solvents with experimental data (R(2) = 0.784), four proteins in aqueous solutions (R(2) = 0.996), and nine organic compounds in nonaqueous solutions (R(2) = 0.834). The temperature dependent behaviors of three solvents, namely, TIP3P water, dimethyl sulfoxide, and cyclohexane have been studied. The major molecular dynamics (MD) settings, such as the sizes of simulation boxes and with/without wrapping the coordinates of MD snapshots into the primary simulation boxes have been explored. We have concluded that our sampling strategy that averaging the mean square displacement collected in multiple short-MD simulations is efficient in predicting diffusion coefficients of solutes at infinite dilution.     

Minimum sequence requirements for selective RNA-ligand binding: a molecular mechanics algorithm using molecular dynamics and free-energy techniques

In vitro evolution techniques allow RNA molecules with unique functions to be developed. However, these techniques do not necessarily identify the simplest RNA structures for performing their functions. Determining the simplest RNA that binds to a particular ligand is currently limited to experimental protocols. Here, we introduce a molecular-mechanics based algorithm employing molecular dynamics simulations and free-energy methods to predict the minimum sequence requirements for selective ligand binding to RNA. The algorithm involves iteratively deleting nucleotides from an experimentally determined structure of an RNA-ligand complex, performing energy minimizations and molecular dynamics on each truncated structure, and assessing which truncations do not prohibit RNA binding to the ligand. The algorithm allows prediction of the effects of sequence modifications on RNA structural stability and ligand-binding energy. We have implemented the algorithm in the AMBER suite of programs, but it could be implemented in any molecular mechanics force field parameterized for nucleic acids. Test cases are presented to show the utility and accuracy of the methodology.     

CHARMM fluctuating charge force field for proteins: II protein/solvent properties from molecular dynamics simulations using a nonadditive electrostatic model

A fluctuating charge (FQ) force field is applied to molecular dynamics simulations for six small proteins in explicit polarizable solvent represented by the TIP4P-FQ potential. The proteins include 1FSV, 1ENH, 1PGB, 1VII, 1H8K, and 1CRN, representing both helical and beta-sheet secondary structural elements. Constant pressure and temperature (NPT) molecular dynamics simulations are performed on time scales of several nanoseconds, the longest simulations yet reported using explicitly polarizable all-atom empirical potentials (for both solvent and protein) in the condensed phase. In terms of structure, the FQ force field allows deviations from native structure up to 2.5 A (with a range of 1.0 to 2.5 A). This is commensurate to the performance of the CHARMM22 nonpolarizable model and other currently existing polarizable models. Importantly, secondary structural elements maintain native structure in general to within 1 A (both helix and beta-strands), again in good agreement with the nonpolarizable case. In qualitative agreement with QM/MM ab initio dynamics on crambin (Liu et al. Proteins 2001, 44, 484), there is a sequence dependence of average condensed phase atomic charge for all proteins, a dependence one would anticipate considering the differing chemical environments around individual atoms; this is a subtle quantum mechanical feature captured in the FQ model but absent in current state-of-the-art nonpolarizable models. Furthermore, there is a mutual polarization of solvent and protein in the condensed phase. Solvent dipole moment distributions within the first and second solvation shells around the protein display a shift towards higher dipole moments (increases on the order of 0.2-0.3 Debye) relative to the bulk; protein polarization is manifested via the enhanced condensed phase charges of typical polar atoms such as backbone carbonyl oxygens, amide nitrogens, and amide hydrogens. Finally, to enlarge the sample set of proteins, gas-phase minimizations and 1 ps constant temperature simulations are performed on various-sized proteins to compare to earlier work by Kaminsky et al. (J Comp Chem 2002, 23, 1515). The present work establishes the feasibility of applying a fully polarizable force field for protein simulations and demonstrates the approach employed in extending the CHARMM force field to include these effects.     

Understanding the adsorption of branched polyamine on surface of gold nanoparticles by molecular dynamics simulations

Molecular dynamics simulations were conducted to characterize the adsorption behavior of branched polyethylenimine (br‐PEI) on the surface of gold nanoparticles (AuNPs). We observed the preferential adsorption of br‐PEI on the [111] surface of AuNPs. Furthermore, br‐PEI maintained a flat arrangement on the surface and wrapped the AuNPs in a stable manner, thereby blocking the adsorption of H₂O molecules and other free br‐PEI molecules. The model and computational results provide theoretical support for relevant experiments. Copyright © 2016 John Wiley & Sons, Ltd.