Nuclephile Substitution CH3Cl - mMD1

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 a nucleophile substitution of a Cl^- by another Cl^- in CH₃Cl is attempted via a meta dynamics calculation.

Input

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 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.

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=50000                                           # 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
##############################################################################

• The INCAR file in this example is the same as in the previous example (Nucleophile Substitution CH3Cl - Standard MD) with the exception of the metadynamics tags HILLS_H and HILLS_W.
• Metadynamics molecular dynamics is formally exact in the limit of infinitesimally small hills (HILLS_H) and infinite update time (HILLS_BIN) for the time-dependent bias potential, hence the parameter [[]] should be as small as possible while HILLS_BIN should be as large as possible.

ICONST

For this example an 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 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

The calculation is executed in one long step. Just simply run the calculation by running your VASP executable.

The time evolution of the collective variable is monitored (and written to timeEvol.dat) by calling the command:

grep fic REPORT |awk '{print $2 }' >> timeEvol.dat

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 ($>3 \AA$) while the transition region (around $0 \AA$) 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.

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

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