This tutorial can be downloaded link.

Intro Tutorial 4: Plotting Density of States and Local Density of States¶

We show how to use the WESTpy Python package to plot the density of states and local density of states from GW quasiparticle energy levels calculated using WEST.

Running DFT and GW calculations¶

We will use the negatively charged nitrogen-vacancy center in diamond as an example. Download the input and pseudopotential files to your current directory:

[1]:

%%bash
wget -N -q https://west-code.org/doc/training/nv_diamond/pw.in
wget -N -q https://west-code.org/doc/training/nv_diamond/wstat.in
wget -N -q https://west-code.org/doc/training/nv_diamond/wfreq.in
wget -N -q http://www.quantum-simulation.org/potentials/sg15_oncv/upf/C_ONCV_PBE-1.2.upf
wget -N -q http://www.quantum-simulation.org/potentials/sg15_oncv/upf/N_ONCV_PBE-1.2.upf

Let’s inspect the pw.in file.

[2]:

%%bash
cat pw.in

&CONTROL
  calculation = 'scf'
  wf_collect = .true.
  prefix = 'nv_diamond'
  outdir = './'
  pseudo_dir = './'
/
&SYSTEM
  ibrav = 0
  ecutwfc = 60
  nspin = 2
  tot_charge = -1
  tot_magnetization = 2.0
  nbnd = 532
  nat = 215
  ntyp = 2
/
&ELECTRONS
  diago_full_acc = .true.
