Molecular Dynamics
Monday 1 – Friday 5 June 2009

Organisers: Mike Allen (University of Warwick), Mark Rodger (University of Warwick), Ben Leimkuhler (University of Edinburgh), Giovanni Ciccotti (University of Rome)

Printed abstracts .

Peter Bolhuis (Amsterdam) Efficient sampling of molecular dynamics trajectories connecting arbitrary metastable states

Conformational changes in proteins (e.g. folding) can exhibit a rough free energy landscape involving multiple (meta)stable states. If there is a separation of time scales between the residence time in these stable states and the molecular crossing time, transitions between stable states can be considered rare events. Markovian state models can describe the long time kinetics of such systems, provided one has access to the hopping rates. For the unbiased study of the dynamics, rate constants, and mechanism of such rare events, transition path sampling has proven to be an effective method. However, in case there are multiple states, only one pair of states can be handled simultaneously. In this work, I present an efficient extension of the path ensemble that includes trajectories connecting two arbitrary stable states within the system. Combining this approach with transition interface sampling one directly obtains an expression for the rate constants of all possible transitions. I will show the efficiency of this approach for several model systems. For the efficiency the path switching behavior turns out to be the key ingredient. If some of the transitions become much more likely than others, one can use a biasing approach. I will show that such a biasing method is effective in some cases, but fails in others.

Stephen Bond (Illinois) Measuring and correcting algorithmic bias in molecular dynamics averages

Trajectories generated by classical molecular dynamics (MD) are best interpreted in a statistical sense due to the weak definition of initial conditions and chaotic nature of the underlying equations of motion. For this reason, accuracy and efficiency of a time integration scheme should be measured with respect to statistical averages, rather than deviations from an "exact trajectory." A practical MD time integrator must be both stable and statistically accurate, allowing for its use over long time intervals with large stepsizes. In this talk, I will survey some results from backward error analysis for geometric integrators and show how (under certain assumptions) these results can be applied to understanding the statistical properties of integrators for MD.

Sara Bonella (Rome 1) Free energy barriers for local hydrogen diffusion in sodium alanates

The single sweep approach, recently introduced by Maragliano and Vanden-Eijnden, combines an accelerated dynamics scheme with a radial basis representation to efficiently reconstruct complex free energy landscapes. In this presentation, single sweep is combined with ab initio molecular dynamics and brute force histogram calculations to investigate the identity of a mobile species observed in experiments that study the first dissociation reaction of sodium alanates - compounds that have attracted considerable interest as potential solid state hydrogen storage materials.

David Coker (Boston) Mixed quantum-classical methods for treating coherent quantum dynamics in model multichromophore light harvesting systems

Recent ultrafast multidimensional nonlinear spectroscopy experiments reveal that the extraordinary energy transfer efficiency of photosynthetic light harvesting antenna systems may arise from long lived quantum coherence of superpositions of multichromophore exciton states. The coherent evolution apparently arises from correlated motion of the chromophores mediated by common modes in their complex protein environments. The slow decoherence of the multichromophore excited state interrupts exciton localization and thus prevents dissipation during energy transport to the reaction center where energy storage is initiated. Understanding the mechanism of preservation of quantum coherence, and the role it plays in these nanoscale energy transmission grids is now amenable to microscopic study using recently developed mixed quantum classical simulation methods. Such simulations may yield important design principles that will be useful as we begin to engineer biomimetic structures for potential photovoltaic applications.

In this presentation we describe a new mixed quantum classical iterative scheme for evolving the full system density matrix that preserves the coherence of quantum superposition states and correctly represents natural dephasing processes in molecular simulations without having to make ad hoc assumptions. The approach is derived from first principles using a path-integral like scheme in which the forward and backward propagators that evolve the full system density operator are combined to give a short-time approximate propagator by linearizing in the difference between forward and backward paths of the environmental variables while keeping track of the interference effects in quantum subsystem variables that describe the quantum subsystem and result in dephasing. We discuss an efficient Monte Carlo approach for importance sampling paths through density matrix state label space, and environmental phase space, as well as the connection of the approach to different approximate schemes. Applications to model multichromophore light harvesting systems are presented.

Gabor Csanyi (Cambridge) Full sampling of atomic configurational spaces

We describe a method to explore the configurational phase space of chemical systems. It is based on the recently proposed Nested Sampling algorithm, and allows us to explore the entire potential energy surface (PES) efficiently in an unbiased way. The algorithm has a simple parameter that directly controls the trade-off between the resolution with which the space is explored and the computational cost. Within this framework, not only does the estimation of expectation values of arbitrary smooth operators at arbitrary temperatures become a simple post-processing step, but by analysing the topology of the samples we are able to visualise the PES in a new and illuminating way. This directly leads to the idea of identifying a discretely valued order parameter with basins and suprabasins of the PES allowing a straightforward and unambiguous definition of macroscopic states of an atomic system and the evaluation of the associated free energies.

