Extra exercise - Vibrational energy

The energy calculated with DFT is the electronic energy at 0 K; however, to model catalytic reactions we are usually intereseted in the energy landscape at finite temperature. In this exercise we will calculate the energy contributions from vibrations and see how they affect the splitting of Ni2 on Ru.

We calculate the vibrational energy in the harmonic approximation using the finite displacement method. For further reading, see for example the following sources:

[1]: Stoffel et al, Angewandte Chemie Int. Ed., 49, 5242-5266 (2010)

[2]: Togo and Tanaka, Scripta Materialia 108, 1-5 (2015)

Vibrational energy of the initial and final states

  1. Calculating vibrations requires tighter convergence than normal energy calculations. Therefore you should first take your already optimised initial and final state geometries from the NEB calculations and relax them further to \(f_\mathrm{max}=0.01\) eV/Å with the QuasiNewton optimiser and an energy cutoff of 500 eV. Converge the eigenstates to 1e-8. (Note that for other systems you might need even tighter convergence!)

Submit the structures to the queue. The optimisation should take 10-15 mins for each structure on 8 cores.

  1. Once you have done this you can calculate the vibrations using the vibrations module in ASE following the template below. We only calculate the vibrations of the adsorbate and assume that the frequencies of the substrate are unchanged - this is a common assumption. Use 4 displacements to fit the frequencies and the same calculator parameters as in a).

Submit the calculations for the initial and the final state to the queue. It will take a while to run, but you can start preparing your analysis (tasks c and d) or read some of the references in the meantime.

You cannot run this code as is - use it as a starting point for your python script!

import numpy as np
from ase.io import read
from ase.vibrations import Vibrations

from gpaw import GPAW, PW

slab = read(structure_file)
slab.calc = GPAW(xc='PBE',
                 mode=PW(500),
                 kpts={'size': kpts_size, 'gamma': True},
                 symmetry={'point_group': False},
                 txt='vib_initial.txt')
energy_initial = slab.get_potential_energy()
np.savetxt('energy_initial.txt', [energy_initial])

vib_initial = Vibrations(slab,
                         indices=N2_indices,
                         nfree=4,
                         name='vib_initial')
vib_initial.run()
vib_initial.summary(log='vib_initial_summary.log')
for i in range(6):
    vib_initial.write_mode(i)

The final three lines write out the frequencies and .traj files with animations of the different modes in order of their energy. Take a look at the vibrational modes in the files. Do you understand why the first mode has a low energy while the last one has a high energy?

  1. Analyse the frequencies in the harmonic approximation:

from ase.thermochemistry import HarmonicThermo

T = 300  # Kelvin
energies = [0.01, 0.05, 0.10]  # An example only - insert your calculated energy in eV levels here
vibs = HarmonicThermo(energies)
Efree = vibs.get_helmholtz_energy(T, verbose=True)
print('Free energy at 300 K: ', Efree)

The verbose keyword gives a detailed description of the different contributions to the free energy. For more information on what the different contributions are, see the ASE background webpage.

Now try to calculate how the different contributions change with temperature. You can for example make a for loop and use the get_entropy() and get_internal_energy() methods (see description here).

  1. Calculate how the vibrational energy affects the overall reaction energy. How does it change with temperature? Which contribution is important for the change in reaction energy?

  2. To make sure that your NEB is converged you should also calculate the vibrational energy of the transition state. Again, this requires tighter convergence than we have used in the NEB exercise. This takes a while to run so to save time, we provide the transition state geometry from a reasonably converged NEB (i.e. \(f_\mathrm{max}=0.01\), a cutoff energy of 800eV and 6x6 k-points) in the file TS.xyz. Calculate the vibrations with these parameters. How many imaginary modes do you get and how do they look? What does this mean?