1	Improved Hybrid Monte Carlo method for conformational sampling
2	Overview
3	Questions related to sampling Sampling Compute equilibrium averages in NVT (or other) ensemble Examples: Equilibrium distribution of solvent molecules in vacancies Free energies Pressure Characteristic conformations
4	Classical molecular dynamics Newton’s equations of motion: Atoms Molecules CHARMM force field (Chemistry at Harvard Molecular Mechanics)
5	Energy Terms Described in the CHARMm forcefield
6	Energy Functions
7	Molecular Dynamics – what does it mean?
8	Hybrid Monte Carlo I
9	Hybrid Monte Carlo II
10	Hybrid Monte Carlo III Hybrid Monte Carlo: Apply stochastic step (e.g., regenerate momenta) Use reversible symplectic integrator for MD to generate the next proposal in MC: Hamiltonian dynamics preserve volume in phase space, and so do symplectic integrators (determinant of Jacobian of map is 1) It is simple to make symplectic integrators time reversible Apply Metropolis MC acceptance criterion
11	Hybrid Monte Carlo IV Advantages of HMC: HMC can propose and accept distant points in phase space Make sure new SHMC has high enough accuracy HMC can move in a biased way, rather than in a random walk like MC (distance n vs sqrt(n)) Make L long enough in SHMC HMC is a rigorous sampling method: systematic sampling errors due to finite step size in MD are eliminated by the Metropolis step of HMC. Make sure bias is eliminated by SHMC
12	Hybrid Monte Carlo V
13	Hybrid Monte Carlo VII The key problem in scaling is the accuracy of the MD integrator Higher order MD integrators could help scaling Creutz and Gocksch (1989) proposed higher order symplectic methods to improve scaling of HMC In MD, however, these methods are more expensive than the gain due to the scaling. They need several force evaluations per step O(N) electrostatic methods may make higher order integrators in HMC feasible for large N
14	Overview
15	Improved HMC Symplectic integrators conserve exactly (within roundoff error) a modified Hamiltonian that for short MD simulations (such as in HMC) stays close to the true Hamiltonian Sanz-Serna & Calvo 94 Our idea is to use highly accurate approximations to the modified Hamiltonian in order to improve the scaling of HMC
16	Shadow Hamiltonian
17	Example Shadow Hamiltonian (partial)
18	SHMC Algorithm
19	SHMC Nearly linear scalability of acceptance rate with system size N Computational cost of SHMC, O(N ^(1+1/2m)) where m is accuracy order of integrator Extra storage (m copies of q and p) Moderate overhead (10% for medium protein such as BPTI)
20	Overview
21	Evaluating MC methods I Is SHMC sampling from desired distribution? Does it preserve detailed balance? Used simple model systems that can be solved analytically. Compared to analytical results and HMC. Examples: Lennard-Jones liquid, butane Is it ergodic? Impossible to prove for realistic problems. Instead, show self-averaging of properties. Computed self-averaging of non-bonded forces and potential energy
22	Evaluating MC methods II Is system equilibrated? Average values of set of properties fluctuate around mean value Convergence to steady state from Different initial conditions Are statistical errors small? Runs about 10 times longer than slowest relaxation in system Estimated statistical errors by block averaging Computed properties (torsion energy, pressure, potential energy) Vary system sizes (4 – 1101 atoms) What are the sampling rates? Cost (in seconds) per new conformation Number of conformations discovered
23	Systems tested
24	ProtoMol: a framework for MD
25	SHMC implementation
26	Experiments: acceptance rates I Numerical experiments confirm the predicted behavior of the acceptance rate with system size. Here, for fixed acceptance rate, the maximum time step for HMC and SHMC is presented
27	Experiments: acceptance rates II
28	Experiments: acceptance rates III
29	Average of observable Average torsion energy for extended atom Butane (CHARMM 28) Each data point is a 114 ns simulation Temperature = 300 K
30	Sampling Metric (or how to count conformations) For each dihedral angle (not including Hydrogen) do this preprocessing: Find local maxima, counting periodicity Label ‘wells’ between maxima During simulation, for each dihedral angle Phi[i]: Determine ‘well’ Phi[i] occupies Assign name of well to a conformation string String determines conformation (extends method by McCammon et al., 1999)
31	Sampling rate for decalanine (dt = 2 fs)
32	Sampling rate for 2mlt
33	Sampling rate comparison C is number of new conformations discovered Cost is total simulation time divided by C Each row found by sweeping through step size (for alanine, between 0.25 and 2 fs; for melittin and bpti between 0.1 and 1 fs) and simulation length L
34	Overview
35	Summary and discussion SHMC has a higher acceptance rate than HMC, particularly as system size and time step increase SHMC discovers new conformations more quickly SHMC requires extra storage and moderate overhead. For large time steps, weights may increase, which harms the variance. This dampens maximum speedup attainable More careful coding is needed for SHMC than HMC For large N, higher order integrators may be competitive with SHMC Instead of reweighting, one may modify the acceptance rule
36	Future work Multiscale problems for rugged energy surface Multiple time stepping algorithms plus constraining Temperature tempering and multicanonical ensemble (e.g., method of Fischer, Cordes, & Schutte 1999) Potential smoothing Combine multiple SHMC runs using method of histograms Include other MC moves (e.g., change essential dihedrals or Chandler’s moves) System size Parallel multigrid or multipole O(N) electrostatics Applications Free energy estimation for drug design Folding and metastable conformation
37	Acknowledgments Graduate student: Scott Hampton Dr. Thierry Matthey, co-developer of ProtoMol, University of Bergen, Norway Students in CSE 598K, “Computational Biology,” Spring 2002 Tamar Schlick for her deligthful new book, Molecular Modeling and Simulation – An Interdisciplinary Guide Dr. Robert Skeel, Dr. Ruhong Zhou, and Dr. Christoph Schutte for valuable discussions Dr. Radford Neal’s presentation “Markov Chain Sampling Using Hamiltonian Dynamics” (http://www.cs.utoronto.ca ) Dr. Klaus Schulten’s presentation “An introduction to molecular dynamics simulations” (http://www.ks.uiuc.edu )