Classical molecular dynamics

This issue of JCP highlights both developments in and applications of classical molecular simulation in 67 articles. A recent issue of JCP focused on electronic structure software, covering a few dozen different packages with different features.¹ The present issue includes discussion of some of the major software packages for classical molecular dynamics (MD), including Amber, CHARMM, GROMACS, LAMMPS, and NAMD. But even more attention is given to development of underlying models and algorithms. These two special issues of JCP together highlight the breadth of research where everything from interactions of light with matter to large biological systems can now be treated if appropriate models and simplifications are applied.

The first MD paper was published in JCP in 1957 by Alder and Wainwright on phase transitions in hard sphere systems,² followed by a paper by the same authors on the methodology of MD.³ Sadly, Bernie Alder passed away just a few months ago in September at the age of 94. Indeed, many of the breakthroughs in the field of classical molecular dynamics have been published in JCP. Analytical theories by the likes of Kirkwood⁴ from the early days of JCP’s history continue to be relevant and still frequently cited in the classical simulation literature. Other breakthroughs that deserve mention are the first simulation of liquid argon using a continuous potential by Rahman⁵ and the first simulations of liquid water by Rahman and Stillinger.⁶ In the sections below, we have organized the contents of the present issue into a number of subtopics even though the categorization in some cases was somewhat arbitrary.

A number of articles in the special issue report on recent algorithmic implementations and advances in performance of commonly used MD engines on modern hardware architectures, most prominently on Graphical Processing Unit (GPU)-accelerated hardware architectures. These include a detailed survey of features and capabilities and recent GPU acceleration of NAMD,⁷ a review of major algorithmic advances in GROMACS focusing on taking full advantage of GPU acceleration,⁸ implementation of GPU-accelerated Martini coarse-grained simulations in ddcMD,⁹ and parallel implementation of hyperdynamics in LAMMPS¹⁰ to perform accelerated MD simulations for solid-state systems. Gecht et al. report on a software tool, MDbenchmark, developed to set up and analyze benchmarks for diverse platforms.¹¹ Of particular interest to faster calculation of electrostatic forces, Predescu et al. report on an optimized charge assignment by the particle mesh calculations,¹² a commonly used method for calculation of electrostatic forces.

Force field development remains an important subject within the molecular simulation community. In this issue, Mohebifar and Rowley report on a new TIP4P-like water model using a Buckingham potential plus an additional r⁻⁸ dispersion term.¹³ Teng and co-workers analyze how different water force field models reproduce the anomalous increase in water diffusivity with increasing (high) pressure.¹⁴ van der Spoel et al. introduce a stable classical potential for linear chemical moieties.¹⁵ Evangelisti et al.¹⁶ as well as O’Hearn et al.¹⁷ developed new models for the Reax reactive force field for gold nanoclusters, respectively, lithium-oxide, whereas Slapikas et al.¹⁸ improved the implementation of the reactive COMB3 potential in the LAMMPS package. A number of extensions to the Amber force field family are introduced. Schepers and Gohlke implement popular dyes and linker molecules¹⁹ and He et al. describe a refinement of the charge calculations for the General Amber force field (GAFF),²⁰ while Bogetti and co-workers introduce an Amber version for protein-mimetic compounds.²¹ Lee et al. introduce support for the Amber force field in the CHARMM graphical user interface (CHARMM-GUI), a tool for simulation preparation.²² Since the concept of sigma-holes and halogen bonding became well known, a steady stream of papers on the topic has been published. Here, Pang et al. describe the implementation of automatic parameterization of halogen bonding to the force field toolkit (ffTk),²³ while Campetella et al. describe an automatic parameterization strategy for compounds involved in halogen bonds based on quantum mechanics.²⁴ Where halogen bonds are typically modeled using a virtual site, Ramasubramani et al. take this one step further and introduce a means to describe anisotropic particles more generally in, e.g., a Lennard-Jones potential.²⁵ Refinement of the Amoeba+ polarizable models to include charge-flux is described by Yang et al.²⁶ Amin et al. evaluate force fields including the polarizable Drude model for peptide–Ca²⁺ interactions.²⁷ Wei et al. introduce an efficient algorithm for computation of polarizable Gaussian multipole electrostatics.²⁸ Finally, Karls and co-workers describe OpenKIM, a cloud-based property calculator based on multiple different interatomic models.²⁹ 

