Mathematical Challenges from Theoretical/Computational Chemistry

The multiple minima problem, discussed elsewhere in this report, can be addressed by approaches other than optimization. Specifically, reasonable sampling strategies can be devised to scan configuration space in the hope of obtaining information about many local minima, maxima, and saddlepoints. Moreover, these strategies can incorporate a variety of chemical information, such as interproton distances from NMR, van der Waals radii (for hard-core exclusion), and secondary-structure elements. Useful sea rch strategies today involve Monte Carlo, molecular dynamics (MD), Brownian and Langevin dynamics, calculations of free-energy perturbations, high-temperature simulations, normal-mode analyses, and various enhanced sampling techniques that involve hybrids of all of the above. As reflected by those many approaches, the problem of inadequate sampling of configurational space is receiving increasing attention as a realization emerges that faster and more powerful computers alone cannot solve this problem in the near future. New methodologies and a hierarchy of approaches at different levels of resolution--in combination with experiment--are needed to attack this sampling problem to advance current capabilities of computational chemistry in connection with biomolecules.

A specific problem involves numerical integration in the context of MD simulations. In this technique, molecular motion is propagated by numerically integrating the classical equations of motion governing a molecular system under the influence of a spe cified force field (McCammon and Harvey, 1987; Allen and Tildesley, 1990). In theory, MD simulations can provide extensive spatial and temporal information. However, inadequate sampling limits the scope of the results that can be obtained in practice. Similar issues arise in other chemical applications, such as quantum mechanics, and the development of improved integration schemes will advance the systems and types of processes that can be simulated on modern computers.

Numerical Methods for Solving Ordinary Differential Equations

Many problems in chemistry can be reduced to the solution of systems of coupled ordinary differential equations (ODEs). Examples include classical and Langevin dynamics, rate equations of kinetic theory, and the time-dependent Schrödinger equation when expanded in a basis set. Thus, numerical integrators used to solve these equations are fundamental tools in computational/theoretical chemistry, and any significant improvement in these integrators (e.g., speedup, long-time stability) results in ad vances throughout the field.

The technology of numerical integrators for solving ODEs has a long history with significant interplay among mathematics, physics, and chemistry. Many of the earliest integrators, such as Runge-Kutta and predictor-corrector integrators, are still in co mmon use, but there have also been recent advances, driven in part by the need for methods that can treat multiple time scales and have greater stability for large-scale coupled nonlinear oscillators commonly found in MD of polymers and biological macromo lecules. The long-time stability of integrators for such systems is a challenging area of mathematical analysis research; perhaps the chemical applications described here will stimulate important developments.

Symplectic Integrators

Symplectic integrators have recently gained attention in the mathematical community and were quickly adapted for use in dynamics calculations in chemistry because of their favorable properties. In applications to Hamiltonian systems, symplectic integra tors have the property of building in Liouville's theorem, whereby areas in phase space are preserved as the system evolves in time. This strong conservation property translates into stability over long-time integrations, an important property in MD calc ulations involving millions and more steps. One consequence of this for constant-energy MD simulations is that except for fluctuations, symplectic integrators with small time steps conserve energy for very long times, whereas nonsymplectic integrators ty pically introduce a systematic drift in the total energy. Time reversibility is another useful practical property of symplectic integrators.

Symplecticness may also be used in determining numerical solutions to the Schrödinger equation. There is an equivalent representation of quantum mechanics in terms of Hamilton's equations (Gray and Verosky, 1994) that makes possible the use of int egrators for the quantum dynamics studies that are used for classical dynamics. Area-preserving mapping are also of interest in their own right in studies of dynamical systems (Meiss, 1992).

A symplectic integrator is any of a class of numerical algorithms for the integration of classical many-body equations of motion that possess favorable properties such as area preservation, energy conservation, and time reversibility.