Peter T. Cummings (Vanderbilt) Direct simulation-based evidence for fluid-solid transition in nanoconfined non-polar fluids

It is well known that the behavior of fluids confined to the order of a few nanometers may differ greatly from that of the corresponding bulk fluid.¹ For example, in the 1980s and early 1990s a large number of studies were reported on a variety of ultrathin liquid films confined between mica surfaces.² In all cases, irrespective of the fluid being studied, a common feature was noted: When the confinement reached the order of several molecular diameters, a rapid many-orders-of-magnitude increase in the viscosity of the confined fluid was observed, together with behavior typical of the stick-slip response of a crystalline solid structure.³ Unfortunately, as attention shifted to the nature of the transition from fluid-like to solid-like behavior a similar consensus was not achieved and, for over a decade, there has been intense debate as to whether it is a first-order (crystallization) or second-order (vitrification) phase transition. Whilst experimental attempts have been made to clarify this issue, and end this debate, they have been severely hampered by the extreme difficulty of observing directly what occurs at this scale.⁴

In contrast to the difficulties involved with experimental observation, molecular simulations techniques have an inherent atomic resolution and, as such, are ideally suited to the study of nanoscale phenomena. Because of this, a number of simulation studies have been undertaken with the aim of elucidating the nature of nanoconfinement-induced phase transitions. Unfortunately, these studies have typically made use of fairly simplistic models (e.g. mica represented by an fcc lattice of Lennard-Jones spheres),^5,6 or have performed the simulations in a manner that may be criticized as biased towards the formation of such structures, leading to concerns that they fail to accurately represent the experimental systems under investigation. Thus, with the aim of elucidating the nature of any change that may occur upon nanoconfinement, here we present unprecedentedly detailed molecular dynamics simulation results for dodecane, and cyclohexane confined between molecularly-smooth, but fully atomistically-detailed mica surfaces.

1. Van Alsten, J. and Granick, S., Langmuir 6, 876 (1990).
2. Alba-Simionesco, C. et al., Journal of Physics-Condensed Matter 18, R15 (2006).
3. Klein, J. and Kumacheva, E., Journal of Chemical Physics 108, 6996 (1998).
4. Bae, S. C. et al., Philosophical Transactions of the Royal Society a-Mathematical Physical and Engineering Sciences 366, 1443 (2008).
5. Cui, S. T., Cummings, P. T., and Cochran, H. D., Fluid Phase Equilibria 183, 381 (2001).
6. Jabbarzadeh, A., Harrowell, P., and Tanner, I., Tribology International 40, 1574 (2007).

Ron Elber (Texas) Milestoning

Milestoning is a combination of space partitioning and short molecular dynamics trajectories that link the partitions using a stochastic model. Kinetics and thermodynamics of the partitions are computed. The theory and the algorithm will be described and applications to peptide folding and conformational transitions in proteins will be discussed.

William Hoover (Ruby Valley) Shockwave Simulation Using Molecular Dynamics
(Joint work with C. G. Hoover)

Shockwaves are characterized for simple fluids. Density, stress, and temperature profiles from molecular dynamics simulations are compared to the predictions of continuum theories. There are many definitions of temperature within the shockwave, both one-dimensional and two-dimensional. Both configurational and kinetic temperatures are considered. Some kinetic definitions require knowing a local instantaneous stream velocity. The application of many of these definitions to the shockwave problem is considered.

Mikko Karttunen (Western Ontario) Multiscale modeling of lipid and surfactant systems

In this talk, I will focus on multiscale modeling¹ of lipid and surfactant systems. The focus will be on techniques based on the so-called Inverse Monte Carlo (IMC) technique and its refinements. The first examples will include single and multicomponent bilayers using different levels of coarse-graining, and the improvement of accuracy of the coarse-grained models by using a thermodynamic constraint.² The motivation for including such a constraint is that the basic IMC does not yield physically meaningful area compressibility or surface tension in coarse-grained bilayers. Then, I will discuss a modification to include some of the internal discrete degrees of freedom of the phospholipid chain to better describe chain ordering.³ Finally, I will discuss some of the advantages and disadvantages of the above approaches, and transferability of potentials. Time permitting, I will also discuss modeling of water using the so-called Mercedes-Benz model in 2- and 3-dimensions.^4,5

