Meta-Dynamics Simulations
In this guide, all necessary information will be given in order to perform
meta-dynamics (MTD) simulations with the xtb program.
For the theory behind our MTD approach please refer to:
S. Grimme, Exploration of Chemical Compound, Conformer, and Reaction Space with Meta-Dynamics Simulations Based on Tight-Binding Quantum Chemical Calculations, J. Chem. Theory Comput., 2019, 155, 2847-2862. DOI: 10.1021/acs.jctc.9b00143 Publication Date (Web): April 3, 2019
From a practical point of view the application of meta-dynamics is quite similar to molecular-dynamic simulations. In MTD simulations a biasing potential given as a sum of Gaussian functions is additionally expressed. The root-mean-square deviation (RMSD) in Cartesian space is chosen as a metric for the collective variables.
All adjustable parameters will be discussed and a guide to how to change them will be given as well as an example.
General MTD Setup
For any MTD calculation a Detailed Input file is necessary to enter
the correct calculation mode. The basic parameters for dynamics are taken
from the $md block as described in the section regarding Molecular Dynamics Simulations.
The $metadyn data group has to be present in the input file.
All available instructions for this data group are shown here:
key
value
description
save
integer
maximal number of structures for rmsd criteria
kpush
real
scaling factor for rmsd criteria
alp
real
width of the gaussian potential used in the rmsd criteria
coord
file
external structure file to initialize rmsd criteria
atoms
list
atoms to include in the rmsd calculation (default: all)
To avoid accidental activation of the bias potential conservative default values
are chosen in the program. So you cannot simply use a commandline-only approach
to perform a MTD calculation. First of all you want to create a metadyn.inp
file with this content
$metadyn
save=10
kpush=1.0
alp=0.2
$end
You can start the metadynamic calculation now by using the --md commandline
flag as
> xtb --md --input metadyn.inp coord
By using the flag --metadyn integer, the number of saved structures may
also be entered via the commandline and need not to be present in the
detailed input:
> xtb --metadyn 10 --input metadyn.inp coord
MTD specific Files
After the xtb program has performed the desired MTD simulation the trajectory of the structures can be found in xtb.trj.
Furthermore, files with the names scoord.* are generated. After every picosecond of simulation the structure at this point
will be written into these files. After a successful completion of the MTD simulation a xtbmdok file will be touched.
The structure and velocities at the end of the simulation will be written into a mdrestart file.
Restart
The mdrestart file can be used to restart an MTD simulation. This can be very helpful for equilibration purposes.
In order to achive this, in the $md block the restart parameter has to be set to true.
> cat xcontrol
$md
restart=true
Example/Case study
To summarize the most important topics of this chapter we will perform an MTD simulation of the water dimer molecule with xTB.
Make sure that xtb is properly set up and you have the following files in your working directory
> cat coord
$coord
-3.41057596145012 0.01397950240733 -2.48858246259228 o
-4.61882519451466 -0.37951683073135 -3.80598534178089 h
-3.93472948277574 1.59474363012965 -1.77520048297923 h
-6.45960866143682 -1.52443579539520 -6.62024481959292 o
-5.26218381176973 -1.40667057754076 -7.97781013203256 h
-6.78373759577982 -3.28799737179945 -6.34039886662289 h
$end
> cat metadyn.inp
$md
time=10
step=1
temp=298
$end
$metadyn
atoms: 1-3
save=10
kpush=0.02
alp=1.2
$end
As you can see, we will run the MTD simulation for 10 ps with a timestep of 1 fs at a temperature of 298 Kelvin. For the meta-dynamics, only the structure of the first water molecule will be taken into account in the rmsd criteria. To start the simulation we call xtb as follows
> xtb --md --input metadyn.inp coord
The output for the example MTD simulation of the water dimer will look like this:
-------------------------------------------------
| Meta Dynamics |
-------------------------------------------------
trajectories on xtb.trj or xtb.trj.<n>
MD time /ps : 10.00
dt /fs : 1.00
SCC accuracy : 1.00
temperature /K : 298.00
max steps : 10000
block length (av.) : 5000
dumpstep(trj) /fs : 100.00 100
dumpstep(coords)/fs: 1000.00 1000
H atoms mass (amu) : 2
# deg. of freedom : 14
SHAKE on. # bonds : 4 all: T
Berendsen THERMOSTAT on
kpush : 0.020
alpha : 1.200
update : 10
time (ps) <Epot> Ekin <T> T Etot
0 0.00 0.00000 0.0198 0. 0. -10.10916
est. speed in wall clock h for 100 ps : 0.01
200 0.20 -10.09118 0.0116 559. 524. -10.12881
400 0.40 -10.11436 0.0105 454. 471. -10.13041
600 0.60 -10.12260 0.0070 431. 316. -10.13157
800 0.80 -10.12671 0.0071 412. 321. -10.13081
... ... ... ... ... ... ...
4800 4.80 -10.13763 0.0084 469. 379. -10.13198
block <Epot> / <T> : -10.13978 465. drift: 0.99D+02 Tbath : 298.
5000 5.00 -10.13775 0.0082 465. 368. -10.13253
5200 5.20 -10.13783 0.0129 469. 582. -10.12808
5400 5.40 -10.13794 0.0105 471. 474. -10.13014
5600 5.60 -10.13804 0.0090 470. 407. -10.13140
... ... ... ... ... ... ...
9800 9.80 -10.13918 0.0083 462. 376. -10.13258
average properties
Epot : -10.1392169717059
Epot (accurate SCC): -10.1402473210558
Ekin : 1.019492766065306E-002
Etot : -10.1290220440452
T : 459.900938472654
thermostating problem
normal exit of md()
In the file xtb.trj we can find our trajectory.
We can analyze the structures now by displaying them in a molecular graphics editor (e.g., MOLDEN, VMD etc. )
or a trajectory analyzer (e.g. TRAVIS).
Constrained MD/MTD simulations
As you may have noticed in the example given above by checking the file xtb.trj, the water dimer dissociates within
the MTD simulation due to the applied bias potential. If you run dynamics for systems that are non-covalently bound,
you may encounter this problem from time to time. To avoid dissociation you can try to confine the simulation in a sphere by
a repulsive potential. For further details check how to confine a cavity in Detailed Input.
To avoid dissociation of the water dimer by a logfermi potential, the input file has to be modified:
> cat metadyn.inp
$md
time=10
step=1
temp=298
$end
$metadyn
atoms: 1-3
save=10
kpush=0.02
alp=1.2
$end
$wall
potential=logfermi
sphere: auto, all
$end
To start the constrained MTD simulation we call xtb as follows:
> xtb --md --input metadyn.inp coord
If you now check the trajectory file, you will see that the water molecules do not separate.
Note
The wall potential does not only work for MD/MTD simulations. It may also be applied in the same manner for single point calculations and geometry optimizations.