Library

Notes: We report results from molecular dynamics simulations for a bistable piecewise-harmonic potential. A new method for molecular dynamics—the Langevin/implicit-Euler scheme—is investigated here and compared to the common Verlet integration algorithm. The implicit scheme introduces new computational and physical features since it (1) does not restrict integration time step to a very small value, and (2) effectively damps vibrational modes ω(very-much-greater-than)ωc, where ωc is a chosen cutoff frequency. The main issue we explore in this study is how different choices of time steps and cutoff frequencies affect computed transition rates. The one-dimensional, double-well model offers a simple visual and computational opportunity for observing the two different damping forces introduced by the scheme—frictional and intrinsic—and for characterizing the dominating force at a given parameter combination. Another question we examine here is the choice of time step below which the Langevin/implicit-Euler scheme produces "correct'' transition rates for a model potential whose energy distribution is "well-described'' classically.

Notes: Results are presented from potential energy minimization of water clusters and from molecular dynamics and Monte Carlo simulations of a liquid water droplet model. A new method for molecular dynamics—the implicit-Euler/Langevin scheme—is used in combination with a truncated Newton minimizer for potential energy functions. Structural and thermodynamic properties are reported for the scheme (with time steps of 5 and 10 fs), compared to a standard explicit formulation (with Δt=1 fs), to a Monte Carlo simulation, and to available experimental data. Results demonstrate that the implicit scheme is computationally feasible for large-scale biomolecular simulations, and that the droplet model can reasonably reproduce general structural features of liquid water. Results also show that the desired behavior is obtained from the implicit formulation: stability over large time steps, and effective damping of the high-frequency vibrational modes. Thus, major "bulk'' properties of the system of interest may be observed more rapidly.

Notes: Abstract Innovative algorithms have been developed during the past decade for simulating Newtonian physics for macromolecules. A major goal is alleviation of the severe requirement that the integration timestep be small enough to resolve the fastest components of the motion and thus guarantee numerical stability. This timestep problem is challenging if strictly faster methods with the same all-atom resolution at small timesteps are sought. Mathematical techniques that have worked well in other multiple-timescale contexts-where the fast motions are rapidly decaying or largely decoupled from others-have not been as successful for biomolecules, where vibrational coupling is strong. This review examines general issues that limit the timestep and describes available methods (constrained, reduced-variable, implicit, symplecttic, multiple-timestep, and normal-mode-based schemes). A section compares results of selected integrators for a model dipeptide, assessing physical and numerical performance. Included is our dual timestep method LN, which relies on an approximate linearization of the equations of motion every Deltat interval (5 fs or less), the solution of which is obtained by explicit integration at the inner timestep Deltatau (e.g., 0.5 fs). LN is computationally competitive, providing 4-5 speedup factors, and results are in good agreement, in comparison to 0.5 fs trajectories. These collective algorithmic efforts help fill the gap between the time range that can be simulated and the timespans of major biological interest (milliseconds and longer). Still, only a hierarchy of models and methods, along with experimentational improvements, will ultimately give theoretical modeling the status of partner with experiment.