1. “Multiscale modeling of emergent materials: biological and soft matter”, T. Murtola, A. Bunker, I. Vattulainen, M. Deserno, and M. Karttunen PCCP 11, 1869-1892 (2009).
2. T. Murtola, E. Falck, M. Patra, M. Karttunen, I. Vattulainen, J. Chem. Phys. 121, 9156-9165.
3. “Systematic coarse-graining from structure using internal states: application to phospholipid / cholesterol bilayer”, T. Murtola, M. Karttunen, I. Vattulainen (submitted).
4. “Microscopic mechanism for cold denaturation”, C.L. Dias, T. Ala-Nissila, M. Karttunen, I. Vattulainen, and M. Grant Phys. Rev. Lett. 100, 118101 (2008).
5. “Three-dimensional "Mercedes-Benz" model for water”, C.L. Dias, T. Ala-Nissila, M. Grant, and M. Karttunen (submitted).

Steven Kenny (Loughborough) Atomistic Modelling of Problems in Materials Science

We will discuss two separate problems that we have been studying, namely the derivation of a multiscale technique which we have applied to modelling nanoindentation and the modelling of long timescale dynamics. In the first problem our main focus has been on modelling nanoindentation whilst in the second we are mainly interested in thin film growth or the modelling of radiation damage in materials.

In the first part of the talk we will present a multiscale model that links finite element methods to classical molecular dynamics. The methodology is novel in that it allows dynamical simulations to be performed in 3D and with arbitrary potential functions. We will illustrate the methodology by demonstrating its use in simulating the nanoindentation process.

Nanoindentation is an experimental technique used to measure the hardness of materials on the nanoscale. Traditionally either molecular dynamics or finite element methods have been used to model this process. The difficulty in applying finite element models to these problems comes from producing a suitable model of the plastic behaviour of the material. Whilst in molecular dynamics models the problem is one of the influence of the boundary conditions on the results.

In the second part of the talk we will present work that we have been performing on modelling long timescale dynamics primarily aimed at either modelling thin film growth or radiation damage. In particular the talk will focus on the calculation of the prefactor for the Arrhenius equation from atomistic simulations, comparing different approaches and discussing some cases where they breakdown.

Karl Kirschner (Fraunhofer Institute SCAI, St Augustin) Connecting quantum mechanics to coarse-grained simulations: the role of force field development

The quality of molecular simulations heavily depends on the reliability of the employed models. Aside from the functional form of the equations and the simulation code itself, force fields can be considered as the foundation of molecular dynamics (MD) simulations. However, there does not exist a well defined route for the development of new parameters when they become needed. In this talk we will present our efforts to derive new force field parameters in a timely manner. Those efforts include in-house tools of various kinds and purpose, ranging from input generating scripts to open source software kits. We will also discuss to what extend some of these tools could contribute towards a well-defined route for force field development. Force fields also represent a convenient link between quantum mechanics calculations, atomistic MD simulations, and coarse-grained (CG) simulations. Towards the latter, we will present our intentions and core features of an up-coming CG code named "Espresso++". Examples will be included at various points of the talk, comprising of ionic liquids, carbohydrates, and PMMA polymer simulations.

Mitchell Luskin (Minnesota) Non-ergodicity of Nosé-Hoover dynamics
(Joint work with Frederic Legoll and Richard Moeckel)

The Nosé-Hoover dynamics is a deterministic method that is commonly used to sample the canonical Gibbs measure. This dynamics extends the physical Hamiltonian dynamics by the addition of a "thermostat" variable that is coupled nonlinearly with the physical variables. The accuracy of the method depends on the dynamics being ergodic. Numerical experiments have been published earlier that are consistent with non-ergodicity of the dynamics for some model problems.

We will present the main ideas of our proof of the non-ergodicity of the Nosé-Hoover dynamics for the one-dimensional non-harmonic oscillator. We will also show that for some multidimensional systems the averaged dynamics for the limit of infinite thermostat "mass" have many invariants, thus giving theoretical support for either non-ergodicity or slow ergodization. Our numerical experiments for a two-dimensional central force problem and the one-dimensional pendulum problem give evidence for non-ergodicity.

We also present numerical experiments for the Nosé-Hoover chain with two thermostats applied to the one-dimensional harmonic oscillator. These experiments seem to support the non-ergodicity of the dynamics if the masses of the reservoirs are large enough and are consistent with ergodicity for more moderate masses.

David Manolopoulos (Oxford) Competing quantum effects in the dynamics of a flexible water model

