<jupyter_start><jupyter_text>Contour deformation In the context of GW method, contour deformation (CD) technique is used in conjunction with resolution of identity (RI) to reduce the formal scaling of the self-energy calculation. Compared to widely used analytic continuation approach it provides a means to evaluate self-energy  directly on the real axis without employing Pade approximants or non-linear least squares fit and potentially offering superior accuracy. Here, we provide a brief outline of the theory behind CD and give an example of the self-energy calculation within CD without invoking RI in order to facilitate comparison with the results prsented above. Detailed discussion of the CD can be found in the following papers:1. <NAME>., <NAME>., <NAME>., & <NAME>. (2018). Core-Level Binding Energies from GW : An Efficient Full-Frequency Approach within a Localized Basis. Journal of Chemical Theory and Computation, 14(9), 4856–4869. https://doi.org/10.1021/acs.jctc.8b004582. <NAME>., <NAME>., <NAME>., <NAME>., <NAME>., <NAME>., & <NAME>. (2011). Electronic properties of interfaces and defects from many-body perturbation theory: Recent developments and applications. Physica Status Solidi (B), 248(2), 275–289. https://doi.org/10.1002/pssb.201046094 CD is used to recast the convolution in the GW expression of self-energy as a difference between two integrals, one which can be performed analytically whereas the other can be evaluated numerically on a relatively small grid. This is achieved by closing the inegration contour as shown below [2]:![Integration contour used to evaluate $\Sigma(\omega)$](CD_scheme.jpg)$$\Sigma(r_1,r_2, \omega) = \frac{i}{2\pi} \int_{-\infty}^{+\infty} e^{i\omega^{\prime} \eta} G(r_1, r_2, \omega + \omega^{\prime}) W(r_1, r_2, \omega^{\prime}) d\omega^{\prime}\\ = \frac{i}{2\pi} \oint_{\Gamma} G(r_1, r_2, \omega + z) W(r_1, r_2, z) dz -  \frac{1}{2\pi} \int_{-\infty}^{+\infty} G(r_1, r_2, \omega + i\omega^{\prime}) W(r_1, r_2, i\omega^{\prime}) d\omega^{\prime}$$Depending on the $\omega$ value the lower-left and the upper-right loops of the contour can enclose one or several poles of the zero-order Green's function whereas the poles of the screened Coulomb interaction never fall within the contour. This allowes to evaluate the countour integral as a sum of corresponding residues with apropriate signs (note that the upper-right loop is traversed counter-clockwise, while the lower-left loop is traversed clockwise). The imaginary axis contribution is calculated using Gauss-Legendre grid. Importantly, the intgrals over the arches vanish iff the screened Coulomb interaction does not contain the exchange contribution.<jupyter_code>import psi4
import numpy as np
import scipy as sp
from matplotlib import pyplot as plt
%matplotlib inline
from IPython.core.display import display, HTML

display(HTML("<style>.container {width:95% !important;}</style>"))


psi4.set_options({'basis' : 'def2-qzvp', 'd_convergence' : 1e-7,'scf_type' : 'out_of_core', 'dft_spherical_points' : 974, 'dft_radial_points' : 150 })

he = psi4.geometry("""He  0.0000 0.0000 0.0000
                       symmetry c1
                       units angstrom
""")



psi4.set_output_file('he_qzvp.out')

scf_e, scf_wfn = psi4.energy('PBE', return_wfn=True)


print("DFT energy is %16.10f" % scf_e)
epsilon = np.asarray(scf_wfn.epsilon_a())
print(epsilon*psi4.constants.hartree2ev)<jupyter_output>DFT energy is    -2.8928225474
[-15.74831893   7.47486413  22.72689859  22.72689859  22.72689859
  51.49043714  95.11061037  95.11061037  95.11061037  95.11061037
  95.11061037 116.71002742 116.71002742 116.71002742 256.25880273
 300.00083801 300.00083801 300.00083801 300.00083801 300.00083801
 300.00083801 300.00083801 420.64172788 420.64172788 420.64172788
 420.64172788 420.64172788 461.58835865 461.58835865 461.58835865]<jupyter_text>``` SCF Total Energy (Ha):      -2.8928225473 (MOLGW) ```<jupyter_code>import GW
