View source for Nuclephile Substitution CH3Cl - mMD1

{{Template:Molecular dynamics - Tutorial}}

== Task ==
In this example a nucleophile substitution of a Cl<sup>-</sup> by another Cl<sup>-</sup> in CH<sub>3</sub>Cl is attempted via a meta dynamics calculation.

== Input ==

=== {{TAG|POSCAR}} ===

    1.00000000000000
      12.0000000000000000    0.0000000000000000    0.0000000000000000
       0.0000000000000000   12.0000000000000000    0.0000000000000000
       0.0000000000000000    0.0000000000000000   12.0000000000000000
 C H Cl
    1   3   2
 cart
          5.91331371  7.11364924  5.78037960
          5.81982231  8.15982106  5.46969017
          4.92222130  6.65954232  5.88978969
          6.47810398  7.03808479  6.71586385
          4.32824726  8.75151396  7.80743202
          6.84157897  6.18713289  4.46842049
*The starting {{TAG|POSCAR}} file for this example can be found under POSCAR.init. It will be needed for the script that runs the job (run.sh).
*A sufficiently large cell is chosen to minimize the interactions between neighbouring cells and hence to simulate an isolated molecular reaction.

=== {{TAG|KPOINTS}} ===
 Automatic
  0
 Gamma
  1  1  1
  0. 0. 0.
*For isolated atoms and molecules interactions between periodic images are negligible (in sufficiently large cells) hence no Brillouin zone sampling is necessary.

=== {{TAG|INCAR}} ===
 {{TAGBL|PREC}}=Low
 {{TAGBL|EDIFF}}=1e-6
 {{TAGBL|LWAVE}}=.FALSE.
 {{TAGBL|LCHARG}}=.FALSE.
 {{TAGBL|NELECT}}=22
 {{TAGBL|NELMIN}}=4
 {{TAGBL|LREAL}}=.FALSE.
 {{TAGBL|ALGO}}=VeryFast
 {{TAGBL|ISMEAR}}=-1
 {{TAGBL|SIGMA}}=0.0516
 
 ############################# MD setting #####################################
 {{TAGBL|IBRION}}=0                                           # MD simulation
 {{TAGBL|NSW}}=50000                                          # number of steps
 {{TAGBL|POTIM}}=1                                            # integration step
 {{TAGBL|TEBEG}}=600                                          # simulation temperature
 {{TAGBL|MDALGO}}=11                                          # metaDynamics with Andersen thermostat
 {{TAGBL|ANDERSEN_PROB}}=0.10                                 # collision probability
 {{TAGBL|HILLS_BIN}}=50                                       # update the time-dependent bias
                                                              # potential every 50 steps
 {{TAGBL|HILLS_H}}=0.005                                      # height of the Gaussian
 {{TAGBL|HILLS_W}}=0.05                                       # width of the Gaussian
 ##############################################################################


*The {{TAG|INCAR}} file in this example is the same as in the previous example ({{TAG|Nucleophile Substitution CH3Cl - Standard MD}}) with the exception of the metadynamics tags {{TAG|HILLS_H}} and {{TAG|HILLS_W}}.
*Metadynamics molecular dynamics is formally exact in the limit of infinitesimally small hills ({{TAG|HILLS_H}}) and infinite update time ({{TAG|HILLS_BIN}}) for the time-dependent bias potential, hence the parameter {{TAG||HILLS_H}} should be as small as possible while {{TAG|HILLS_BIN}} should be as large as possible.
*A random seed for reproducibility is used in this calculations which is added later by the ''run'' script.
*Meta dynamics is formally exact in the limit of infinitesimally small hills and infinite update times for the time-dependent bias potential. Hence the parameter {{TAG|HILLS_H}} should be as small as possible while {{TAG|HILLS_BIN}} should be as large as possible.

=== {{TAG|ICONST}} ===
For this example an {{TAG|ICONST}} file is used which looks like:
 R 1 5 0
 R 1 6 0
 S 1 -1 5
