Back to the main Index

Postprocessing tools for GW calculations

This notebook explains how to use AbiPy and matplotlib to visualize the results produced by the GW code. The self-energy code (optdriver 4) saves the final results in the SIGRES.nc file while the screening code (optdriver 3) stores the inverse dielectric matrix in the SCR.nc file.

Table of Contents

[back to top]

Let's start by importing the basic modules we will need for this tutorial.

In [1]:
# Use this at the beginning of your script so that your code will be compatible with python3
from __future__ import print_function, division, unicode_literals

import warnings 
warnings.filterwarnings("ignore")  # Ignore warnings

from abipy import abilab
abilab.enable_notebook() # This line tells AbiPy we are running inside a notebook
import abipy.data as abidata

# This line configures matplotlib to show figures embedded in the notebook.
# Replace `inline` with `notebook` in classic notebook
%matplotlib inline   

# Option available in jupyterlab. See https://github.com/matplotlib/jupyter-matplotlib
#%matplotlib widget 

How to visualize QP corrections

[back to top]

As usual, we start by opening the netcdf file with abiopen:

In [2]:
sigres = abilab.abiopen(abidata.ref_file("tgw1_9o_DS4_SIGRES.nc"))
print(sigres)
================================= File Info =================================
Name: tgw1_9o_DS4_SIGRES.nc
Directory: /Users/gmatteo/git_repos/abipy/abipy/data/refs
Size: 1005.34 kb
Access Time: Tue Oct 15 01:07:46 2019
Modification Time: Wed Mar 20 16:53:35 2019
Change Time: Wed Mar 20 16:53:35 2019

================================= Structure =================================
Full Formula (Si2)
Reduced Formula: Si
abc   :   3.823046   3.823046   3.823046
angles:  60.000000  60.000000  60.000000
Sites (2)
  #  SP       a     b     c
---  ----  ----  ----  ----
  0  Si    0     0     0
  1  Si    0.25  0.25  0.25

Abinit Spacegroup: spgid: 0, num_spatial_symmetries: 48, has_timerev: True, symmorphic: True

============================== Kohn-Sham bands ==============================
Number of electrons: 8.0, Fermi level: 6.236 (eV)
nsppol: 1, nkpt: 6, mband: 100, nspinor: 1, nspden: 1
smearing scheme: , tsmear_eV: 1.088, occopt: 1
Direct gap:
    Energy: 2.513 (eV)
    Initial state: spin=0, kpt=[+0.000, +0.000, +0.000], weight: 0.031, band=3, eig=5.951, occ=2.000
    Final state:   spin=0, kpt=[+0.000, +0.000, +0.000], weight: 0.031, band=4, eig=8.464, occ=0.000
Fundamental gap:
    Energy: 0.570 (eV)
    Initial state: spin=0, kpt=[+0.000, +0.000, +0.000], weight: 0.031, band=3, eig=5.951, occ=2.000
    Final state:   spin=0, kpt=[+0.500, +0.500, +0.000], weight: 0.094, band=4, eig=6.521, occ=0.000
Bandwidth: 12.096 (eV)
Valence maximum located at:
    spin=0, kpt=[+0.000, +0.000, +0.000], weight: 0.031, band=3, eig=5.951, occ=2.000
Conduction minimum located at:
    spin=0, kpt=[+0.500, +0.500, +0.000], weight: 0.094, band=4, eig=6.521, occ=0.000


=============================== QP direct gaps ===============================
QP_dirgap: 3.537 (eV) for K-point: [-0.250, -0.250, +0.000], spin: 0
QP_dirgap: 4.357 (eV) for K-point: [-0.250, +0.250, +0.000], spin: 0
QP_dirgap: 4.117 (eV) for K-point: [+0.500, +0.500, +0.000] X, spin: 0
QP_dirgap: 8.711 (eV) for K-point: [-0.250, +0.500, +0.250] W, spin: 0
QP_dirgap: 3.297 (eV) for K-point: [+0.500, +0.000, +0.000] L, spin: 0
QP_dirgap: 3.126 (eV) for K-point: [+0.000, +0.000, +0.000] $\Gamma$, spin: 0


============== QP results for each k-point and spin (all in eV) ==============
K-point: [-0.250, -0.250, +0.000], spin: 0
   band      e0     qpe  qpe_diago   vxcme  sigxme  sigcmee0        vUme  \