Notes: As molecular dynamics simulations continue to provide important insights into biomolecular structure and function, a growing demand for increasing the time span of the simulations is emerging. Our focus here is developing a new algorithm, LIN (Langevin/implicit-Euler/normal mode), that combines normal-mode and implicit-integration techniques, for large time step biomolecular applications. In the normal-mode phase of LIN, we solve an approximate linearized Langevin formulation to resolve the rapidly varying components of the motion. In the implicit phase, we resolve the remaining components of the motion by numerical integration with the implicit-Euler scheme. Developments of the normal-mode phase of LIN are discussed in this paper. Specifically, we solve two crucial issues of the method. The first involves how to choose and how often to update the Hessian approximation for the linearized Langevin equation. This approximation must be computationally feasible and physically reasonable to capture the motion in the higher end of the vibrational spectrum. Three such general Hessian approximations are discussed. The related issue—the frequency of the Hessian update—is analyzed by projecting the motion onto the different vibrational modes. This analysis demonstrates that a one-picosecond interval is reasonable for updating the Hessian in the model system examined here.In this connection, we illustrate that the high-frequency motions are highly localized while the low-frequency motions are delocalized. We also show rigorously that the mode amplitudes are inversely proportional to the frequency (consistent with the equipartition theorem), with 90% of the displacement fluctuations coming from a very small group of low-frequency modes. Anharmonic effects essentially influence the low-frequency modes. The second issue involves how to solve the linearized Langevin equation at large timesteps correctly, where the usual discretized formulation of the random force is invalid. This is accomplished by using analytic expressions for the distributions associated with positions and velocities of the individual oscillators as a function of frequency, obtained as the solution of the corresponding Fokker–Planck equation. We apply LIN with these developments to the nucleic acid component deoxycytidine with timesteps ranging from 100 to 1000 fs. We demonstrate that LIN is stable in these simulations, with energies fluctuating about the same values—and possessing overall similar dynamical features—in comparison to 1 fs explicit simulations, though the fluctuations are significantly larger at larger timesteps. Moreover, continuous dynamics is maintained, and pathway information can be obtained. Computational performance is competitive only at very large time steps: a gain factor of 3–4 is obtained for runs with 1000 fs time steps. Larger gains may be achieved for biomolecules, where sparsity and parallelization can be exploited significantly.

Notes: As a balanced response to the Wu and Watts critique of two novel methods for molecular dynamics, the record is set straight with respect to LIs damping, the mixing of two different aspects of LI, the effects of quantum mechanics, and the performance of LIN. © 1995 American Institute of Physics.

