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