Symplectic integrators may be implicit or explicit. In explicit methods, the solution at the end of the time step is obtained by performing operations on the variables at the beginning of each time step. Symbolically, we write y = f(y, t, ...), where f is some nonlinear f unction, t is the time step, y is the approximation to the solution y at time nt, and the dots indicate other parameters or pre vious solutions (e.g., y, y). With implicit integrators, the final solutions are functions of both the initial and the final variables (y = f(y, y, t, ...)), and so coupled nonlinear equations must generally be solved at each time step to propagate the trajectory. The explicit versions generally involve simple algorithms that (for propagation only) use modest memory, while implicit methods involve more comple x algorithms but are often more powerful for treating systems with disparate time scale dynamics, as discussed below.

The development of symplectic integrators has involved significant interplay among mathematicians, physicists, and chemists. Seminal work on symplectic integrators was done by both physicists and mathematicians (Ruth, 1983; Feng, 1986; Candy and Rozmus , 1991; McLachlan and Atela, 1992; Okunbor and Skeel, 1992; Calvo and Sanz-Serna, 1993) based on second- and third-order explicit approaches and Runge-Kutta methods. Implicit approaches were developed in parallel (Channell and Scovel, 1990; De Frutos and Sanz-Serna, 1992). Recently, these ideas have found their way into the chemistry community (Gray et al., 1994) with promising r esults. The Verlet integrator (Verlet, 1967), already in common use, was found to be symplectic, thereby explaining the good associated stability observed in practice. However, symplectic integrators that improve on previously available methods have als o been developed (Gray et al., 1994). Initial applications using these methods suggest that they may become favored for simulations of polymer dynamics and related problems with small time steps.

The Time Step Problem in Molecular Dynamics

Although standard explicit schemes, such as the Verlet and related methods, are simple to formulate and fast to propagate, they impose a severe constraint on the maximum time step possible. Instability--uncontrollable growth of coordinates and velociti es--occurs for step sizes much larger than 1 femtosecond (10 second). This step size is determined by the period associated w ith high-frequency modes present in all macromolecules, and it contrasts with the much longer time scales (up to 10 seconds) that govern key conformational changes (e.g., folding) in macromolecules. This dispari ty in time scales urges the development of methods that increase the time step for biomolecular simulations. Even if the stability of the numerical formulation can be ensured, an important issue concerning the reliability of the results arises, since vib rational modes in molecular systems are intricately coupled.

Standard techniques of effectively freezing the fast vibrational modes by a constrained formulation (Ryckaert et al., 1977; van Gunsteren and Berendsen, 1977; van Gunsteren, 1980; Miyamoto and Kollman, 1992) increase the time step by a small factor such as two, still with added complexity at each step. The multiple time step approaches for updating the slow and fast forces provide additional speedup (Streett et al., 1978; Grubmüller et al., 1991; Tuckerman and Berne, 1992; Watanabe and Karplus, 19 93), although some stability issues are also involved (Biesiadecki and Skeel, 1993).

Implicit Integration Schemes

There are well-known numerical techniques for solving differential equations describing physical processes with multiple time scales (Gear, 1971; Dahlquist and Björck, 1974). Various implicit formulations are available that balance stability, accu racy, and complexity. However, the standard implicit techniques used by numerical analysts (Kreiss, 1991) have not been directly applicable to MD simulations of macromolecules, for the following reasons.

First, such implicit schemes are often designed to suppress the rapidly decaying component of the motion. This is a valid approach when the contribution of these components becomes negligible for sufficiently long times, as is the case for the second t erm in y(t) = exp (- t) + exp (- 100 t). However, this situation does not hold for biomolecular systems because of the intricate vibrational coupling. It is well recognized that concerted conformational transitions (e.g., in hinge-bending proteins) requ ire a cooperative mechanism driven by small- scale fluctuations to make possible a large-scale collective displacement. Thus, although the absence of the positional fluctuations associated with these high-frequency modes may not by itself be a severe problem, the absence of the energies associated with these modes may be undesirable for proteins and nucleic acids, si nce cooperative motions among the correlated vibrational modes may rely on energy transfer from these high-frequency modes.

