Search

This article considers the application of Langevin dynamics to sampling and investigates how to choose the damping parameter in Langevin dynamics for the purpose of maximizing thoroughness of sampling. Also, it considers the computation of measures of sampling thoroughness., Comment: 21 pages, 6 figures

International audience; NAMD is a molecular dynamics program designed for high-performance simulations of very large biological objects on central processing unit (CPU)-and graphics processing unit (GPU)-based architectures. NAMD offers scalable performance on petascale parallel supercomputers consisting of hundreds of thousands of cores, as well as on inexpensive commodity clusters commonly found in academic environments. It is written in C++ and leans on Charm++ parallel objects for optimal performance on low-latency architectures. NAMD is a versatile, multipurpose code that gathers state-of-the-art algorithms to carry out simulations in apt thermodynamic ensembles, using the widely popular CHARMM, AMBER, OPLS and GROMOS biomolecular force fields. Here, we review the main features of NAMD that allow both equilibrium and enhanced-sampling molecular dynamics simulations with numerical efficiency. We describe the underlying concepts utilized by NAMD and their implementation, most notably for handling long-range electrostatics, controlling the temperature, pressure and pH, applying external potentials on tailored grids, leveraging massively parallel resources in multiple-copy simulations, as well as hybrid QM/MM descriptions. We detail the variety of options offered by NAMD for enhanced-sampling simulations aimed at determining free-energy differences of either alchemical or geometrical transformations, and outline their applicability to specific problems. Last, we discuss the roadmap for the development of NAMD and our current efforts towards achieving optimal performance on GPUbased architectures, for pushing back the limitations that have prevented biologically realistic billion-atom objects to be fruitfully simulated, and for making large-scale simulations less expensive and easier to set up, run and analyze. NAMD is distributed free of charge with its source code at www.ks.uiuc.edu.

Driving molecular dynamics simulations with data-guided collective variables offer a promising strategy to recover thermodynamic information from structure-centric experiments. Here, the 3-dimensional electron density of a protein, as it would be determined by cryo-EM or X-ray crystallography, is used to achieve simultaneously free-energy costs of conformational transitions and refined atomic structures. Unlike previous density-driven molecular dynamics methodologies that determine only the best map-model fits, our work uses the recently developed Multi-Map methodology to monitor concerted movements within equilibrium, non-equilibrium, and enhanced sampling simulations. Construction of all-atom ensembles along chosen values of the Multi-Map variable enables simultaneous estimation of average properties, as well as real-space refinement of the structures contributing to such averages. Using three proteins of increasing size, we demonstrate that biased simulation along reaction coordinates derived from electron densities can serve to induce conformational transitions between known intermediates. The simulated pathways appear reversible, with minimal hysteresis and require only low-resolution density information to guide the transition. The induced transitions also produce estimates for free energy differences that can be directly compared to experimental observables and population distributions. The refined model quality is superior compared to those found in the Protein DataBank. We find that the best quantitative agreement with experimental free-energy differences is obtained using medium resolution (~5 Å) density information coupled to comparatively large structural transitions. Practical considerations for generating transitions with multiple intermediate atomic density distributions are also discussed.

Fast methods for calculating two-body interactions have many applications, and for molecular science and cosmology, it is common to employ periodic boundary conditions. However, for the 1/r potential, the energy and forces are ill-defined. Adopted here is the model given by the classic Ewald sum. For the fast calculation of two-body forces, the most celebrated method is the fast multipole method and its tree-code predecessor. However, molecular simulations typically employ mesh-based approximations and the fast Fourier transform. Both types of methods have significant drawbacks, which, in most respects, are overcome by the less well-known multilevel summation method (MSM). Presented here is a realization of the MSM, which can be regarded as a multilevel extension of the (smoothed) particle mesh Ewald (PME) method, but with the Ewald softening replaced by one having a finite range. The two-level (single-grid) version of MSM requires fewer tuning parameters than PME and is marginally faster. Additionally, higher-level versions of MSM scale well to large numbers of processors, whereas PME and other two-level methods do not. Although higher-level versions of MSM are less efficient on a single processor than the two-level version, evidence suggests that they are more efficient than other methods that scale well, such as the fast multipole method and tree codes.

