Metadynamics Simulation: How to define collective variables

[psriv]

Dear Admin and VASP users,

If we consider the dissociation process, say H2CO3 into CO2+H2O (J. Chem Phy 126, 204315 2007). Suppose molecules is

O2
|
H6--O4--C1--O3--H5

In the paper they took C–O and O–H coordination numbers as CV. Since VASP manual does not explain how to set CVs in ICONST file very clearly it would be really good if someone can explain:

1) How to set the coordination number of, say C-O or C-H bond in ICONST file?
2) Or if I want to take the distance between O2 and center of bond (CB) between C1--O3 atoms. Then how ICONST will look like?
3) How to define q_i and c_i in general in ICONST file as given on http://cms.mpi.univie.ac.at/vasp/vasp/I ... sec_iconst

TIA,

Best regards,
PS

[admin]

This is too much to be explained in one message.
Let me give you a simple example in terms of interatomic distances how the collective variable (CV) works.
In the mentioned paper a conversion is considered from the configuration in Fig. 1b to that in Fig. 1d.
The change of bonding is the following. H5 goes away from O4 and connects to O3, and O3 goes away from C1.
During the reaction the CV should change continuously, should be smooth (no jumps) and should include
all changes of the bonding. The following definition fulfils the criteria:
R 4 5 0
R 1 3 0
R 3 5 0
S 1. 1. -1. 0
When running a MD all atoms are free to move with the constraint CV=r1+r2-r3.
For starting geometry (SG) and the geometry of the product (GP) you will get the variation
of the gradient averaged out to zero (gradient: grep b_m REPORT; the value of CV: grep cc REPORT).
For any configuration with the CV between SG and PG you will get a point for sampling
the free-energy barrier. 8-10 points are usually enough for the integration.
Note that this is example to explain how the CV works, not an example of the metadynamics.

Here is a quick test if the coordinate is good to describe the dissociation H2CO3 --> H2O + CO2 .
Putting H2CO3 into a box 8x9x10A, using just Gamma point and running the slow growth MD
(J. Phys. Chem. B 1997, 101, 7877) at 300K, with 300 eV cutoff, the timestep = 0.5 fs and a growth of
the coordinate by 0.001 A in each step. The calculation finished in app. 1 hour (16 CPUs) and
the result is the following:

http://snag.gy/HcMAE.jpg

The curve of gradients starts and ends at the zero line. Zero-gradient points are defined
by minima at the line of integrated gradients (red line). The transition state is defined
by the maximum at the red line. The change of the free energy of the dissociation (integrated
between the starting configuration and the transition state) is 137.2 kJ/mol, nicely comparable
with 155.2 kJ/mol from the metadynamics (JCP 2007, 126, 204315). At increased precision
(2x2x2 kpoints, finished in 2.8 hours) the barrier is 141.0 kJ/mol. For the reversed reaction
the integration between the final coordination and the TS gives 164.3 kJ/mol. Both
the correct shape of the calculated curve of gradients and the reasonable value of the activation
energy show that the collective variable r1+r2-r3 is correct.

[psriv]

Dear VASP Admin,

Thank you so much for explaining it so nicely. I have slightly different case I am studying exfoliation of a 2D material from its 3D counterpart. Since the sheet has some defects so there are some bonds formed between the two layers. I want to study the dissociation of those particular bonds. In my ICONST file I took (say 1-5 are in 1st layer and 6-10 atoms are in 2nd layer).
R 1 10 5
R 2 9 5
R 3 8 5
R 4 7 5

I am using following metadynamics tags:
MDALGO = 21
HILLS_BIN = 10
HILLS_H = 0.01
HILLS_W = 0.01
STATUS = 5

I apply bias potential only on those particular bonds. Simply MD shows dissociation of those bonds, but if I do metadynamics simulation, some bond shrink and some expand. I was expecting metadynamics to run much faster than simple MD but in my case it's not. It would really a great help for me, if you can suggest and comment whether my choice for CV is right or not. Or if there can be better CV's (I checked with fixing the bond-lengths of in-plane bonds but that also didn't work).

Thanking you for your suggestions.

Best regards,
PS

[admin]

To understand the reaction path one needs to see the starting configuration
and the configuration of the product.

[mhl980]

Dear admin,

As I know the activation energy (energy barrier) does not depends on temperature. Is that correct?.
If the activation energy (barrier energy) does not depends on temperature do we necessary to perform the simulation at 300K ?.

Regards

[gnayak]

This is too much to be explained in one message.
Let me give you a simple example in terms of interatomic distances how the collective variable (CV) works.
In the mentioned paper a conversion is considered from the configuration in Fig. 1b to that in Fig. 1d.
The change of bonding is the following. H5 goes away from O4 and connects to O3, and O3 goes away from C1.
During the reaction the CV should change continuously, should be smooth (no jumps) and should include
all changes of the bonding. The following definition fulfils the criteria:
R 4 5 0
R 1 3 0
R 3 5 0
S 1. 1. -1. 0
When running a MD all atoms are free to move with the constraint CV=r1+r2-r3.
For starting geometry (SG) and the geometry of the product (GP) you will get the variation
of the gradient averaged out to zero (gradient: grep b_m REPORT; the value of CV: grep cc REPORT).
For any configuration with the CV between SG and PG you will get a point for sampling
the free-energy barrier. 8-10 points are usually enough for the integration.
Note that this is example to explain how the CV works, not an example of the metadynamics.

Here is a quick test if the coordinate is good to describe the dissociation H2CO3 --> H2O + CO2 .
Putting H2CO3 into a box 8x9x10A, using just Gamma point and running the slow growth MD
(J. Phys. Chem. B 1997, 101, 7877) at 300K, with 300 eV cutoff, the timestep = 0.5 fs and a growth of
the coordinate by 0.001 A in each step. The calculation finished in app. 1 hour (16 CPUs) and
the result is the following:

http://snag.gy/HcMAE.jpg

The curve of gradients starts and ends at the zero line. Zero-gradient points are defined
by minima at the line of integrated gradients (red line). The transition state is defined
by the maximum at the red line. The change of the free energy of the dissociation (integrated
between the starting configuration and the transition state) is 137.2 kJ/mol, nicely comparable
with 155.2 kJ/mol from the metadynamics (JCP 2007, 126, 204315). At increased precision
(2x2x2 kpoints, finished in 2.8 hours) the barrier is 141.0 kJ/mol. For the reversed reaction
the integration between the final coordination and the TS gives 164.3 kJ/mol. Both
the correct shape of the calculated curve of gradients and the reasonable value of the activation
energy show that the collective variable r1+r2-r3 is correct.

Dear VASP Admin,
This post is really helpful.

I am trying to migrate of impurity atom from the initial to final position and trying to look at the free-energy barrier.

In the following case, I am not able to assign my CV in terms of cartesian coordinates in ICONST file. There is no example available "how to assign the cartesian coordinates".

If you could help regarding this, will be appreciated.

Thank you and regards,
Ganesh K Nayak

[ankitma]

I am trying to do something similar with an aim to use cartasian coordinates as CV but not able to get it write. Is there an example?

[fabien_tran1]

Hi,

Can you please upload the files INCAR, POSCAR and ICONST to see what you have tried. And also explain why it does not work.

[fabien_tran1]

Have you tried the flags cX, cY and cZ, as explained at the page describing the ICONST file?