Differences

This shows you the differences between two versions of the page.

--- exercises:2016_uzh_cmest:path_optimization_neb [2016/10/20 18:24] – [Visualize the trajectory] tmueller
+++ exercises:2016_uzh_cmest:path_optimization_neb [2020/08/21 10:15] (current) – external edit 127.0.0.1
@@ Line 177: / Line 177: @@
 This may take a couple of hours. Continue with the exercises below once the calculation finishes.
-====== Visualize the trajectory ======
+====== 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:
+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'':
 <code>
@@ Line 215: / Line 215: @@
 </note>
+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'':
+<code>
+awk '/E =/ {i=i+1; printf "%s %16.8f\n", i, $6}' ethane_neb_aba_8r.xyz
+</code>
+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.
+<code - 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
+</code>
+Once this run completes, you should find a file ''ethane_TS_va-VIBRATIONS-1.mol''.
+Now we are going to use the application //molden// (which you can load using ''module load molden'') to visualize the vibrational modes:
+<code>
+$ molden ethane_TS_va-VIBRATIONS-1.mol
+</code>
+Click the //Norm. Mode// checkbox in the //Molden Control// window to list all the modes. What is the lowest frequency you get? By clicking on it you can visualize it.
+The presence of a  negative (imaginary) mode means that it is actually a transition state (and not stable).
+Now repeat the same steps presented here for the bead with the lowest energy. What is now the first frequency you get in the list? Is this geometry stable?
+Please note: while you should get only 18 different frequencies you get 21 instead. That means that 3 frequencies are global rotations instead of modes in the molecule and should be ignored when looking for negative frequencies to identify whether a conformer is stable or not.