Periodicity is a popular way to handle boundaries in atomistic simulations in the bulk. However, slowly decaying Coulomb interactions lead to ill-defined energies and forces. The Ewald sum is a common way of giving meaning to these terms, but its derivation, based on Fourier series, is not transparent. Described in this article is a derivation based on the limiting behaviour as the cutoff goes to infinity of a truncated and shifted Coulomb potential. This observation strengthens the case for the use of the original Ewald sum. And such a construction may offer advantages for computation and/or modelling.

The choice of the parameter value for regularized inverse problems is critical to the results and remains a topic of interest. This article explores a criterion for selecting a good parameter value by maximizing the probability of the data, with no prior knowledge of the noise variance. These concepts are developed for l 2 and consequently l 1 regularization models by way of their Bayesian interpretations. Based on these concepts, an iterative scheme is proposed and demonstrated to converge accurately, and analytical convergence results are provided that substantiate these empirical observations. For some of the most common inverse problems, including MRI, SAR, denoising, and deconvolution, an extremely efficient algorithm is derived, making the iterative scheme very attractive for real case use. The computational concerns associated with the general case for any inverse problem are also carefully addressed. A robust set of 1D and 2D numerical simulations confirm the effectiveness of the proposed approach.

The article by Sarich and Schutte takes a broad view of the use of Markov state models, with biomolecules as the intended application, and proposes novel approaches based on mathematical insights.

Considered is the construction of transition paths of conformational changes for proteins and other macromolecules, using methods that do not require the generation of dynamics trajectories. Special attention is given to the use of a reduced set of collective variables for describing such paths. A favored way to define transition paths is to seek channels through the transition state having cross sections with a high reactive flux (density of last hitting points of reactive trajectories). Given here is a formula for reactive flux that is independent of the parameterization of “collective variable space.” This formula is needed for the principal curve of the reactive flux (as in the revised finite temperature string method) and for the maximum flux transition (MaxFlux) path. Additionally, a resistance functional is derived for narrow tubes, which when minimized yields a MaxFlux path. A strategy for minimization is outlined in the spirit of the string method. Finally, alternative approaches based on determining trajectories of high probability are considered, and it is observed that they yield paths that depend on the parameterization of collective variable space, except in the case of zero temperature, where such a path coincides with a MaxFlux path.

The multilevel summation method for calculating electrostatic interactions in molecular dynamics simulations constructs an approximation to a pairwise interaction kernel and its gradient, which can be evaluated at a cost that scales linearly with the number of atoms. The method smoothly splits the kernel into a sum of partial kernels of increasing range and decreasing variability with the longer-range parts interpolated from grids of increasing coarseness. Multilevel summation is especially appropriate in the context of dynamics and minimization, because it can produce continuous gradients. This article explores the use of B-splines to increase the accuracy of the multilevel summation method (for nonperiodic boundaries) without incurring additional computation other than a preprocessing step (whose cost also scales linearly). To obtain accurate results efficiently involves technical difficulties, which are overcome by a novel preprocessing algorithm. Numerical experiments demonstrate that the resulting method offers substantial improvements in accuracy and that its performance is competitive with an implementation of the fast multipole method in general and markedly better for Hamiltonian formulations of molecular dynamics. The improvement is great enough to establish multilevel summation as a serious contender for calculating pairwise interactions in molecular dynamics simulations. In particular, the method appears to be uniquely capable for molecular dynamics in two situations, nonperiodic boundary conditions and massively parallel computation, where the fast Fourier transform employed in the particle–mesh Ewald method falls short.

Spline interpolation is a splendid tool for multiscale approximation on unbounded domains. In particular, it is well suited for use by the multilevel summation method (MSM) for calculating a sum of pairwise interactions for a large set of particles in linear time. Outlined here is an algorithm for spline interpolation on unbounded domains that is efficient and elegant though not so simple. Further gains in efficiency are possible via quasi-interpolation, which compromises collocation but with minimal loss of accuracy. The MSM, which may also be of value for continuum models, embodies most of the best features of both hierarchical clustering methods (tree methods, fast multipole methods, hierarchical matrix methods) and FFT-based 2-level methods (particle–particle particle–mesh methods, particle–mesh Ewald methods).