The development of approaches to enhance sampling of complex systems and to accurately determine free energies associated with the status of complex systems continues to be a particularly active area. Enhanced sampling remains of paramount importance, with a large number of contributions devoted to this issue. Takaba et al. propose an edge expansion parallel cascade selection MD (eePaCS-MD) to investigate the large-amplitude motions of proteins without prior knowledge of the conformational transitions.³⁰ eePaCS-MD belongs to the class of adaptive enhanced sampling methods that rely on iteratively initiating independent unbiased MD simulations from structures selected from previous runs according to specific rules. A different class of enhanced sampling methods rely on biasing factors that alter the normal dynamical propagation. Peter et al. present the path correlated MD (CORE-MD) method, which biases the propagation of the system using information from the autocorrelation of the path integral over the reduced action.³¹ Kokh et al. propose the Random Acceleration MD Simulation (RAMD) to capture the slow dissociation of ligands from proteins and other macromolecules by the introduction of a small randomly oriented force applied to the center of mass (COM) of the ligand during the simulations.³² A related method is the Peptide Gaussian accelerated molecular dynamics (Pep-GaMD) proposed by Wang and Miao to enhance the sampling of peptide binding and calculate binding free energies and kinetic rate constants.³³ Hartmann et al. propose the Infinite Switch Simulated Tempering in Force (FISST), an enhanced sampling tool that is based on introducing a linear biasing force on a pre-determined collective variable and treating the bias as a replica-exchange simulated tempering in the infinite swapping limit.³⁴ Karabin and Stuart explore simulated annealing with adaptive cooling rates, designed to make use of the heat capacity of the system for optimizing atomistic structures.³⁵ Wu and Brooks propose a reformulation of the Self-Guided Langevin Dynamics (SGLD) simulation methods to enhance conformational sampling through a non-Hamiltonian propagation promoting low frequency motions.³⁶ Ovchinnikov et al. present the Restrained Locally Enhanced Sampling (RLES) method.³⁷ This type of strategy introduces multiple replicas along selected degrees of freedom to boost the statistical sampling. Here, a restraint potential is introduced to drive the many-replica system to the canonical ensemble corresponding to the physical, single-replica system. RLES is found to explore the space of configurations with an efficiency comparable to that of temperature replica exchange.

Several contributions focus on fundamental issues regarding dynamical propagation. Monmarché et al. describe the BOUNCE integrator, a velocity jump process allowing for an adaptive timestep specific to each atom–atom pair as an alternative to multi-timestep algorithms for fast and accurate MD simulations.³⁸ Bernetti and Bussi describe a first-order barostat, a stochastic cell-rescaling constant pressure control algorithm that samples volume fluctuations correctly, which can be easily implemented in existing codes.³⁹ Finkelstein et al. explore the necessary conditions for a statistically accurate discrete-time integration of Langevin equations reproducing correct configurational Boltzmann sampling and free-particle diffusion.⁴⁰ Rupakheti et al. examine the fundamental significance of the dual-thermostat extended Lagrangian propagation that is used to treat the inducible degrees of freedom dynamically in the context of a polarizable force field based on classical Drude oscillators by formulating the statistical mechanics of an extended canonical ensemble.⁴¹ Cherniavskyi et al. explore the ability of hybrid Monte Carlo (MC)-MD enhanced sampling techniques called MDAS (MD with Alchemical Steps) for characterizing the phase diagram of a coarse-grained model (Martini) of heterogeneous phospholipid membranes.⁴² The MC-MD combination is shown to be a useful approach to enhance the study of phase separation and microdomain formation in biological membranes.

Two contributions consider hybrid schemes combining MC sampling of the discrete state in MD trajectories. Grünewald et al. describe a titratable Martini model for constant pH MD simulations,⁴³ and Michael et al. describe a hybrid MC-MD for handling random amino acid mutations and substitutions in protein design—the latter being also handled to treat protonation events.⁴⁴ Quantitative determination of free energies in complex systems continues to be important. Urano et al. explore the treatment of long-range Coulombic energy calculation for net charged systems using the fast multipole method in free energy calculations,⁴⁵ and Ito and Cui introduce a novel staged transformation protocol for multi-level free energy simulations⁴⁶ used primarily to connect a low-level molecular mechanical force field model to a high-level quantum mechanical representation.

A critically important area of active research concerns the development of methods to determine the transition pathways in slow processes. Wang and Elber assess the impact of a biasing potential introduced in exact calculations of kinetics based on the algorithm WARM (Wind-Assisted Reweighted Milestoning).⁴⁷ The use of a biasing force (wind) significantly enhances the sampling of otherwise rare trajectories but also introduces an exponential weight to the trajectories that significantly impacts the value of the statistics. It is suggested that the algorithm is most useful to observe rare events occurring on steep landscapes. Wu et al. introduce a method to optimize a conformational pathway through a space of well-chosen reduced variables to characterize large conformational transitions in proteins.⁴⁸ The adaptively biased path optimization strategy utilizes unrestricted, enhanced sampling in the region of a path in the reduced variable space to identify a broad path between two stable end states. The pathway method is applied to the inactivation transition of Src tyrosine kinase catalytic domain, providing new insights into this well-studied conformational equilibrium. Ribeiro and co-workers apply a combination of machine learning and infrequent metadynamics to study kinetics of drug dissociation.⁴⁹ 

