Tutorial 1: the ionization potential and electron affinity of ozone

In this tutorial, we will calculate the ionisation potential and electron affinity of ozone.

The input

The input file for this calculation can be downloaded here. Just to briefly highlight the most important details of the workflow block

2  "workflow": {
3    "functional": "ki",
4    "method": "dscf",

here we select the KI functional (as opposed to KIPZ),

3    "functional": "ki",
4    "method": "dscf",
5    "init_orbitals": "kohn-sham",

specifies that we are going to calculate the screening parameters via a ΔSCF procedure, whereby we compute the energies of various \(N\), \(N-1\), and \(N+1\)-electron systems (see the theory section for details), and

4    "method": "dscf",
5    "init_orbitals": "kohn-sham",
6    "alpha_numsteps": 1,

specifies that we have chosen to use the Kohn-Sham orbitals as our variational orbitals. This is common practice for molecules.

Meanwhile, the atoms block describes the both the cell and the atoms it contains. If you are familiar with Quantum ESPRESSO input files then most of this should look very familiar to you (albeit in JSON format).

Running the calculation

In order to run the calculation, simply run

koopmans run ozone.json | tee ozone.md

Tip

In order to run in parallel, set the PARA_PREFIX environment variable to mpirun -np 4 (or similar)

The outline

First, let us inspect the contents of ozone.md. This provides us with an outline of the workflow that we ran. It is in markdown format, which is rendered nicely by modern code editors such as vscode.

After the header we can see there are a list of calculations that have been run by koopmans. These come under three headings.

Initialization

The first step in any Koopmans calculation is the initialization step. In this step we initialize the density and the variational orbitals.

  • Initialization

    • ✅ 01-dft_init_nspin1 completed

    • ✅ 02-dft_init_nspin2_dummy completed

    • ✅ 03-convert_files_from_spin1to2 completed

    • ✅ 04-dft_init_nspin2 completed

For this calculation we can see that koopmans has run four calculations. These initialize the density with the PBE density. Indeed, from this point onwards in the calculation the density will never change, because the KI functional yields the same density as the base functional. (N.B. This is not true of KIPZ.)

These PBE calculations also have provided us with our variational orbitals – we can see that the calculation has selected the PBE Kohn-Sham orbitals as the variational orbitals (because earlier we set "init_orbitals": "kohn-sham").

But why four calculations? The reason for this is that the calculations we will perform later involve the addition/removal of a single electron, which means the density we need to generate here must correspond to a nspin = 2 calculation. However, we know ozone is a closed-shell molecule and simply performing a nspin = 2 PBE calculation risks introducing spin contamination (i.e. falling into a local minimum where \(n^\uparrow(\mathbf{r}) \neq n^\downarrow(\mathbf{r})\)).

This sequence of four calculations is designed to avoid this; we first optimize the density constrained to be spin-unpolarized, and only once that density has been minimized do we lift this restriction. This additional step can be disabled by adding "fix_spin_contamination": false to the workflow block of ozone.json.

The input and output Quantum ESPRESSO files can all be found in the corresponding directory 01-koopmans-dscf/01-initialization/.

Calculating the screening parameters

The second step in the calculation involves the calculation of the screening parameters:

  • Calculate Screening Via DSCF

    • Iteration 1

      • ✅ 01-ki completed

      • Orbital 1

        • ✅ 01-dft_n-1 completed

      • Orbital 2

        • ✅ 01-dft_n-1 completed

      • Orbital 3

        • ✅ 01-dft_n-1 completed

etc. Here, we are calculating the screening parameters using the ΔSCF method. For filled orbitals (orbitals 1-9 of ozone) this requires an \(N-1\)-electron PBE calculation where we freeze the ith orbital, empty it, and allow the rest of the density to relax. This yields \(E_i(N-1)\). Meanwhile, \(E(N)\), \(\varepsilon_{i}^\alpha(1)\), and \(\varepsilon_{i}^0(1)\) are all obtained during the trial KI calculation calc_alpha/ki. Together with the value of our guess for the screening parameters (\(\alpha^0_i = 0.6\)), this is sufficient to update our guess for \(\alpha_i\) (see the theory section for details).

The procedure for empty orbitals is slightly different, as we can see when it comes to orbital 10:

  • Orbital 10

    • ✅ 01-dft_n+1_dummy completed

    • ✅ 02-pz_print completed

    • ✅ 03-dft_n+1 completed

where now we must call Quantum ESPRESSO several times in order to obtain \(E_i(N+1)\).

Click here for detailed descriptions of each calculation
01-dft_n+1_dummy and 02-pz_print

preliminary calculations that generate files required by the subsequent constrained DFT calculation

03-dft_n+1

a \(N+1\)-electron PBE calculation where we freeze the 10th orbital, fill it, and allow the rest of the density to relax. This yields \(E_i(N+1)\)

At the end of this section we can see a couple of tables:

α

1

2

3

4

5

6

7

8

9

10

0

0.6

0.6

0.6

0.6

0.6

0.6

0.6

0.6

0.6

0.6

1

0.655691

0.727571

0.78386

0.663859

0.772354

0.726848

0.729968

0.741899

0.779264

0.717389

ΔEi - λii (eV)

1

2

3

4

5

6

7

8

9

10

0

-0.47736

-0.920499

-1.12409

-0.452149

-1.05055

-0.830893

-0.785139

-0.888597

-1.05016

0.711209

The first table lists the screening parameters \(\alpha_i\) that we obtained – we can see from row 0 we started with a guess of \(\alpha_i = 0.6\) for every orbital i, and row 1 shows the alpha values.

The second table lists \(\Delta E_i - \lambda_{ii}^\alpha\). This is a measure of how well converged the alpha values are: if this value is close to zero, then the alpha values are well-converged. Note that the values listed above correspond to our starting guess of \(\alpha_i = 0.6\); this table does not show how well-converged the final alpha values are.

Note

In order to see how well-converged our new screening parameters are, try increasing alpha_numsteps in the input file from 1 to 2. Can you make sense of the contents of the resulting tables?

The final calculation

Having determined the screening parameters, the final KI calculation is now run:

  • ✅ 03-ki_final completed

Workflow complete 🎉

The files

The markdown output provides an outline of the workflow that we ran, but for more details we can inspect the Quantum ESPRESSO input and output files themselves. These can be found in a nested directory structure that mirrors the outline we just followed. For example, the input and output files for the final KI calculation can be found in 01-koopmans-dscf/03-ki_final/. To get a sense of the entire directory structure, you can run the command

tree -I TMP-CP 01-koopmans-dscf/

Extracting the ionization potential and electron affinity

Let’s now extract the KI ionization potential and electron affinity for ozone from our calculation.

The ionization potential (IP) corresponds to the negative of the energy of the HOMO (highest occupied molecular orbital). If you open 01-koopmans-dscf/03-ki_final/ki_final.cpo and search near the bottom of the file you will see a section something like

...
HOMO Eigenvalue (eV)

-12.5199

LUMO Eigenvalue (eV)

-1.8218

Eigenvalues (eV), kp =   1 , spin =  1

-40.1869  -32.9130  -24.2288  -19.6841  -19.4902  -19.2696  -13.6037  -12.7618  -12.5199

Empty States Eigenvalues (eV), kp =   1 , spin =  1

-1.8218

Electronic Gap (eV) =    10.6981


Eigenvalues (eV), kp =   1 , spin =  2

-40.1869  -32.9130  -24.2288  -19.6841  -19.4902  -19.2696  -13.6037  -12.7618  -12.5199

Empty States Eigenvalues (eV), kp =   1 , spin =  2

-1.8218

Electronic Gap (eV) =    10.6981
...

Very clearly we can see the HOMO eigenvalue of -12.52 eV. Thus we have a KI IP of 12.52 eV. This compares extremely well to the experimental value of ~ 12.5 eV, and is a marked improvement on the PBE result of 7.95 eV (which we can obtain from the HOMO Eigenvalue in 01-initialization/01-dft_init_nspin2/04-dft_init_nspin2.cpo).

Meanwhile, the electron affinity (EA) corresponds to the negative of the energy of the LUMO (lowest unoccupied molecular orbital). From the same files we can see that the KI EA is 1.82 eV (cf. ~ 2.1 eV experiment, 6.17 eV PBE)

Tip

If you prefer working within python, you need not write a script to parse the contents of .cpo files in order to extract the IP and EA. Instead, koopmans will have generated a python-readable file ozone.pkl containing all of the important calculation data.

You can read these files like so:

"""Parse the .pkl output of a koopmans workflow."""

from koopmans import io

# Load the workflow object
wf = io.read('ozone.pkl')

# Access the results from the very last calculation
results = wf.calculations[-1].results

# Calculate the IP and EA
ip = -results['homo_energy']
ea = -results['lumo_energy']

# Print
print(f' IP = {ip:.2f} eV')
print(f' EA = {ea:.2f} eV')

Indeed, it is also possible to run the workflow from within python (rather than calling koopmans from the command line)

"""Run the ozone workflow with a python script."""

from koopmans.io import read

wf = read('ozone.json')
wf.run()

in which case you have immediate access to the workflow object wf rather than having to load a .pkl file.