FES convergence issue on metadynamics simulation

[Lee Sangmin]

Hi. I have some problems while running OpenMM 7.5 well-tempered metadynamics module. I'm trying to calculate the two-dimensional FES utilizing metadynamics module. Two bias variables are defined accordingly:

biasCA = BiasVariable(cvCA,0.0,2.0,0.1,False)
biastheta = BiasVariable(cvtheta, -rad/2, rad/2, distrad, False)
(rad is a value of 3.14159265 and disrad is a value of 0.174532925, which is equal to rad/18)

Now the metadynamics object is defined as follows:

meta = Metadynamics(system, [biasCA,biastheta], temperature, 10, 0.25*kilojoules_per_mole, 500, 500, biasDir='./')

Simulation was running pretty well without any issue, but I've found that the free energy results got highly fluctuated depending on the running time. I tried to check the convergence of FES by running multiple simulations with different runtime (20ns, 40ns, 60ns, ...). Although short simulation may represent fluctuating behavior due to several FES regions unexplored, even FES with more than 200ns indicate same issue (like contour lines are irregularly positioned and free energy value is changed quite widely up to ~10kT), which makes me hard to figure out whether FES is well-converged or not.

I wonder is there any way to rationalize the convergence of FES results. Do I need to change the biasFactor and Gaussian height value to search for the proper setting? To add, I'd like to understand the way to interpret the bias npy file written during the simulation. Is it possible to analyze the FES exploration using such file?

[Peter Eastman]

So far as I know, there's no good way to determine convergence except for seeing when the results stop changing. You only know about the regions it has explored. You don't know whether there are still regions it hasn't explored yet. Maybe someone with more experience running metadynamics simulations can offer more guidance.

The free energy can be directly calculated from the biases. The getFreeEnergy() function performs the calculation for you.

[Emilio Galllicchio]

One way to test convergence is to monitor the number of independent visits to a set of chosen conformational states. If a state is visited only once during the simulation, its free energy is unlikely to be converged.

Testing convergence of free energy landscapes for bounded regions of conformational space, such as one defined by periodic torsional angles, is easier than for unbounded regions, such as those defined by distances. In the former case, one would expect to eventually visit all combinations of angles to some degree. In the latter case, as Peter pointed out, there isn't a general way to know if there are important unexplored regions.

The well-tempered metadynamics algorithm implemented in OpenMM has a temperature boost factor parameter that can be used to overcome free energy barriers. It's a double-edged sword, however: the higher the temperature factor, the more metadynamics explores uninteresting high free energy regions of conformational space where it could become lost.

The .npy file written by Metadynamics holds the biasing potential. It can be independently loaded into a Continuous{1,2,3}DFunction() to run biased MD with the biasing potential built by Metadynamics. See: https://github.com/openmm/openmm/blob/m ... ynamics.py