thermo support for metadynamics

[liyigerry]

Originally, transition-based reweighting analysis method (tram) has been developed to deal with enhanced sampling, such as umbrella sampling or replica exchange. However, metadynamics has not been well supported with tram in documents. Can anyone provide detailed examples?

For the sake of simplicity, we have only one trajectory generated from metadynamics. In this trajectory, each snapshot corresponds a Hamiltonian due to adding bias more frequently than step of snapshot output. So we can generate ttrajs with linear index of length of trajectory and get dtrajs by clustering on collective variables of metadynamics.

In documents of pyemma, the btrajs[i][t, k] is the reduced bias energy in trajectory i and time step t for the k’th thermostate. Here, we have only one trajectory, so i is 0.

How to understand the bias at step t, or snapshot t, for the thermostate of k? For example, What is the difference between bias for btrajs[0, 0] (step 0 and thermostate 0) and btrajs[0, 1] (step 0 and thermostate 1).

Waiting for your help.

Thanks.

[fabian-paul]

Hi Yi Li,

this question has come up before, see issue https://github.com/markovmodel/PyEMMA/issues/1014 and the example .

While it is possible to use TRAM for metadynamics simulations, there are a couple of practical challenges that have to be addressed for metadynamics data. These are mainly:

• Memory requirements: the RAM usage for the biases scales as T*K*K where T is the aggregated trajectory length and K is the number of themodynamic states (i. e. bias potentials). As the number K is very large in metadynamics compared to other methods, the RAM requirements should be checked beforehand. In principle it would be possible to compress the biases in RAM (e. g. by not storing them directly but by re-implementing the bias potential function inside Pyemma. However this functionality is currently not provided by PyEmma.)
• Ergodicity: TRAM is based on MSMs and therefore required ergodic (i. e. reversible) sampling in each of the thermodynamic states. In every thermodynamic state, there must be a trajectory that goes back to the every that was visited (see issue https://github.com/markovmodel/PyEMMA/issues/1052 and slides). This is a much stronger requirement than any other methods. (It's the price you pay for using the advantages of MSM and in my opinion it is the more correct way of doing things. If your data is not ergodic, you better should know about this.) Pyemma provides tools for checking the ergodicity of the sampling. They are integrated in TRAM.

If it is clear to you that these challenges can be addressed, you can dive more into technical information that I linked to above (issue and example).

Fabian

[fabian-paul]

In the context of metadynamics, k is the number of deposited Gaussians. If you fill the btrajs[0] array, you would need to write a loop over all time-steps t and all numbers of Gaussians k and set:

btrajs[0][t, k] = bias energy of configuration x[t] evaluated with the bias potential where k Gaussians are active

Note that this means that you have to "look into the past and the future". k is not the number of Gaussians that are active at time t but there needs to be a loop over all k from zero to the maximum number of Gaussians.

I think this information is not directly provided by Plumed and you need to reconstruct this information from the hills file and possibly from the colvars file.

Filling the btrajs array can be done using the np.cumsum function, which returns the partial sum of its input array. The line should roughly look like this (see also the example)

btrajs[0][t, :] = np.cumsum([G(m[k], sigma, colvar[0][t]) for k in range(0, K)])

where G(mean, sigma, x) is the Gaussian function and colvar[0][t] is the value of the colvars in simulation 0 and at time t.