====== Nudged Elastic Band ======
In this exercise you will compute the energy profile for a simple reaction in a planar cluster of 7 Ar atoms with the NEB method.
The NEB method requires at least the start- and end-configuration between which the reaction path should be computed. In addition it is advantageous to add a guess for an intermediate configuration. This is useful in particular when multiple reaction paths are possible and one is interested in a specific one. By adding a suitable intermediate configuration, the optimization can be guided towards the path of interest. However, keep in mind that the intermediate configuration is merely an initial guess and it's further optimized by the NEB algorithm.
The aim of this exercise is to compute the activation energy required for moving atom 2 (dark blue) into the center of the cluster.
In 2D there are at least two possible paths:
^ Path 1 ^ Path 2 ^
| Direct exchange with the central atom | Coordinated rotation of three atoms |
| {{neb_path1.gif|}} | {{neb_path2.gif|}} |
===== Path 1: direct exchange =====
=== 1. Step ===
Save the following commented CP2K input file to a file named ''neb1.inp''
&GLOBAL
RUN_TYPE BAND
PROJECT_NAME neb1 ! the neb calculation will produce a few output files, that will be labeled with this name
&END GLOBAL
&MOTION
&CONSTRAINT ! we are working in 2d, and to keep things simple, we fix the z position of all the atoms (atom 1 to atom7)
&FIXED_ATOMS
COMPONENTS_TO_FIX Z
LIST 1..7
&END
&END
&BAND
NPROC_REP 1 ! numebr of processors needed for each replica
BAND_TYPE IT-NEB ! type of NEB to be performed
NUMBER_OF_REPLICA 10 ! number of replicas required
&OPTIMIZE_BAND ! parameters for the optimizer. This section is required and should not be changed
OPT_TYPE DIIS
&DIIS
MAX_STEPS 1000
N_DIIS 3
&END
&END
&REPLICA ! coordinates of the satting configuration (REQUIRED: has to be the first one in the input)
&COORD
-0.0000000000 0.0000000000 -0.0000000000
3.8030201671 -0.0000003430 -0.0000000000
-3.8030201671 0.0000003430 0.0000000000
1.9019125593 3.2944696295 -0.0000000000
1.9019119654 -3.2944699726 0.0000000000
-1.9019119654 3.2944699726 0.0000000000
-1.9019125593 -3.2944696295 0.0000000000
&END
&END REPLICA
&REPLICA ! intermediate configuration
&COORD
2.215467 -0.734540 0.000000
2.144383 0.900403 0.000000
-3.803020 0.000000 0.000000
1.901913 3.294470 -0.000000
1.901912 -3.294470 0.000000
-1.901912 3.294470 0.000000
-1.901913 -3.294470 0.000000
&END
&END REPLICA
&REPLICA ! ending configuration (REQUIRED: has to be the last one in the input)
&COORD
3.8030201671 -0.0000003430 -0.0000000000
0.0000000000 0.0000000000 -0.0000000000
-3.8030201671 0.0000003430 0.0000000000
1.9019125593 3.2944696295 -0.0000000000
1.9019119654 -3.2944699726 0.0000000000
-1.9019119654 3.2944699726 0.0000000000
-1.9019125593 -3.2944696295 0.0000000000
&END
&END REPLICA
&END BAND
&END MOTION
&FORCE_EVAL ! this section is the same as usual
METHOD FIST
&MM
&FORCEFIELD
&SPLINE
EMAX_SPLINE 10000
&END
&NONBONDED
&LENNARD-JONES
atoms Ar Ar
EPSILON [K_e] 119.8
SIGMA [angstrom] 3.401
RCUT [angstrom] 25.0
&END LENNARD-JONES
&END NONBONDED
&CHARGE
ATOM Ar
CHARGE 0.0
&END CHARGE
&END FORCEFIELD
&POISSON
PERIODIC NONE
&EWALD
EWALD_TYPE none
&END EWALD
&END POISSON
&END MM
&SUBSYS
&CELL
ABC [angstrom] 20 20 20
PERIODIC NONE
&END CELL
&COORD !starting coordinates, to be repeated here.
Ar -0.0000000000 0.0000000000 -0.0000000000
Ar 3.8030201671 -0.0000003430 -0.0000000000
Ar -3.8030201671 0.0000003430 0.0000000000
Ar 1.9019125593 3.2944696295 -0.0000000000
Ar 1.9019119654 -3.2944699726 0.0000000000
Ar -1.9019119654 3.2944699726 0.0000000000
Ar -1.9019125593 -3.2944696295 0.0000000000
&END COORD
&END SUBSYS
&END FORCE_EVAL
=== 2. Step: Run CP2K ===
$ cp2k.popt -i neb1.inp -o neb1.out
=== 3. Step ===
For the NEB calcualtions, CP2k produces a few output files. The most important are:
* neb1.out : standard CP2K output file. It tells you whether that the calculation is completed.\\ (See:[[single_point_calculation|Computation of the Lennard Jones curve for two Ar atoms]].Part I, Step 3)
* neb1-pos-Replica_nr_XXX-1.xyz : those are the replica optimization trajectories. You get a trajectory for each replica.
* neb1-BANDXXX.out : geometry optimization output for each replica.
=== 4. Step: Checking the trajectory ===
Here is a short script to create a ''movie.xyz'' file showing the trajectory from starting to ending point.
Use this procedure to make sure that the trajectory you obtain is the one you actaully want to study. The movie.xyz can be read by VMD.
$ for a in 01 02 03 04 05 06 07 08 09 10 ; do tail -n 9 neb1-pos-Replica_nr_${a}-1.xyz >> movie.xyz ; done
=== 5. Step: Generating the energy profile ===
Here is a short script to create an energy profile as a function of the replica number.
$ for a in 1 2 3 4 5 6 7 8 9 10 ; do grep ENERGY neb1-BAND${a}.out | tail -n 1 | awk '{print '${a}',$9}' ; done > neb1_profile
The energy profile will be stored in the file ''neb1_profile''. Any plotting program should be able to visualize it, e.g. gnuplot:
$ echo "plot 'neb1_profile' w l; pause mouse" | gnuplot
===== Path 2: Coordinated rotation =====
Here is the input file for the PATH2.
Following the same procedure as above, you can obtain a trajectory and an energy profile.
&GLOBAL
RUN_TYPE BAND
PROJECT_NAME neb2
&END GLOBAL
&MOTION
&CONSTRAINT
&FIXED_ATOMS
COMPONENTS_TO_FIX Z
LIST 1..7
&END
&END
&BAND
NPROC_REP 1
BAND_TYPE IT-NEB
NUMBER_OF_REPLICA 20
&OPTIMIZE_BAND
OPT_TYPE DIIS
&DIIS
MAX_STEPS 1000
N_DIIS 3
&END
&END
&REPLICA
&COORD
-0.0000000000 0.0000000000 -0.0000000000
3.8030201671 -0.0000003430 -0.0000000000
-3.8030201671 0.0000003430 0.0000000000
1.9019125593 3.2944696295 -0.0000000000
1.9019119654 -3.2944699726 0.0000000000
-1.9019119654 3.2944699726 0.0000000000
-1.9019125593 -3.2944696295 0.0000000000
&END
&END REPLICA
&REPLICA
&COORD
3.8030201671 -0.0000003430 -0.0000000000
1.9019125593 3.2944696295 -0.0000000000
-3.8030201671 0.0000003430 0.0000000000
0 0 0
1.9019119654 -3.2944699726 0.0000000000
-1.9019119654 3.2944699726 0.0000000000
-1.9019125593 -3.2944696295 0.0000000000
&END
&END REPLICA
&REPLICA
&COORD
1.9019125593 3.2944696295 -0.0000000000
0 0 0
-3.8030201671 0.0000003430 0.0000000000
3.8030201671 -0.0000003430 -0.0000000000
1.9019119654 -3.2944699726 0.0000000000
-1.9019119654 3.2944699726 0.0000000000
-1.9019125593 -3.2944696295 0.0000000000
&END
&END REPLICA
&END BAND
&END MOTION
&FORCE_EVAL
METHOD FIST
&MM
&FORCEFIELD
&SPLINE
EMAX_SPLINE 10000
&END
&NONBONDED
&LENNARD-JONES
atoms Ar Ar
EPSILON [K_e] 119.8
SIGMA [angstrom] 3.401
RCUT [angstrom] 25.0
&END LENNARD-JONES
&END NONBONDED
&CHARGE
ATOM Ar
CHARGE 0.0
&END CHARGE
&END FORCEFIELD
&POISSON
PERIODIC NONE
&EWALD
EWALD_TYPE none
&END EWALD
&END POISSON
&END MM
&SUBSYS
&CELL
ABC [angstrom] 20 20 20
PERIODIC NONE
&END CELL
&COORD
Ar -0.0000000000 0.0000000000 -0.0000000000
Ar 3.8030201671 -0.0000003430 -0.0000000000
Ar -3.8030201671 0.0000003430 0.0000000000
Ar 1.9019125593 3.2944696295 -0.0000000000
Ar 1.9019119654 -3.2944699726 0.0000000000
Ar -1.9019119654 3.2944699726 0.0000000000
Ar -1.9019125593 -3.2944696295 0.0000000000
&END COORD
&END SUBSYS
&END FORCE_EVAL
===== Questions =====
* Sketch the reaction profile for the two paths presented, reporting shape and activation energy.
* What are the major differences between the profiles? How do you think the result will change by adopting different number of replicas?