0     0  -5.046  -4.760     -4.561 -10.493 -16.549     6.542  1.344e-321   
1     1   2.411   2.281      2.238 -11.089 -13.671     2.410  1.344e-321   
2     2   4.008   3.867      3.825 -10.689 -12.566     1.694  1.344e-319   
3     3   4.008   3.867      3.825 -10.689 -12.566     1.695  1.334e-322   
4     4   6.963   7.405      7.523  -9.434  -5.330    -3.544  8.053e-322   
5     5   8.979   9.462      9.594  -9.212  -4.504    -4.093  1.636e-152   
6     6  11.702  12.132     12.276 -10.569  -5.276    -4.719  1.344e-321   
7     7  11.702  12.133     12.278 -10.569  -5.276    -4.717  1.344e-321   

     ze0  
0  0.589  
1  0.751  
2  0.768  
3  0.769  
4  0.789  
5  0.786  
6  0.748  
7  0.748  
 
K-point: [-0.250, +0.250, +0.000], spin: 0
   band      e0     qpe  qpe_diago   vxcme  sigxme  sigcmee0        vUme  \
0     0  -4.064  -3.960     -3.904 -10.617 -16.264     5.808  1.340e-319   
1     1   0.418   0.303      0.258 -10.811 -14.348     3.376   0.000e+00   
2     2   2.088   1.979      1.944 -10.067 -12.427     2.216   0.000e+00   
3     3   4.590   4.437      4.391 -10.912 -12.687     1.576   0.000e+00   
4     4   8.311   8.793      8.925  -9.604  -5.066    -3.924  5.774e-313   
5     5  10.711  11.188     11.335 -10.115  -4.841    -4.649  1.334e-322   
6     6  11.052  11.455     11.573 -10.233  -5.401    -4.310  2.668e-322   
7     7  11.990  12.395     12.526  -9.894  -4.357    -5.001   0.000e+00   

     ze0  
0  0.652  
1  0.723  
2  0.758  
3  0.770  
4  0.786  
5  0.764  
6  0.773  
7  0.756  
 
K-point: [+0.500, +0.500, +0.000] X, spin: 0
   band      e0     qpe  qpe_diago   vxcme  sigxme  sigcmee0        vUme  \
0     0  -1.936  -1.940     -1.947 -10.762 -15.448     4.680  1.732e-312   
1     1  -1.936  -1.944     -1.942 -10.762 -15.448     4.675  6.719e-322   
2     2   3.030   2.851      2.795 -10.568 -12.966     2.163   0.000e+00   
3     3   3.030   2.851      2.795 -10.568 -12.966     2.163  5.774e-313   
4     4   6.521   6.968      7.084  -9.078  -5.038    -3.478  5.774e-313   
5     5   6.521   6.969      7.084  -9.078  -5.038    -3.477  1.334e-322   
6     6  15.976  16.184     16.280 -10.528  -3.786    -6.438   0.000e+00   
7     7  15.976  16.185     16.282 -10.528  -3.786    -6.437  5.385e-322   

     ze0  
0  0.700  
1  0.701  
2  0.762  
3  0.761  
4  0.795  
5  0.795  
6  0.685  
7  0.684  
 
K-point: [-0.250, +0.500, +0.250] W, spin: 0
   band      e0     qpe  qpe_diago   vxcme  sigxme  sigcmee0        vUme  \
0     0  -1.760  -1.773     -1.808 -10.873 -15.503     4.611  6.722e-320   
1     1  -1.760  -1.793     -1.778 -10.873 -15.503     4.582  1.334e-322   
2     2   1.982   1.867      1.830  -9.991 -12.434     2.290   0.000e+00   
3     3   1.982   1.867      1.830  -9.991 -12.434     2.290   0.000e+00   
4     4  10.159  10.578     10.696  -9.647  -4.967    -4.142  6.722e-320   
5     5  10.159  10.577     10.697  -9.647  -4.967    -4.143   0.000e+00   
6     6  10.902  11.366     11.509 -10.031  -4.701    -4.722   0.000e+00   
7     7  10.902  11.366     11.509 -10.031  -4.701    -4.722   0.000e+00   

     ze0  
0  0.691  
1  0.691  
2  0.757  
3  0.757  
4  0.779  
5  0.779  
6  0.764  
7  0.764  
 
K-point: [+0.500, +0.000, +0.000] L, spin: 0
   band      e0     qpe  qpe_diago   vxcme  sigxme  sigcmee0        vUme  \
0     0  -3.760  -3.647     -3.578 -10.759 -16.296     5.719   0.000e+00   
1     1  -1.137  -1.182     -1.198 -10.178 -14.330     4.091  5.543e-311   
2     2   4.726   4.559      4.509 -11.006 -12.812     1.588   0.000e+00   
3     3   4.726   4.559      4.509 -11.006 -12.812     1.588   0.000e+00   
4     4   7.408   7.856      7.978 -10.054  -5.736    -3.748   0.000e+00   
5     5   9.251   9.772      9.915  -9.698  -4.959    -4.074   0.000e+00   
6     6   9.251   9.771      9.916  -9.698  -4.959    -4.075   0.000e+00   
7     7  13.458  13.754     13.847  -8.011  -2.413    -5.210   0.000e+00   

     ze0  