Notes: The notion of error in practical molecular and Langevin dynamics simulations of large biomolecules is far from understood because of the relatively large value of the timestep used, the short simulation length, and the low-order methods employed. We begin to examine this issue with respect to equilibrium and dynamic time-correlation functions by analyzing the behavior of selected implicit and explicit finite-difference algorithms for the Langevin equation. We derive: local stability criteria for these integrators; analytical expressions for the averages of the potential, kinetic, and total energy; and various limiting cases (e.g., timestep and damping constant approaching zero), for a system of coupled harmonic oscillators. These results are then compared to the corresponding exact solutions for the continuous problem, and their implications to molecular dynamics simulations are discussed. New concepts of practical and theoretical importance are introduced: scheme-dependent perturbative damping and perturbative frequency functions. Interesting differences in the asymptotic behavior among the algorithms become apparent through this analysis, and two symplectic algorithms, "LIM2'' (implicit) and "BBK'' (explicit), appear most promising on theoretical grounds. One result of theoretical interest is that for the Langevin/implicit-Euler algorithm ("LI'') there exist timesteps for which there is neither numerical damping nor shift in frequency for a harmonic oscillator. However, this idea is not practical for more complex systems because these special timesteps can account only for one frequency of the system, and a large damping constant is required. We therefore devise a more practical, delay-function approach to remove the artificial damping and frequency perturbation from LI. Indeed, a simple MD implementation for a system of coupled harmonic oscillators demonstrates very satisfactory results in comparison with the velocity-Verlet scheme. We also define a probability measure to estimate individual trajectory error. This framework might be useful in practice for estimating rare events, such as barrier crossing. To illustrate, this concept is applied to a transition-rate calculation, and transmission coefficients for the five schemes are derived. © 1996 American Institute of Physics.

Notes: Two algorithms are presented for integrating the Langevin dynamics equation with long numerical time steps while treating the mass terms as finite. The development of these methods is motivated by the need for accurate methods for simulating slow processes in polymer systems such as two-site intermolecular distances in supercoiled DNA, which evolve over the time scale of milliseconds. Our new approaches refine the common Brownian dynamics (BD) scheme, which approximates the Langevin equation in the highly damped diffusive limit. Our LTID ("long-time-step inertial dynamics") method is based on an eigenmode decomposition of the friction tensor. The less costly integrator IBD ("inertial Brownian dynamics") modifies the usual BD algorithm by the addition of a mass-dependent correction term. To validate the methods, we evaluate the accuracy of LTID and IBD and compare their behavior to that of BD for the simple example of a harmonic oscillator. We find that the LTID method produces the expected correlation structure for Langevin dynamics regardless of the level of damping. In fact, LTID is the only consistent method among the three, with error vanishing as the time step approaches zero. In contrast, BD is accurate only for highly overdamped systems. For cases of moderate overdamping, and for the appropriate choice of time step, IBD is significantly more accurate than BD. IBD is also less computationally expensive than LTID (though both are the same order of complexity as BD), and thus can be applied to simulate systems of size and time scale ranges previously accessible to only the usual BD approach. Such simulations are discussed in our companion paper, for long DNA molecules modeled as wormlike chains. © 2000 American Institute of Physics.

Notes: We apply our new algorithms presented in the companion paper (LTID: long-time-step inertial dynamics, IBD: inertial Brownian dynamics) for mass-dependent Langevin dynamics (LD) with hydrodynamics, as well as the standard Brownian dynamical (BD) propagator, to study the thermal fluctuations of supercoiled DNA minicircles. Our DNA model accounts for twisting, bending, and salt-screened electrostatic interactions. Though inertial relaxation times are on the order of picoseconds, much slower kinetic processes are affected by the Brownian (noninertial) approximation typically employed. By comparing results of LTID and IBD to those generated using the standard (BD) algorithm, we find that the equilibrium fluctuations in writhing number, Wr, and radius of gyration, Rg, are influenced by mass-dependent terms. The autocorrelation functions for these quantities differ between the BD simulations and the inertial LD simulations by as much as 10%. In contrast, when the nonequilibrium process of relaxation from a perturbed state is examined, all methods (inertial and diffusive) yield similar results with no detectable statistical differences between the mean folding pathways. Thus, while the evolution of an ensemble toward equilibrium is equally well described by the inertial and the noninertial methods, thermal fluctuations are influenced by inertia. Examination of such equilibrium fluctuations in a biologically relevant macroscopic property—namely the two-site intermolecular distance—reveals mass-dependent behavior: The rate of juxtaposition of linearly distant sites along a 1500-base pair DNA plasmid, occurring over time scales of milliseconds and longer, is increased by about 8% when results from IBD are compared to those from BD. Since inertial modes that decay on the picosecond time scale in the absence of thermal forces exert an influence on slower fluctuations in macroscopic properties, we advocate that IBD be used for generating long-time trajectories of supercoiled DNA systems. IBD is a practical alternative since it requires modest computational overhead with respect to the standard BD method. © 2000 American Institute of Physics.

Notes: We present an efficient new method termed LN for propagating biomolecular dynamics according to the Langevin equation that arose fortuitously upon analysis of the range of harmonic validity of our normal-mode scheme LIN. LN combines force linearization with force splitting techniques and disposes of LIN's computationally intensive minimization (anharmonic correction) component. Unlike the competitive multiple-timestepping (MTS) schemes today—formulated to be symplectic and time-reversible—LN merges the slow and fast forces via extrapolation rather than "impulses;" the Langevin heat bath prevents systematic energy drifts. This combination succeeds in achieving more significant speedups than these MTS methods which are limited by resonance artifacts to an outer timestep less than some integer multiple of half the period of the fastest motion (around 4–5 fs for biomolecules). We show that LN achieves very good agreement with small-timestep solutions of the Langevin equation in terms of thermodynamics (energy means and variances), geometry, and dynamics (spectral densities) for two proteins in vacuum and a large water system. Significantly, the frequency of updating the slow forces extends to 48 fs or more, resulting in speedup factors exceeding 10. The implementation of LN in any program that employs force-splitting computations is straightforward, with only partial second-derivative information required, as well as sparse Hessian/vector multiplication routines. The linearization part of LN could even be replaced by direct evaluation of the fast components. The application of LN to biomolecular dynamics is well suited for configurational sampling, thermodynamic, and structural questions. © 1998 American Institute of Physics.