My Community

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

Hi,

Could you please also provide the INCAR and OUTCAR files?

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

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.

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

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.

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

Dear mo_salha,

I checked your POSCAR file with the VESTA viewer. The output of VASP seems to be correct.
In the appended figure I colored the Mg atom (atom number 1469) in pink and the oxygen atom (atom number 1341) in blue.
Without taking into account periodic boundary conditions I get a distance of roughly 15 Angstroem. When taking into account periodic boundary conditions I get a distance of 11.11Angstroem. See the bottom of the figure for the exact distance value. To me it seems like you are selecting the wrong atom numbers and therefore get unexpected distances. I would recommend using VESTA for atom selections in the ICONST file since it definitely uses the same atom ordering as VASP.

All the best Jonathan

Hi again guys

Metadynamics seems to be working well so far, yet almost too well
I'm interested in the Mg-O distance but only up to around 6 Angstroms, however bias keeps being added to high energy regions of the Mg-O distance, leading to a deep exploration of distances > 10 Ang. (See attached plot)

I'm assuming this is happening because I havn't specified any constraints on the maximum Mg-O distance and I'm getting runaway bias
Any tips on controlling this?

Best,
Mo

Dear mo_salha,

In VASP it is not possible to set a maximum bond length. But you could run VASP from within the ASE environment in python3.
There you could use the constrain module to apply a maximum bond length. Check this site for more information:
https://wiki.fysik.dtu.dk/ase/ase/constraints.html

All the best Jonathan

Thanks Johnathan I'll give it a go

Mo

Hi again on a seperate note I managed to compile lammps with plumed and vasp support using this site
https://www.vasp.at/wiki/index.php/Runn ... _in_LAMMPS

Everything seems to be in working order except that my lammps output seems not to be recognising my ML_FF after I specified:
units metal
dimension 3
boundary p p p
atom_style atomic

# Read LAMMPS data file from ASE conversion
read_data geometry.data

# MLFF potential (adjust to your MLFF library/path)
pair_style vasp
pair_coeff * * ML_FF H O Mg Cl

# Setup neighbor and timestep
neighbor 2.0 bin
neigh_modify delay 0 every 1 check yes

timestep 0.001 # 1 fs

# Enable PLUMED
fix plumed all plumed plumedfile plumed.dat

# Thermo output
thermo 100
thermo_style custom step temp pe etotal press

# Run
run 50000

It's normal for my ML_FF made in vasp to be a binary right?
It seems lammps is expecting some text file format

Mo

Attached is the error:

Loading intel-mkl/2024.1.0
  Loading requirement: intel-tbb/2021.9.0-gcc-12.2.0
LAMMPS (19 Nov 2024 - Development - )
Reading data file ...
  triclinic box = (0 0 0) to (23.657023 22.85021 31.353957) with tilt (10.535818 30.265145 -1.3143321)
WARNING: Triclinic box skew is large. LAMMPS will run inefficiently. (../domain.cpp:221)
  1 by 2 by 2 MPI processor grid
  reading atoms ...
  1488 atoms
  read_data CPU = 0.036 seconds
pair_coeff vasp:
  VASP ML force field file   = ML_FF
  Setting up force field .....
EEEEEEE  RRRRRR   RRRRRR   OOOOOOO  RRRRRR      ###     ###     ###
E        R     R  R     R  O     O  R     R     ###     ###     ###
E        R     R  R     R  O     O  R     R     ###     ###     ###
EEEEE    RRRRRR   RRRRRR   O     O  RRRRRR       #       #       #
E        R   R    R   R    O     O  R   R
E        R    R   R    R   O     O  R    R      ###     ###     ###
EEEEEEE  R     R  R     R  OOOOOOO  R     R     ###     ###     ###


EEEEEEE  RRRRRR   RRRRRR   OOOOOOO  RRRRRR      ###     ###     ###
E        R     R  R     R  O     O  R     R     ###     ###     ###
E        R     R  R     R  O     O  R     R     ###     ###     ###
EEEEE    RRRRRR   RRRRRR   O     O  RRRRRR       #       #       #
E        R   R    R   R    O     O  R   R
E        R    R   R    R   O     O  R    R      ###     ###     ###
EEEEEEE  R     R  R     R  OOOOOOO  R     R     ###     ###     ###

ERROR: Unknown file format, input file does not seem to contain a valid header.

  ---->  I REFUSE TO CONTINUE WITH THIS SICK JOB ... BYE!!! <----