Numerous studies have identified large quantum mechanical effects in the dynamics of liquid water. In this paper, we suggest that these effects may have been overestimated due to the use of rigid water models and flexible models in which the intramolecular interactions were described using simple harmonic functions. To demonstrate this, we introduce a new simple point charge model for liquid water, q-TIP4P/F, in which the O-H stretches are described by Morse-type functions. We have parameterized this model to give the correct liquid structure, diffusion coefficient, and infra-red absorption frequencies in quantum (path integral-based) simulations. The model also reproduces the experimental temperature-variation of the liquid density and affords reasonable agreement with the experimental melting temperature of hexagonal ice at atmospheric pressure. By comparing classical and quantum simulations of the liquid, we find that quantum mechanical fluctuations increase the rates of translational diffusion and orientational relaxation in our model by a factor of around 1.15. This effect is much smaller than that observed in all previous simulations of simple empirical water models, which have found a quantum effect of at least 1.4 regardless of the quantum simulation method or the water model employed. The small quantum effect in our model is a result of two competing phenomena. Intermolecular zero point energy and tunneling effects destabilize the hydrogen bonding network, leading to a less viscous liquid with a larger diffusion coefficient. However this is offset by intramolecular zero point motion, which changes the average water monomer geometry resulting in a larger dipole moment, stronger intermolecular interactions, and slower diffusion. We end by suggesting, on the basis of simulations of other potential energy models, that the small quantum effect we find in the diffusion coefficient is associated with the ability of our model to produce a single broad O-H stretching band in the infra-red absorption spectrum.

Glenn Martyna (Watson Research Centre, IBM) Treating Manybody Polarization and Manybody Dispersion in Complex Systems: The Quantum Drude Oscillator Formalism

There are many physical systems where the nonpairwise additive nature of polarization and dispersion interactions becomes very important, in particular, the complex heterogeneous systems of interest in chemistry, biology and physics. For example, the dipole moment of water changes from 1.85 Debye in the gas phase to approximately 2.6 Debye in the bulk liquid and attains intermediate values at hydrophobic interfaces due to manybody polarization. Similarly, although the bulk properties of hydrophobic fluids can be modeled using a pair potential, this underestimates the surface tension by 30% due to manybody dispersion interactions. In order to model both the full manybody polarization and dispersion interactions in atomic and molecule systems, a system of quanized Drude oscillators is introduced and a O(N) simulation method based on quantum path integrals is described using diagrammatic expansions of the propagator. Applications to simple systems, xenon, and more complex systems, water, are given.

Florian Muller-Plathe (Darmstadt) Trends in materials simulation - a molecular dynamics miscellany

Molecular Dynamics simulations are at present arriving in applied science and engineering as a tool to help solve real problems of real materials. As a consequence, simulators have to face new challenges, in addition to finding realistic simulation models: composite materials; materials at complex interfaces; complicated transport processes; materials properties being inherently of multi-scale nature; interaction, for which there exist no trusted physical descriptions; the bridge to continuum models; and so on. We will sketch a few recent developments in the fields of interfaces, transport and multi-scale methods.

Adrian Mulholland (Bristol) Protein dynamics and enzyme catalysis: insights from simulations

Enzymes exhibit complex dynamics and conformational changes associated with their catalytic cycles. The question of the possible functional roles of these motions is at the heart of current debates in structural biology.

Simulations have an important part to play in resolving these challenging issues. Classical molecular dynamics simulations can identify functionally relevant motions on the nanosecond timescale: for example, simulations of human scavenger decapping enzyme have identified a cooperative periodic opening and closing of the dimer. For modelling enzyme-catalysed reaction mechanisms, combined quantum mechanics/molecular mechanics (QM/MM) methods are a good approach. QM/MM calculations can now be carried out with high-level ab initio methods capable of high accuracy. Calculated activation energies for some enzyme reactions, such as those catalysed by chorismate mutase and p-hydroxybenzoate hydroxylase, agree well with experiment, indicating that dynamical effects on the rate are relatively small. Quantum tunnelling can be a significant factor, e.g. in the key proton transfer step of aromatic amine dehydrogenase. Semiempirical QM/MM VTST/SCT calculations give results (e.g. kinetic isotope effects) in good agreement with experiment, and show that quantum tunnelling is dominant in the reaction. Simulations indicate that long-range coupled motion of the protein is apparently not involved in 'driving' tunnelling in this case.

Complex conformational effects are observed in some enzymes: for example, hydrolysis of the 'sleep inducer' oleamide in fatty acid amide hydrolase appears to involve a minor conformation of the protein. QM/MM modelling of citrate synthase has suggested an unusual mechanism involving arginine acting as an acid; this mechanism also suggests a means of coupling chemical and conformational changes in this enzyme.