/
K_POINTS gamma
CELL_PARAMETERS angstrom
10.70400  0.000000  0.000000
0.000000  10.70400  0.000000
0.000000  0.000000  10.70400
ATOMIC_SPECIES
C  12.01099968  C_ONCV_PBE-1.2.upf
N  14.00699997  N_ONCV_PBE-1.2.upf
ATOMIC_POSITIONS crystal
C        0.000019794   0.000070305   0.000070305
C        0.083283093   0.083408574   0.083408574
C       -0.000180249   0.166803354   0.166803354
C        0.082881581   0.250195244   0.250195244
C        0.166624862  -0.000041790   0.166735686
C        0.249868283   0.083201633   0.250010134
C        0.166624862   0.166735686  -0.000041790
C        0.249868283   0.250010134   0.083201633
C       -0.000091814   0.000127381   0.333479766
C        0.083066258   0.083584948   0.416725766
C       -0.000487477   0.167194483   0.500175895
C        0.083392434   0.250002399   0.583473949
C        0.166924803   0.000258154   0.500150151
C        0.250013786   0.083347137   0.583306800
C        0.166311809   0.166704028   0.333362376
C        0.249570188   0.249881859   0.416552187
C        0.000115132   0.000238438   0.666764846
C        0.083394701   0.083489806   0.750066588
C        0.000026260   0.166833904   0.833359608
C        0.083323351   0.250184829   0.916656698
C        0.166736956   0.000070304   0.833353141
C        0.250075226   0.083408574   0.916616438
C        0.166806534   0.166683497   0.666597859
C        0.250134179   0.250159636   0.749908003
C       -0.000091814   0.333479766   0.000127381
C        0.083066258   0.416725766   0.083584948
C       -0.000487477   0.500175895   0.167194483
C        0.083392434   0.583473949   0.250002399
C        0.166311809   0.333362376   0.166704028
C        0.249570188   0.416552187   0.249881859
C        0.166924803   0.500150151   0.000258154
C        0.250013786   0.583306800   0.083347137
C       -0.000658632   0.333776112   0.333776112
C        0.082664024   0.417142882   0.417142882
C       -0.000820418   0.500581677   0.500581677
C        0.082964776   0.583567558   0.583567558
C        0.166378380   0.333444811   0.500019126
C        0.249908129   0.416758715   0.583268171
C        0.166378380   0.500019126   0.333444811
C        0.249908129   0.583268171   0.416758715
C        0.000384359   0.333113519   0.666472256
C        0.083422600   0.416888861   0.750099921
C        0.001008133   0.500569049   0.834341482
C        0.083928842   0.583364715   0.917262190
C        0.166794030   0.333479766   0.833241534
C        0.250251595   0.416725765   0.916399608
C        0.166656573   0.500161197   0.666656696
C        0.250240644   0.583409895   0.749872737
C        0.000115132   0.666764846   0.000238438
C        0.083394701   0.750066588   0.083489806
C        0.000026260   0.833359608   0.166833904
C        0.083323351   0.916656698   0.250184829
C        0.166806534   0.666597859   0.166683497
C        0.250134179   0.749908003   0.250159636
C        0.166736956   0.833353141   0.000070304
C        0.250075226   0.916616438   0.083408574
C        0.000384359   0.666472256   0.333113519
C        0.083422600   0.750099921   0.416888861
C        0.001008133   0.834341482   0.500569049
C        0.083928842   0.917262190   0.583364715
C        0.166656573   0.666656696   0.500161197
C        0.250240644   0.749872737   0.583409895
C        0.166794030   0.833241534   0.333479766
C        0.250251595   0.916399608   0.416725765
C       -0.000213111   0.666782588   0.666782588
C        0.083311937   0.749952076   0.749952076
C        0.000052516   0.833385863   0.833385863
C        0.083382825   0.916716171   0.916716171
C        0.166905085   0.666764845   0.833448478
C        0.250156456   0.750066587   0.916728048
C        0.166905085   0.833448478   0.666764845
C        0.250156456   0.916728048   0.750066587
C        0.333402340  -0.000041790  -0.000041790
C        0.416676784   0.083201633   0.083201633
C        0.333159833   0.166493178   0.166493178
C        0.416349440   0.249682793   0.249682793
C        0.500029025  -0.000354839   0.166704028
C        0.583218841   0.082903541   0.249881860
C        0.500029025   0.166704028  -0.000354839
C        0.583218841   0.249881860   0.082903541
C        0.333370681  -0.000354839   0.333362375
C        0.416548508   0.082903541   0.416552188
C        0.332959588   0.166292937   0.499855389
C        0.416063422   0.249396771   0.582921731
C        0.500111458  -0.000288269   0.500019126
C        0.583425371   0.083241482   0.583268171
C        0.499651689   0.166054115   0.332985041
C        0.582748739   0.249135425   0.416082085
C        0.333350150   0.000139883   0.666597859
C        0.416826287   0.083467528   0.749908001
C        0.333470008   0.166803351   0.833153102
C        0.416861892   0.250195244   0.916214932
C        0.500146416   0.000127380   0.833241534
C        0.583392421   0.083584946   0.916399607
C        0.500065835   0.166814168   0.666338186
C        0.583619715   0.251058000   0.748936005
C        0.333370681   0.333362375  -0.000354839
C        0.416548508   0.416552188   0.082903541
C        0.332959588   0.499855389   0.166292937
C        0.416063422   0.582921731   0.249396771
C        0.499651689   0.332985041   0.166054115
C        0.582748739   0.416082085   0.249135425
C        0.500111458   0.500019126  -0.000288269
C        0.583425371   0.583268171   0.083241482
C        0.332720767   0.332985042   0.332985042
C        0.415802073   0.416082083   0.416082083
C        0.333250645   0.499852804   0.499852804
C        0.417222928   0.583187872   0.583187872
C        0.498591763   0.331925116   0.499441582
C        0.583396493   0.416729832   0.579206565
C        0.498591763   0.499441582   0.331925116
C        0.583396493   0.579206565   0.416729832
C        0.333480822   0.333399188   0.666338185
C        0.417724652   0.416953059   0.748936006
C        0.333861137   0.500175896   0.832845871
C        0.416669048   0.583473950   0.916725784
C        0.500442762   0.333776113   0.832674716
C        0.583809537   0.417142883   0.915997374
C        0.501709296   0.500779318   0.665000232
C        0.585243702   0.584379787   0.748657890
C        0.333350150   0.666597859   0.000139883
C        0.416826287   0.749908001   0.083467528
C        0.333470008   0.833153102   0.166803351
C        0.416861892   0.916214932   0.250195244
C        0.500065835   0.666338186   0.166814168
C        0.583619715   0.748936005   0.251058000
C        0.500146416   0.833241534   0.000127380
C        0.583392421   0.916399607   0.083584946
C        0.333480822   0.666338185   0.333399188
C        0.417724652   0.748936006   0.416953059
C        0.333861137   0.832845871   0.500175896
C        0.416669048   0.916725784   0.583473950
C        0.501709296   0.665000232   0.500779318
C        0.585243702   0.748657890   0.584379787
C        0.500442762   0.832674716   0.333776113
C        0.583809537   0.915997374   0.417142883
C        0.333666749   0.666621439   0.666621439
C        0.416918318   0.750005875   0.750005875
C        0.333500557   0.833359606   0.833359606
C        0.416851480   0.916656697   0.916656697
C        0.499780168   0.666472252   0.833717707
C        0.583555519   0.750099920   0.916755948
C        0.499780168   0.833717707   0.666472252
C        0.583555519   0.916755948   0.750099920
C        0.666816800   0.000258154   0.000258154
C        0.749973452   0.083347138   0.083347138
C        0.666522036   0.166292937   0.166292937
C        0.749588381   0.249396774   0.249396774
C        0.833264514   0.000139885   0.166683496
C        0.916574652   0.083467530   0.250159637
C        0.833264514   0.166683496   0.000139885
C        0.916574652   0.250159637   0.083467530
C        0.666685774  -0.000288268   0.333444809
C        0.749934823   0.083241482   0.416758717
C        0.666519454   0.166583992   0.499852804
C        0.749854524   0.250556280   0.583187872
C        0.833323351  -0.000010077   0.500161196
C        0.916539384   0.083573998   0.583409896
C        0.833004838   0.166814169   0.333399186
C        0.915602654   0.251058001   0.416953059
C        0.666827846  -0.000010077   0.666656695
C        0.750076548   0.083573997   0.749872735
C        0.666842545   0.167194483   0.832845873
C        0.750140602   0.250002397   0.916725783
C        0.833431502   0.000238437   0.833448480
C        0.916733238   0.083489806   0.916728049
C        0.833288094   0.167000099   0.666621440
C        0.916672524   0.250251670   0.750005875
C        0.666685774   0.333444809  -0.000288268
C        0.749934823   0.416758717   0.083241482
C        0.666519454   0.499852804   0.166583992
C        0.749854524   0.583187872   0.250556280
C        0.833004838   0.333399186   0.166814169
C        0.915602654   0.416953059   0.251058001
C        0.833323351   0.500161196  -0.000010077
C        0.916539384   0.583409896   0.083573998
C        0.666108231   0.331925113   0.331925113
C        0.745873217   0.416729834   0.416729834
C        0.831666888   0.335042646   0.500779317
C        0.915324539   0.418577041   0.584379782
C        0.831666888   0.500779317   0.335042646
C        0.915324539   0.584379782   0.418577041
C        0.667445965   0.335042648   0.665000235
C        0.751046436   0.418577041   0.748657890
C        0.667248326   0.500581676   0.832512930
C        0.750234210   0.583567558   0.916298128
C        0.833138909   0.333113518   0.833717710
C        0.916766570   0.416888861   0.916755950
C        0.839493491   0.496368856   0.672826841
C        0.916577565   0.582927166   0.749910917
C        0.666827846   0.666656695  -0.000010077
C        0.750076548   0.749872735   0.083573997
C        0.666842545   0.832845873   0.167194483
C        0.750140602   0.916725783   0.250002397
C        0.833288094   0.666621440   0.167000099
C        0.916672524   0.750005875   0.250251670
C        0.833431502   0.833448480   0.000238437
C        0.916733238   0.916728049   0.083489806
C        0.667445965   0.665000235   0.335042648
C        0.751046436   0.748657890   0.418577041
C        0.667248326   0.832512930   0.500581676
C        0.750234210   0.916298128   0.583567558
C        0.839493491   0.672826841   0.496368856
C        0.916577565   0.749910917   0.582927166
C        0.833138909   0.833717710   0.333113518
C        0.916766570   0.916755950   0.416888861
C        0.663035502   0.672826834   0.672826834
C        0.749593820   0.749910920   0.749910920
C        0.667235700   0.834341481   0.834341481
C        0.750031368   0.917262190   0.917262190
C        0.833449242   0.666782586   0.833120239
C        0.916618725   0.749952077   0.916645283
C        0.833449242   0.833120239   0.666782586
C        0.916618725   0.916645283   0.749952077
N        0.658425487   0.491758845   0.491758845

