unable to analyse the metadynamics output file of the cp2k

[Deewan Singh Teja]

Hey Good afternoon everyone

Can anybody known how to analysis the cp2k metadynamics file? Please help if possible

Thank you

Deewan

[Marcella Iannuzzi]

Please refer to some tutorial like

https://www.cp2k.org/exercises:2015_cecam_tutorial:mtd1

Regards

Marcella

[Karl Irikura]

I have the same question.  That tutorial does not explain how to extract free energy information from the output files.  

Thanks for any tips,

Karl

[David M. Sherman]

You're right: the tutorials don't see to explain this bit at all...

You need to use the "graph" program which comes with cp2k. When you do a metadynamics run, the resulting *.restart file will contain a record of all of the HILLS  added during the metadynamics simulation. Note that if you do a multiple-walker run, the *.restart file in any one of the walker directories will also have all the hills from all the other walkers. Let's assume you did a metadynamics simulation with 2 COLVARS and this generated a restart file named "mysim.restart".  To turn that record of HILLS in the restart file into a one-dimensional free-energy surface with respect to the first COLVAR, you run the program graph.psmp using the following syntax:

graph.psmp -cp2k -ndim 2 -ndw 1 -file mysim.restart

If you want  a one-dimensional free energy surface with respect to the second COLVAR, you run

graph.psmp -cp2k -ndim 2 -ndw 2  -file mysim.restart

(Note: BE VERY CAREFUL HERE!  Even though you only want a one dimensional fes, the argument to -ndim MUST be the actual number of COLVARS that were used in the metadynamics run.  The argument to -ndw is which COLVAR(S) you want to generate the free-energy surface with respect to.  If you accidently use -ndim 1, you will get what looks like a perfectly sensible fes, but it will be completely erroneous).

Running the graph program will generate a file called fes.dat that you can plot with gnuplot to see the free energy surface.  Here's an example of a gnuplot script to

generate the plot (note: I'm converting my Cl coordination number COLVAR into a sensible number by multiplying it by the number of Cl atoms in the cell  and converting the free energy from AU to kJ/mole).

set xlabel "Cl CN"
#set terminal svg
set xrange[-1:4]
set yrange[-230:-120]
plot "fes.dat" using ($1*4):($2*2625.5) with lines
pause -1 "Hit any key to continue"

Now,  if you want a 2-dimensional plot of the the free energy surface (i.e., over both

COLVARS) then you run:

graph.psmp -cp2k -ndim 2 -ndw 1 2 -file mysim.restart

and plot that using gnuplot.  Here's a gnuplot input script to generate a contour map:

set view map
unset surface
set contour
set cntrparam levels 30
set dgrid3d 100,100,1
set xlabel "Cl CN"
set ylabel "O CN"
set style line 1 lc rgb "black" lw 1
set style increment user
#set terminal svg
splot [-1:4] [-1:5] "fes.dat" using ($1*4.):($2*55):($3*2625.5) with lines
pause -1 "Hit any key to continue"

Have fun!

[Karl Irikura]

Thanks, David, the "graph" program is working for me.

What are the energy units in the file "fes.dat"?  The well depths are < 0.1 units.  

More generally, is it better to run one series of simulations for each collective variable, or a single series with two collective variables?  In my case there are two titratable groups and I'm trying to calculate pKa. 

Thanks,

Karl