MICCoM School 2022 Ex#1 : The vertical ionization potential of SiH4¶
We are going to compute the vertical ionization potential (VIP) of the silane molecule using two levels of theory, namely density functional theory (DFT) with the semilocal exchange-correlation functional PBE, and full-frequency G\(_0\)W\(_0\).
The VIP is the minimum amount of energy that is needed to remove an electron from a molecule. When the removal does not take into account ion relaxation processes (change of molecular geometry), the ionization is said to be vertical.
The experimental value for the VIP of silane is tabulated in the NIST chemistry webbook: http://webbook.nist.gov/chemistry/form-ser/.
VIP of silane = \(12.3\) eV (experimental value)
~ We would like to compute the VIP of silane from first principles using DFT and GW. ~
Let’s see which level of theory performs better.
1.1 : DFT calculation¶
We are going to numerically solve the Kohn-Sham equations for valence electrons only. Electronic wavefunctions and single-particle energies are obtained by solving the secular equation:
\begin{equation} H^{\text{KS}}[\rho] \, \psi_n = \varepsilon_n^{\text{KS}} \psi_n \end{equation}
The Kohn-Sham Hamiltonian of \(N\) interacting electrons includes the following terms:
\begin{equation} H^{\text{KS}}[\rho] = T + V_{\text{ion}} + V_{\text{Hartree}}[\rho] + V_{\text{xc}}[\rho] \end{equation}
where \(T\) is the kinetic energy, \(V_{\text{ion}}\) is the electron-ion potential, \(V_{\text{Hartree}}\) is the Hartree potential, and \(V_{\text{xc}}\) is the exchange-correlation potential.
Note that the secular equation needs to be solved self-consistenly since the operator that we are diagonalizing depends on the eigenvectors (\(\rho = \sum_n |\psi_n|^2\)).
We are going to use plane-waves as basis set, which means that an electronic wavefunction \(\psi_n\) is obtained by computing the coefficients \(C_n\) of the plane-waves expansion:
\begin{equation} \psi_n(\mathbf{r}) = \sum_{\mathbf{G}} C_n(\mathbf{G}) e^{i\mathbf{G}\cdot\mathbf{r}} \end{equation}
In this way the kinetic energy of a wavefunction is \(\sum_{\mathbf{G}} \left| C_n(\mathbf{G}) \right|^2 \mathbf{G}^2\). The total number of plane-waves is determined once we specify the size of the simulation box and the kinetic energy cutoff.
The electron-ion potetial \(V_{\text{ion}}\) is solved separately for each given atomic species. Pseudopotential files containing this information (previously computed) are typically loaded at the beginning of the calculation. Pseudopotentials can be downloaded from several websites, for instance:
SG15 : http://www.quantum-simulation.org/potentials/sg15_oncv/
Quantum ESPRESSO : https://www.quantum-espresso.org/pseudopotentials/
In this tutorial we will be using only norm-conserving pseudopotentials, i.e., no ultrasoft or projector augmented wave (PAW).
The Hartree potential is determined self-consistely as \(V_{\text{Hartree}}(\mathbf{r}) = \int d\mathbf{r^\prime} \rho(\mathbf{r^\prime}) \frac{1}{|\mathbf{r}-\mathbf{r^\prime}|}\).
The exact form for the exchange and correlation potential is unknown. Many approximations are available. The user needs to make this decision. In this tutorial we will be using the PBE exchange-correlation functional, which belongs to the family of semilocal generalized gradient approximation (GGA).
If we interpret \(\varepsilon^{\text{KS}}_n\) as electronic excitation energies (which is physically incorrect but people still do this as an approximation), the VIP will simply be the energy of the highest occupied molecular orbital (HOMO) obtained in a DFT calculation.
\begin{equation} \text{VIP}^{\text{DFT}} = 0 - \varepsilon^{\text{KS}}_{\text{HOMO}} \end{equation}
Question: Where did the core electrons go?
Question: How many valence electrons do we have in SiH4?
Question: Plane-waves satisfy periodic boundary conditions, how can we simulate an isolated molecule?
The following files are needed in order to compute the electronic structure of SiH\(_4\) with DFT:
Si_ONCV_PBE-1.2.upf: pseudopotential file for SiH_ONCV_PBE-1.2.upf: pseudopotential file for Hpw.in: input file for the DFT calculation (executablepw.x)
[1]:
# download the pseudopotential files
!wget -N -q http://www.quantum-simulation.org/potentials/sg15_oncv/upf/Si_ONCV_PBE-1.2.upf
!wget -N -q http://www.quantum-simulation.org/potentials/sg15_oncv/upf/H_ONCV_PBE-1.2.upf
# download the input file
!wget -N -q https://west-code.org/doc/training/silane/pw.in
Let’s give a quick look at the input for DFT (description of the input variables for pw.x can be found here: https://www.quantum-espresso.org/Doc/INPUT_PW.html).
[2]:
!cat pw.in
&control
calculation = 'scf'
restart_mode = 'from_scratch'
pseudo_dir = './'
outdir = './'
prefix = 'silane'
wf_collect = .TRUE.
/
&system
ibrav = 1
celldm(1) = 20
nat = 5
ntyp = 2
ecutwfc = 25.0
nbnd = 10
assume_isolated ='mp'
/
&electrons
diago_full_acc = .TRUE.
/
ATOMIC_SPECIES
Si 28.0855 Si_ONCV_PBE-1.2.upf
H 1.00794 H_ONCV_PBE-1.2.upf
ATOMIC_POSITIONS bohr
Si 10.000000 10.000000 10.000000
H 11.614581 11.614581 11.614581
H 8.385418 8.385418 11.614581
H 8.385418 11.614581 8.385418
H 11.614581 8.385418 8.385418
K_POINTS {gamma}
We have placed the silane molecule in a cubic box (edge = 20 bohr), and specified a kinetic energy cutoff of 25 Ry.
We can run the DFT calculation invoking the executable pw.x on 4 cores. This calculation should finish in just a few seconds.
[3]:
!mpirun -n 4 pw.x -i pw.in > pw.out
Let’s extract the energy levels of the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) predicted by DFT.
[4]:
!grep 'highest occupied, lowest unoccupied level' pw.out
highest occupied, lowest unoccupied level (ev): -8.2251 -0.4647
HOMO : \(-8.2251\) eV
LUMO : \(-0.4647\) eV
Let’s plot this result.
[5]:
import matplotlib.pyplot as plt
# VIP
y = {}
y['exp'] = [12.3] # experimental
y['dft'] = [8.2251] # DFT
# colors of the plot
c = {}
c['exp'] = 'black'
c['dft'] = 'red'
# plot
x = list(range(1, len(y)+1))
labels = y.keys()
fig, ax = plt.subplots(1, 1)
for counter, label in enumerate(labels) :
for a in y[label] :
ax.hlines(a, x[counter]-0.25, x[counter]+0.25, color=c[label])
plt.xticks(x, labels)
plt.ylabel('Energy (eV)')
plt.title('VIP')
plt.show()
Question: Is this level of agreement ok?
1.2 : Full-frequency G0W0 calculation¶
The GW method in the framework of many-body perturbation theory (MBPT) can improve the result of DFT by computing single-particle excitations as the eigenvalues of the quasiparticle (QP) Hamiltonian:
\begin{equation} H^{\text{QP}}[\rho] = T + V_{\text{ion}} + V_{\text{Hartree}}[\rho] + \Sigma \end{equation}
Note that \(H^{\text{QP}}\) can be obtained from \(H^{\text{KS}}\) by replacing the exchange and correlation potential with the self-energy \(\Sigma\). In the GW approximation the electron self-energy can be written as :
\begin{equation} \Sigma(\mathbf{r},\mathbf{r^\prime};E) = i \int \frac{d\omega}{2\pi} G(\mathbf{r},\mathbf{r^\prime};E-\omega) W(\mathbf{r},\mathbf{r^\prime};\omega) \end{equation}
\(G\) is the Green’s function and \(W=\epsilon^{-1}v\) is the screened Coulomb potential. Because both \(G\) and \(W\) depend self-consistently on the solution of the QP Hamiltonian, typically the eigenvalues of the operator of \(H^{\text{QP}}\) are obtained using the eigenvalues and eigenvectors of \(H^{\text{KS}}\) (computed previously in DFT) and using first order perturbation theory:
\begin{equation} E^{\text{QP}}_n = \varepsilon^{\text{KS}}_n + \left\langle \psi_n \right | \Sigma(E^{\text{QP}}_n) - V_{\text{xc}}\left| \psi_n \right\rangle \end{equation}
Following the workflow discussed in J. Chem. Theory Comput. 11, 2680 (2015), we use the DFT output to:
Compute the dielectric screening : executable
wstat.xCompute the quasiparticle energies : executable
wfreq.x
1.2.1 : Calculation of dielectric screening¶
The wstat executable will read the output of DFT and will compute the static dielectric screening adopting the projective dielectric eigendecomposition (PDEP) technique. An iterative diagonalization of the dielectric matrix is performed:
\begin{equation} \epsilon_{\mathbf{G},\mathbf{G^\prime}}(\omega=0) = \sum_{i=1}^{\texttt{n_pdep_eigen}} \, \phi_i(\mathbf{G}) \, \lambda_i \phi^\ast_i(\mathbf{G^\prime}) \end{equation}
where \(\phi_i\) and \(\lambda_i\) are the eigenvectors and eigenvalues of the dielectric matrix (extressed here in terms of plane-waves). The eigenvectors \(\phi_i\) are used as the basis set for the subsequent GW calculation.
The diagonalization is started with random potentials and then iterated using the Davidson algorithm:
Let’s download and give a quick look at the input for wstat.x (description of the input variables for wstat.x can be found here: https://west-code.org/doc/West/latest/manual.html#wstat-control).
[6]:
!wget -N -q https://west-code.org/doc/training/silane/wstat.in
!cat wstat.in
input_west:
qe_prefix: silane
west_prefix: silane
outdir: ./
wstat_control:
wstat_calculation: S
n_pdep_eigen: 50
The variable silane instructs the code to find the output of DFT in the folder silane.save. We are computing 50 eigenpotentials \(\phi_i\), i.e., 50 eigenvectors of the static dielectric matrix. This step uses the occupied single-particle states and energies obtained in the DFT calculation. Unoccupied states are not needed.
We can run the PDEP calculation invoking the executable wstat.x on 16 cores. This calculation should finish in 1 minute.
[7]:
!mpirun -n 16 wstat.x -i wstat.in > wstat.out
We load the output (in JSON format) of the wstat.x calculation, which can be found in the silane.wstat.save folder.
[8]:
import json
# read data from JSON file
with open('silane.wstat.save/wstat.json', 'r') as file :
data = json.load(file)
# print the data
print(json.dumps(data, indent=2))
{
"runjob": {
"startdate": "11Oct2022",
"starttime": " 9:32:59",
"completed": true,
"endtime": " 9:33:15",
"enddate": "11Oct2022"
},
"software": {
"package": "WEST",
"program": "WSTAT",
"version": "5.2.0",
"westgit": "v5.2.0-10-g918d532",
"website": "https://west-code.org",
"citation": "M. Govoni et al., J. Chem. Theory Comput. 11, 2680 (2015).",
"qeversion": "7.1"
},
"config": {
"io": {
"islittleendian": true
}
},
"parallel": {
"nranks": 16,
"nimage": 1,
"npool": 1,
"nbgrp": 1,
"nrg": 16,
"nproc": 16,
"nthreads": 1
},
"input": {
"input_west": {
"qe_prefix": "silane",
"west_prefix": "silane",
"outdir": "./"
},
"wstat_control": {
"wstat_calculation": "S",
"n_pdep_eigen": 50,
"n_pdep_times": 4,
"n_pdep_maxiter": 100,
"n_dfpt_maxiter": 250,
"n_pdep_read_from_file": 0,
"n_steps_write_restart": 1,
"trev_pdep": 0.001,
"trev_pdep_rel": 0.1,
"tr2_dfpt": 1e-12,
"l_kinetic_only": false,
"l_minimize_exx_if_active": false,
"l_use_ecutrho": false,
"qlist": [
1
]
},
"server_control": {
"document": "{}"
}
},
"system": {
"basis": {
"npw": {
"proc": [
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0
],
"min": 0,
"max": 0,
"sum": 0
},
"ngm": {
"proc": [
4220,
4221,
4219,
4221,
4221,
4221,
4220,
4221,
4221,
4220,
4219,
4219,
4219,
4220,
4220,
4220
],
"min": 4219,
"max": 4221,
"sum": 67522
},
"gamma_only": true,
"ecutwfc:ry": 25.0,
"ecutrho:ry": 100.0
},
"cell": {
"units": "a.u.",
"omega": 8000.0,
"a1": [
20.0,
0.0,
0.0
],
"a2": [
0.0,
20.0,
0.0
],
"a3": [
0.0,
0.0,
20.0
],
"b1": [
0.3141592653589793,
0.0,
0.0
],
"b2": [
0.0,
0.3141592653589793,
0.0
],
"b3": [
0.0,
0.0,
0.3141592653589793
],
"alat": 20.0,
"tpiba": 0.3141592653589793
},
"electron": {
"nbnd": 10,
"nkstot": 1,
"nspin": 1,
"nelec": 8.0,
"npol": 1,
"lsda": false,
"noncolin": false,
"lspinorb": false
},
"3dfft": {
"s": [
64,
64,
64
],
"p": [
64,
64,
64
]
},
"bzsamp": {
"k": [
{
"id": 1,
"crystcoord": [
0.0,
0.0,
0.0
]
}
]
}
},
"memory": {
"units": "Mb",
"evc": 0.08026123046875,
"nlpp": 0.12841796875,
"rhor": 0.25,
"rhog": 0.032196044921875,
"gshells": 0.005401611328125,
"dvg": 1.605224609375,
"dng": 1.605224609375,
"hr_distr": 0.30517578125,
"vr_distr": 0.30517578125,
"dvpsi": 0.0321044921875,
"dpsi": 0.0321044921875
},
"exec": {
"ndav": 4,
"davitr": [
{
"dav_iter": -1,
"ev": [
-1.2605591279261816,
-1.1812479933930013,
-1.176354788237097,
-1.1712738214218055,
-0.8162645071961094,
-0.8100056551237361,
-0.8037713586928801,
-0.6207256375128,
-0.6139258457189167,
-0.4504239969889595,
-0.44258236373617965,
-0.42285739065633154,
-0.4161871352235839,
-0.4086761256499228,
-0.391144269933381,
-0.35154650952224475,
-0.22034435672551075,
-0.21296061444937364,
-0.2042616451193836,
-0.1661200380207801,
-0.15969273612456536,
-0.1553415063345141,
-0.15237493663785076,
-0.14389848531187638,
-0.1365968664788392,
-0.12511600591857588,
-0.11539384948929213,
-0.10459429582946164,
-0.10042992216715237,
-0.09870244795049053,
-0.0933238874905026,
-0.09167788000062861,
-0.08518312554603533,
-0.07958471780926124,
-0.0776708939312429,
-0.07307036957745469,
-0.06891873838631551,
-0.060148023225904235,
-0.05778630984761288,
-0.05592401268009498,
-0.05092294690647727,
-0.04951946848860927,
-0.0476517799892079,
-0.04435874746873535,
-0.04317218464886519,
-0.040665473458473905,
-0.03778191195635325,
-0.03576602983411223,
-0.033769213117953605,
-0.015542819711099503
],
"conv": [
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false
],
"notcnv": 50,
"time_elap:sec": 3.51153302192688,
"time_elap:hum": "03.5s",
"time_iter:sec": 3.5115301609039307,
"time_iter:hum": "03.5s",
"sternop_ncalls": 1011,
"dfpt_tr2": 1e-06,
"dfpt_dim": 50,
"diago_dim": 50
},
{
"dav_iter": 1,
"ev": [
-1.2754684424024114,
-1.191995341184158,
-1.19198786300694,
-1.1918721110320232,
-0.82442753885941,
-0.8243639321874088,
-0.8242997840716799,
-0.6361356957558751,
-0.6296568274548767,
-0.6296413545568273,
-0.5006636696332897,
-0.5005795189579708,
-0.500474475279866,
-0.43005011856959735,
-0.4300320948062763,
-0.4299825294011507,
-0.23235168084292807,
-0.2322517154513132,
-0.23209974647635992,
-0.1829326559805843,
-0.1827335781160475,
-0.18229299877133223,
-0.17784527910811446,
-0.17672947835944156,
-0.17667590632711627,
-0.14560918460051153,
-0.14528730945197754,
-0.14486852686919827,
-0.12075118485366688,
-0.11910453951818724,
-0.1187348760281049,
-0.11809073390359559,
-0.11578278244178404,
-0.1151722894313686,
-0.11430674937951284,
-0.11350388815667996,
-0.11006967778786642,
-0.09208669672138343,
-0.09144853305507969,
-0.0897120173120519,
-0.07827072092186087,
-0.07558536284116134,
-0.07372706056817617,
-0.07142441109651232,
-0.06943098454262646,
-0.06809364629565738,
-0.06268888490486112,
-0.062021508729807016,
-0.06072976596578535,
-0.060135374072292344
],
"conv": [
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false
],
"notcnv": 50,
"time_elap:sec": 7.870642900466919,
"time_elap:hum": "07.8s",
"time_iter:sec": 4.215182065963745,
"time_iter:hum": "04.2s",
"sternop_ncalls": 1724,
"dfpt_tr2": 1e-12,
"dfpt_dim": 50,
"diago_dim": 100
},
{
"dav_iter": 2,
"ev": [
-1.2754720666724284,
-1.191998889630293,
-1.1919914251534018,
-1.1918750476969548,
-0.8244331949249563,
-0.824371989387179,
-0.824307057750549,
-0.6363480461593288,
-0.629663885414068,
-0.6296546797640005,
-0.5007531759617246,
-0.5007191828843528,
-0.5006632705684765,
-0.4301010435507772,
-0.43008785306001956,
-0.4300639248433781,
-0.23245848670842764,
-0.232454446880072,
-0.23244933014889632,
-0.18328587917955322,
-0.1832547553729702,
-0.18324633004836113,
-0.17843114260158843,
-0.17754830176814526,
-0.17754103307355587,
-0.14595932159715894,
-0.145955923078163,
-0.1459411471756543,
-0.1226019592550607,
-0.12013930359811296,
-0.12013399805600997,
-0.12012057383845803,
-0.11636716556326777,
-0.11636290017226038,
-0.11529987441574399,
-0.11529799856558509,
-0.11525428280816115,
-0.09406113606695948,
-0.09404985091462528,
-0.09401167062537157,
-0.07992389685332407,
-0.07987094865761875,
-0.07975914970189521,
-0.07468254818337605,
-0.07301849155829937,
-0.07296297688645226,
-0.06562513611180669,
-0.06555066994951836,
-0.06528330422704633,
-0.06306269344176937
],
"conv": [
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
false,
false,
false,
false,
false,
true,
false,
true,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false,
false
],
"notcnv": 21,
"time_elap:sec": 12.536376953125,
"time_elap:hum": "12.5s",
"time_iter:sec": 4.1379029750823975,
"time_iter:hum": "04.1s",
"sternop_ncalls": 1692,
"dfpt_tr2": 1e-12,
"dfpt_dim": 50,
"diago_dim": 150
},
{
"dav_iter": 3,
"ev": [
-1.2754720668029091,
-1.1919988897734295,
-1.191991425327886,
-1.1918750479082716,
-0.8244331954209921,
-0.8243719902168194,
-0.8243070584764758,
-0.6363480751501617,
-0.6296638864684181,
-0.62965468160695,
-0.5007531991331637,
-0.5007192212407726,
-0.5006633284291045,
-0.4301010610680996,
-0.4300878999594423,
-0.43006394104361023,
-0.2324590857953143,
-0.23245478344153164,
-0.2324498222228378,
-0.18328843644430254,
-0.18325666138108107,
-0.1832473382383492,
-0.1784332377577128,
-0.1775504004743457,
-0.17754325566955023,
-0.14596095608969584,
-0.14595895912881743,
-0.1459466107908167,
-0.12261797473339454,
-0.12014634880115369,
-0.12014593756057011,
-0.12014040871072855,
-0.11637378892529171,
-0.11636989698973867,
-0.11531375801416825,
-0.11530875850124506,
-0.11530541512826781,
-0.09409890231649157,
-0.09409787245061647,
-0.0940955088831779,
-0.07996894456596207,
-0.07996593163736887,
-0.07996545801532576,
-0.07480835355986729,
-0.07311589574225359,
-0.07311484170817655,
-0.06579115177015801,
-0.06578527255722108,
-0.06578041128186739,
-0.06314475074029666
],
"conv": [
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true,
true
],
"notcnv": 0,
"time_elap:sec": 14.89906096458435,
"time_elap:hum": "14.8s",
"time_iter:sec": 1.7204980850219727,
"time_iter:hum": "01.7s",
"sternop_ncalls": 708,
"dfpt_tr2": 1e-12,
"dfpt_dim": 21,
"diago_dim": 171
}
]
},
"timing": {
"WSTAT": {
"cpu:sec": 13.765601,
"cpu:hum": "13.7s",
"wall:sec": 16.18815803527832,
"wall:hum": "16.1s",
"nocalls": 2
},
"wstat_readin": {
"cpu:sec": 0.4568550000000001,
"cpu:hum": "00.4s",
"wall:sec": 0.5870640277862549,
"wall:hum": "00.5s",
"nocalls": 1
},
"fetch_input": {
"cpu:sec": 0.029162000000000132,
"cpu:hum": "< 00.1s",
"wall:sec": 0.0364842414855957,
"wall:hum": "< 00.1s",
"nocalls": 2
},
"fft": {
"cpu:sec": 0.05137000000000036,
"cpu:hum": "< 00.1s",
"wall:sec": 0.09902167320251465,
"wall:hum": "< 00.1s",
"nocalls": 9
},
"fft_scatter": {
"cpu:sec": 3.6312919999999504,
"cpu:hum": "03.6s",
"wall:sec": 3.8297791481018066,
"wall:hum": "03.8s",
"nocalls": 20889
},
"init_vloc": {
"cpu:sec": 0.05437700000000012,
"cpu:hum": "< 00.1s",
"wall:sec": 0.05521893501281738,
"wall:hum": "< 00.1s",
"nocalls": 2
},
"init_us_1": {
"cpu:sec": 0.009448999999999819,
"cpu:hum": "< 00.1s",
"wall:sec": 0.013067007064819336,
"wall:hum": "< 00.1s",
"nocalls": 2
},
"v_of_rho": {
"cpu:sec": 0.020774999999999988,
"cpu:hum": "< 00.1s",
"wall:sec": 0.032330989837646484,
"wall:hum": "< 00.1s",
"nocalls": 1
},
"v_xc": {
"cpu:sec": 0.018656999999999924,
"cpu:hum": "< 00.1s",
"wall:sec": 0.027471065521240234,
"wall:hum": "< 00.1s",
"nocalls": 1
},
"v_h": {
"cpu:sec": 0.0016580000000001593,
"cpu:hum": "< 00.1s",
"wall:sec": 0.002373933792114258,
"wall:hum": "< 00.1s",
"nocalls": 1
},
"do_setup": {
"cpu:sec": 0.05048600000000003,
"cpu:hum": "< 00.1s",
"wall:sec": 0.0683588981628418,
"wall:hum": "< 00.1s",
"nocalls": 1
},
"init_pw_ar": {
"cpu:sec": 0.04259199999999996,
"cpu:hum": "< 00.1s",
"wall:sec": 0.05210399627685547,
"wall:hum": "< 00.1s",
"nocalls": 1
},
"hinit0": {
"cpu:sec": 0.033872000000000124,
"cpu:hum": "< 00.1s",
"wall:sec": 0.034713029861450195,
"wall:hum": "< 00.1s",
"nocalls": 1
},
"chidiago": {
"cpu:sec": 13.244492,
"cpu:hum": "13.2s",
"wall:sec": 15.499743938446045,
"wall:hum": "15.4s",
"nocalls": 1
},
"sqvc_init": {
"cpu:sec": 0.016258000000000106,
"cpu:hum": "< 00.1s",
"wall:sec": 0.016334056854248047,
"wall:hum": "< 00.1s",
"nocalls": 1
},
"randomize": {
"cpu:sec": 0.2202679999999999,
"cpu:hum": "00.2s",
"wall:sec": 0.22018694877624512,
"wall:hum": "00.2s",
"nocalls": 1
},
"paramgs": {
"cpu:sec": 0.059492000000000544,
"cpu:hum": "< 00.1s",
"wall:sec": 0.06182503700256348,
"wall:hum": "< 00.1s",
"nocalls": 5
},
"dfpt": {
"cpu:sec": 12.170046,
"cpu:hum": "12.1s",
"wall:sec": 12.78049111366272,
"wall:hum": "12.7s",
"nocalls": 5
},
"init_us_2": {
"cpu:sec": 0.0030950000000018463,
"cpu:hum": "< 00.1s",
"wall:sec": 0.004007816314697266,
"wall:hum": "< 00.1s",
"nocalls": 5
},
"init_us_2:cp": {
"cpu:sec": 0.003003000000001421,
"cpu:hum": "< 00.1s",
"wall:sec": 0.003938198089599609,
"wall:hum": "< 00.1s",
"nocalls": 5
},
"fftw": {
"cpu:sec": 7.157574000000035,
"cpu:hum": "07.1s",
"wall:sec": 7.495359897613525,
"wall:hum": "07.4s",
"nocalls": 20880
},
"alphapc": {
"cpu:sec": 0.04244199999999232,
"cpu:hum": "< 00.1s",
"wall:sec": 0.04595136642456055,
"wall:hum": "< 00.1s",
"nocalls": 221
},
"wstat_restar": {
"cpu:sec": 0.5095869999999998,
"cpu:hum": "00.5s",
"wall:sec": 1.5846247673034668,
"wall:hum": "01.5s",
"nocalls": 4
},
"pdep_write": {
"cpu:sec": 0.48690699999999687,
"cpu:hum": "00.4s",
"wall:sec": 1.6841545104980469,
"wall:hum": "01.6s",
"nocalls": 750
},
"linstern": {
"cpu:sec": 11.062415000000001,
"cpu:hum": "11.0s",
"wall:sec": 11.60523533821106,
"wall:hum": "11.6s",
"nocalls": 171
},
"stern": {
"cpu:sec": 9.671373999999979,
"cpu:hum": "09.6s",
"wall:sec": 10.131323099136353,
"wall:hum": "10.1s",
"nocalls": 5135
},
"h_psi": {
"cpu:sec": 9.006320000000056,
"cpu:hum": "09.0s",
"wall:sec": 9.43653130531311,
"wall:hum": "09.4s",
"nocalls": 5135
},
"h_psi:pot": {
"cpu:sec": 8.84657099999999,
"cpu:hum": "08.8s",
"wall:sec": 9.262455224990845,
"wall:hum": "09.2s",
"nocalls": 5135
},
"vloc_psi": {
"cpu:sec": 7.466788999999979,
"cpu:hum": "07.4s",
"wall:sec": 7.8167266845703125,
"wall:hum": "07.8s",
"nocalls": 5135
},
"h_psi:calbec": {
"cpu:sec": 0.8840789999999679,
"cpu:hum": "00.8s",
"wall:sec": 0.9279425144195557,
"wall:hum": "00.9s",
"nocalls": 5135
},
"calbec": {
"cpu:sec": 0.8083199999999717,
"cpu:hum": "00.8s",
"wall:sec": 0.8619239330291748,
"wall:hum": "00.8s",
"nocalls": 5135
},
"add_vuspsi": {
"cpu:sec": 0.36807999999999197,
"cpu:hum": "00.3s",
"wall:sec": 0.3934800624847412,
"wall:hum": "00.3s",
"nocalls": 5135
},
"alphapv": {
"cpu:sec": 0.5223740000000312,
"cpu:hum": "00.5s",
"wall:sec": 0.5605552196502686,
"wall:hum": "00.5s",
"nocalls": 5135
},
"build_hr": {
"cpu:sec": 0.015184000000001419,
"cpu:hum": "< 00.1s",
"wall:sec": 0.015445947647094727,
"wall:hum": "< 00.1s",
"nocalls": 4
},
"diagox": {
"cpu:sec": 0.03941000000000017,
"cpu:hum": "< 00.1s",
"wall:sec": 0.05038881301879883,
"wall:hum": "< 00.1s",
"nocalls": 4
},
"redistr_vr": {
"cpu:sec": 0.0027069999999991268,
"cpu:hum": "< 00.1s",
"wall:sec": 0.0043947696685791016,
"wall:hum": "< 00.1s",
"nocalls": 3
},
"update_vr": {
"cpu:sec": 0.029883000000001658,
"cpu:hum": "< 00.1s",
"wall:sec": 0.032776832580566406,
"wall:hum": "< 00.1s",
"nocalls": 3
},
"chidiago:las": {
"cpu:sec": 1.000000000139778e-05,
"cpu:hum": "< 00.1s",
"wall:sec": 7.152557373046875e-06,
"wall:hum": "< 00.1s",
"nocalls": 1
},
"refresh_vr": {
"cpu:sec": 0.008898999999999546,
"cpu:hum": "< 00.1s",
"wall:sec": 0.00884103775024414,
"wall:hum": "< 00.1s",
"nocalls": 1
},
"pdep_db": {
"cpu:sec": 0.033900999999998405,
"cpu:hum": "< 00.1s",
"wall:sec": 0.22748708724975586,
"wall:hum": "00.2s",
"nocalls": 1
}
}
}
Now we can plot the 50 eigenvalues \(\lambda_i\) of the dielectric matrix.
[9]:
import numpy as np
import matplotlib.pyplot as plt
# y : eigenvalues of the dielectic matrix
y = data['exec']['davitr'][-1]['ev']
# x : index of eigenvalue
x = np.arange(1, len(y)+1, 1)
print('y : ', y)
print('x : ', x)
# plot
fig, ax = plt.subplots(1, 1)
ax.plot(x, y, 'ro-')
plt.xlabel('i')
plt.ylabel('Eigenvalue')
plt.title('PDEP')
plt.show()
y : [-1.2754720668029091, -1.1919988897734295, -1.191991425327886, -1.1918750479082716, -0.8244331954209921, -0.8243719902168194, -0.8243070584764758, -0.6363480751501617, -0.6296638864684181, -0.62965468160695, -0.5007531991331637, -0.5007192212407726, -0.5006633284291045, -0.4301010610680996, -0.4300878999594423, -0.43006394104361023, -0.2324590857953143, -0.23245478344153164, -0.2324498222228378, -0.18328843644430254, -0.18325666138108107, -0.1832473382383492, -0.1784332377577128, -0.1775504004743457, -0.17754325566955023, -0.14596095608969584, -0.14595895912881743, -0.1459466107908167, -0.12261797473339454, -0.12014634880115369, -0.12014593756057011, -0.12014040871072855, -0.11637378892529171, -0.11636989698973867, -0.11531375801416825, -0.11530875850124506, -0.11530541512826781, -0.09409890231649157, -0.09409787245061647, -0.0940955088831779, -0.07996894456596207, -0.07996593163736887, -0.07996545801532576, -0.07480835355986729, -0.07311589574225359, -0.07311484170817655, -0.06579115177015801, -0.06578527255722108, -0.06578041128186739, -0.06314475074029666]
x : [ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
49 50]
As we can see, the eigenvalues of the dielectric matrix decay to zero. The number of PDEPs, n_pdep_eigen, controls the size of the basis set used in the GW calculation. n_pdep_eigen is a parameter of the simulation that needs to be converged. Typically, n_pdep_eigen is set to multiple of the number of electrons in the system.
Question: How do we compute 70 eigenpotentials?
1.2.2 : Calculation of GW electronic structure¶
The GW electronic structure is computed by computing the dielectric screening at multipole frequencies with Lanczos iterations and by treating the frequency integration of the correlation part of the self energy with the contour deformation techinique.
Let’s download and give a quick look at the input for wfreq.x (description of the input variables for wfreq.x can be found here: https://west-code.org/doc/West/latest/manual.html#wfreq-control).
[10]:
!wget -N -q https://west-code.org/doc/training/silane/wfreq.in
!cat wfreq.in
input_west:
qe_prefix: silane
west_prefix: silane
outdir: ./
wstat_control:
wstat_calculation: S
n_pdep_eigen: 50
wfreq_control:
wfreq_calculation: XWGQ
n_pdep_eigen_to_use: 50
qp_bandrange: [1,5]
n_refreq: 300
ecut_refreq: 2.0
The keyword XWGQ instructs the code to compute the exact exchange operator (X), the screened Coulomb operator (W), the Green’s function (G), and the quasiparticle energies (Q).
In this case we are using the 50 eigenpotentials \(\phi_i\) obtained by wstat.x.
We are going to compute the quasiparticle energies of the states 1, 2, 3, 4 (HOMO), and 5 (LUMO). Remember that in this example (spin-unpolarized) each state can host two electrons.
The frequency parameters control the contour deformation integration. The frequency dependence of the dielectric screening is sampled in the energy window \([0,2]\) Ry using 300 points.
We can run the GW calculation invoking the executable wfreq.x on 16 cores. This calculation should finish in 1 minute.
[11]:
!mpirun -n 16 wfreq.x -i wfreq.in > wfreq.out
We load the output (in JSON format) of the wfreq.x calculation, which can be found in the silane.wfreq.save folder.
[12]:
import json
# read data from JSON file
with open('silane.wfreq.save/wfreq.json', 'r') as file :
data = json.load(file)
# print the data
print(json.dumps(data, indent=2))
{
"runjob": {
"startdate": "11Oct2022",
"starttime": " 9:33:34",
"completed": true,
"endtime": " 9:33:43",
"enddate": "11Oct2022"
},
"software": {
"package": "WEST",
"program": "WFREQ",
"version": "5.2.0",
"westgit": "v5.2.0-10-g918d532",
"website": "https://west-code.org",
"citation": "M. Govoni et al., J. Chem. Theory Comput. 11, 2680 (2015).",
"qeversion": "7.1"
},
"config": {
"io": {
"islittleendian": true
}
},
"parallel": {
"nranks": 16,
"nimage": 1,
"npool": 1,
"nbgrp": 1,
"nrg": 16,
"nproc": 16,
"nthreads": 1
},
"input": {
"input_west": {
"qe_prefix": "silane",
"west_prefix": "silane",
"outdir": "./"
},
"wstat_control": {
"wstat_calculation": "S",
"n_pdep_eigen": 50,
"n_pdep_times": 4,
"n_pdep_maxiter": 100,
"n_dfpt_maxiter": 250,
"n_pdep_read_from_file": 0,
"n_steps_write_restart": 1,
"trev_pdep": 0.001,
"trev_pdep_rel": 0.1,
"tr2_dfpt": 1e-12,
"l_kinetic_only": false,
"l_minimize_exx_if_active": false,
"l_use_ecutrho": false,
"qlist": [
1
]
},
"wfreq_control": {
"wfreq_calculation": "XWGQ",
"n_pdep_eigen_to_use": 50,
"qp_bandrange": [
1,
5
],
"qp_bands": [
0
],
"macropol_calculation": "N",
"n_lanczos": 30,
"n_imfreq": 128,
"n_refreq": 300,
"ecut_imfreq": 100.0,
"ecut_refreq": 2.0,
"wfreq_eta": 0.0036749308824745378,
"n_secant_maxiter": 21,
"trev_secant": 0.0036749308824745378,
"l_enable_lanczos": true,
"l_enable_off_diagonal": false,
"o_restart_time": 0.0,
"ecut_spectralf": [
-2.0,
1.0
],
"n_spectralf": 204
}
},
"system": {
"basis": {
"npw": {
"proc": [
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0
],
"min": 0,
"max": 0,
"sum": 0
},
"ngm": {
"proc": [
4220,
4221,
4219,
4221,
4221,
4221,
4220,
4221,
4221,
4220,
4219,
4219,
4219,
4220,
4220,
4220
],
"min": 4219,
"max": 4221,
"sum": 67522
},
"gamma_only": true,
"ecutwfc:ry": 25.0,
"ecutrho:ry": 100.0
},
"cell": {
"units": "a.u.",
"omega": 8000.0,
"a1": [
20.0,
0.0,
0.0
],
"a2": [
0.0,
20.0,
0.0
],
"a3": [
0.0,
0.0,
20.0
],
"b1": [
0.3141592653589793,
0.0,
0.0
],
"b2": [
0.0,
0.3141592653589793,
0.0
],
"b3": [
0.0,
0.0,
0.3141592653589793
],
"alat": 20.0,
"tpiba": 0.3141592653589793
},
"electron": {
"nbnd": 10,
"nkstot": 1,
"nspin": 1,
"nelec": 8.0,
"npol": 1,
"lsda": false,
"noncolin": false,
"lspinorb": false
},
"3dfft": {
"s": [
64,
64,
64
],
"p": [
64,
64,
64
]
},
"bzsamp": {
"k": [
{
"id": 1,
"crystcoord": [
0.0,
0.0,
0.0
]
}
]
}
},
"exec": {
"Q": {
"secitr": 3,
"en": [
{
"ksb": [
1,
1,
1
],
"ein": [
-13.230981270352434,
-16.26615213933269,
-16.103461398940283
],
"eout": [
-16.26615213933269,
-16.103461398940283,
-16.126442757174104
],
"sc_ein": [
[
1.7659827273250182,
0.014583989091543714
],
[
3.581007935231201,
0.07514251063272293
],
[
3.4409541840918707,
0.11069082478558719
]
]
},
{
"ksb": [
1,
1,
2
],
"ein": [
-8.226552436785038,
-12.145989533047741,
-12.038568029539881
],
"eout": [
-12.145989533047741,
-12.038568029539881,
-12.041375786118818
],
"sc_ein": [
[
0.010466906176382348,
0.003272239246325086
],
[
0.7305227660066568,
0.00892881625011567
],
[
0.7073751375953274,
0.008794218426373911
]
]
},
{
"ksb": [
1,
1,
3
],
"ein": [
-8.22537802992713,
-12.141823660074568,
-12.034279625593381
],
"eout": [
-12.141823660074568,
-12.034279625593381,
-12.037079367870797
],
"sc_ein": [
[
0.013781153568594457,
0.003272405433181099
],
[
0.7333805411639928,
0.008931013424447288
],
[
0.710217882884513,
0.008796293848555962
]
]
},
{
"ksb": [
1,
1,
4
],
"ein": [
-8.225073262212328,
-12.141850045890173,
-12.034301905490201
],
"eout": [
-12.141850045890173,
-12.034301905490201,
-12.037101945463592
],
"sc_ein": [
[
0.013237978473109088,
0.0032729463236466855
],
[
0.7329844903483206,
0.008932371065598663
],
[
0.7098182772432421,
0.00879763962964248
]
]
},
{
"ksb": [
1,
1,
5
],
"ein": [
-0.4646702125746965,
0.664173458065683,
0.663067500001202
],
"eout": [
0.664173458065683,
0.663067500001202,
0.663067500001202
],
"sc_ein": [
[
-0.3481319149046877,
-1.058286063240667e-05
],
[
-0.3718350341861773,
-0.00011661816924339793
],
[
-0.37181078103599796,
-0.00011661816924339793
]
]
}
]
}
},
"output": {
"Q": {
"bandmap": [
1,
2,
3,
4,
5
],
"indexmap": [
[
1,
1
],
[
2,
2
],
[
3,
3
],
[
4,
4
],
[
5,
5
]
],
"K000001": {
"eks": [
-13.230981270352434,
-8.226552436785038,
-8.22537802992713,
-8.225073262212328,
-0.4646702125746965
],
"z": [
0.6612260601413072,
0.868604768882124,
0.8686339641874008,
0.868616662913504,
0.9803946557176213
],
"eqpLin": [
-16.26615213933269,
-12.145989533047741,
-12.141823660074568,
-12.141850045890173,
0.664173458065683
],
"eqpSec": [
-16.126442757174104,
-12.041375786118818,
-12.037079367870797,
-12.037101945463592,
0.663067500001202
],
"sigma_diff": [
-0.022981358233819825,
-0.002807756578937437,
-0.0027997422774165845,
-0.0028000399733918565,
0.0
],
"occupation": [
2.0,
2.0,
2.0,
2.0,
0.0
],
"sigmax": [
-17.605224493338945,
-15.765070382904447,
-15.764310528252162,
-15.76413620844094,
-0.5854584095249137
],
"vxcl": [
-11.249025118576247,
-11.242266869012985,
-11.241788303510699,
-11.241686110388825,
-2.0850079337127707
],
"vxcnl": [
0.0,
0.0,
0.0,
0.0,
0.0
],
"hf": [
-6.3561993747626975,
-4.5228035138914615,
-4.522522224741464,
-4.522450098052116,
1.4995495241878571
],
"sigmac_eks": {
"re": [
1.7659827273250182,
0.010466906176382348,
0.013781153568594457,
0.013237978473109088,
-0.3481319149046877
],
"im": [
0.014583989091543714,
0.003272239246325086,
0.003272405433181099,
0.0032729463236466855,
-1.058286063240667e-05
]
},
"sigmac_eqpLin": {
"re": [
3.581007935231201,
0.7305227660066568,
0.7333805411639928,
0.7329844903483206,
-0.3718350341861773
],
"im": [
0.07514251063272293,
0.00892881625011567,
0.008931013424447288,
0.008932371065598663,
-0.00011661816924339793
]
},
"sigmac_eqpSec": {
"re": [
3.4409541840918707,
0.7073751375953274,
0.710217882884513,
0.7098182772432421,
-0.37181078103599796
],
"im": [
0.11069082478558719,
0.008794218426373911,
0.008796293848555962,
0.00879763962964248,
-0.00011661816924339793
]
}
}
}
},
"timing": {
"WFREQ": {
"cpu:sec": 8.173569,
"cpu:hum": "08.1s",
"wall:sec": 8.818508863449097,
"wall:hum": "08.8s",
"nocalls": 2
},
"wfreq_readin": {
"cpu:sec": 0.4345310000000001,
"cpu:hum": "00.4s",
"wall:sec": 0.5477490425109863,
"wall:hum": "00.5s",
"nocalls": 1
},
"fetch_input": {
"cpu:sec": 0.0393199999999998,
"cpu:hum": "< 00.1s",
"wall:sec": 0.05291390419006348,
"wall:hum": "< 00.1s",
"nocalls": 2
},
"fft": {
"cpu:sec": 0.036698999999999815,
"cpu:hum": "< 00.1s",
"wall:sec": 0.08657431602478027,
"wall:hum": "< 00.1s",
"nocalls": 14
},
"fft_scatter": {
"cpu:sec": 2.313853999999959,
"cpu:hum": "02.3s",
"wall:sec": 2.429741621017456,
"wall:hum": "02.4s",
"nocalls": 14473
},
"init_vloc": {
"cpu:sec": 0.054491999999999985,
"cpu:hum": "< 00.1s",
"wall:sec": 0.05463600158691406,
"wall:hum": "< 00.1s",
"nocalls": 2
},
"init_us_1": {
"cpu:sec": 0.00997300000000001,
"cpu:hum": "< 00.1s",
"wall:sec": 0.011365890502929688,
"wall:hum": "< 00.1s",
"nocalls": 2
},
"v_of_rho": {
"cpu:sec": 0.022467999999999932,
"cpu:hum": "< 00.1s",
"wall:sec": 0.02729201316833496,
"wall:hum": "< 00.1s",
"nocalls": 1
},
"v_xc": {
"cpu:sec": 0.0381419999999999,
"cpu:hum": "< 00.1s",
"wall:sec": 0.04444599151611328,
"wall:hum": "< 00.1s",
"nocalls": 2
},
"v_h": {
"cpu:sec": 0.0015659999999999563,
"cpu:hum": "< 00.1s",
"wall:sec": 0.0016570091247558594,
"wall:hum": "< 00.1s",
"nocalls": 1
},
"do_setup": {
"cpu:sec": 0.048203999999999914,
"cpu:hum": "< 00.1s",
"wall:sec": 0.0633389949798584,
"wall:hum": "< 00.1s",
"nocalls": 1
},
"init_pw_ar": {
"cpu:sec": 0.040395999999999876,
"cpu:hum": "< 00.1s",
"wall:sec": 0.0461888313293457,
"wall:hum": "< 00.1s",
"nocalls": 1
},
"hinit0": {
"cpu:sec": 0.033833,
"cpu:hum": "< 00.1s",
"wall:sec": 0.03384113311767578,
"wall:hum": "< 00.1s",
"nocalls": 1
},
"solve_hf": {
"cpu:sec": 0.06213999999999986,
"cpu:hum": "< 00.1s",
"wall:sec": 0.07319092750549316,
"wall:hum": "< 00.1s",
"nocalls": 1
},
"sigmavxc": {
"cpu:sec": 0.0049640000000001905,
"cpu:hum": "< 00.1s",
"wall:sec": 0.006000041961669922,
"wall:hum": "< 00.1s",
"nocalls": 1
},
"init_us_2": {
"cpu:sec": 0.0016629999999988598,
"cpu:hum": "< 00.1s",
"wall:sec": 0.0018558502197265625,
"wall:hum": "< 00.1s",
"nocalls": 3
},
"init_us_2:cp": {
"cpu:sec": 0.001620000000000843,
"cpu:hum": "< 00.1s",
"wall:sec": 0.0018200874328613281,
"wall:hum": "< 00.1s",
"nocalls": 3
},
"fftw": {
"cpu:sec": 4.50848900000002,
"cpu:hum": "04.5s",
"wall:sec": 4.702910661697388,
"wall:hum": "04.7s",
"nocalls": 14439
},
"sigmax": {
"cpu:sec": 0.03873700000000002,
"cpu:hum": "< 00.1s",
"wall:sec": 0.041181087493896484,
"wall:hum": "< 00.1s",
"nocalls": 1
},
"sqvc_init": {
"cpu:sec": 0.03898700000000055,
"cpu:hum": "< 00.1s",
"wall:sec": 0.039562225341796875,
"wall:hum": "< 00.1s",
"nocalls": 7
},
"ffts": {
"cpu:sec": 0.014574000000000087,
"cpu:hum": "< 00.1s",
"wall:sec": 0.01490020751953125,
"wall:hum": "< 00.1s",
"nocalls": 20
},
"write_hf": {
"cpu:sec": 0.0,
"cpu:hum": "< 00.1s",
"wall:sec": 0.0015702247619628906,
"wall:hum": "< 00.1s",
"nocalls": 1
},
"wlanczos": {
"cpu:sec": 3.395701,
"cpu:hum": "03.3s",
"wall:sec": 3.589892864227295,
"wall:hum": "03.5s",
"nocalls": 1
},
"pdep_read": {
"cpu:sec": 0.06580900000000334,
"cpu:hum": "< 00.1s",
"wall:sec": 0.08992171287536621,
"wall:hum": "< 00.1s",
"nocalls": 100
},
"alphapc": {
"cpu:sec": 0.19857400000001313,
"cpu:hum": "00.1s",
"wall:sec": 0.20917320251464844,
"wall:hum": "00.2s",
"nocalls": 283
},
"lan_H": {
"cpu:sec": 6.412891999999999,
"cpu:hum": "06.4s",
"wall:sec": 6.658470630645752,
"wall:hum": "06.6s",
"nocalls": 9
},
"h_psi": {
"cpu:sec": 5.180463000000021,
"cpu:hum": "05.1s",
"wall:sec": 5.406877517700195,
"wall:hum": "05.4s",
"nocalls": 270
},
"h_psi:pot": {
"cpu:sec": 5.140214000000003,
"cpu:hum": "05.1s",
"wall:sec": 5.367469787597656,
"wall:hum": "05.3s",
"nocalls": 270
},
"vloc_psi": {
"cpu:sec": 4.892422000000003,
"cpu:hum": "04.8s",
"wall:sec": 5.10514497756958,
"wall:hum": "05.1s",
"nocalls": 270
},
"h_psi:calbec": {
"cpu:sec": 0.1674900000000079,
"cpu:hum": "00.1s",
"wall:sec": 0.1779468059539795,
"wall:hum": "00.1s",
"nocalls": 270
},
"calbec": {
"cpu:sec": 0.16089799999999777,
"cpu:hum": "00.1s",
"wall:sec": 0.17319059371948242,
"wall:hum": "00.1s",
"nocalls": 270
},
"add_vuspsi": {
"cpu:sec": 0.0712160000000015,
"cpu:hum": "< 00.1s",
"wall:sec": 0.07684707641601562,
"wall:hum": "< 00.1s",
"nocalls": 270
},
"brak": {
"cpu:sec": 0.09926699999999933,
"cpu:hum": "< 00.1s",
"wall:sec": 0.10172724723815918,
"wall:hum": "00.1s",
"nocalls": 9
},
"sw_restart": {
"cpu:sec": 0.1092929999999992,
"cpu:hum": "00.1s",
"wall:sec": 0.17614316940307617,
"wall:hum": "00.1s",
"nocalls": 4
},
"chi_invert": {
"cpu:sec": 0.01360599999999934,
"cpu:hum": "< 00.1s",
"wall:sec": 0.017594099044799805,
"wall:hum": "< 00.1s",
"nocalls": 1
},
"glanczos": {
"cpu:sec": 4.018922,
"cpu:hum": "04.0s",
"wall:sec": 4.241312026977539,
"wall:hum": "04.2s",
"nocalls": 1
},
"write_over": {
"cpu:sec": 0.001440000000000552,
"cpu:hum": "< 00.1s",
"wall:sec": 0.011359930038452148,
"wall:hum": "< 00.1s",
"nocalls": 5
},
"write_gfreq": {
"cpu:sec": 0.003713999999998663,
"cpu:hum": "< 00.1s",
"wall:sec": 0.029255390167236328,
"wall:hum": "< 00.1s",
"nocalls": 5
},
"sg_restart": {
"cpu:sec": 0.003068000000000737,
"cpu:hum": "< 00.1s",
"wall:sec": 0.01123499870300293,
"wall:hum": "< 00.1s",
"nocalls": 5
},
"solve_qp": {
"cpu:sec": 0.16948000000000008,
"cpu:hum": "00.1s",
"wall:sec": 0.24360990524291992,
"wall:hum": "00.2s",
"nocalls": 1
},
"coll_gw": {
"cpu:sec": 0.020823000000000036,
"cpu:hum": "< 00.1s",
"wall:sec": 0.03419899940490723,
"wall:hum": "< 00.1s",
"nocalls": 1
},
"read_over": {
"cpu:sec": 0.0018779999999996022,
"cpu:hum": "< 00.1s",
"wall:sec": 0.0029571056365966797,
"wall:hum": "< 00.1s",
"nocalls": 5
},
"read_gfreq": {
"cpu:sec": 0.003434999999999633,
"cpu:hum": "< 00.1s",
"wall:sec": 0.01312708854675293,
"wall:hum": "< 00.1s",
"nocalls": 5
},
"read_hf": {
"cpu:sec": 0.00034199999999984243,
"cpu:hum": "< 00.1s",
"wall:sec": 0.00039505958557128906,
"wall:hum": "< 00.1s",
"nocalls": 1
},
"sigmac_i": {
"cpu:sec": 0.020918999999997467,
"cpu:hum": "< 00.1s",
"wall:sec": 0.021049976348876953,
"wall:hum": "< 00.1s",
"nocalls": 4
},
"sigmac_r": {
"cpu:sec": 0.0012529999999983943,
"cpu:hum": "< 00.1s",
"wall:sec": 0.0013022422790527344,
"wall:hum": "< 00.1s",
"nocalls": 4
},
"wfreq_db": {
"cpu:sec": 0.0202650000000002,
"cpu:hum": "< 00.1s",
"wall:sec": 0.038127899169921875,
"wall:hum": "< 00.1s",
"nocalls": 1
}
}
}
We can print the quasiparticle energy levels using the WESTpy Python package, which offers useful pre- and post-processing tools for WEST calculations.
[13]:
from westpy import wfreq2df
df, data = wfreq2df('silane.wfreq.save/wfreq.json')
display(df)
_ _ _____ _____ _____
| | | | ___/ ___|_ _|
| | | | |__ \ `--. | |_ __ _ _
| |/\| | __| `--. \ | | '_ \| | | |
\ /\ / |___/\__/ / | | |_) | |_| |
\/ \/\____/\____/ \_/ .__/ \__, |
| | __/ |
|_| |___/
WEST version : 5.2.0
Today : 2022-10-11 09:34:04.514495
| k | s | n | eks | eqpLin | eqpSec | sigmax | sigmac_eks | sigmac_eqpLin | sigmac_eqpSec | vxcl | vxcnl | hf | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | 1 | -13.230981 | -16.266152 | -16.126443 | -17.605224 | 1.765983 | 3.581008 | 3.440954 | -11.249025 | 0.0 | -6.356199 |
| 1 | 1 | 1 | 2 | -8.226552 | -12.145990 | -12.041376 | -15.765070 | 0.010467 | 0.730523 | 0.707375 | -11.242267 | 0.0 | -4.522804 |
| 2 | 1 | 1 | 3 | -8.225378 | -12.141824 | -12.037079 | -15.764311 | 0.013781 | 0.733381 | 0.710218 | -11.241788 | 0.0 | -4.522522 |
| 3 | 1 | 1 | 4 | -8.225073 | -12.141850 | -12.037102 | -15.764136 | 0.013238 | 0.732984 | 0.709818 | -11.241686 | 0.0 | -4.522450 |
| 4 | 1 | 1 | 5 | -0.464670 | 0.664173 | 0.663068 | -0.585458 | -0.348132 | -0.371835 | -0.371811 | -2.085008 | 0.0 | 1.499550 |
All energies in the table are given in eV. The full-frequency G\(_0\)W\(_0\) energies correspond to the eqpSec column.
Full explanation of the keys:
k: k-point indexs: spin indexn: state indexeks: Kohn-Sham energy \(\varepsilon^{\text{KS}}_n\), obtained in the DFT calculationeqpLin: quasiparticle energy \(E^{\text{QP}}_n\), obtained by approximating the self-energy to first order in the frequencyeqpSec: quasiparticle energy \(E^{\text{QP}}_n\), obtained by using the secant method to solve the frequency-dependency in the quasiparticle equationsigmax: exchange self-energysigmac_eks: correlation self-energy, evaluated ateks(reandimidentify the real and imaginary parts respectively)sigmac_eqpLin: correlation self-energy, evaluated ateqpLinsigmac_eqpSec: correlation self-energy, evaluated ateqpSecvxcl: contribution to the energy given by the semilocal exchange-correlation functionalvxcnl: contribution to the energy given by the excact exchange (when using hybrid functional)hf: quasiparticle energy according to perturbative Hartree-Fock, i.e., no correlation self-energy
Let’s use WESTpy to plot the density of states (DOS), obtained as a sum of Gaussian functions
\begin{equation} \text{DOS}(E) = \sum_{i=1}^{N_{\text{states}}} \frac{1}{\sqrt{2\pi\sigma^2}} e^{-\frac{(E-E_i)^2}{2\sigma^2}} \end{equation}
We use \(\sigma\) = 0.2 eV.
[14]:
from westpy import ElectronicStructure
# extracting energy levels from the data
y = {}
y['dft'] = data['output']['Q']['K000001']['eks']
y['gw'] = data['output']['Q']['K000001']['eqpSec']
es = ElectronicStructure()
for key in ['dft','gw'] :
es.addKey(key,key)
for b, energy in enumerate(y[key]) :
es.addDataPoint([1,1,b+1],key,y[key][b])
es.plotDOS(kk=[1],ss=[1],energyKeys=['dft','gw'],energyRange=[-18.,2.,0.01],fname='dos.png')
Requested (emin,emax) : -18.0 2.0
Detected (emin,emax) : -16.126442757174104 0.663067500001202
output written in : dos.png
waiting for user to close image preview...
Now going back to the problem of computing the vertical ionization potential of SiH\(_4\), in the GW case we have:
\begin{equation} \text{VIP}^{\text{GW}} = 0 - E^{\text{QP}}_{\text{HOMO}} \end{equation}
The HOMO (n=4 in the table) predicted by GW is located at \(-12.037102\) eV. Let’s update the VIP plot to compare DFT and GW results.
[15]:
import matplotlib.pyplot as plt
# VIP
y = {}
y['exp'] = [12.3] # experimental
y['dft'] = [8.2251] # dft
y['gw'] = [12.037102] # gw
# colors
c = {}
c['exp'] = 'black'
c['dft'] = 'red'
c['gw'] = 'blue'
# plot
x = list(range(1, len(y)+1))
labels = y.keys()
fig, ax = plt.subplots(1, 1)
counter = 0
for counter, label in enumerate(labels) :
for a in y[label] :
ax.hlines(a, x[counter]-0.25, x[counter]+0.25, color=c[label])
plt.xticks(x, labels)
plt.ylabel('Energy (eV)')
plt.title('VIP')
plt.show()
If we take into account finite size effects, we have:
\begin{equation} \text{VIP}^{\text{GW}} = E_{\text{vacuum}} - E^{\text{QP}}_{\text{HOMO}} \end{equation}
We extract the vacuum energy from the output of the DFT calculation.
[16]:
!grep 'Corrected vacuum level' pw.out
Corrected vacuum level = 0.24581357 eV
The vacuum energy is \(0.24581357\) eV. Let’s update the VIP plot again.
[17]:
import matplotlib.pyplot as plt
# vacumm
vacuum = 0.24581357
# VIP
y = {}
y['exp'] = [12.3] # experimental
y['dft'] = [8.2251] # dft
y['dft_vac'] = [8.2251+vacuum] # dft
y['gw'] = [12.037102] # gw
y['gw_vac'] = [12.037102+vacuum] # gw
# colors
c = {}
c['exp'] = 'black'
c['dft'] = 'red'
c['dft_vac'] = 'red'
c['gw'] = 'blue'
c['gw_vac'] = 'blue'
# plot
x = list(range(1, len(y)+1))
labels = y.keys()
fig, ax = plt.subplots(1, 1)
for counter, label in enumerate(labels) :
for a in y[label] :
ax.hlines(a, x[counter]-0.25, x[counter]+0.25, color=c[label])
plt.xticks(x, labels)
plt.ylabel('Energy (eV)')
plt.title('VIP')
plt.show()
Question: Which level of theory better describes electronic excitations, DFT or GW?
Finally, we inspect the convergence of GW with respect to the size of the basis set used to represent the dielectric matrices. We have computed 50 eigenpotentials using wstat.x. We can now use fewer eigenpotentials in wfreq.x and test the convergece.
Let’s prepare the input for 5 calculations done with 10, 20, 30, 40, 50 eigenpotentials. This is achieved by setting the n_pdep_eigen_to_use keyword in the input file.
[18]:
import yaml
# create 5 files, each with a different value of the variable n_pdep_eigen_to_use
for npdep in [10, 20, 30, 40, 50] :
# read data wfreq.in
with open('wfreq.in', 'r') as file :
input_data = yaml.load(file, Loader=yaml.FullLoader)
# modify the input_data
input_data['wfreq_control']['n_pdep_eigen_to_use'] = npdep
# write data to YAML file
with open(f'wfreq_{npdep}.in', 'w') as file :
yaml.dump(input_data, file, sort_keys=False)
[19]:
# give a quick look
!cat wfreq_10.in
!cat wfreq_20.in
!cat wfreq_30.in
!cat wfreq_40.in
!cat wfreq_50.in
input_west:
qe_prefix: silane
west_prefix: silane
outdir: ./
wstat_control:
wstat_calculation: S
n_pdep_eigen: 50
wfreq_control:
wfreq_calculation: XWGQ
n_pdep_eigen_to_use: 10
qp_bandrange:
- 1
- 5
n_refreq: 300
ecut_refreq: 2.0
input_west:
qe_prefix: silane
west_prefix: silane
outdir: ./
wstat_control:
wstat_calculation: S
n_pdep_eigen: 50
wfreq_control:
wfreq_calculation: XWGQ
n_pdep_eigen_to_use: 20
qp_bandrange:
- 1
- 5
n_refreq: 300
ecut_refreq: 2.0
input_west:
qe_prefix: silane
west_prefix: silane
outdir: ./
wstat_control:
wstat_calculation: S
n_pdep_eigen: 50
wfreq_control:
wfreq_calculation: XWGQ
n_pdep_eigen_to_use: 30
qp_bandrange:
- 1
- 5
n_refreq: 300
ecut_refreq: 2.0
input_west:
qe_prefix: silane
west_prefix: silane
outdir: ./
wstat_control:
wstat_calculation: S
n_pdep_eigen: 50
wfreq_control:
wfreq_calculation: XWGQ
n_pdep_eigen_to_use: 40
qp_bandrange:
- 1
- 5
n_refreq: 300
ecut_refreq: 2.0
input_west:
qe_prefix: silane
west_prefix: silane
outdir: ./
wstat_control:
wstat_calculation: S
n_pdep_eigen: 50
wfreq_control:
wfreq_calculation: XWGQ
n_pdep_eigen_to_use: 50
qp_bandrange:
- 1
- 5
n_refreq: 300
ecut_refreq: 2.0
We run them all and stash the output files. These calculations should finish in 3 minutes.
[20]:
!rm -f silane.wfreq.save/*.json
!mpirun -n 16 wfreq.x -i wfreq_10.in > wfreq_10.out
!mv silane.wfreq.save/wfreq.json wfreq_10.json
!mpirun -n 16 wfreq.x -i wfreq_20.in > wfreq_20.out
!mv silane.wfreq.save/wfreq.json wfreq_20.json
!mpirun -n 16 wfreq.x -i wfreq_30.in > wfreq_30.out
!mv silane.wfreq.save/wfreq.json wfreq_30.json
!mpirun -n 16 wfreq.x -i wfreq_40.in > wfreq_40.out
!mv silane.wfreq.save/wfreq.json wfreq_40.json
!mpirun -n 16 wfreq.x -i wfreq_50.in > wfreq_50.out
!mv silane.wfreq.save/wfreq.json wfreq_50.json
We plot the quasiparticle energy levels as a function of n_pdep_eigen_to_use.
[21]:
import json
import matplotlib.pyplot as plt
# load data from wfreq_XX.json
data = {}
for npdep in [10, 20, 30, 40, 50] :
with open(f'wfreq_{npdep}.json') as file :
data[npdep] = json.load(file)
# energy levels
y = {}
for npdep in [10, 20, 30, 40, 50] :
y[npdep] = data[npdep]['output']['Q']['K000001']['eqpSec']
# colors
c = {}
c[10] = 'red'
c[20] = 'blue'
c[30] = 'green'
c[40] = 'purple'
c[50] = 'black'
# plot
x = list(range(1, len(y)+1))
labels = y.keys()
fig, ax = plt.subplots(1, 1)
for counter, label in enumerate(labels) :
for a in y[label] :
ax.hlines(a, x[counter]-0.25, x[counter]+0.25, color=c[label])
plt.xticks(x, labels)
plt.ylabel('Energy (eV)')
plt.title('Energy structure')
plt.show()
Question: How many PDEPs are needed to converge the quasiparticle energies?