Nuclephile Substitution CH3Cl - mMD2

Overview >Liquid Si - Standard MD > Liquid Si - Freezing > Nucleophile Substitution CH3Cl - Standard MD > Nuclephile Substitution CH3Cl - mMD1 > Nuclephile Substitution CH3Cl - mMD2 > Nuclephile Substitution CH3Cl - mMD3 > Nuclephile Substitution CH3Cl - SG > Nuclephile Substitution CH3Cl - BM > List of tutorials

Task

In this example the nucleophile substitution of a Cl^- by another Cl^- in CH₃Cl via meta dynamics is correctly simulated by using extra "repulsive potential walls."

Input

POSCAR

   1.00000000000000
     12.0000000000000000    0.0000000000000000    0.0000000000000000
     0.0000000000000000    12.0000000000000000    0.0000000000000000
     0.0000000000000000    0.0000000000000000    12.0000000000000000
   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

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.

INCAR

PREC=Low
EDIFF=1e-6
LWAVE=.FALSE.
LCHARG=.FALSE.
NELECT=22
NELMIN=4
LREAL=.FALSE.
ALGO=VeryFast
ISMEAR=-1
SIGMA=0.0516

############################# MD setting #####################################
IBRION=0                                           # MD simulation
NSW=1000                                           # number of steps
POTIM=1                                            # integration step
TEBEG=600                                          # simulation temperature
MDALGO=11                                          # metaDynamics with Andersen thermostat
ANDERSEN_PROB=0.10                                 # collision probability
HILLS_BIN=50                                       # update the time-dependent bias
                                                           # potential every 50 steps
HILLS_H=0.005                                      # height of the Gaussian
HILLS_W=0.05                                       # width of the Gaussian
##############################################################################

• Same INCAR file as in the previous example (Nuclephile Substitution CH3Cl - mMD).

ICONST

R 1 5 0
R 1 6 0
S 1 -1 5

Calculation

In principle, 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 (which is linked with a lower barrier than the SN2 reaction). In order to avoid this undesired process, an extra bias potential ("repulsive walls") is used whose role is to restrict our sampling to a relevant region (approx. $-3 \AA <$ collective variable $< 3 A$). In fact, the positions of walls can be chosen arbitrarily - we only require that the region between the walls contains all the information we are interested in (in this case we want to see free-energy minima for both "reactant" and "product" as well as the "transition state"). In order for the walls to be effective, we also require that they are significantly higher than the expected reaction barrier (otherwise the likelihood to cross the wall during meta dynamics would be higher than that for the barrier). From the potential energy profile (static calculations not reported here) we obtained a reasonable guess for the reaction barrier - it is about 0.4 eV - hence the height for the wall of 1 eV should be sufficient.

For practical reasons, we split our (presumably long) meta dynamics calculation into shorter runs of length of 1000 fs (NSW=1000; POTIM=1). At the end of each run the HILLSPOT file (containing bias potential generated in previous run) must be copied to the PENALTYPOT fiel (the input file with bias potential) - this is done automatically in the script run which looks as follows:

#!/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

Please run this script by typing:

bash ./run

• Importantly the user has to adjust the runvasp variable, which holds the command for the executable command.
• The calculation is run for 50 ps. For convenience the calculation is split into 50 calculations each of length 1 ps.

Download

Overview >Liquid Si - Standard MD > Liquid Si - Freezing > Nucleophile Substitution CH3Cl - Standard MD > Nuclephile Substitution CH3Cl - mMD1 > Nuclephile Substitution CH3Cl - mMD2 > Nuclephile Substitution CH3Cl - mMD3 > Nuclephile Substitution CH3Cl - SG > Nuclephile Substitution CH3Cl - BM > List of tutorials