Predicting the vibrational lifetime of CO on a copper surface using electronic friction

It is recommended to read the theory section (Part 1 of this tutorial) before attempting these calculations. In this section, we will calculate the phonon lifetime of the internal stretch of CO on a copper surface, both for an aperiodic vibration and a full overlayer vibration.

Calculation expense and runtime

In the following tutorial, progressively larger supercells of CO/Cu(100) will be investigated. The larger cells require significant computational resources and runtimes in order to calculate. If these are unavailable, it is recommended to perform the calculations for the smallest cells only, as the tutorial will provide the results for larger cells.

Workflow

We are going to carry out:

1. Geometry Optimization: Obtain a stable structure of CO overlayer on a Cu(100) surface
2. Generate supercells: Replicate the primitive cell provided to generate supercells of CO/Cu(100)
3. Calculate normal modes: Carry out a finite-difference phonon calculation using ASE to evaluate the normal modes and Hessian for each cell.
4. Calculate electronic friction: Calculate electronic friction tensors for each cell using FHI-aims.
5. Plot vibrational lifetimes: Project electronic friction tensor along normal modes to generate lifetimes.

Materials

Input files and scripts can be found here

Related starting files

file or dir name	description
control.in	FHI-aims control file, used as a template for all DFPT calculations
geometry.in	FHI-aims geometry file of benzene.
submit.sh	Example script for submission of calculations to HPC (only an example, all calculations can be run locally)
generate_supercells.py	Python script to periodically repeat unit cell using ASE.
vib.py	Python script to perform normal mode / Hessian calculation using ASE.
project_tensor.py	Python script to project the Cartesian electronic friction tensor onto normal modes
collect_lifetimes.py	Python script to just collect the lifetime data from different supercells
plot_q0_lifetime.py	Python script to plot the convergence of the \(\mathbf{q}=0\) lifetime with respect to cell size
plot_qave_lifetime.py	Python script to plot the convergence of the \(\bar{q}\) lifetime with respect to cell size

Generate supercells using ASE

python generate_supercells.py

This should generate subdirectories and geometry files for each investigated unit cell.

Calculate normal modes

Unfortunately, normal mode analysis is not yet available during the finite difference calculation to evaluate electronic friction. Thus, we perform an additional finite difference calculation to evaluate the Hessian and perform normal mode analysis. Note: this script uses start.in for the geometry file input name, so please rename geometry.in when employing this script.

With increasing unit cell size, we can reduce the k_grid size, to ensure that the calculations are performed at a constant \(\mathbf{k}\)-mesh density. Suggested mesh sizes for this system are: 32 32 1 for primitive unitcell, 16 16 1 for the \(2\times2\) cell, 10 10 1 for the \(3\times3\) cell, and 8 8 1 for the \(4\times4\) cell. Please refer to the recommended reference for details on determining a converged \(\mathbf{k}\)-mesh size for electronic friction calculations.

The run command is:

python vib.py > normal_modes.out

Output files

file or dir name	description
normal_modes.out	Contains list of frequencies and Hessian.
vib.X.pckl	Pickle files for each displacement.
parameters.ase	Calculator settings from ASE.

Calculate electronic friction for each cell

The provided control.in provides several options to control the calculation of electronic friction:

calculate_friction: either numerical_friction (finite difference) or DFPT (density functional perturbation theory) are options.

friction_broadening_width: controls the width of the broadening function in eV, should be carefully tested by the user.

friction_perturbation: controls the \(\hbar\omega\) term that enters the electronic friction evaluation (in eV). From the normal mode analysis, \(0.252\) eV is the suggested perturbation energy for the internal stretch of CO. This should be adjusted to the perturbing energy of the mode investigated.

friction_delta_type: controls the type of broadening function. The developers recommend gaussian, though little difference in electronic friction has been observed for different broadening functions.

Since the expression for electronic friction includes a sum over occupied and unoccupied states, it is important to include unoccupied eigenstates in the calculation, we can do so using the empty_states keyword, the sensitivity of the electronic friction to the number of states should be tested by the user.

First, copy control.in to each subdirectory

find . -type d -exec cp control.in {} ';'

With increasing unit cell size, we can reduce the k_grid size. Suggested mesh sizes are: 32 32 1 for primitive unitcell, 16 16 1 for the \(2\times2\) cell, 10 10 1 for the \(3\times3\) cell, and 8 8 1 for the \(4\times4\) cell.

Secondly, we need to add the keyword calculate_friction .true. to each C and O atom in each geometry.in. This enables the calculation of the friction tensor to just the C and O atoms.

sed -i '' '/C$/s/$/\ncalculate_friction .true./' */geometry.in

and

sed -i '' '/O$/s/$/\ncalculate_friction .true./' */geometry.in

The evaluation of electronic friction requires careful checking of convergence with respect to k-mesh size

Copper in particular is difficult to obtain converged electronic friction for, due to its low density of states at the Fermi level. See the following reference for a thorough investigation of the convergence behaviour of CO/Cu(100) for electron-phonon coupling calculations.

We are now ready to submit the calculations, either by copying a submit script to each directory and submitting to a HPC scheduler or by running FHI-aims locally.

Once the calculations are finished successfully, each subdirectory should be populated by a friction_tensor.out file.

Output files

file or dir name	description
friction_tensor.out	Mass-weighted friction tensor in ps\(^{-1}\).
friction_eigenvectors.jmol	Eigenvectors of the friction tensor, the "friction modes" can be visualized with Jmol.

Calculate vibrational lifetimes

For each cell, we will now calculate the vibrational lifetimes by projecting the electronic friction tensor onto the normal modes:

python project_tensor.py <friction_tensor.out> <normal_modes.out> > lifetimes.out

Plot vibrational lifetimes

We will collect the vibrational lifetimes in a single file:

python collect_lifetimes.py <calculation_dirs>

Output files

file or dir name	description
lifetimes.txt	Contains list of calculated lifetimes for different unit cell sizes.

We can now plot the convergence with respect to supercell size of the in-phase CO stretch lifetime on the surface with \(\mathbf{q}=0\).

python plot_q0_lifetime.py

Producing:

As can be seen in the plot above, for CO/Cu(100), the in-phase coherent stretch of CO does not reach convergence with respect to unit cell size (\(N_\mathrm{CO}\) is the number of CO molecules in the unit cell) with these settings. The low density-of-states around the Fermi level for copper makes this difficult to converge. If the reader is interested, they can reattempt this calculation with more layers, a larger \(\mathbf{k}\)-point mesh, or a different metal surface (hint: platinum has a much higher density of states at the Fermi level.)

We can do the same for the aperiodic vibrational lifetime for the CO stretch by summing over \(\mathbf{q}\) (see theory section):

python plot_qave_lifetime.py

Producing:

We can see the aperiodic stretch lifetime is easier to converge, with a relatively stable lifetime reached for cells with 8 or more CO molecules. The converged lifetime of approximately \(2\) ps is comparable to the experimental measure of the lifetime after excitation by infra-red (measured by Morin et al, J. Chem. Phys., 1992), suggesting that this system may dephase quickly after light excitation.