This is a relatively big system with 215 atoms and 862 valence electrons, and requires an explicit treatment of spin polarization. As detailed in Yu et al., J. Chem. Theory Comput. 18, 4690-4707 (2022), WEST can use GPU accelerators to carry out large-scale GW calculations. For instance, we can use four GPUs for the ground state DFT calculation, and 64 GPUs for the GW calculation.

[ ]:

%%bash
mpirun -n 4 pw.x -i pw.in > pw.out
mpirun -n 64 wstat.x -nimage 16 -npool 2 -nbgrp 2 -i wstat.in > wstat.out
mpirun -n 64 wfreq.x -nimage 8 -npool 2 -nbgrp 4 -i wfreq.in > wfreq.out

Note that we have used the image, pool, and band group parallelization levels, as discussed in Tutorial 3. By using all the parallelization levels, we avoid carrying out fast Fourier transforms (FFTs) on multiple GPUs, which is inefficient because of the all-to-all communication involved in parallel FFTs.

If the reader does NOT have the computational resources to run the calculation, the WEST output file needed for the next step can be directly downloaded as:

[3]:

%%bash
mkdir -p nv_diamond.wfreq.save
wget -N -q https://west-code.org/doc/training/nv_diamond/wfreq.json -O nv_diamond.wfreq.save/wfreq.json

