1#!/usr/local/bin/python3.8 2r""" 3GW with scissors operator 4========================= 5 6This example shows how to generate an energy-dependent scissors operator 7by fitting the GW corrections as function of the KS eigenvalues. 8We then use the scissors operator to correct the KS band structure 9computed on a high symmetry k-path. Finally, the LDA and the QPState band 10structure are plotted with matplotlib. 11""" 12import abipy.data as abidata 13from abipy.abilab import abiopen, ElectronBandsPlotter 14 15# Get the QP results from the SIGRES.nc database. 16sigma_file = abiopen(abidata.ref_file("si_g0w0ppm_nband30_SIGRES.nc")) 17 18# Let's have a look at the QP correction as function of the KS energy. 19# Don't shift KS eigenvalues to have zero energy at the Fermi energy. 20# because we need the absolute values for the fit. 21# The qpeme0(e0) curve consists of two branches: 22# the one in the [-6, 5.7] eV interval associated to valence states 23# and the one in the [6.1, 15] eV interval associated to conduction bands. 24# We will fit these results with 2 functions defined in these two domains. 25sigma_file.plot_qps_vs_e0(e0=None) 26sigma_file.plot_qps_vs_e0() 27 28qplist_spin = sigma_file.qplist_spin 29 30# Define the two domains and construct the scissors operator 31domains = [[-10, 6.1], [6.1, 18]] 32scissors = qplist_spin[0].build_scissors(domains, bounds=None) 33 34# Read the KS band energies computed on the k-path 35ks_bands = abiopen(abidata.ref_file("si_nscf_GSR.nc")).ebands 36 37# Read the KS band energies computed on the Monkhorst-Pack (MP) mesh 38# and compute the DOS with the Gaussian method 39ks_mpbands = abiopen(abidata.ref_file("si_scf_GSR.nc")).ebands 40ks_edos = ks_mpbands.get_edos() 41 42# Apply the scissors operator first on the KS band structure 43# along the k-path then on the energies computed with the MP mesh. 44qp_bands = ks_bands.apply_scissors(scissors) 45 46qp_mpbands = ks_mpbands.apply_scissors(scissors) 47 48# Compute the DOS with the modified QPState energies. 49qp_edos = qp_mpbands.get_edos() 50 51# Plot the LDA and the QPState band structure with matplotlib. 52plotter = ElectronBandsPlotter() 53 54plotter.add_ebands("LDA", ks_bands, edos=ks_edos) 55plotter.add_ebands("LDA+scissors(e)", qp_bands, edos=qp_edos) 56 57# By default, the two band energies are shifted wrt to *their* fermi level. 58# Use e=0 if you don't want to shift the eigenvalus 59# so that it's possible to visualize the QP corrections. 60plotter.combiplot(title="Silicon band structure") 61 62plotter.gridplot(title="Silicon band structure") 63sigma_file.close() 64