Molecular simulation:
Algorithmic and Mathematical aspects

Institut Henri Poincaré

December 2004, 1-3

The main aim of this workshop is to mix chemists, physicits, biologists and mathematicans to exchange informations about molecular dynamic, long time integration and possible common problems. It follows the Prestissimo Workshop 2003 on the same topic.

This workshop will be held in the Amphithéâtre Darboux of the  Institut Henri Poincaré, in Paris.

If you wish to participate, please send an email to for information.

This workshop is organized by the INRIA Cooperative Research Action PRESTISSIMO and supported by the DFG-Schwerpunktprogramm 1095: Analysis, Modeling and Simulation of Multiscale Problems.

December 1:

Morning session (Chairman: Claude Le Bris)
09:00-10:00 Welcome
10:00-11:00 Andrew Stuart (Warwick)
Infinite Dimensional MCMC methods (Abstract,Slides)
11:30-12:30 Brian Laird (Kansas)
Direct Calculation of Crystal-Melt Interfacial Free Energies using Molecular-Dynamics Simulation (Abstract,Slides)

Afternoon session (Chairman: Olivier Pironneau)
14:00-15:00 Christof Schütte  (Berlin) 
Metastability in molecular dynamics: Mathematical Backgroud and Application to DNA (Abstract)
15:30-16:30 Mark Tuckerman  (NYU - Courant Institute)
Efficient sampling methods for molecular dynamics: Long time steps, variable transformations and adiabatic dynamics (Abstract,Slides)

December 2:

Morning session (Chairman: Gilles Zerah)
10:00-11:00 Arthur Voter (Los Alamos)
Accelerated Molecular Dynamics Methods (Abstract)
11:30-12:30 John Straub (Boston)
Atomistic and Reduced Models of Protein Folding and Aggregation (Abstract,Slides)

Afternoon session (Chairman: Christian Lubich)
14:00-14:30 Chris Sweet (Leicester)
Hamiltonian multiple thermostat techniques for thermostatting (Abstract)
14:30-15:00 Caroline Lasser (München)
Surface hopping for conical intersections (Abstract)
15:00-15:30 Erwan Faou (INRIA Rennes)
A Poisson integrator for Gaussian wavepacket dynamics (Abstract,Slides)
16:00-17:00 Martin Kröger (ETH Zürich)
Mesoscale simulation of complex fluids (Abstract)

December 3:

Morning session (Chairman: Andrew Stuart)
10:00-11:00 Tony Shardlow (Manchester)
Dissipative Paricle Dynamics (Abstract)
11:30-12:30 Robert Skeel (Purdue)
What Makes Molecular Dynamics Work? (Abstract,Slides)

Afternoon session (Chairman: Ben Leimkuhler)
14:00-14:30 Frédéric Legoll (IMA, University of Minnesota)
High-order averaging schemes for molecular dynamics simulations (Abstract,Slides)
15:00-16:00 Aiichiro Nakano (USC, Los Angeles)
Large Space-Time Multiscale Simulations of Nanosystems (Abstract,Slides)

Scientific Committee

Poster Session

During the conference, a poster session will be organized. For more information please contact