Plotting density of states¶

The density of states (DOS) is defined as follows: \begin{equation} \text{DOS}(E) = \sum_{\sigma=1}^{N_{\text{spin}}} \sum_{i=1}^{N_{\text{band}}} \delta (E-E_i^{\sigma}) \end{equation} where \(E_i^{\sigma}\) is the energy of the \(i\)th wavefunction in the \(\sigma\) spin channel obtained at the PBE or GW level of theory, and \(\delta\) is the Dirac delta function. In practice, we model the delta function by a Gaussian function with a chosen width. Here we are going to use 0.1 eV.

To plot the DOS, we first collect the output of WEST into a pandas dataframe (see also Tutorial 1).

[4]:

from westpy import *

df, data = wfreq2df('nv_diamond.wfreq.save/wfreq.json')
display(df)


 _    _ _____ _____ _____
| |  | |  ___/  ___|_   _|
| |  | | |__ \ `--.  | |_ __  _   _
| |/\| |  __| `--. \ | | '_ \| | | |
\  /\  / |___/\__/ / | | |_) | |_| |
 \/  \/\____/\____/  \_/ .__/ \__, |
                       | |     __/ |
                       |_|    |___/

WEST version     :  5.1.0
Today            :  2022-08-22 16:47:58.485276

	k	s	n	eks	eqpLin	eqpSec	sigmax	sigmac_eks	sigmac_eqpLin	sigmac_eqpSec	vxcl	vxcnl	hf
