Nudged elastic band

Regtest see:

cp2k/tests/NEB/

Manual see:

%MOTION%BAND

Reference:

J. Chem. Phys. 113, 9978 (2000) for IT-NEB;

J. Chem. Phys. 113, 9901 (2000) for CI-NEB

Introduction

Consider a potential energy surface (PES) as below, the goal is to find the saddle point, transition structure between reactant and prodcut. The commonly-used algorithms include NEB,DIMER etc

Maxima or Saddle Points (maximum in one direction but minimum in other directions) have one or more negative (imaginary) frequencies.

Mathematically, the saddle point should fulfill the requirements:

1. The gradient should be zero, $\frac{dE}{dr} = 0 $

2. The sign of Hessian should be negative, $\frac{d^2E}{dr^2} < 0 $

3. Only ONE negative (imaginary) frequency

Exercies

In this exercise, we will study the $S_N2$ nucleophilic substitution reaction

$$ \text{Cl}^- + \text{CH}_3\text{Cl} \longleftrightarrow \text{Cl}\text{CH}_3 + \text{Cl}^-. $$

For energy and force evaluation, we will use the Parameterization Method 6 (PM6), which is a relatively inexpensive semiempirical model for electronic structure evaluation. For accurate reaction characterizations, ab initio methods, such as DFT or higher level, should be used.

In order to find the activation barrier for the reaction we will use the NEB method. The NEB method provides a way to find the minimum energy path (MEP) between two different conformations of the system corresponding local potential energy minima. The MEP is a one-dimensional path on the $3N$-dimensional potential energy landscape and every point on the MEP is a potential energy minimum in all directions perpendicular to the path. MEP passes at least one saddle point and the energy of the highest saddle point is the peak of the activation barrier of the reaction.

Geometry optimizations

To start the NEB calculation, we first need two different geometries at local energy minima, which we can obtain by running geometry optimization calculations. Use the following CP2K input script to obtain the first geometry:

geo.inp
&GLOBAL
  PRINT_LEVEL LOW
  PROJECT ch3cl
  RUN_TYPE GEO_OPT                # Geometry optimization calculation
&END GLOBAL

&MOTION
  &GEO_OPT                        # Parameters for GEO_OPT convergence
    MAX_FORCE 1.0E-4
    MAX_ITER 2000
    OPTIMIZER BFGS
    &BFGS
      TRUST_RADIUS [bohr] 0.1
    &END
  &END GEO_OPT
&END MOTION

&FORCE_EVAL
  METHOD Quickstep                # Quickstep - Electronic structure methods
  &DFT
   CHARGE -1                      # There is a negatively charged anion
    &QS
      METHOD PM6                  # Parametrization Method 6
      &SE
      &END SE
    &END QS
    &SCF                          # Convergence parameters for force evaluation
      SCF_GUESS ATOMIC
      EPS_SCF 1.0E-5
      MAX_SCF 50
       &OUTER_SCF
          EPS_SCF 1.0E-7
          MAX_SCF 500
       &END
    &END SCF
    &POISSON                      # POISSON solver for non-periodic calculation
      PERIODIC NONE
      PSOLVER WAVELET
    &END
  &END DFT

  &SUBSYS
    &CELL
      ABC 10.0 10.0 10.0
      PERIODIC NONE
    &END CELL
    &COORD
C      -4.03963494       0.97427857      -0.29785096
H      -3.97152206       2.11080568      -0.35500445
H      -3.95108814       0.19729141       0.53163699
H      -3.95810737       0.47922964      -1.32151084
Cl     -5.77885435       1.07089312      -0.04609427
Cl     -1.77974793       0.99422072      -0.29785096
    &END COORD
   &END SUBSYS
&END FORCE_EVAL
TASK 1
  • Run the geometry optimization calculation
  • Modify the initial geometry guess in the input and run a second optimization to obtain the other local minimum geometry (The other Cl needs to make a covalent bond with the C).

NEB

Now we are ready to start the NEB. Obtain the two optimized geometries from the GEO_OPT trajectories (last step in the ch3cl-pos-1.xyz file) to a separate folder and name them init.xyz and final.xyz. Now the NEB calculation can be run by the following input script:

neb.inp
&GLOBAL
  PRINT_LEVEL LOW
  PROJECT ch3cl
  RUN_TYPE BAND                      # Nudged elastic band calculation
&END GLOBAL

&MOTION
  &BAND
    NUMBER_OF_REPLICA 10             # Number of "replica" geometries along the path
    K_SPRING 0.05
    &OPTIMIZE_BAND
      OPT_TYPE DIIS
      &DIIS
        MAX_STEPS 1000
      &END
    &END
    BAND_TYPE CI-NEB                 # Climbing-image NEB
    &CI_NEB
      NSTEPS_IT  5                   # First take 5 normal steps, then start CI
    &END
    &REPLICA
      COORD_FILE_NAME init.xyz
    &END
    &REPLICA
      COORD_FILE_NAME final.xyz
    &END
    &PROGRAM_RUN_INFO
     INITIAL_CONFIGURATION_INFO
    &END
  &END BAND
&END MOTION

&FORCE_EVAL
  METHOD Quickstep
  &DFT
    ... same as GEO_OPT ...
  &END DFT
  &SUBSYS
    &CELL
      ABC 10.0 10.0 10.0
      PERIODIC NONE
    &END CELL
    &TOPOLOGY
      COORD_FILE_NAME init.xyz
      COORDINATE xyz
    &END TOPOLOGY
  &END SUBSYS
&END FORCE_EVAL

In the main output of the NEB calculation, you will find a number of the following sections:

*******************************************************************************
 BAND TYPE                     =                                          CI-NEB
 BAND TYPE OPTIMIZATION        =                                              SD
 STEP NUMBER                   =                                               1
 NUMBER OF NEB REPLICA         =                                              10
 DISTANCES REP =        0.252266        0.229810        0.225933        0.227113
                        0.201788        0.138731        0.138676        0.204817
                        0.252328
 ENERGIES [au] =      -24.815004      -24.813368      -24.808500      -24.803002
                      -24.800607      -24.802583      -24.805589      -24.809114
                      -24.813372      -24.815004
 BAND TOTAL ENERGY [au]        =                             -248.08548199708440
 *******************************************************************************

These sections show for every replica geometry along the NEB trajectory, the distance to its neighbors and its energy. The final section corresponds to the converged NEB trajectory.