Second, implicit schemes with known high stability (e.g., implicit Euler) can introduce numerical damping (Zhang and Schlick, 1993). This has prompted the application of such implicit schemes to the Langevin dynamics formulation, which involves frictio nal and Gaussian random forces in addition to the systematic force to mimic molecular collisions and therefore a thermal reservoir. This stabilizes implicit discretizations and can be used to quench high-frequency vibrational modes (Peskin and Schlick, 1 989; Schlick and Peskin, 1989), but unphysical increased rigidity can result (Zhang and Schlick, 1993). Therefore, more rigorous approaches are required to resolve the subdynamics correctly, such as by combining normal-mode techniques with implicit integrat ion (Zhang and Schlick, 1994); significant linear algebra work in the spectral decomposition is necessary for feasibility for macromolecular systems. For example, banded structures for the Hessian approximation (see related discussion in the section on m ultivariate minimization beginning on page 68) can be exploited in the linearized equations of motion. There has also been some work on implicit schemes that do not have inherent damping, but preliminary experience suggests that for nonlinear systems, de sirable energy conservation properties can be obtained only up to moderate time steps (Simo et al., 1992; Zhang and Schlick, 1995). In particular, serious resonance problems have been noted (Mandziuk and Schlick, 1995).

Third, implicit schemes for multiple time scale problems increase complexity, since they involve solution of a nonlinear system or minimization of a nonlinear function at each time step. Therefore, very efficient implementations of these additional com putations are necessary, and even then, computational gain (with respect to standard "brute-force" integrations at small time ste ps) can be realized only at very large time steps.

Future Prospects

The preceding subsections have described several recent accomplishments in the development of integration methods in MD simulations and have also outlined several important challenges for the future. What makes these integration problems particularly c hallenging is the fact that solutions demand much more than straightforward application of standard mathematical techniques. At this point it appears that the optimal algorithms for MD will require a combination of methods and strategies discussed above, including symplectic and implicit numerical integration schemes that have minimal intrinsic damping, and correct resolution of the subdynamics of the system by some other technique (e.g., normal-mode analysis). Undoubtedly, high-performance implementati ons will make possible a gain of several orders of magnitude in the simulation times, and there are certainly additional gains to be achieved by clever programming strategies.

References

Allen, M.P., and D.J. Tildesley, 1990, Computer Simulation of Liquids, Oxford University Press, New York.

Biesiadecki, J.J., and R.D. Skeel, 1993, Dangers of multiple-time-step methods, J. Comput. Phys. 109:318-328.

Calvo, M.P., and J.M. Sanz- Serna, 1993, The development of variable- step symplectic integrators, with application to the two- body p roblem, SIAM J. Sci. Comput. 14:936.

Candy, J., and W. Rozmus, 1991, A symplectic integration algorithm for separable Hamiltonian functions, J. Comput. Phys. 92:230.

Channell, P.J., and J.C. Scovel, 1990, A symplectic integration of Hamiltonian systems, Nonlinearity 3:231.

Dahlquist, G., and Å. Björck, 1974, Numerical Methods, Prentice Hall, Englewood Cliffs, N.J.

De Frutos, J., and J.M. Sanz- Serna, 1992, An easily implementable fourth- order method for the time integration of wave problems, J. Comput. Phys. 103:160.

Feng, K., 1986, Difference schemes for Hamiltonian formalism and symplectic geometry, J. Comput. Math. 4:279-289.

Gear, C.W., 1971, Numerical Initial Value Problems in Ordinary Differential Equations, Prentice Hall, Englewood Cliffs, N.J.

Gray, S.K., and J.M. Verosky, 1994, Classical Hamiltonian structures in wave packet dynamics, J. Chem. Phys. 100:5011-5022.

Gray, S.K., D.W. Noid, and B.G. Sumpter, 1994, Symplectic integrators for large scale molecular dynamics simulations: A comparison of several explicit methods, J. Chem. Phys. 101:4062-4072.

