Path optimization using NEB

In the last exercise you have 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 three geometry files you should put/create in a new exercise directory:

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.

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:

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

We replaced the CP2K executable cp2k.sopt with cp2k.popt, which is a parallel version of CP2K. By prefixing the command with mpirun -np 8 we tell it to run it using the MPI system using 8 cores. And finally to have the command continue to run even if you log out, we prefixed everything with nohup. The ampersand & at the end is to run everything in the background.

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:

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:

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.

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:

$ molden ethane_TS_va-VIBRATIONS-1.mol

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.