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 "from_scratch": true,
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 ozone.json | tee ozone.out
Tip
In order to run in parallel, set the PARA_PREFIX environment variable to mpirun -np 4 (or similar)
The output
First, let us inspect the contents of ozone.out: after the header we can see there are a list of Quantum ESPRESSO calculations that have been run by koopmans. These come under three main headings.
Initialization
The first step in any Koopmans calculation is the initialization step. In this step we initialize the density and the variational orbitals.
14 Initialisation of density and variational orbitals
15 ==================================================
16 Running init/dft_init_nspin1... done
17 Running init/dft_init_nspin2_dummy... done
18 Running init/dft_init_nspin2... done
19 Overwriting the variational orbitals with Kohn-Sham orbitals
20 Copying the spin-up variational orbitals over to the spin-down channel
For this calculation we can see that koopmans has run three PBE calculations. These initialize the density with the PBE density. Indeed, from this point onwards in the calcuation 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 three PBE 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 three calculations is designed to avoid this; we first optimise the density constrained to be spin-unpolarized, and only once that density has been minimised 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 for this first step can all be found in the directory init/.
Calculating the screening parameters
The second step in the calculation involves the calculation of the screening parameters:
22 Calculating screening parameters
23 ================================
24 Running calc_alpha/ki... done
25
26 Orbital 1
27 ---------
28 Running calc_alpha/orbital_1/dft_n-1... done
29
30 Orbital 2
31 ---------
32 Running calc_alpha/orbital_2/dft_n-1... done
33
34 Orbital 3
35 ---------
36 Running calc_alpha/orbital_3/dft_n-1... done
37
38 Orbital 4
39 ---------
40 Running calc_alpha/orbital_4/dft_n-1... done
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:
61 Orbital 10
62 ----------
63 Running calc_alpha/orbital_10/pz_print... done
64 Running calc_alpha/orbital_10/dft_n+1_dummy... done
65 Running calc_alpha/orbital_10/dft_n+1... done
66
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
- pz_print and dft_n+1_dummy
preliminary calculations that generate files required by the subsequent constrained DFT calculation
- 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:
67
68 alpha
69 1 2 3 4 ... 7 8 9 10
70 0 0.60000 0.600000 0.600000 0.600000 ... 0.600000 0.600000 0.600000 0.60000
71 1 0.65568 0.727565 0.783855 0.663858 ... 0.729884 0.741895 0.779259 0.71739
72
73 [2 rows x 10 columns]
74
75 Delta E_i - epsilon_i (eV)
76 1 2 3 ... 8 9 10
77 0 -0.477272 -0.920449 -1.124074 ... -0.888568 -1.05013 0.711214
78
79 [1 rows x 10 columns]
80
81 Screening parameters have been determined but are not necessarily converged
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 - \varepsilon_{i}^\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 input and output Quantum ESPRESSO files for this step can be found in the directory calc_alpha/.
The final calculation
Having determined the screening parameters, the final KI calculation is now run:
83 Final KI calculation
84 ====================
85 Running final/ki_final... done
86
87 Workflow complete
The input and output Quantum ESPRESSO files for this step can be found in the directory final/.
Extracting the ionisation potential and electron affinity
Let’s now extract the KI ionisation potential and electron affinity for ozone from our calculation.
The ionisation potential (IP) corresponds to the negative of the energy of the HOMO (highest occupied molecular orbital). If you open 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 init/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 section in final/ki_final.cpo 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 final/ki_final.cpo in order to extract the IP and EA. Instead, koopmans will have generated a python-readable file ozone.kwf containing all of the important calculation data.
You can read these files like so:
from koopmans import io
# Load the workflow object
wf = io.read('ozone.kwf')
# 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)
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 in the .kwf file.