1. “Biomolecular simulation and modelling: status, progress and prospects”, M.W. van der Kamp et al. Royal Soc. Interface 5 Suppl 3, 173-190 (2008).
2. “High-accuracy computation of reaction barriers in enzymes”, F. Claeyssens et al. Angewandte Chemie Intl. Edn. 45, 6856-6859 (2006).
3. “Cooperative symmetric to asymmetric conformational transition of the apo-form of scavenger decapping enzyme revealed by simulations”, U. Pentikäinen et al. Proteins: Struct. Funct. Bioinf. 70, 498-508 (2008).
4. “Conformational effects in enzyme catalysis: reaction via a high energy conformation in fatty acid amide hydrolase”, A. Lodola et al. Biophys. J. 92, L20-L22 (2007).
5. “Atomic Description of an Enzyme Reaction Dominated by Proton Tunneling”, L. Masgrau et al. Science 312, 237-241 (2006).

Massimo Noro (Unilever Discover Port Sunlight) Modelling (skin) lipid bilayers
(joint work with Chinmay Das, Peter Olmsted, Dominic J Tildesley)

Skin lipids models are essential to understand the interaction of personal care products with the outer layers of skin. We have studied model systems at the atomistic level, to answer questions about the structure and dynamics of these layers. In doing so we have identified the role of the three main lipid components (ceramides, cholesterol, and free fatty acids) in building the typical lamellar structure and generating the local mechanical and barrier properties, which are essential for a healthy skin.¹ We will demonstrate how atomistic simulations of skin lipids are used to calculate accurately the free energy cost for crossing the skin protective barrier.²

1. Notman, R., DenOtter W., Noro, M.G., Anwar, J. and Briels, W., “Simulations of skin barrier function: Free energies of hydrophobic and hydrophilic transmembrane pores in ceramide bilayers”, Biophysics J., 95, 4763, 2008.
2. C. Das, P. Olmsted, M.G. Noro “Simulations studies of Stratum Corneum lipid mixtures”, Biophysics J. (submitted)

Rebecca Notman (University of Warwick) Molecular dynamics simulations of peptides interacting with quartz surfaces

The specific binding of peptide sequences to inorganic surfaces is of considerable interest for its applications in biotechnology and materials science. Whilst there has been rapid experimental growth in this area, a fundamental understanding of the interactions that regulate binding is still lacking. In this work we aim to characterise the interactions of strong-binding peptides with quartz surfaces and elucidate the key mechanisms of binding. An increased understanding of peptide-inorganic interactions would be a significant step towards the control of molecular recognition and self assembly processes.

We have used molecular dynamics (MD) simulations to build up a detailed picture of the peptide-quartz interface system. Replica exchange MD simulations explored the solution structures of quartz binding peptides to identify key structural properties of the peptides that may be important for binding. We have also carried out simulations of quartz surfaces in water in order to investigate the water structure at these surfaces and the role that this may play in the peptide-surface interactions. Finally, we have carried out simulations of the peptides binding to the quartz surface where we gain further insights into the relationship between peptide structure and strong-binding function.

Frank Pinski (Cincinnati) Transition paths for molecular motion: What are the characteristics of the most probable paths?
(Work done with Andrew Stuart, Mathematics Institute, University of Warwick.)

We explore transition paths of particles moving across energy barriers. The underlying motion is described by Brownian Dynamics, the overdamped limit of Langevin’s equations. The density of paths is governed by a well-known probability measure. For a variety of models, we use gradient descent to find the Most Probable Paths (MPPs), paths that exist at a probability maximum. In the talk, we will describe characteristics of such paths. Some are unexpected; others can be anticipated. Many of the features of paths in Lennard-Jones clusters are also seen in simple one- and two-dimensional models, giving a framework for a discussion of our findings.

Matej Praprotnik (ETH Zürich) Adaptive Resolution Molecular Dynamics Simulation

The Adaptive Resolution Scheme (AdResS) for hybrid atomistic/mesoscale molecular dynamics (MD) simulations will be presented. The key feature of this particle-based method is that it allows for a dynamical change of molecular resolution by changing the number of molecular degrees of freedom on-the-fly during the course of an MD simulation. So far, we have applied AdResS to molecular liquids. The system is modeled in different domains at different levels of detail while the liquid molecules move freely between the regions. In order to cover even wider range of length scales we have also included the continuum description of a liquid in our model and performed a concurrent triple-scale simulation of the liquid. The triple-scale scheme, which is shown to correctly describe the hydrodynamics, successfully sorts out the problem of large molecule insertion in the hybrid particle-continuum simulations of molecular liquids. This approach allows for efficient grand-canonical MD simulations of open systems involving an explicit solvent, e.g., liquid water.

David Quigley (Warwick) Simulating orientational specificity in the growth of calcite on self-assembled monolayers

Crystalline mineral phases appear in nature as a major component of biocomposite materials such as shells, coral, teeth and bone. The organisms which exploit these phases exert a remarkable degree of control over the morphology and orientation of the growing crystal, sometimes employing multiple polymorphs within a single structure. There is considerable interest in understanding and reproducing the mechanisms behind this control, with the aim of developing biomimetic processes for materials engineering.

