Path optimization using NEB
In one of the last exercises you calculated the energy for Ethane for two slightly different geometries and noticed that the geometry optimization was not able to change one structure into the other with lower energy. As presented in the lecture, it may happen quiet often that a minimization algorithm gets stuck in a local minimum, respectively it is not guaranteed to find the global minimum.
In this exercise, we will therefore perform Nudged Elastic Band (NEB) calculations using the same molecule as before and investigate the energy path between the two geometries.
Following are four geometry files you should put/create in a new exercise directory:
- ethane_1_opt.xyz
8 i = 12, E = -14.9518242480 C 6.7709731556 5.9999999991 6.2147000005 C 5.2290258859 6.0000000007 6.2147000001 H 4.8193442852 6.0000000005 5.1955521086 H 4.8183461262 6.8823861870 6.7231540475 H 4.8183461248 5.1176138164 6.7231540484 H 7.1806551657 5.9999999994 5.1955522543 H 7.1816533264 6.8823860579 6.7231539750 H 7.1816533250 5.1176139384 6.7231539741
- ethane_s1.xyz
8 i = 76, E = -14.9559722838 C 6.7640435612 6.0000000003 5.9997401503 C 5.2359564388 5.9999999997 6.0002598497 H 4.8332682152 6.0000000001 7.0232923722 H 4.8330445407 5.1142018851 5.4887808109 H 4.8330445407 6.8857981148 5.4887808113 H 7.1667317848 5.9999999999 4.9767076278 H 7.1669554593 6.8857981149 6.5112191891 H 7.1669554593 5.1142018852 6.5112191887
- ethane_ts.xyz
8 i = 76, E = -14.9518421887 C 6.7713019119 5.9963482236 5.9999305287 C 5.2296433443 6.0033677689 6.0065640271 H 4.8226439717 6.0198987997 7.0258669051 H 4.8144385670 5.1135492476 5.5115353199 H 4.8205805610 6.8779314573 5.4838302216 H 7.1765327350 5.1402836770 5.4459664211 H 7.1832404413 6.9024073310 5.5339472993 H 7.1852801526 5.9447026992 7.0166786010
- ethane_s2.xyz
8 i = 76, E = -14.9559544815 C 6.7635882192 5.9976047386 6.0012623335 C 5.2356460989 6.0047571469 6.0004006617 H 4.8320044755 6.0067336493 7.0232527262 H 4.8290963082 5.1202363024 5.4891970991 H 4.8366693774 6.8919269855 5.4881514124 H 7.1628913027 5.1334052558 6.5507178317 H 7.1674882265 5.9517120047 4.9798028986 H 7.1709320051 6.9026197936 6.4738574951
You could in principle start from the geometries you already optimized. In fact, the files ethane_1_opt.xyz and ethane_s1.xyz are geometry optimizations resulting from the previous ethane1.xyz and ethane2.xyz with slightly different settings (more details below) and some modifications.
The input file you need then looks as follows:
- ethane_neb_aba.inp
&GLOBAL PROJECT ethane_neb_aba RUN_TYPE BAND PRINT_LEVEL MEDIUM &END GLOBAL &FORCE_EVAL METHOD Quickstep ! Electronic structure method (DFT,...) &DFT BASIS_SET_FILE_NAME BASIS_MOLOPT POTENTIAL_FILE_NAME POTENTIAL &POISSON ! Solver requested for non periodic calculations PERIODIC XYZ &END POISSON &SCF ! Parameters controlling the convergence of the scf. This section should not be changed. SCF_GUESS ATOMIC EPS_SCF 1.0E-6 MAX_SCF 300 &END SCF &XC ! Parameters needed to compute the electronic exchange potential &XC_FUNCTIONAL PBE &END XC_FUNCTIONAL &END XC &END DFT &SUBSYS &CELL ABC 12. 12. 12. PERIODIC XYZ &END CELL &TOPOLOGY ! Section used to center the atomic coordinates in the given box. Useful for big molecules &CENTER_COORDINATES &END COORD_FILE_FORMAT xyz COORD_FILE_NAME ./ethane_1_opt.xyz &END &KIND H ELEMENT H BASIS_SET DZVP-MOLOPT-GTH POTENTIAL GTH-PBE-q1 &END KIND &KIND C ELEMENT C BASIS_SET DZVP-MOLOPT-GTH POTENTIAL GTH-PBE-q4 &END KIND &END SUBSYS &END FORCE_EVAL &MOTION &BAND BAND_TYPE CI-NEB NUMBER_OF_REPLICA 8 K_SPRING 0.05 &CONVERGENCE_CONTROL MAX_FORCE 0.0010 RMS_FORCE 0.0050 &END ROTATE_FRAMES TRUE ALIGN_FRAMES TRUE &CI_NEB NSTEPS_IT 2 &END &OPTIMIZE_BAND OPT_TYPE DIIS OPTIMIZE_END_POINTS FALSE &DIIS MAX_STEPS 1000 &END &END &PROGRAM_RUN_INFO &END &CONVERGENCE_INFO &END &REPLICA COORD_FILE_NAME ./ethane_s1.xyz &END &REPLICA COORD_FILE_NAME ./ethane_ts.xyz &END &REPLICA COORD_FILE_NAME ./ethane_s2.xyz &END &END BAND &END MOTION
One notable difference to the previous input files is the specification of the periodic boundary conditions in the &POISSON section and in the &CELL (and the increased size of the box – now 12 Å instead of 10 Å). This is to use a different solver configuration which behaves better in combination with NEB and if the box is big enough it will not change the physics.
To run this simulation, we use a slightly different command which will return right away:
$ nohup mpirun -np 4 cp2k.popt -i ethane_neb_aba.inp -o ethane_neb_aba.out &
In the background, the process is still running, which you can verify by either watching the changes to the output file (exit this command with CTRL+c) using
$ tail -f ethane_neb_aba.out
or by looking at the list of your processes:
$ ps uxf
Note the use of mpirun -np 4 together with the cp2k.popt executable. This tells the machine to run it using the MPI parallelization system on 4 cores, which should speed up your calculations. The prefix nohup lets the command continue even if you log out and the ampersand & makes it run in the background. Remember that tcopt2 is shared among you and that there is only a limited number of cores. Before launching, you can use the command top to check if there are enough cores available for your calculation.
This may take a couple of hours. Continue with the exercises below once the calculation finishes.
Visualize the trajectory and plot the energy curve
When you take another peek at the input file you used to run the calculation, you will notice that we specified NUMBER_OF_REPLICA 8, which means that CP2K will generate in total 8 beads (3 we specified in the &REPLICA sections, 5 will be generated automatically by interpolation).
You should therefore find 8 files named ethane_neb_aba-pos-Replica_nr_1-1.xyz..ethane_neb_aba-pos-Replica_nr_8-1.xyz in your exercise directory, containing the optimization of each bead. To get the trajectory over the band, we extract the last frame (see the tip below) and write it into a separate file named ethane_neb_aba_8r.xyz:
$ for i in {1..8} ; do tail -n 10 "ethane_neb_aba-pos-Replica_nr_${i}-1.xyz" >> ethane_neb_aba_8r.xyz ; done
Look at the movement again using VMD.
A XYZ file consists of one or multiple blocks (frames) of the following:
8 <-- the number of atoms i = 0, E = -14.9559722838 <-- some comment, CP2K writes the number of the iteration and the resulting energy here C 6.7640435612 6.0000000003 5.9997401503 <-- the atomic symbol and position C 5.2359564388 5.9999999997 6.0002598497 <-- .. for the specified number of atoms H 4.8332682152 6.0000000001 7.0232923722 H 4.8330445407 5.1142018851 5.4887808109 H 4.8330445407 6.8857981148 5.4887808113 H 7.1667317848 5.9999999999 4.9767076278 H 7.1669554593 6.8857981149 6.5112191891 H 7.1669554593 5.1142018852 6.5112191887
Extracting the last frame (the optimized geometry in case of a geometry optimization) is therefore simple once you know the number of atoms using the command tail, which can be used to get the last N lines of a file using the switch -n N.
In our case of 8 atoms this is:
$ tail -n 10 geo_opt_output.xyz
Since CP2K writes the energy in the comment line of each frame in the XYZ file (see tip above), we can extract the energy values for each bead directly from the newly generated ethane_neb_aba_8r.xyz:
awk '/E =/ {i=i+1; printf "%s %16.8f\n", i, $6}' ethane_neb_aba_8r.xyz
Create a plot for this energy curve.
Vibrational analysis
To verify whether the point at the highest energy is actually a transition state, we will be doing a vibrational analysis.
First identify the bead with the highest energy (see exercise above) and create a new XYZ file named ethane_neb_aba_TS.xyz with the respective coordinates (extracted from either the correct ethane_neb_aba-pos-Replica_nr_N-1.xyz file or the ethane_neb_aba_8r.xyz).
Use the following input file and the same command as above (with different input and output file names of course) to generate the analysis.
- ethane_TS_va.inp
&GLOBAL PROJECT ethane_TS_va RUN_TYPE NORMAL_MODES PRINT_LEVEL MEDIUM &END GLOBAL &FORCE_EVAL METHOD Quickstep ! Electronic structure method (DFT,...) &DFT BASIS_SET_FILE_NAME BASIS_MOLOPT POTENTIAL_FILE_NAME POTENTIAL &POISSON ! Solver requested for non periodic calculations PERIODIC XYZ &END POISSON &SCF ! Parameters controlling the convergence of the scf. This section should not be changed. SCF_GUESS ATOMIC EPS_SCF 1.0E-7 MAX_SCF 300 &END SCF &XC ! Parametes needed to compute the electronic exchange potential &XC_FUNCTIONAL PBE &END XC_FUNCTIONAL &END XC &END DFT &SUBSYS &CELL ABC 12. 12. 12. PERIODIC XYZ &END CELL &TOPOLOGY ! Section used to center the atomic coordinates in the given box. Useful for big molecules &CENTER_COORDINATES &END COORD_FILE_FORMAT xyz COORD_FILE_NAME ./ethane_neb_aba_TS.xyz &END &KIND H ELEMENT H BASIS_SET DZVP-MOLOPT-GTH POTENTIAL GTH-PBE-q1 &END KIND &KIND C ELEMENT C BASIS_SET DZVP-MOLOPT-GTH POTENTIAL GTH-PBE-q4 &END KIND &END SUBSYS &END FORCE_EVAL &VIBRATIONAL_ANALYSIS NPROC_REP 1 DX 0.01 FULLY_PERIODIC &PRINT &MOLDEN_VIB &END &CARTESIAN_EIGS &END &PROGRAM_RUN_INFO &EACH REPLICA_EVAL 1 &END &END &END &END
Once this run completes, you can grep for the modes frequency with:
$ grep "VIB|Frequency" ethane_TS_va.out
The presence of a negative (imaginary) mode means that it is actually a transition state (and not stable). Are there any of that kind for this structure ? If so, what are the atomic displacements associated with it ?
Note that cp2k also produces a ethane_TS_va-VIBRATIONS-1.mol file. It is intended to be read by the Molden program, which makes the visualization of the modes possible.