0	1	1	321	8.367174	7.400018	7.433568	-20.102276	3.000145	3.222744	3.214358	-15.953647	0.0	-4.148628
1	1	1	322	8.367273	7.400794	7.433616	-20.101597	2.999942	3.222335	3.214098	-15.953158	0.0	-4.148439
2	1	1	323	8.376845	7.430813	7.438766	-20.127620	3.018200	3.234410	3.232714	-15.956949	0.0	-4.170671
3	1	1	324	8.377936	7.432131	7.442610	-20.138275	3.013117	3.232256	3.229776	-15.973121	0.0	-4.165155
4	1	1	325	8.381440	7.390587	7.456276	-20.008882	2.970645	3.196673	3.181567	-15.902148	0.0	-4.106734
...	...	...	...	...	...	...	...	...	...	...	...	...	...
395	1	2	516	20.940372	22.386984	23.649553	-8.283733	-4.933001	-3.719844	-3.964094	-14.956837	0.0	6.673104
396	1	2	517	20.950590	22.336122	23.597504	-8.047657	-4.891736	-3.703635	-3.921624	-14.616189	0.0	6.568532
397	1	2	518	20.950620	22.335941	23.597143	-8.047741	-4.891770	-3.703932	-3.922213	-14.616398	0.0	6.568656
398	1	2	519	20.964816	22.421303	23.723918	-8.304636	-4.927318	-3.661473	-3.916795	-14.987825	0.0	6.683188
399	1	2	520	20.964843	22.420867	23.729650	-8.303952	-4.927378	-3.661139	-3.925628	-14.986806	0.0	6.682855

400 rows × 13 columns

We create an ElectronicStructure object and attach to it the DFT and GW energy levels.

[5]:

es = ElectronicStructure()

dos_keys = ['eks','eqpSec']

for key in dos_keys :
    es.addKey(key,key)
    for i, row in df.iterrows() :
        k, s, n = int(row['k']), int(row['s']), int(row['n'])
        es.addDataPoint([k,s,n],key,row[key])

We call the plotDOS method to plot the DOS.

[6]:

es.plotDOS(kk=[1],ss=[1,2],energyKeys=['eks'],energyRange=[8,20,0.01],fname='dos_dft.jpg')
es.plotDOS(kk=[1],ss=[1,2],energyKeys=['eqpSec'],energyRange=[8,20,0.01],fname='dos_gw.jpg')
es.plotDOS(kk=[1],ss=[1],energyKeys=['eks','eqpSec'],energyRange=[11,20,0.01],fname='dos_gw2.jpg')
es.plotDOS(kk=[1],ss=[2],energyKeys=['eks','eqpSec'],energyRange=[11,20,0.01],fname='dos_gw3.jpg')

Requested (emin,emax) :  8 20
Detected  (emin,emax) :  8.36717413576688 20.964843493350457
output written in :  dos_dft.jpg
waiting for user to close image preview...

../../_images/tutorials_basic_basic_004_20_1.png

Requested (emin,emax) :  8 20
Detected  (emin,emax) :  7.433567942407069 23.951547524074822
output written in :  dos_gw.jpg
waiting for user to close image preview...

../../_images/tutorials_basic_basic_004_20_3.png

Requested (emin,emax) :  11 20
Detected  (emin,emax) :  7.433567942407069 23.951547524074822
output written in :  dos_gw2.jpg
waiting for user to close image preview...

../../_images/tutorials_basic_basic_004_20_5.png

Requested (emin,emax) :  11 20
Detected  (emin,emax) :  7.437083758901655 23.93715508366057
output written in :  dos_gw3.jpg
waiting for user to close image preview...

../../_images/tutorials_basic_basic_004_20_7.png

Plotting local density of states¶

The local density of states (LDOS) is defined as follows: \begin{equation} \text{LDOS}(z,E) = \sum_{\sigma=1}^{N_{\text{spin}}} \sum_{i=1}^{N_{\text{band}}} \int \frac{\text{d}x}{L_x} \int \frac{\text{d}y}{L_y} \vert \psi_i^{\sigma} (x,y,z)\vert ^2 \delta(E-E_i^{\sigma}) \end{equation} where \(\psi_i^{\sigma}\) is the \(i\)-th wavefunction in the \(\sigma\) spin channel obtained at the PBE level of theory, with energy \(E_i^{\sigma}\) obtained at the PBE or GW level of theory. \(L_x\) and \(L_y\) are the lengths of the x and y axes of the simulation box, respectively, whereas \(z\) corresponds to the z axis of the simulation box. \(\delta\) is the Dirac delta function, modeled by a Gaussian function with a width of 0.05 eV in our case.

