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:

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 8 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

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.