My Research in Applied Mathematics:

Scientific Computing & Numerical Analysis

I am particulary interested in numerical algorithms that preserve and/or exploit special structure which mathematical problems inherit from the natural systems they model.

Bibliography

Background

The fields of structural biology and computational chemistry provide many mathematical problems with a common feature---Multiple Scales:
The most interesting things occur on large scales in space or time, but numerical methods for these problems can get stuck treating phenomena at scales which are perhaps many orders of magnitude smaller!

A problem with multiple scales in time is the computer simulation of molecular dynamics. For a large biomolecule, the motion of interest can occur over time intervals as long as one second. However, due to the very fast vibration of atoms in the molecule, the basic timestep for computer simulation is the femtosecond---10^{-15} seconds. The cost of each timestep is dominated by evaluation of the interatomic forces.
Practical numerical methods for such a problem should possess two properties:

1. Frugal consumption of force evaluations
We would like to perform the costly force evaluations as infrequently as possible, but two difficulties present themselves:
  • events such as collisions might occasionally require frequent force evaluation for correct numerical resolution
  • the fastest forces might always require frequent force evalutions.
2. Good long-time behavior
Most systems of interest are modelled by equations of motion for which analytic solutions are impossible to obtain, but we often know that solutions should possess certain properties (conserved quantities, geometric structure, time reversibility, thermodynamic properties,...). We would like to formulate, understand and apply numerical schemes which mimic these attributes.

My Recent Work

Long-time Integration Methods in Molecular Dynamics Simulation
Here in Tamar Schlick's group we develop novel approaches for the numerical treatment of the stochastic differential equations (Langevin Dynamics).

  • Implicit time discretizations, which can be numerically stable at very large timesteps, are considered in our paper on the LIN method [3] . In LIN, the implicit midpoint method is combined with linearized dynamics. Timesteps can be taken 30-60 times larger than traditional integration schemes, with resulting trajectories which are faithful to the modelled system. In practice, the extra work stemming from the implict method yields a computational speedup of less than double-despite the large timesteps.

  • The LN method, (replacement of the implicit component with a simpler treatment) is proposed [3] and further studied in [6] . LN timesteps are 6-10 times larger than traditional methods, with resulting computational speedup factor of 4-5. Our most recent LN work separates the systematic forces according to their spatial scale, and treats each with an appropriate (and efficient) timestep. Recent results [7] show that Langevin simulations of proteins with the LN method can be carried out in one-tenth (or less) the CPU time required by traditional methods; while reproducing dynamic, structural and energetic features of small-timestep trajectories. LN's ability to treat the slowest systematic components at long timesteps stems from the extrapolation treatment by which the slow forces are introduced into the numerical Langevin dynamics. We have studied this approach [8] in comparison with another widely-used slow force treatment: the impulse. Impulse methods are most suitable for moderately increasing the efficiency of constant-energy simulations, but stability considerations related to the fastest motion impose a severe limit on the size of timestep for the slow forces. Extrapolation methods do not suffer from the dramatic stability limitations, but fail to possess favorable energy-conservation properties. However, we show that extrapolation methods are quite suitable for efficient long-timestep simulations stochastic Langevin dynamics.

    An interesting review of current research in molecular dynamics algorithms is by Ron Elber, "Novel methods for molecular dynamics simulations",
    Current Opinion in Structural Biology, 6,232-235(1996).

    Constrained Dynamics
    One possible remedy is to constrain some of the fastest motion, enabling a longer timestep in the numerical simulation.

  • We studied structure preserving (symplectic) time discretizations for constrained molecular dynamics in [1] . In that work, we show that the structure of the molecule can be exploited to characterize the convergence of traditional methods for solving the nonlinear equations associated with the constraints and formulate new, more efficient, solution methods.

  • We have applied efficient constrained techniques from molecular dynamics to a particle model for conservative multibody dynamics problems in [2]. In such systems, constraints can be used to model rigid bodies and to restrict the motions of the bodies relative to one another.

  • The need to accurately resolve occasional strains and collisions in multibody systems while maintaining overall efficiency led us to seek adaptive variable timestep methods which are faithful to the long-time dynamics of conservative multibody systems. We presented a structure preserving time-reversible variable timestep time discretization for constrained dynamics in [5], along with a model for the dynamics of elastic rods.

    Structured Linear Algebra Problems

  • In constrained molecular and multibody dynamics, the structure of the modelled system can be used to determine the topology of the constraints. This topology in turn can be used to formulate and analyze efficient methods for solving the nonlinear equations of constraint. (See [1] and [2].)

  • In [4] we consider the implementation of eigenvalue algorithms for Hamiltonian matrices. These methods preserve special spectral properties which are destroyed by off-the-shelf eigenvalue packages.

  • In molecular modelling, Normal Mode Analysis is an important analytical tool. NM analysis amounts to finding eigenvalues and eigenvectors of the mass weighted Hessian (second spatial derivative of the potential energy). When the potential is restricted to include only local interactions, the resulting Hessian can be viewed as a collection of small matrices, coupled in a known way. I have written an efficient parallel code for this problem. The algorithm is described in a working paper you can get from me by mail.
    Back to My Home Page.

    Computational Mathematics and Chemistry Group's Home Page

    barth@kzoo.edu