Atomic oscillations present in classical molecular dynamics restrict the step size that can be used. Multiple time stepping schemes offer only modest improvements, and implicit integrators are costly and inaccurate. The best approach may be to actually remove the highest frequency oscillations by constraining bond lengths and bond angles, thus permitting perhaps a 4-fold increase in the step size. However, omitting degrees of freedom produces errors in statistical averages, and rigid angles do not bend for strong excluded volume forces. These difficulties can be addressed by an enhanced treatment of holonomic constrained dynamics using ideas from papers of Fixman (1974) and Reich (1995, 1999). In particular, the 1995 paper proposes the use of “flexible” constraints, and the 1999 paper uses a modified potential energy function with rigid constraints to emulate flexible constraints. Presented here is a more direct and rigorous derivation of the latter approach, together with justification for the use of constraints in molecular modeling. With rigor comes limitations, so practical compromises are proposed: simplifications of the equations and their judicious application when assumptions are violated. Included are suggestions for new approaches.

Binding free energies are calculated for the phosphorylated and unphosphorylated complexes between the kinase inducible domain (KID) of the DNA transcriptional activator cAMP response element binding (CREB) protein and the KIX domain of its coactivator, CREB-binding protein (CBP). To our knowledge, this is the first application of a method based on a potential of mean force (PMF) with restraining potentials to compute the binding free energy of protein-protein complexes. The KID:KIX complexes are chosen here because of their biological relevance to the DNA transcription process and their relatively small size (81 residues for the KIX domain of CBP, and 28 residues for KID). The results for pKID:KIX and KID:KIX are −9.55 and −4.96 kcal/mol, respectively, in good agreement with experimental estimates (−8.8 and −5.8 kcal/mol, respectively). A comparison between specific contributions to protein-protein binding for the phosphorylated and unphosphorylated complexes reveals a dual role for the phosphorylation of KID at Ser-133 in effecting a more favorable free energy of the bound system: 1), stabilization of the unbound conformation of phosphorylated KID due to favorable intramolecular interactions of the phosphate group of Ser-133 with the charged groups of an arginine-rich region spanning both α-helices, which lowers the configurational entropy; and 2), more favorable intermolecular electrostatic interactions between pSer-133 and Arg-131 of KID, and Lys-662, Tyr-658, and Glu-666 of KIX. Charge reduction through ligand phosphorylation emerges as a possible mechanism for controlling the unbound state conformation of KID and, ultimately, gene expression. This work also demonstrates that the PMF-based method with restraining potentials provides an added benefit in that important elements of the binding pathway are evidenced. Furthermore, the practicality of the PMF-based method for larger systems is validated by agreement with experiment. In addition, we provide a somewhat differently structured exposition of the PMF-based method with restraining potentials and outline its generalization to systems in which both protein and ligand may adopt unbound conformations that are different from those of the bound state.

In long-time numerical integration of Hamiltonian systems, and especially in molecular dynamics simulation, it is important that the energy is well conserved. For symplectic integrators applied with sufficiently small step size, this is guaranteed by the existence of a modified Hamiltonian that is exactly conserved up to exponentially small terms. This article is concerned with the simplified Takahashi--Imada method, which is a modification of the St"ormer--Verlet method that is as easy to implement but has improved accuracy. This integrator is symmetric and volume-preserving, but no longer symplectic. We study its long-time energy conservation and give theoretical arguments, supported by numerical experiments, which show the possibility of a drift in the energy (linear or like a random walk). With respect to energy conservation, this article provides empirical and theoretical data concerning the importance of using a symplectic integrator.

The multilevel summation method (MSM) offers an efficient algorithm utilizing convolution for evaluating long-range forces arising in molecular dynamics simulations. Shifting the balance of computation and communication, MSM provides key advantages over the ubiquitous particle–mesh Ewald (PME) method, offering better scaling on parallel computers and permitting more modeling flexibility, with support for periodic systems as does PME but also for semiperiodic and nonperiodic systems. The version of MSM available in the simulation program NAMD is described, and its performance and accuracy are compared with the PME method. The accuracy feasible for MSM in practical applications reproduces PME results for water property calculations of density, diffusion constant, dielectric constant, surface tension, radial distribution function, and distance-dependent Kirkwood factor, even though the numerical accuracy of PME is higher than that of MSM. Excellent agreement between MSM and PME is found also for interface potentials of air–water and membrane–water interfaces, where long-range Coulombic interactions are crucial. Applications demonstrate also the suitability of MSM for systems with semiperiodic and nonperiodic boundaries. For this purpose, simulations have been performed with periodic boundaries along directions parallel to a membrane surface but not along the surface normal, yielding membrane pore formation induced by an imbalance of charge across the membrane. Using a similar semiperiodic boundary condition, ion conduction through a graphene nanopore driven by an ion gradient has been simulated. Furthermore, proteins have been simulated inside a single spherical water droplet. Finally, parallel scalability results show the ability of MSM to outperform PME when scaling a system of modest size (less than 100 K atoms) to over a thousand processors, demonstrating the suitability of MSM for large-scale parallel simulation.