0  0.623  
1  0.737  
2  0.772  
3  0.772  
4  0.785  
5  0.783  
6  0.783  
7  0.760  
 
K-point: [+0.000, +0.000, +0.000] $\Gamma$, spin: 0
   band     e0    qpe  qpe_diago   vxcme  sigxme  sigcmee0  vUme    ze0
0     0 -6.145 -5.719     -5.412 -10.388 -16.893     7.239   0.0  0.581
1     1  5.951  5.806      5.764 -11.253 -12.628     1.188   0.0  0.777
2     2  5.951  5.806      5.764 -11.253 -12.628     1.188   0.0  0.777
3     3  5.951  5.806      5.764 -11.253 -12.628     1.187   0.0  0.777
4     4  8.464  8.931      9.065 -10.042  -5.563    -3.878   0.0  0.777
5     5  8.464  8.932      9.066 -10.042  -5.563    -3.877   0.0  0.777
6     6  8.464  8.933      9.067 -10.042  -5.563    -3.876   0.0  0.778
7     7  9.205  9.664      9.798 -10.748  -5.734    -4.421   0.0  0.773
 


Let's have a look at the KS energies used to compute the Green's function $G_0$, the RPA screening $W_0$ and the $G_0W_0$ self-energy:

In [3]:
sigres.ebands.plot();

The SIGRES file contains the KS as well as the QP direct gaps for all the k-points included in the calculation (kptgw). To plot the differente QP - KS, use:

In [4]:
sigres.plot_qpgaps();

For the absolute QP gaps:

In [5]:
sigres.plot_qpgaps(plot_qpmks=False);

To plot the QP results as a function of the initial KS energy:

In [6]:
sigres.plot_qps_vs_e0(); #tight_layout=True);

and use with_fields to filter the quantity of interest:

In [7]:
sigres.plot_qps_vs_e0(with_fields=["vxcme", "sigxme", "sigcmee0"], sharey=True);

To plot the QP energies on top of the KS energies used in the SIGMA run:

In [8]:
sigres.plot_qpbands_ibz();

To plot the KS band structure with markers whose size is proportional to the QP correction and whose direction gives the sign of the correction:

In [9]:
sigres.plot_ksbands_with_qpmarkers(fact=1000);

We can also plot the $<\Psi^{KS}_{mk}\,|\,\Psi^{QP}_{nk}>$ coefficients for given spin and k-point:

In [10]:
sigres.plot_eigvec_qp(spin=0, kpoint=0);

In this case, we have a diagonal matrix because the wavefunctions are not updated ($G_0W_0$). The scenario is completely different if you start to perform self-consistent calculations with update of the QP amplitudes.

Plotting the spectral function

[back to top]

This example shows how to plot the $G_0W_0$ spectral functions $A(\omega)$ at the $\Gamma$ point. See also lesson tgw2_4

In [11]:
with abilab.abiopen(abidata.ref_file("al_g0w0_sigmaw_SIGRES.nc")) as al_sigres:
    # Plot A(w) for the first spin, the gamma point, and bands in [0,1,2,3]
    al_sigres.plot_spectral_functions();

Analyzing multiple SIGRES files with robots

[back to top]

To analyze the convergence of the QP results, we can use the SigresRobot. Let's build our robot from a list of SIGRES files.

In [12]:
# List of SIGRES files computed with different values of nband.
filenames = [
    "si_g0w0ppm_nband10_SIGRES.nc",
    "si_g0w0ppm_nband20_SIGRES.nc",
    "si_g0w0ppm_nband30_SIGRES.nc",
]

filepaths = [abidata.ref_file(fname) for fname in filenames]

robot = abilab.SigresRobot.from_files(filepaths)

Then we plot the convergence of the QP direct gap as a function of the number of bands in the self-energy for all the k-points available in the netcdf files:

In [13]:
robot.plot_qpgaps_convergence(sortby="sigma_nband", sharey=False);

If we are interested in the convergence of the real/imaginary part of the self-energy and of the renormalization factor ...

In [14]:
robot.plot_qpdata_conv_skb(spin=0, kpoint=(0, 0, 0), band=3, sortby="sigma_nband");

We can also plot the QP data as a function of the KS energies on the same figure with:

In [15]:
robot.plot_qpfield_vs_e0("qpeme0", sortby="sigma_nband");
In [16]:
#robot.get_qpgaps_dataframe(spin=0, kpoint=(0, 0, 0))

Back to the main Index