Hybrid Monte Carlo Methods

For a simulation to equilibrate quickly it is extremely important (if not sufficient) to use methods of exploring a representive selection of configurations as quickly as possible. The power of a MD algorithm is that it takes relatively large steps in phase space and allows the positions and velocities of all the degrees of freedom of the system to be updated simultaneously. However, one disadvantage is that the calculated trajectories are not necessarily ergodic, particularly when the system contains harmonic degrees of freedom. MC simulations do not have this problem because they can improve ergodicity by making random changes to the configurations. The drawback to traditional Metropolis MC simulations is the rate at which they explore phase space, particularly for dense systems or those with long-range interactions [28,63,89].

Given that the two techniques complement each other in their ability to explore the phase space, a variety of hybrid methods have been devised, in which the simulation algorithm alternates between MD and MC [7,39,51,53,57,63,73,110]. These HMC algorithms combine the large steps taken in phase space by MD with the ability of MC to change direction of the trajectory randomly. Thus the MC part of the simulation ensures ergodicity and eliminates inaccuracies in the energy and the MD part speeds up the simulation by allowing large steps to be taken in phase space. The HMC algorithm was originally developed for use in quantum chromodynamics and traditionally has been used for condensed-matter systems [39,53,57,73,110]. Recently, different variants have been successfully developed for biomolecular simulations [51,63,159]. Authors of [63] applied HMC to protein folding and authors of [51] used HMC to sample multi-minima energy landscapes of triribonucleotide, a small RNA segment. Zhang [159] showed that it can also be used to test protein potential functions and used it to refine protein structures. Gromov and Pablo [57] implemented HMC algorithms to RESPA and used it to efficiently simulate mixtures of Lennard-Jones chains.

We plan to develop a MOLLY-based HMC algorithm for MTS simulations. For a HMC method to be valid, it has to satisfy the property of detailed balance. It has been shown that the property of detailed balance needed is satisfied if the integrator is both reversible and symplectic [110,136]. The multiscale MOLLY integrator, that we are proposing and is described in Sections 2.1 and 2.2, is both symplectic and reversible. This MOLLY based MTS-HMC technique, in particular, will be used to compute free energies of the systems. The free energy is one of the most important quantities in thermodynamics as free energy difference gives us binding constants of different ligands to a binding site and the partition coefficient of a solute in different solvents, two quantities widely used in drug design. This will be attempted in collaboration with Prof. Edward Maginn of the department of chemical engineering, who has developed configurational biased MC methods [99]; HMC will be a way of generalizing his methods to arbitrary systems, particularly biomolecules. We also intend to combine it with computational geometry algorithms, to produce flexible docking programs, cf. [55,85,86,93,94,95,127,128]. This part will be done in collaboration with Prof. Danny Chen in our department. In these simulations geometric methods will produce geometrically feasible dockings and the HMC free energy computation will serve as a score, cf. [133].