The application of a symplectic integrator to a Hamiltonian system formally conserves the value of a modified, or shadow, Hamiltonian defined by some asymptotic expansion in powers of the step size. An earlier article describes how it is possible to construct highly accurate shadow Hamiltonian approximations using information readily available from the numerical integration. This article improves on this construction by giving formulas of order up to 24 (not just up to 8) and by greatly reducing both storage requirements and roundoff error. More significantly, these high order formulas yield remarkable results not evident for 8th order formulas, even for systems as complex as the molecular dynamics of water. These numerical experiments not only illuminate theoretical properties of shadow Hamiltonians but also give practical information about the accuracy of a simulation. By removing systematic energy fluctuations, they reveal the rate of energy drift for a given step size and uncover the ill effects of using switching functions that do not have enough smoothness.

NAMD is a parallel molecular dynamics code designed for high-performance simulation of large biomolecular systems. NAMD scales to hundreds of processors on high-end parallel platforms, as well as tens of processors on low-cost commodity clusters, and also runs on individual desktop and laptop computers. NAMD works with AMBER and CHARMM potential functions, parameters, and file formats. This article, directed to novices as well as experts, first introduces concepts and methods used in the NAMD program, describing the classical molecular dynamics force field, equations of motion, and integration methods along with the efficient electrostatics evaluation algorithms employed and temperature and pressure controls used. Features for steering the simulation across barriers and for calculating both alchemical and conformational free energy differences are presented. The motivations for and a roadmap to the internal design of NAMD, implemented in C++ and based on Charm++ parallel objects, are outlined. The factors affecting the serial and parallel performance of a simulation are discussed. Finally, typical NAMD use is illustrated with representative applications to a small, a medium, and a large biomolecular system, highlighting particular features of NAMD, for example, the Tcl scripting language. The article also provides a list of the key features of NAMD and discusses the benefits of combining NAMD with the molecular graphics/sequence analysis software VMD and the grid computing/collaboratory software BioCoRE. NAMD is distributed free of charge with source code at www.ks.uiuc.edu.

A reaction probability is required to calculate the rate constant of a diffusion-dominated reaction. Due to the complicated geometry and potentially high dimension of the reaction probability problem, it is usually solved by a Brownian dynamics simulation, also known as a random walk or path integral method, instead of solving the equivalent partial differential equation by a discretization method. Building on earlier work, this article completes the development of a robust importance sampling algorithm for Brownian dynamics—i.e., biased Brownian dynamics with weight control—to overcome the high energy and entropy barriers in biomolecular association reactions. The biased Brownian dynamics steers sampling by a bias force, and the weight control algorithm controls sampling by a target weight. This algorithm is optimal if the bias force and the target weight are constructed from the solution of the reaction probability problem. In reality, an approximate reaction probability has to be used to construct the bias force and the target weight. Thus, the performance of the algorithm depends on the quality of the approximation. Given here is a method to calculate a good approximation, which is based on the selection of a reaction coordinate and the variational formulation of the reaction probability problem. The numerically approximated reaction probability is shown by computer experiments to give a factor-of-two speedup over the use of a purely heuristic approximation. Also, the fully developed method is compared to unbiased Brownian dynamics. The tests for human superoxide dismutase, Escherichia coli superoxide dismutase, and antisweetener antibody NC6.8, show speedups of 17, 35, and 39, respectively. The test for reactions between two model proteins with orientations shows speedups of 2578 for one set of configurations and 3341 for another set of configurations.

We study the position recurrence relation of several existing numerical integrators for the Langevin equation and use the modified equation approach to analyse their accuracy. We show that for the harmonic oscillator, the BBK integrator converges weakly with order 1 while the vGB82 and Langevin impulse (LI)‡ integrator converge weakly with order 2. We also study a restricted class of velocity definitions—those that lead to explicit starting procedures. We show that some recurrence relations exact for constant force, can achieve the exact virial relation by a proper definition of velocity, extending the result of Pastor et al. on the analysis of BBK integrators in 1988.