*This file is the same as in the previous example {{TAG|Nucleophile Substitution CH3Cl - Standard MD}} with the exception that the 5 at the fourth entry in the third row specifies that a bias potential is applied to the special coordinate.

== Calculation ==



=== Running the calculation ===
The mass for hydrogen in this example is set 3.016 a.u. corresponding to the tritium isotope. This way larger timesteps can be chosen for the MD. 
For practical reasons, we split our (pressumably long) meta dynamics calculation into shorter runs of length of 1000 fs for each ({{TAG|NSW}}=1000 and {{TAG|POTIM}}=1). At the end of each run, the {{TAG|HILLSPOT}} file (containing the bias potential generated in the previous run) must be copied to the {{TAG|PENALTYPOT}} file (the input file with bias potential). Meta dynamics is a fully stochastic simulation, the results can depend on the random numbers used in the calculations. To ensure best reproducibility the calculation is started from the same random seed (''rseed="RANDOM_SEED =         311137787                0                0"''). After each step the last random seed is copied to the {{TAG|INCAR}} file for the new step. This way reproducability is ensured between each step.

The calculation is executed using the "run" script:
 #!/bin/bash
 
 runvasp="mpirun -np x executable_path"
 
 # ensure that this sequence of MD runs is reproducible
 cp POSCAR.init POSCAR
 cp INCAR.init INCAR
 rseed="RANDOM_SEED =         311137787                0                0"
 echo $rseed >> INCAR
 
 
 i=1
 while [ $i -le 50 ] 
 do
 
   # start vasp
   $runvasp
 
   # ensure that this sequence of MD runs is reproducible
   rseed=$(grep RANDOM_SEED REPORT |tail -1)
   cp INCAR.init INCAR
   echo $rseed >> INCAR
 
   # use bias potential generated in previous mMD run
   cp HILLSPOT PENALTYPOT
 
   # use the last configuration generated in the previous
   # run as initial configuration for the next run
   cp CONTCAR POSCAR
 
   # backup some important files
   cp REPORT REPORT.$i
   cp vasprun.xml vasprun.xml.$i
 
   let i=i+1
 done
*The user has to adjust the ''runvasp'' variable, which holds the command for the executable command. 
*Please run this script by typing:
 bash ./run

===Expectation===
The meta dynamics simulation pushes the system against the reaction barrier. The amplitude of oscillation of the collective variable increases (as a larger and larger region of configuration space is visited) and at some point the collective variable switches from a positive value (corresponding to the reactant) to a negative value (corresponding to the product).


===Reality===
Most likely the collective variable increases to unexpectedly large values (<math>>3 \AA</math>) while the transition region (around <math>0 \AA</math>) is either not sampled or it is smpled only in the final part of the simulation. In other words, the simulation spent most of its time sampling an uninteresting part of the configuration space related to the dissociation of the vdW complex instead pushing the configuration over the barrier. This is because meta dynamics always seeks for the path of least resistance. In the case of our model system this corresponds to the dissociation of the vdW complex because the barrier this process is much lower than that for the SN2 reaction.

To verify this the time evolution of the collective variable is monitored (and written to timeEvol.dat) by calling the command:
 bash ../timeEv.sh
The time evolution is visualized using the command:
 gnuplot -e "set terminal jpeg; set xlabel 'Timestep'; set ylabel 'CV (Ang)'; set style data lines; plot 'timeEvol.dat'" > timeEvol.jpg
It should look like the following:

[[File:TimeEvol mMD1 CH3Cl.jpg|300px]]

We indeed see that the transition region is never sampled and that the meta dynamics simulation takes the path of least resistance, which is the dissociation of the vdW complex.


===Solution===
In order to restrict the sampling into the part of configuration space that is relevant for the process of interest (say between <math>-2 \AA</math>  and <math>2 \AA</math>) we must use some trick which is explained in the next example.

== Download == 

[http://www.vasp.at/vasp-workshop/examples/CH3Cl_mMD1.tgz CH3Cl_mMD1.tgz]

{{Template:Molecular dynamics - Tutorial}}

[[Category:Examples]]