One such biomimetic process is the crystallisation of calcite on a monolayer of self-assembled carboxylic acid molecules.¹ This has been studied in a range of experiments, demonstrating orientational specificity with a preference for nucleation at the (0112) plane. In contrast, purely epitaxial arguments suggest the (0001) plane as preferable. Furthermore, the specificity is sensitive to the length of the chain, exhibiting an odd-even effect.² Simulation of the amorphous to crystalline transition can playa valuable role in understanding this phenomenon at the nucleation stage.

Previous simulations have been hampered by use of elevated temperatures to overcome the timescale problem associated with simulating crystallisation.³ This precludes the possibility of treating the monolayer dynamically, and the inclusion of explicit water. In contrast, the metadynamics method⁴ has proven useful in the context of freezing from the melt without the need to impose unrealistic thermal conditions.⁵ We have adapted our implementation of this method⁶ to the crystallisation of amorphous calcium carbonate⁷ and applied to heterogeneous nucleation on self-assembled monolayers (SAMs).

By including the dynamics of the organic chains, we are able to simulate the orientational selectivity in this system without imposing any structural defects. Furthermore, we are able to investigate the role of head-group ionisation, chain flexibility and chain length in the selection process, generating a number of crystal orientations which agree directly with experiment. As a result of this investigation we propose a "complimentary nucleation" process in which both the organic and mineral phases adapt to achieve a local charge epitaxy and trigger a nucleation event.

1. A. Travaille. L. Kaptijn, P. Verwer, B. Hulsken, J. Elemans, R. Nolte, and H. van Kempen, J. Am. Chem. Soc. 125, 11571 (2003).
2. Y. Han and J. Aizenberg, Angew. Chem. Int. Edit. 42, 3668 (2003).
3. C. L. Freeman, J. H. Harding, and D. M. Duffy, Langmuir 24, 9607 (2008).
4. A. Laio and M. Parrinello, Proc. Nat. Acad. Sci. 99, 12562 (2002).
5. D. Quigley and P. M. Rodger, J. Chem. Phys. 128 (2008), ISSN 0021-9606.
6. D. Quigley and P. M. Rodger, Mol. Simul. 35, 613 (2009).
7. D. Quigley and P. M. Rodger, J. Chem. Phys. 128 (2008), ISSN 0021-9606.

Weiqing Ren (New York) A seamless multiscale method and its application to complex fluids

I will present a seamless multiscale method for the study of multiscale problems. The multiscale method aims at capturing the macroscale behavior of the system which is modeled by a (incomplete) macroscale model. This is done by coupling the macro model with a micro model: The macro model provides the necessary constraint for the micro model and the micro model supplies the missing data (e.g. the constitutive relation or the boundary conditions) needed in the macro model. In the multiscale method, the macro and micro models evolve simultaneously using different time steps, and they exchange data at every step. The micro model uses its own appropriate (micro) time step. The macro model uses a macro time step but runs at a slower pace than required by accuracy and stability considerations in order for the micro dynamics to have sufficient time to adapt to the environment provided by the macro state. The method has the advantage that it does not require the reinitialization of the micro model at each macro time step or each macro iteration step. The data computed from the micro model is implicitly averaged over time. In this talk, I will discuss the algorithm of the multiscale method, the error analysis, and its application to complex fluids. If time allows, I will also discuss the stability of different coupling schemes in domain-decomposition type of multiscale methods.

Debra Searles (Griffith) The dissipation function and the properties of highly nonequilibrium systems

The role of the dissipation function in the study of highly nonequilibrium systems will be discussed. The dissipation function was first identified in development of fluctuation theorems for nonequilibrium systems. It can be shown to be equivalent to the rate of spontaneous entropy production under some circumstances. The fluctuation theorem can then be used to show that, on average, the time integral of the dissipation function is greater than zero for nonequilibrium systems, demonstrating the emergence of irreversibility and the second law of thermodynamics. However, in the last couple of years it has become apparent that the dissipation function has a more central role in the study of nonequilibrium behaviour. It is the fundamental argument of linear and nonlinear response theory and controls the relaxation of systems to equilibrium.

One of the practical implications of these results is in the study of complex, confined systems that are driven far from equilibrium and thermostatted using realistic thermostats. In the past, expressions for nonlinear response only applied to homogeneously driven and thermostatted systems and now they can be extended to treat these more complicated systems. The dissipation function will be discussed and results of simulations on confined systems presented.

Christof Schuette (Free University Berlin) Markov state models: Theory and applications

