graph.popt tool for metadynamics

[Laura Bonometti]

Hi everyone!

I am using CP2K for running metadynamics simulations. I used the graph.popt tool for obtaining the free energy of my system, but actually I can't understand how it is calculated. I read that the final energy results are in hartree but if I convert them in eV they do not make any sense in terms of order of magnitude, even if the obtained curve is in agreement with what I expected. I know that the free energy should be calculated considering the probability distribution of the CVs but if I try to reconstruct the curve in this way, I don't obtain the same results obtained by exploiting the graph.popt tool. I also tried to use the python script found at  https://www.cp2k.org/exercises:2021_uzh_acpc2:ex03 and the range of free energy is reasonable in this case, but completely different respect to what obtained by using graph.popt, so I think that I'm not understanding what is printed by this program.

Can you help me to solve this question?

Kind regards,

Laura

[Marcella Iannuzzi]

Dear Laura,

Difficult to guess what is going wrong with graph.popt. Maybe there is some normalisation factor that is wrongly set for the calculation of a multi-dimensional FES.

Maybe more information would help.

The python script you are indicating, though, it is not generating a FES from a metadynamics  potential. It builds a probability histogram from a standard molecular dynamics sampling.

This cannot be directly applied to an enhanced sampling, because you are not taking into account the bias potential in collecting the probability distribution. 

Regards

Marcella

[Laura Bonometti]

Dear Marcella,

Thank you for your help!

I attach here the restart file used for obtaining the FES from graph.popt if this can help.

For running the calculation with graph.popt I wrote: graph.popt -cp2k -ndim 2 -ndw 1 2 -file metadyn_Ni_6Li_new-1_59500.restart -out fes.dat .

Do I have to add some other information?

Maybe I let the calculation run for too long, depositing too much Gaussians, so maybe I have to stop the calculation/consider another restart file related to a lower number of computed steps, obviously considering that the particle reached the expected position. 

Regards,

Laura

[Marcella Iannuzzi]

Dear Laura, 

The run you sent me is done for two collective variables, two distances, that behave in a very similar way (isn't this choice redundant). 

The oscillations of these two variables are quite large, mostly between 1 and 8 Å.

 I would say that the changes are too fast due to too large hills. But I cannot be sure because I do not know what these distances are.

The height of the hills used for this run is ~2.7 and they are spawned every 30 fs. 

It is possible that the space has not been properly sampled and the FES resulting from the analysis done by graph is not reliable. 

This is obviously only a guess, since I do not know the system and the type of topography that you are expecting. 

As a rule of thumb, the height of the hills should be one order of magnitude smaller than the smaller energy barrier that you want to resolve.

If you are aiming at barriers of 27 eV, 2.7 eV is OK

Regards

Marcella