Those who develop MD engines understandably put significant effort into optimizing simulation parameters and performance, but what comes before and after a simulation contributes a lot to obtaining useful results. To paraphrase Tolstoy, simulation methods are (largely) alike, but each type of analysis is used in its own way.⁵⁰ On the very general side, Roe and Brooks describe a ten-step protocol for preparation and equilibration of explicit solvent simulations,⁵¹ and Brehm et al. give an overview of TRAVIS, a package with many novel methods for post-processing and analyzing MD trajectories.⁵² 

Other contributions describe specific types of analyses: Bey et al. discuss how to measure the line tension at a solid–liquid–vapor contact line.⁵³ Hsu et al. present methods to bias MD simulations to agree with x-ray solution scattering data.⁵⁴ Melo et al. use dynamic network analysis to identify important residues and information pathways in allosteric enzymes.⁵⁵ Monroe et al. use water as an example to explore methods of thermodynamic extrapolation of structural properties to temperatures and densities that differ from those chosen for simulation.⁵⁶ Mori et al. present a cross-entropy minimization method for finding the reaction coordinate from a large number of collective variables in complex molecular systems and show applications to the alanine dipeptide.⁵⁷ Rice et al. present a new approach to identify chemical species from MD simulations of reacting materials under extreme temperatures and pressures.⁵⁸ Naullage et al. discuss how to identify ice-binding surfaces.⁵⁹ A further paper opines on solutions to challenges facing the biomolecular (and materials) modeling community in validation and reproducibility, data portability, and wider use of coarse-graining and multi-scale tools.⁶⁰ 

The application of MD to biological systems has a long and rich history. The Nobel Prize in Chemistry 2013 was awarded jointly to Martin Karplus, Michael Levitt, and Arieh Warshel “for the development of multiscale models for complex chemical systems,” best exemplified by biomolecular systems and processes. Several papers investigating biomolecular systems are included in the special issue.

Two papers report on the timely topic of virus modeling: a study of the catalytic and conformational dynamics of the protease PLpro in several coronaviruses, to inform design of inhibitors,⁶¹ and simulation of the helical assembly of nucleocapsid proteins in the Ebola virus, to understand what molecules make it stable and can thus be potential therapeutic targets.⁶² Another paper describes the use of MD simulations in characterizing the interactions between signaling pathway protein SH3 domains and proline-rich segments to construct free energy landscapes and validate binding sites.⁶³ 

Given their partially ordered structure, biological membranes have been consistently among the best suited biomolecular systems characterized by classical MD simulations, particularly since their fluid structures might not be fully probed experimentally. The special issue includes a number of articles focusing on investigation of various molecular properties of membranes with MD simulations. Some articles focus primarily on permeation of molecular species across the lipid bilayer^64,65 or diffusion of ions across a sugar layer on the surface of endothelial cells induced by applying an external electric field.⁶⁶ These studies establish that the calculated permeability from MD simulations can be largely affected by the nature and conditions of the simulations, by the technique used for quantifying permeation events and even by the lipid composition.^64,65 In another study,⁶⁷ MD simulations of asymmetric lipid bilayers are used to show how the resulting imbalance in counterion distribution might contribute to the membrane potential. Free energy perturbation of formation of controlled cylindrical defects in lipid bilayers has been employed by Kumar et al. to quantify the energetics associated with the formation of defects and pores in lipid bilayers and their dependence on the lipid composition.⁶⁸ Dixit and Lazaridis return to the old problem of pore formation in lipid bilayers using free energy calculations.⁶⁹ Finally, Martini coarse-grained MD simulations are used to somewhat address the challenging case of lipopolysaccharides in the outer membrane and their slow mixing.⁷⁰ 

The popularity of MD as a method for modeling materials more broadly has grown dramatically. Studies of four such systems are included in the special issue: (1) use of reactive MD to model the precipitation of colloidal gels, which act as a binding phase in concrete;⁷¹ (2) analyses of the structure and ordering of water at TiO₂ surfaces, relevant to photoelectrochemical dissociation of water for photocatalysis;⁷² (3) identification of creep mechanisms in oriented, crystalline polyethylene fibers occurring via dislocation nucleation;⁷³ and (4) simulation of phase transitions in multi-component simple liquids, with application to microphase separation in block copolymer melts.⁷⁴ 

The present special issue highlights that development of models and tools remains an important part of the field of classical MD simulations. JCP is one of the strongholds for these kinds of articles, and many such articles have historically been very well-cited. It should be mentioned that a large fraction of the community focuses rather on applications of these tools, in broad areas ranging from characterization and design of materials, to modeling and simulations of complex biomolecular systems and assemblies, to drug design. All these applications rely on fundamental algorithms and physical models, including efficient and faster implementations on modern hardware, such as those described here. The future of MD simulation as a complementary method to experiment is highly promising, and we can only expect more studies in which molecular simulations of increasing accuracy, realism, and quantification will provide novel insight into the mechanisms underlying diverse atomic-scale processes.