The efficient determination of reliable models for the long-term dynamical behavior is one of the grand challenges in molecular dynamics. Direct numerical simulation of dynamical effects that are crucial for biomolecular function in many cases requires infeasibly long simulation times. Therefore the question of how to construct Markov state models (MSM) for the effective dynamics attracted much attention recently. The talk will present some new approaches to this problem. In particular we will discuss how accurate such MSM can be and how to algorithmically construct them from parallel molecular dynamics simulations.

Robert D. Skeel (Purdue) Maximum flux transition paths of conformational change

Given two metastable states A and B of a biomolecular system, the problem is to calculate the likely paths of the transition from A to B. Such a calculation is more informative and more manageable if done for a reduced set of collective variables chosen so that paths cluster in collective variable space. The computational task becomes that of computing the "center" of such a cluster. A good way to define the center employs the concept of a committor, whose value at a point in collective variable space is the probability that a trajectory at that point will reach B before A. The committor "foliates" the transition region into a collection of isocommittors. The maximum flux transition path (MFTP) is defined as a path that crosses each isocommittor at a point which (locally) has the highest crossing rate of distinct reactive trajectories. (The MFTP bears no close relation to the MaxFlux method of Huo and Straub.) To make the calculation tractable, three uncontrolled approximations are introduced - one fewer than what is needed to derive a minimum free energy path (MFEP). It is shown that an MFTP is qualitatively superior to an MFEP (though quantatively not too different). These qualities aid in the construction of simple and robust algorithms. Such an algorithm and its performance is discussed.

Christopher Sweet (Notre Dame) Simulation of Long Timescale Dynamics of Proteins Using Normal-Mode Langevin

The Normal Mode Langevin dynamics integrator (NML) is a novel method that approximates the kinetics or thermodynamics of a biomolecule by a reduced model based on a normal mode decomposition of the dynamical space. The basis set uses the eigenvectors of a mass re-weighted Hessian matrix, calculated with a biomolecular force field. This set is ordered according to the eigenvalues, which has a physical meaning (square-root of the mode frequency). The low frequency eigenvalues correspond to more collective motions, whereas the highest frequency eigenvalues are the limiting factor for the stability of the integrator.

In this method the higher frequency modes are minimized while respecting the subspace of low frequency dynamical modes, which are accurately propagated. To provide accurate partitioning of the system throughout the simulation requires expensive re-diagonalization at regular intervals, which generally scales as O(N³) time and O(N²) memory. To overcome this limitation an inexpensive method of re-diagonalization has been derived that captures the essential features of the low frequency modes, while meeting the stability criterion of the integrator with respect to the largest included eigenvalue. This method scales as O(N log N) and O(N) time, thus retaining the potential efficiency increase of the underlying NML method. The combined method we call Coarse-grained Normal Mode Langevin (CNML). CNML is capable of accommodating large conformational change of proteins, as shown in simulation of WW domain folding.

Karl Travis (Sheffield University) A study of radiation-induced disorder in Pu-doped ceramics using molecular dynamics and topological analysis

Zirconolite (CaZrTi₂O₇), is one of the phases of SYNROC-C, a ceramic based wasteform proposed for the immobilization of plutonium. Naturally occurring zirconolite, present in mineral deposits, is known to immobilize Pu for geological time periods. However, a full understanding of why it has such a strong radiation resistance and details of the mechanism which lies behind the radiation-induced crystalline to metamict phase transition, is currently lacking.

Molecular Dynamics (MD) is the most appropriate numerical tool to study the formation and migration of atomic defects created when a heavy recoil nucleus is released following α-decay by a Pu atom. By tracking the motions of ca. 3/4 million atoms it is possible to study in real time the growth of damage cascades and the subsequent recovery (or lack thereof), by annealing processes, using realistic recoil energies.

We report the first ever MD simulations of radiation damage in Pu-doped crystalline zirconolite. The damage is characterized using both defect analysis and topological modeling. We find Zirconolite is very radiation resistant in agreement with experimental observations. Topological analysis of damaged structures yields significant information that could provide an explanation for this phenomenon; the final structures are determined largely by the re-arrangement of TiO_x polyhedra into stable configurations of an edge sharing nature.

Mark Tuckerman (New York) Developments in conformational sampling using molecular dynamics

Sampling conformational equilibria on rough energy landscapes is one of the computational grand challenge problems. If met, many important problems, most notably biomolecular structure prediction, could be significantly impacted. In this talk, I will discuss a variety of approaches that can accelerate conformational sampling molecular dynamics from simple mass scaling to adiabatic dynamics and spatial warping transformations. The performance of these approaches will be demonstrated for simple toy examples, polymer chains, where comparisons to replica-exchange Monte Carlo will be made, and small peptides in solution.

Eric Vanden-Eijnden (New York) Theory and Modeling of Reactive Events