This paper shows that in molecular dynamics (MD) when constant- energy (NVE) simulations of Newton's equations of motion are attempted using the multiple time stepping (MTS) integrator Verlet-I/r-RESPA/Impulse, there are nonlinear instabilities when the longest step size is a third or possibly a fourth of the period(s) of the fastest motion(s) in the system. This is demonstrated both through a thorough set of computer experiments and through the analysis of a nonlinear model problem. The numerical experiments include not only the unconstrained dynamics simulation of a droplet of flexible water and a flexible protein, but also the constrained dynamics simulation of a solvated protein, representing a range of simulation protocols commonly in use by biomolecular modelers. The observed and predicted instabilities match exactly. Previous work has identified and explained a linear instability for Verlet-I/r-RESPA/Impulse at around half the period of the fastest motion. Mandziuk and Schlick discovered nonlinear resonances in single time stepping MD integrators, but unstable nonlinear resonances for MTS integrators are reported here for the first time. This paper also offers an explanation on the instability of MTS constrained molecular dynamics simulations of explicitly solvated proteins. More aggressive multiple step sizes are possible with mild Langevin coupling or targeted Langevin coupling, and its combination with the mollified Impulse method permits step sizes 3 to 4 times larger thanVerlet-I/r-RESPA/Impulse while still retaining some accuracy.

The best simple method for Newtonian molecular dynamics is indisputably the leapfrog Stormer-Verlet method. The appropriate generalization to simple Langevin dynamics is unclear. An analysis is presented comparing an ‘impulse method’ (kick; fluctuate; kick), the 1982 method of van Gunsteren and Berendsen, and the Brunger-Brooks-Karplus (BBK) method. It is shown how the impulse method and the van Gunsteren-Berendsen methods can be implemented as efficiently as the BBK method. Other considerations suggest that the impulse method is the best basic method for simple Langevin dynamics, with the van Gunsteren-Berendsen method a close contender.

In this paper we show the possibility of using very mild stochastic damping to stabilize long time step integrators for Newtonian molecular dynamics. More specifically, stable and accurate integrations are obtained for damping coefficients that are only a few percent of the natural decay rate of processes of interest, such as the velocity autocorrelation function. Two new multiple time stepping integrators, Langevin Molly (LM) and Brunger–Brooks–Karplus–Molly (BBK–M), are introduced in this paper. Both use the mollified impulse method for the Newtonian term. LM uses a discretization of the Langevin equation that is exact for the constant force, and BBK–M uses the popular Brunger–Brooks–Karplus integrator (BBK). These integrators, along with an extrapolative method called LN, are evaluated across a wide range of damping coefficient values. When large damping coefficients are used, as one would for the implicit modeling of solvent molecules, the method LN is superior, with LM closely following. However, with mild damping of 0.2 ps−1, LM produces the best results, allowing long time steps of 14 fs in simulations containing explicitly modeled flexible water. With BBK–M and the same damping coefficient, time steps of 12 fs are possible for the same system. Similar results are obtained for a solvated protein–DNA simulation of estrogen receptor ER with estrogen response element ERE. A parallel version of BBK–M runs nearly three times faster than the Verlet-I/r-RESPA (reversible reference system propagator algorithm) when using the largest stable time step on each one, and it also parallelizes well. The computation of diffusion coefficients for flexible water and ER/ERE shows that when mild damping of up to 0.2 ps−1 is used the dynamics are not significantly distorted.

One of the most fruitful ways to analyze the effects of discretization error in the numerical solution of a system of differential equations is to examine the "modified equations," which are equations that are exactly satisfied by the (approximate) discrete solution. These do not actually exist in general but rather are defined by an asymptotic expansion in powers of the discretization parameter. Nonetheless, if the expansion is suitably truncated, the resulting modified equations have a solution which is remarkably close to the discrete solution. In the case of a Hamiltonian system of ordinary differential equations, the modified equations are also Hamiltonian if and only if the integrator is symplectic. Evidence for the existence of a Hamiltonian for a particular calculation is obtained by calculating modified Hamiltonians and monitoring how well they are conserved. Also, energy drifts caused by numerical instability are better revealed by evaluating modified Hamiltonians. Doing this calculation would normally be complicated and highly dependent on the details of the method, even if differences are used to approximate derivatives. A relatively simple procedure is presented here, nearly independent of the internal structure of the integrator, for obtaining highly accurate estimates for modified Hamiltonians. As a bonus of the method of construction, the modified Hamiltonians are exactly conserved by a numerical solution in the case of a quadratic Hamiltonian.