import time 
start = time.time()
# Analytic calculation for reference purposes
gw_par = {'no_qp' : 1, 'nv_qp' : 1, 'nomega_sigma' : 501, 'step_sigma' : 0.01, 'gl_npoint' : 200, 'analytic_W' : True, 'low_mem' : True }
gw_he_test = GW.GW_DFT(scf_wfn, he, gw_par)
gw_he_test.print_summary()
end = time.time()
print("GW calculation took %10.2f s" % (end - start))<jupyter_output>Number of basis functions:  30
occ/virt: 1/29
Attempting to create RI basis set for DEF2-QZVP (RIFIT)... 
Auxiliary basis set has been generated!
Number of auxiliary basis functions:  47
Fraction of HF exchange is  0.000
Running in production mode!
Shape of the omega_grid_all is  (2, 501)
Caculating GW self-energy via contour deformation
Analytic W has been requested; performing RPA calculation
Shape of omega tensor is  (30, 30, 29)
Calculation of the integral term requires    0.007 Gb
Calculation of the residue term requires     0.018 Gb
Using low-memory algorithm
Finished calculating self-energy
Performing one-shot G0W0
SigX - Vxc
[-0.34438168  0.13933937]
Perfoming graphic solution of the inverse Dyson equation
Done!
E^lin, eV  E^graph, eV  Z 
   -23.672423     -23.476796       0.913907
    11.007629      11.006144       0.989647
GW calculation took       0.29 s<jupyter_text>``` GW eigenvalues (eV)           E0         SigX-Vxc       SigC          Z        E_qp^lin     E_qp^graph   1    -15.748327    -9.371124     0.703643     0.913908   -23.669607   -23.477018   2      7.474860     3.791494    -0.221863     0.989658    11.007575    11.006253```<jupyter_code>start = time.time()
# Analytic calculation for reference purposes
gw_par = {'no_qp' : 1, 'nv_qp' : 1, 'nomega_sigma' : 501, 'step_sigma' : 0.01, 'gl_npoint' : 200, 'analytic_W' : True, 'low_mem' : True, 'evgw_iter' : 5 }
gw_he_test_ev = GW.GW_DFT(scf_wfn, he, gw_par)
gw_he_test_ev.print_summary()
end = time.time()
print("GW calculation took %10.2f s" % (end - start))
# What if I request more qp-s and converge their energies together?
start = time.time()
# Analytic calculation for reference purposes
gw_par = {'no_qp' : 1, 'nv_qp' : 29, 'nomega_sigma' : 501, 'step_sigma' : 0.01, 'gl_npoint' : 200, 'analytic_W' : True, 'low_mem' : True, 'evgw_iter' : 10 }
gw_he_test1_ev = GW.GW_DFT(scf_wfn, he, gw_par)
gw_he_test1_ev.print_summary()
end = time.time()
print("GW calculation took %10.2f s" % (end - start))<jupyter_output>Number of basis functions:  30
occ/virt: 1/29
Attempting to create RI basis set for DEF2-QZVP (RIFIT)... 
Auxiliary basis set has been generated!
Number of auxiliary basis functions:  47
Fraction of HF exchange is  0.000
Running in production mode!
Shape of the omega_grid_all is  (30, 501)
Caculating GW self-energy via contour deformation
Analytic W has been requested; performing RPA calculation
Shape of omega tensor is  (30, 30, 29)
Calculation of the integral term requires    0.007 Gb
Calculation of the residue term requires     0.018 Gb
Using low-memory algorithm
Finished calculating self-energy
Performing one-shot G0W0
SigX - Vxc
[-0.34438168  0.13933937  0.25006238  0.25006238  0.25006238  0.3270321
  0.38471001  0.38471001  0.38471001  0.38471001  0.38471001  0.47563179
  0.47563179  0.47563179  0.78111999  0.55286242  0.55286242  0.55286242
  0.55286242  0.55286242  0.55286242  0.55286242  0.72129897  0.72129897
  0.72129897  0.72129897  0.72129897  0.83403239  0.83403239  0.834[...]
StarcoderdataJupyterStructuredCleanDedup