Erwan Faou and Christian Lubich
A Poisson integrator for Gaussian wavepacket dynamic
We consider the variational approximation of the time-dependent Schrödinger equation by Gaussian wavepackets.  The corresponding finite-dimensional dynamical system inherits a Poisson (or non-canonically symplectic) structure from the Schrödinger equation by its construction via  the Dirac-Frenkel-McLachlan variational principle. The variational splitting  between kinetic and potential energy turns out to yield an explicit,  easily implemented numerical scheme. This method is a time-reversible Poisson integrator, which also preserves the L² norm and linear and angular momentum. Using backward error analysis, we show long-time energy conservation for this splitting scheme. In the semi-classical limit the numerical approximations to position and momentum converge to those obtained by applying the Störmer-Verlet method to the classical limit system.
Martin Kröger
Mesoscale simulation of complex fluids
In this talk we present simple, mesoscale models and results obtained via nonequilibrium molecular and brownian dynamics for fluids such as polymers, ferrofluids, micellar systems, liquid crystals. We then give some details about a number of algorithmic and methodical challenging aspects related with these simulations:
i) Construction of the entanglement network,
ii) Chebyshev-polynomials for the calculation of hydrodynamic interactions,
iii) Multiplostats and Beyond-equilibrium molecular dynamics.
Brian Laird and Ruslan L. Davidchack
Direct Calculation of Crystal-Melt Interfacial Free Energies using Molecular-Dynamics Simulation 
The crystal-melt interfacial free energy, the work required to create a unit area of interface between a crystal and its own melt, is  a controlling property in the kinetics and morphology of crystal growth and nucleation, especially in the case of dendritic growth. Despite the technological importance of this quantity, accurate experimental data is difficult to obtain. The paucity of experimental measurements has motivated the development of a variety of novel computational methods to determine the interfacial free energy via molecular simulation. After a short tutorial on thermodynamic integration techniques for free energy calculation, I will introduce our method of cleaving walls for the calculation of the crystal-melt interfacial free energy, and a competing method based on fluctuation spectra. Results for a variety of simple systems will be presented to give a broad picture of the interaction and crystal structure dependence of the interfacial free energy. The results will be discussed in relation to popular empirical theories of the interfacial free energy.
Caroline Lasser
Surface hopping for conical intersections
The nucleonic evolution near electron energy level crossings plays an important  role in laser-controlled femtosecond chemistry. The notoriously high number of degrees of freedom (three times the number of nuclei) suggests computational schemes, which are based on the adiabatic decoupling provided by the time-dependent Born-Oppenheimer approximation. Such approaches, however, break down in the presence of conical level crossings. We present a novel algorithm, which implements a Markov process combining classical transport and non-adiabatic transitions. Opposed to the well-established trajectory surface hopping algorithms of chemical physics, our rigorously derived approach is based on precise criteria for hopping-time and position as well as an explicit Landau-Zener type transition rate for the non-adiabatic transfer.
Frédéric Legoll
High-order averaging schemes for molecular dynamics simulations
Many thermodynamical properties of chemical systems are defined as phase space averages of observables. A way to compute these averages is to use Molecular Dynamics and to compute time averages on long trajectories.
In this work, we propose high-order averaging formulae to compute time averages, in order to improve the convergence rate with respect to the simulation time. Many numerical examples will be provided, as well as error bounds.
We will also discuss the application of these averaging formulae for the numerical integration of highly oscillatory problems.
This is joint work with Eric Cancès, Claude Le Bris and Gabriel Turinici (CERMICS and INRIA Rocquencourt), and with François Castella, Philippe Chartier and Erwan Faou (INRIA Rennes).
Aiichiro Nakano, Rajiv K. Kalia, and Priya Vashishta
Large Space-Time Multiscale Simulations of Nanosystems
We are developing a hierarchical simulation framework, in which accurate but compute-intensive calculations are invoked only where and when they are needed. Our multiscale simulation embeds atomistic simulation based on reactive molecular dynamics (MD) methods and quantum mechanical (QM) calculation based on the density functional theory (DFT) in a continuum model.
Our scalable computational framework combines data locality principles and divide-and-conquer strategies: Linear-scaling simulation algorithms based on space-time multiresolution techniques; topology-preserving computational-space decomposition with wavelet-based adaptive load balancing; spacefilling-curve-based data compression for compressed software pipelines; and interactive and immersive visualization of massive simulation data. Large-scale atomistic simulations have been applied to various nanotechnology areas. These include the design of tough nanophase ceramics and nanocomposites, in silico mechanical testing based on nanoindentation, hypervelocity impact, nano-energetic materials, colloidal and epitaxial quantum dots (QDs), and surface switching of self-assembled monolayers.
Christof Schütte
Metastability in molecular dynamics: Mathematical Backgroud and Application to DNA
The talk will be concerned with the reliable analysis of the main conformations of a biomolecule under consideration. In a conformation, the large scale geometric
structure of the molecule is understood to be conserved, whereas on smaller scales the system may well rotate, oscillate or fluctuate. Furthermore, transitions between main conformations are rare events or, in other words, a typical trajectory of a molecular system stays for long periods of time within the conformation, while exits are long-term events. Hence, the term conformation includes both geometric and dynamical aspects.  From the geometrical point of view, conformations are understood to represent all molecules with the same large scale geometric structure and may thus be identified with a subset of the state space. From the dynamical point of view, a conformation typically persists for long periods of time (compared to the fastest molecular motions) such that the associated subset of the state space is metastable and the resulting macroscopic dynamical behavior can be described as a flipping process between the
metastable subsets. Understanding of conformation dynamics---that is the statistics of the flipping process and its mechanisms---is crucial to the understanding of biomolecular flexibility and activity.
We present a novel method for the identification of the most important conformations of a biomolecular system from molecular dynamics time series by means of Hidden Markov Models (HMMs) with output behavior given by stochastic differential equations.
The performance of the resulting method is illustrated by numerical tests and by application to molecular dynamics time series of a DNA segment.
Tony Shardlow
Dissipative Paricle Dynamics
I will introduce Dissipative Particle Dynamics (DPD), a system of stochastic differential equation (SDEs) used to model fluid flows. I shall discuss the numerical analysis of SDEs, focusing on DPD, and shall prove that DPD is geometrically ergodic in one spatial dimension and explain the difficulties in extending this proof to two and three dimensions.
Robert Skeel
What Makes Molecular Dynamics Work?
The numerical treatment of molecular dynamics (MD) is problematic due to the exponential growth of error with time. The technique of shadowing can show that the computed solution is very nearly the exact trajectory for a slightly different initial value. This type of result would be adequate, and it has been shown to work for fairly complicated problems. However, it is highly unlikely that a (useful) shadowing result is possible for highly elliptic dynamical systems such as MD. Rather, it is suggested that a general treatment for MD be based on randomness in the initial values and applying the concept of weak convergence from stochastic differential equations. Weak convergence requires that expectations be accurately computed for smooth distributions of initial values. In this setting it is plausible that accurate solutions can be obtained for very long intervals of time. There remain questions concerning the accuracy of different numerical integrators in this weak sense, and these questions are explored. In the case of ergodic Hamiltonian systems, evidence is presented suggesting that consistent integrators give convergent results on very long intervals if the integrator nearly conserves energy on very long intervals and conserves volume in phase space. However, for certain practical reasons it seems that the stronger property of being symplectic is needed, and this is explained.
John Straub
Atomistic and Reduced Models of Protein Folding and Aggregation
The thermodynamics and dynamics of atomic-level and coarse-grained models of protein folding and aggregation are explored. Two coarse-grained protein models are developed and examined. The interaction potentials depend on the radial distance between interaction sites, as well as the relative orientation of the sites. In one case, the potentials are fit to statistical data on protein structural preferences; in the other the potentials are more empirical in nature. Varying chain lengths are studied and the thermodynamics of the coil-to-helix transition are analyzed in terms of the Zimm-Bragg model. Results agree qualitatively and quantitatively with all-atom Monte Carlo simulations and other reduced-model Monte Carlo simulations. The kinetics of the coil-to-helix transition is examined in order to explore the possible importance of multiple-nucleation sites for helix formation that may be important to the mechanism of helix formation, even for relatively short chain lengths. Comparison is made with results derived from variational calculations of diffusional transition pathways and analytic models. Prospects for the use of coarse-grained potentials in exploring protein aggregation are discussed.
Andrew Stuart
Infinite Dimensional MCMC methods
There are a wide variety of sampling problems which are infinite dimensional in character. Examples include transtion path sampling in chemistry and nonlinear filtering in signal processing. In both these examples the basic object to sample is a path in time, and is hence infinite dimensional. We describe an abstract MCMC method for sampling such problems, based  on generalizing the Langevin sampler to infinite dimensions. This leads naturally to stochastic PDEs which, in their invariant measure, sample from the required distribution.
Chris Sweet
Hamiltonian multiple thermostat techniques for thermostatting
Molecular dynamics trajectories that sample from the Gibbs distribution can be generated by introducing a modified Hamiltonian with additional degrees of freedom as described by Nosé. To achieve the ergodicity required for canonical sampling, a number of techniques have been proposed based on incorporating additional thermostatting variables, such as Nosé-Hoover chains. For Nosé dynamics, it is often stated that the system is driven to equilibrium through a resonant interaction between the self-oscillation frequency of the thermostat variable and a natural frequency of the underlying system. By the introduction of multiple thermostat Hamiltonian formulations, which are not restricted to chains, it has been possible to clarify this perspective, using harmonic models, and exhibit practical deficiencies of the standard Nosé-chain approach. As a consequence of this it has been possible to propose both a Hamiltonian chains method, Nosé-Poincaré chains,  and a new "recursive thermostatting" procedure which both obtain canonical sampling. These methods also have the advantage that the choice of Nosé mass is essentially independent of the system to be thermostatted.
Mark Tuckerman
Efficient sampling methods for molecular dynamics: Long time steps, variable transformations and adiabatic dynamics
One of the computational grand challenge problems is the development of methodology capable of sampling conformational equilibria in systems with rough energy landscapes. If met, many important problems, most notably protein folding, could be significantly impacted. In this talk, I will discuss several new approaches for enhancing conformational sampling in molecular dynamics calculations. First, a new method will be presented for overcoming the resonance problem that plagues multiple time scale molecular dynamics methods and imposes limits on the size of the outer time step. It will be shown that the new method allows time steps as large as 100 fs to be achieved for difficult problems such as flexible models of water and biological macromolecules. Next, I will discuss a new approach that combines molecular dynamics with a novel variable transformation designed to warp configuration space in such a way that barriers are reduced and attractive basins stretched. The new method rigorously preserves equilibrium properties while leading to very large enhancements in sampling efficiency. The method will be demonstrated on the Honeycutt-Thirumalai model and compared to parallel tempering. Finally, it will be shown that, when variable transformations are combined with an adiabatic separation between a reactive subspace of coordinates and the remaining degrees of freedom, an efficient method for exploring free energy surfaces emerges. Here, studies of small peptides both in the gas-phase and in solution will be used to demonstrate the method and comparisons will be made to multi-dimensional umbrella sampling.
Arthur Voter
Accelerated Molecular Dynamics Methods
A significant problem in the atomistic simulation of materials is that molecular dynamics simulations are limited to nanoseconds, while important reactions and diffusive events often occur on time scales of microseconds and longer.  Although rate constants for slow events can be computed directly using transition state theory (with dynamical corrections, if desired, to give exact rates), this requires first knowing the transition state. Often, however, we cannot even guess what events will occur.  For example, in vapor-deposited metallic surface growth, surprisingly complicated exchange events are pervasive.  I will discuss recently developed methods (hyperdynamics, parallel replica dynamics, and temperature accelerated dynamics) for treating this problem of complex, infrequent-event processes.  The idea is to directly accelerate the dynamics to achieve longer times without prior knowledge of the available reaction paths.  I will present our latest method developments and some recent applications, including metallic surface growth, deformation and dynamics of carbon nanotubes, and annealing after radiation damage events in MgO.