Grubmüller, H., H. Heller, A. Windemuth, and K. Schulten, 1991, Generalized verlet algorithm for efficient molecular dynamics simulations with long- range interactions, Mol. Simul. 6:121-142.

Kreiss, H.-O., 1991, Problems with different time scales, Acta Numerica 1:101-139.

Mandziuk, M., and T. Schlick, 1995, Resonance in chemical- system dynamics simulated by the implicit- midpoint scheme, Chem. Phys. Lett., in press.

McCammon, J.A., and S.C. Harvey, 1987, Dynamics of Proteins and Nucleic Acids, Cambridge University Press, Cambridge.

McLachlan R.I., and P. Atela, 1992, The accuracy of symplectic integrators, Nonlinearity 5:541.

Meiss, J.D., 1992, Symplectic maps, variational principles, and transport, Rev. Mod. Phys. 64:795.

Miyamoto, S., and P.A. Kollman, 1992, SETTLE: An analytical version of the SHAKE and RATTLE algorithm for rigid water models, J. Comput. Chem. 13:952-962.

Okunbor, D.I., and R.D. Skeel, 1992, Canonical numerical methods for molecular dynamics simulations, Math. Com. 59:439.

Peskin, C.S., and T. Schlick, 1989, Molecular dynamics by the backward- Euler method, Commun. Pure Appl. Math. 42:1001-1 031.

Ruth, R.D., 1983, A canonical integration technique, IEEE Trans. Nucl. Sci. NS-30:2669.

Ryckaert, J.P., G. Ciccotti, and H.J.C. Berendsen, 1977, Numerical integration of the cartesian equations of motion of a system with constraints: Molecular dynamics of n- alkanes, J. Comput. Phys. 23:327-341.

Schlick, T., and C.S. Peskin, 1989, Can classical equations simulate quantum- mechanical behavior? A molecular dynamics investiga tion of a diatomic molecule with a Morse potential, Commun. Pure Appl. Math. 42:1141-1163.

Simo, J.C., N. Tarnow, and K.K Wong, 1992, Exact energy- momentum conserving algorithms and symplectic schemes for nonlinear dyna mics, Computer Methods in Applied Mechanics and Engineering 100:63-116.

Streett, W.B., D.J. Tildesley, and G. Saville, 1978, Multiple time step methods in molecular dynamics, Mol. Phys. 35:639-648.

Stuart, A.M., and A.R. Humphries, 1994, Model problems in numerical stability theory for initial value problems, SIAM Review 36:226-257.

Tuckerman, M.E., and B.J. Berne, 1992, Molecular dynamics in systems with multiple time scales: systems with stiff and soft degrees of freedom and with short and long range forces, J. Comput. Chem. 95:8362-8364.

van Gunsteren, W.F., 1980, Constrained dynamics of flexible molecules, Mol. Phys. 40:1015-1019.

van Gunsteren, W.F., and H.J.C. Berendsen, 1977, Algorithms for macromolecular dynamics and constraint dynamics, Mol. Phys. 34:1311-1327.

Verlet, L., 1967, Computer "experiments" on classical fluids. I. Thermodynamical properties of Lennard- Jones molecules, Phys . Rev. 159:98.

Watanabe, M., and M. Karplus, 1993, Dynamics of molecules with internal degrees of freedom by multiple time- step methods, J. Chem. Phys. 99:8063-8074.

Zhang, G., and T. Schlick, 1993, LIN: A new algorithm to simulate the dynamics of biomolecules by combining implicit- integration and normal mode techniques, J. Comput. Chem. 14:1212-1233.

Zhang, G., and T. Schlick, 1994, The Langevin/implicit-Euler/normal-mode scheme (LIN) for molecular dynamics at large time steps, J. Chem. Phys. 101:4995-5012.

Zhang, G., and T. Schlick, 1995, Implicit discretization schemes for Langevin dynamics, Mol. Phys., in press.

NAS Home Page | NAP Home Page | Reading Room | Report Home Page