To plot the LDOS along the z axis of the simulation box, in addition to the Kohn-Sham and quasiparticle energy levels, we need the wavefunctions averaged over z. Check out the westpp.x input description and generate an input file for WEST called westpp.in.

Download this file to your current working directory:

[7]:

%%bash
wget -N -q https://west-code.org/doc/training/nv_diamond/westpp.in

Let’s inspect the westpp.in file, input for westpp.x.

[8]:

%%bash
cat westpp.in

input_west:
    qe_prefix: nv_diamond
    west_prefix: nv_diamond
    outdir: ./

wstat_control:
    wstat_calculation: S
    n_pdep_eigen: 862
    tr2_dfpt: 0.00000001
    n_steps_write_restart: 0

wfreq_control:
    wfreq_calculation: XWGQ
    n_pdep_eigen_to_use: 862
    macropol_calculation: C
    o_restart_time: -1.0
    qp_bandrange: [321,480]

westpp_control:
    westpp_calculation: W
    westpp_range: [321,480]
    westpp_format: z
    westpp_sign: False

Run westpp.x.

[ ]:

%%bash
mpirun -n 4 westpp.x -i westpp.in > westpp.out

This computes the wavefunctions averaged over the z axis, and stores the results in several wfcKxxxxxxBxxxxxx.plavz files in nv_diamond.westpp.save. Again, if the reader does NOT have the computational resources to run the calculation, the output files needed for the next step can be downloaded as:

[9]:

%%bash
wget -N -q https://west-code.org/doc/training/nv_diamond/nv_diamond.westpp.save.tar.gz
tar zxf nv_diamond.westpp.save.tar.gz

We create an ElectronicStructure object and register the energy levels and their associate wavefunction files. The energy levels are sorted into three categories, namely the bulk states, the spin up defect states, and the spin down defect states. States strongly localized around the nitrogen-vacancy center are considered to be the defect states. We call the plotLDOS method to plot the LDOS.

[10]:

es2 = ElectronicStructure()

# defect states
defect_states = [430,431,432]

dos_keys = ['bulk','up','down']
for key in dos_keys :
    es2.addKey(key,key)

for i, row in df.iterrows() :
    k, s, n = int(row['k']), int(row['s']), int(row['n'])
    kindex = f"K{k+(s-1)*len(data['system']['bzsamp']['k']):06d}"
    if n in defect_states :
        if s == 1 :
            es2.addDataPoint([k,s,n],'up',row['eqpSec'])
            es2.addDataPoint([k,s,n],'wfcFile',f'nv_diamond.westpp.save/wfcK{s:06d}B{n:06d}.plavz')
        else :
            es2.addDataPoint([k,s,n],'down',row['eqpSec'])
            es2.addDataPoint([k,s,n],'wfcFile',f'nv_diamond.westpp.save/wfcK{s:06d}B{n:06d}.plavz')
    else :
        if n in range(321,481) :
            es2.addDataPoint([k,s,n],'bulk',row['eqpSec'])
            es2.addDataPoint([k,s,n],'wfcFile',f'nv_diamond.westpp.save/wfcK{s:06d}B{n:06d}.plavz')

es2.plotLDOS(kk=[1],ss=[1,2],energyKeys=['bulk','up','down'],sigma=0.05,energyRange=[11,20,0.01],fname='ldos_gw.jpg')

Requested (emin,emax) :  11 20
Detected  (emin,emax) :  7.433567942407069 22.976032527957337
output written in :  ldos_gw.jpg
waiting for user to close image preview...

../../_images/tutorials_basic_basic_004_33_1.png

Where the bulk states, the spin up defect states, and the spin down defect states are plotted in black, red, and blue, respectively. The reader is encouraged to adapt the script to plot the LDOS at the PBE level of theory, and compare it to the GW result. For more information on how LDOS is computed see: Yu et al., J. Chem. Theory Comput. 18, 4690-4707 (2022).