In the first part of the talk, I will explain why we may need to go beyond the standard framework of transition state theory (TST) to describe activated processes and reactive events, and I will present another framework, termed transition path theory (TPT), that permits to do that. Unlike TST, which gives mainly an expression for the rate of the reactive event, TPT describes more fully the statistical properties of the reactive trajectories (i.e. those trajectory by which the reactive event occurs), in particular in terms of their probability density function and their probability current.

In the second part of the talk, I will describe how TPT can be use to design and/or improve numerical methods for computing the pathways and rate of reactive events. I will focus in particular on the string method and extensions upon milestoning.

Jiri Vanicek (EPFL Lausanne) Efficient evaluation of accuracy of quantum molecular dynamics using dephasing representation

Ab initio methods for electronic structure of molecules have reached a satisfactory accuracy for calculation of static properties, but remain too expensive for quantum dynamical calculations. We propose an efficient semiclassical method for evaluating the accuracy of a lower level quantum dynamics, as compared to a higher level quantum dynamics, without having to perform any quantum dynamics. The method is based on dephasing representation of quantum fidelity and its feasibility is demonstrated on the photodissociation dynamics of CO₂. We suggest how to implement the method in existing molecular dynamics codes and describe a simple test of its applicability.

Preprint: B. Li, C. Mollica, and J. Vanicek: “Efficient evaluation of accuracy of molecular quantum dynamics using dephasing representation,” arXiv:0905.1630v1 [physics.chem-ph] (2009), http://arxiv.org/abs/0905.1630v1

Rodolphe Vuilleumier (Université Pierre et Marie Curie, Paris) Microscopic flow around a diffusing particle

The fundamental model for diffusion is the Stokes-Einstein approach. It relates the diffusion coefficient to the friction coefficient which is in turn extracted from the hydrodynamic friction felt by a moving sphere in a liquid environment. The hydrodynamic flow around the sphere, commonly called Stokes flow, depends on the boundary conditions on the sphere. Two boundary conditions are usually considered: stick boundary conditions (no velocity at contact), as in the orginal approach, or slip boundary conditions (finite tangential velocity). As these two boundary conditions lead to different diffusion coefficient as a function of particle size, the usual way to investigate these from numerical simulations has been to compute the diffusion coefficient as a function of particle size. This leads for diffusion of Lennard-Jones particles to the suggestion of slip boundary conditions. Here, we want to present a direct calculation of the microscopic velocity field around a diffusing particle from numerical simulations. This allow for comparison between the atomic model and the hydrodynamics approach. It is first demonstrated that the hydrodynamics flow is well recovered after only a few atomic radius from the tagged particle. However, two effects of the velocity fluctuations of the diffusing particle are evidenced. The boundary conditions are shown to include a finite normal velocity at contact, which would seem at first in contradiction with the non-penetrability of the particles but is an effect of fluctuations. Then, the flux of momentum in the flow is not determined solely by viscosity but also from the diffusion of the tagged particle.

This gives new insight in the diffusion process in liquids and in particular on the role of fluctuations at the atomic level. We will also try to show how this method could be also used in order to provide a more mechanistic description of the diffusion processes, in terms similar to reaction paths and transition states.

David Wales (Cambridge) Energy landscapes and dynamics: from clusters to biomolecules

Coarse-graining the potential energy surface into the basins of attraction of local minima, provides a computational framework for investigating structure, dynamics and thermodynamics in molecular science. Steps between local minima form the basis for global optimisation via basin-hopping and for calculating thermodynamic properties using the superposition approach and basin-sampling. To treat global dynamics we must include transition states of the potential energy surface, which link local minima via steepest-descent paths. We may then apply the discrete path sampling method, which provides access to rate constants for rare events. In large systems the paths between minima with unrelated structures may involve hundreds of stationary points of the potential energy surface. New algorithms have been developed for both geometry optimisation and making connections between distant local minima, which allow us to treat such systems. Applications will be presented for a wide variety of atomic and molecular clusters, glass-formers and biomolecules. Results for folding, misfolding and aggregation of peptides and proteins illustrate how experimental time and length scales can now be addressed for such systems.

Selected Publications:
1. B. Strodel and D.J. Wales, Chem. Phys. Lett., 466, 105-115 (2008).
2. D.J. Wales and T.V. Bogdan, J. Phys. Chem. B, 110, 20765-20776 (2006).
3. D.J. Wales, Int. Rev. Phys. Chem., 25, 237-282 (2006)
4. D.J. Wales, “Energy Landscapes”, Cambridge University Press, Cambridge, 2003.
5. D.J. Wales and H.A. Scheraga, Science, 285, 1368-1372 (1999)

See also:
Mathematics Research Centre
Mathematical Interdisciplinary Research at Warwick (MIR@W)
Past Events 
Past Symposia