Scalable hybrid Monte Carlo statistical sampling

For a simulation to equilibrate quickly it is extremely important, if not sufficient, to use methods that explore a representative selection of configurations as quickly as possible. The power of an 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. Monte Carlo (MC) simulations do not have this problem because they can improve ergodicity by making random changes to the configurations. The drawback to traditional MC simulations is the slow rate at which they explore phase space, particularly for dense systems or those with long-range interactions [11,39,59].

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 [4,22,31,33,36,39,43,68]. These hybrid Monte Carlo (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 while the MD part speeds up the simulation by allowing large steps to be taken in phase space. The great difficulty of HMC is that it does not scale well with system size, and thus it becomes less attractive for large biomolecules. We propose a biased HMC algorithm that scales almost linearly, and would give a speedup of an order of magnitude over conventional HMC for a system with 100,000 atoms, at the cost of moderately increased storage, cf. Sections 2.2 and 3.3.