Metadynamics does not work with small time step

[susilehtola]

The biasing potential in metadynamics is controlled by the MD step number. This means that the biasing potential gets switched on instantaneously when a small time step is used, resulting in the system exploding. For instance, when running systems with hydrogens with shake=0 and the usual 0.5 fs time step, one experiences an explosion already in the initial dynamics run since the starting configuration is included in the metadynamics potential.

The correct way to switch on the biasing potential is to do it as a function of actual time.

[susilehtola]

We encountered this issue in metadynamics on the vanillin molecule with the coordinates

$coord
       -3.14393748772232        3.33813352345031       -0.00014119336956      o
        7.86002766715000       -1.65701448465138        0.00040621434150      o
        5.74049590710396       -2.49527551648547       -0.00050479923313      c
        3.42304695453588       -0.98774709392443       -0.00010485496334      c
        1.06591048798478       -2.15083142745214       -0.00014844660978      c
       -1.17270680517221       -0.77415535108070       -0.00013844779300      c
       -3.30174527376127       -2.20634525593637       -0.00025807238808      o
       -5.72053722271011       -1.08375306666073       -0.00010984349014      c
       -1.01555970570177        1.88182884434437       -0.00005730429937      c
        1.35264373756029        3.03267846294878        0.00013226089640      c
        3.55184020306940        1.64453015125235        0.00021056287264      c
       -2.66751301607659        5.09916430174068       -0.00002324857945      h
        5.40066653268587       -4.56215905836274        0.00078445601463      h
        0.91214108164020       -4.18420996464937       -0.00021662002928      h
       -7.04162340875582       -2.66213089601776       -0.00034302738473      h
       -6.02711939534645        0.06913218178066        1.68326359513705      h
       -6.02710355274802        0.06966531366495       -1.68313322383375      h
        1.43463686418473        5.07677471557364        0.00015134556396      h
        5.37644021153172        2.55171839991760        0.00024009577804      h
$end

and the input

$md
   temp=300
   time=10.0
   dump=5.0
   step=0.5
   shake=0
$end
$metadyn
   save=10
   kpush=1.0
   alp=0.2
$end

The result does not appear to be fully reproducible, though. I wonder if another issue is that atoms can be placed on top of each other in the random initialization.

[susilehtola]

For instance, running

 xtb --omd --input metadyn.inp coord 

on my laptop with xtb-6.5.1-1.fc37.x86_64 I got the following result

         time (ps)    <Epot>      Ekin   <T>   T     Etot
      0    0.00      0.00000   0.0613    0.    0.   -33.29787
est. speed in wall clock h for 100 ps :  0.36
    200    0.10    -32.60311   0.3225 1967. 3573.   -32.42622
    400    0.20    -32.71420   0.3797 3016. 4207.   -32.48540
    600    0.30    -26.03888   0.3612 3399. 4002.     0.36117
    800    0.40    -19.53728   0.3006 3463. 3330.     0.30058
   1000    0.50    -15.63373   0.2510 3379. 2781.     0.25098
   1200    0.60    -13.03028   0.2104 3241. 2331.     0.21038
   1400    0.70    -11.17014   0.1771 3084. 1963.     0.17714
   1600    0.80     -9.77474   0.1499 2924. 1661.     0.14992
   1800    0.90     -8.68926   0.1276 2770. 1414.     0.12764
   2000    1.00     -7.82077   0.1094 2624. 1212.     0.10941

after which xtb crashed. However, other runs worked fine.