How to implement Metadynamics

[mo_salha]

Hi all
I'm trying to follow the example given here ( wiki/index.php/Nuclephile_Substitution_CH3Cl_-_SG ) to do some metadynamics runs with a machine learning forcefield we developed. The manual assures me that the langevin thermostat is compatible with metadynamics but I find very little documentation on implementing it. Attached is a .zip file containing my inputs as I try to replicate the slow growth approach seen in the linked example.

It's stated that "The values of all collective variables for each MD step are listed in REPORT-file, check the lines after the string
Metadynamics" yet I find no string related to metadynamics in my REPORT output, nor in the OUTCAR

========================================
MD step No. 1
========================================

SHAKE converged in 12 steps

>Const_coord
cc> S -2.93121 -2.93121 -0.389181E-05

>Blue_moon
lambda |z|^(-1/2) GkT |z|^(-1/2)*(lambda+GkT)
b_m> 0.244390E+01 0.172014E+01 0.589344E-03 0.420488E+01

>Energies
E_tot E_pot E_kin EPS ES
e_b> -0.26417879E+02 -0.26930592E+02 0.53219318E+00 0.00000000E+00 0.00000000E+00

>Temperature
T_sim T_inst
tmprt> 600.000 882.258

>Thermostat, num. of collisions: 1

RANDOM_SEED = 311137787 352078 2

========================================
MD step No. 2
========================================
(etc.)

Any idea if I'm missing any lines that might 'activate' metadynamics?

Thanks
Mo

[fabien_tran1]

Hi,

Could you please also provide the INCAR and OUTCAR files?

[mo_salha]

Hi
I've attached the INCAR but could only show a small part of the OUTCAR since it's too big to attach

Hope this helps

[fabien_tran1]

Hi,

In ICONST, the collective variables are those whose status (last value at each line) is 5. So, you need to replace 0 by 5 where appropriate. Then, the values will be printed in REPORT.

[mo_salha]

Thanks Fabian I got it working!

I adjusted my ICONST file to try and have the collective variable set to the Mg O distance in my model, which are atoms 3 and 2 in my POSCAR

ICONST:
R 2 3 5
~
~
~

Is there a way I can view this collective variable changing with the energy? I tried looking into the HILLSPOT file but the values for Mg-O distance should be > 2A and not the ~1.57 I keep seeing in the HILLSPOT.

How am I supposed to read the HILLSPOT file:
1.57080 1.57080 1.57080 0.00004 0.00003 -0.00001 0.00500 0.05000
1.57080 1.57079 1.57080 0.00014 0.00009 -0.00004 0.00500 0.05000
1.57079 1.57079 1.57080 0.00024 0.00014 -0.00009 0.00500 0.05000
1.57078 1.57079 1.57080 0.00035 0.00022 -0.00013 0.00500 0.05000
1.57078 1.57078 1.57080 0.00048 0.00035 -0.00013 0.00500 0.05000
1.57077 1.57078 1.57080 0.00060 0.00054 -0.00006 0.00500 0.05000
1.57076 1.57077 1.57080 0.00073 0.00078 0.00005 0.00500 0.05000
1.57075 1.57077 1.57080 0.00084 0.00101 0.00017 0.00500 0.05000
1.57074 1.57076 1.57080 0.00095 0.00126 0.00030 0.00500 0.05000
1.57073 1.57076 1.57080 0.00107 0.00158 0.00051 0.00500 0.05000
etc etc
~
~

Many thanks for your help
Mo

[fabien_tran1]

Hi,

I got help from a colleague (Jonathan Lahnsteiner) specialist in this topic. His answer is:

The HILLSPOT file has the same format as the PENALTIPOT file. This is described on these two wiki pages:
https://www.vasp.at/wiki/index.php/HILLSPOT
https://www.vasp.at/wiki/index.php/PENALTYPOT
So for the last line of the example which the user gives this would be
1.57073 1.57076 1.57080 0.00107 0.00158 0.00051 0.00500 0.05000
centers of gaussian
x1= 1.57073, x2 = 1.57076, x3=1.57080, x4=0.00107, x5=0.00158, x6=0.00051 and
h=0.00500 the height oof the gaussian (h*exp(-x2/w)) and the width = 0.00051.

The user get's such a line for every time step in the MD simulations. To see the total bias potential along a collective variable one has to add up all the lines from the HILLSPOT file.
For example if the user has a ICONST file for a single bond bwteen atom 1 and 4 it would look like
R 1 4 5 (where 5 tells vasp use bias potential).

Then for 3 time steps the HILLSPOT file could look like
1.567 0.00500 0.05000
1.555 0.00500 0.05000
1.543 0.00500 0.05000

To visualize the total bias potential the user would have to add
0.05exp(-(1.567-x)2/0.0052)+0.05exp(-(1.555-x)2/0.0052)+0.05*exp(-(1.543-x)2/0.0052)
Then he could plot the bias potential for this collective variable. If there are more collective variables add the gaussians up for every of those.

[mo_salha]

Amazing thanks again to you both.
I'm now able to plot the data from the HILLSPOT and corrected it by changing my ICONST to:
R 1341 1469 5

Where atom 1341 is an O and atom 1469 is Mg. Found these two atoms using ovito and the distance between them is about ~2A
Why then are my HILLSPOT and the fic_p in my REPORT file, printing values for the collective variable in the range of ~11?
HILLSPOT:
11.12033 0.00500 0.05000
11.12436 0.00500 0.05000
11.12915 0.00500 0.05000
11.13882 0.00500 0.05000
11.15184 0.00500 0.05000

REPORT:
========================================
MD step No. 1
========================================

Atomic velocities initialized by STEP_tb

Lattice velocities initialized by STEP_tb

p_b> 0.000 -3.275 3.516 0.241

t_b> 300.000 300.895 348.208 300.990

>Energies
E_tot E_pot E_kin EPS ES
e_b> -0.72636906E+04 -0.73216998E+04 0.58009204E+02 0.00000000E+00 0.00000000E+00

>Metadynamics
fic_p> 11.11752

RANDOM_SEED = 311137787 17892 0

If the collective variable is the Mg-O distance between two atoms that are coordinated they cant be 11A away from eachother?
I can't imagine OVITO is listing the atoms differently to how they appear in the POSCAR (which I've attached below)

Thanks again for all your help on this

Mo