File: theory_mbt.html

package info (click to toggle)
abinit 7.8.2-2
links: PTS, VCS
area: main
in suites: jessie, jessie-kfreebsd
size: 278,292 kB
ctags: 19,095
sloc: f90: 463,759; python: 50,419; xml: 32,095; perl: 6,968; sh: 6,209; ansic: 4,705; fortran: 951; objc: 323; makefile: 43; csh: 42; pascal: 31
file content (1107 lines) | stat: -rw-r--r-- 39,075 bytes
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head><title>Many-Body Theory in ABINIT</title></head><body bgcolor="#ffffff">
<hr>
<h1>Many-Body Theory in ABINIT</h1>

<hr>
<p>
The aim of this section is to introduce the Green's function formalism, the concept of self-energy
and the set of coupled equations proposed by Hedin whose self-consistent solution, in principle, 
gives the exact Green's function of the interacting system.
We mainly focus on the aspects of the theory that are important for understanding the 
different steps of the calculation and the role played by the input variables used to control the run.
A much more consistent and rigorous introduction to the physical concept of 
Green's function and self-energy can be found in any standard textbook on Many-Body theory, see for example [1,2,3].

</p>
<h5>Copyright (C) 2000-2014 ABINIT group (MG,MS)
<br> This file is distributed under the terms of the GNU General Public License, see
~abinit/COPYING or <a href="http://www.gnu.org/copyleft/gpl.txt">
http://www.gnu.org/copyleft/gpl.txt </a>.
<br> For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
</h5>

<script type="text/javascript" src="list_internal_links.js"> </script>

<h3><b>Contents of the lesson:</b></h3>


<ul>
<li><a href="theory_mbt.html#green_self_intro">1</a> The Green's function and the self-energy 
  </li><li><a href="theory_mbt.html#hedins_equations">2</a> Hedin's equations
  </li><li><a href="theory_mbt.html#gw_approximation">3</a> The GW approximation 
  </li><li><a href="theory_mbt.html#perturbative">4</a> Perturbative approach 
  </li><li><a href="theory_mbt.html#RPA_Fourier_space">5</a> The RPA polarizability in Fourier space
  </li><li><a href="theory_mbt.html#oscillator_notes">6</a>  Notes on the calculation of oscilaltor matrix elements
  </li><li><a href="theory_mbt.html#hilbert_transform">7</a>  Hilbert transform method 
  </li><li><a href="theory_mbt.html#evaluation_gw_sigma">8</a> Evaluation of the GW self-energy
  </li><li><a href="theory_mbt.html#contour_deformation">9</a> The contour deformation technique
  </li><li><a href="theory_mbt.html#references">10</a> References

</li>
</ul>

<hr>


<h4><a name="green_self_intro"></a></h4>

<h3><b>1. Green's function and self-energy.</b></h3>


<p>
The time-ordered Green's function <i>G</i>(12), also called the propagator, 
defines the probability amplitude for the propagation of an added or removed electron in a many-body system. 
Since the probability amplitude is simply given by the overlap between the final and the initial state, 
<i>G</i>(12) can be expressed as 

</p>
<p align="center">

    <img src="./theory_mbt/G_def.png">
    
</p>


where the matrix element is taken in the Heisenberg representation, T is the time-ordering operator 
and the creation and annihilation field operators act on the the ground state of the <i>N</i>-electron many-body Hamiltonian
(the conventions used in the equations are explained in the section on
<a href="theory_mbt.html#notations">notation</a>).
The propagator contains only part of the full information carried by 
the many-body wave function, but it includes the relevant part
for the study of charged excitations. Also, any single-particle operator
acting on the system can be evaluated once the Green's function is known.
<p>
Useful physical information about the charged excitation energies 
of the many-body system can be obtained by expressing the propagator in the so-called Lehmann representation [1-2-3].
To this purpose it is useful to introduce the following notation to denote the 
charged excitation energies of the <i>N</i>-electron system[5]:

</p>
<p align="center">

    <img src="./theory_mbt/charged_exc_ene.png">
    
</p>


where <i>E<sub>0</sub><sup>N</sup></i> is the ground state energy of the electron system with <i>N</i>
electrons, and <i>i</i> is the set of quantum numbers labeling the excited states with <i>N</i>1 electrons. 
Finally, <i>&#956;</i> is the chemical potential of the system.
Other important quantities that will be used in the following are
the so-called Lehmann amplitudes defined, in the Schrdinger representation, by

<p align="center">

    <img src="./theory_mbt/lehmann_ampl.png">
    
</p>


The Lehmann representation of the Green's function

<p align="center">

    <img src="./theory_mbt/Lehmann.png">
    

</p>


makes it clear that, in the frequency domain, the time-ordered Green's function 
contains the complete excitation spectrum corresponding to excitations of an (<i>N</i>-1)-particle and an (<i>N</i>+1)-particle
system. Hence, locating the poles of the Green's function in the
complex plane provides the information needed to interpret those
processes measured in experiments in which a single electron is
inserted to or removed from the system. The figure below gives a
schematical representation of the location of the poles of the
time-ordered Green's function.
<p align="center">

    <img src="./theory_mbt/poles_G.png">
    
</p>


<p>

Where the ionisation potential is the energy required to remove an electron from 
the system, the electron affinity to add an electron, and the chemical potential 
is typically taken to be in the middle of the gap. For a metallic system these 
energies are all equal to each other, and there is no gap.
</p>


<p>
<!--
The equation of motion technique tries to solve the problem by differentiating the Green's function a number of times.
This approach leads to a series of coupled differential equations involving functions with a progressively increasing number of particles
whose solution has the same complexity as the original Schr\"odinger problem.
The calculation of the fully interacting Green's function can be thus reduced 
to the problem of finding the self-energy of the system.
However, at this stage, the question of how G can be calculated remains.
-->

The Dyson equation 

</p>
<p align="center">

    <img src="./theory_mbt/G_Dyson.png">
    
</p>


establishes a connection between the fully interacting <i>G</i> and the propagator, <i>G<sub>0</sub></i>, 
of an approximate non-interacting system
through a (non-local, non-Hermitian and time dependent) potential called the self-energy, &#931;.
Since <i>G<sub>0</sub></i> is supposed to be known exactly, the problem of calculating <i>G</i>(12) 
has now been reduced to the calculation of the self-energy.

<p>
The self-energy is not a mere mathematical device used in a roundabout way to obtain <i>G</i> but is has a direct physical meaning.
The knowledge of the self-energy operator, allows one to 
describe the quantum-mechanical state of a renormalized electron in the many-body system
by solving the quasiparticle (QP) equation[5]:

</p>
<p align="center">

    <img src="./theory_mbt/qp_equation.png">
    
</p>


<!--
Although the above equation looks similar to a mean-field approach,
it does not constitute a mean-field formulation since the self-energy takes  all dynamic many-electron processes into account.
-->
The QP eigenstates so obtained can be used to construct <i>G</i><!--
Unlike the KS approach, the QP equation involves the non Hermitian Sigma and therefore
the eigenvalues are complex-valued.
-->
according to the Lehmann representation. Note that the QP equation
departs from the Kohn Sham equation since the QP eigenvectors and
eigenvalues do have a direct physical meaning: they can be used to
obtain both the charge density of the interacting system and to
describe the properties of charged excitations.
<br>
<br>
<br>
<hr><br>


<h3><a name="hedins_equations"></a></h3>

<h3><b> 2. Hedin's equations.</b></h3>


<p>
In 1965 Hedin[6] showed how to derive a set of coupled integro-differential equations
whose self-consistent solution, in principle, gives the exact self-energy of the system and
therefore the exact <i>G</i>.
<!--
W(12) as fundamental building block.
The dynamically screened interaction generated at position $(1)$ by an electron at $(2)$
is expressed in terms of the bare Coulomb term and the inverse dielectric matrix as
-->
</p>
<p>
The fundamental building blocks employed in the formalism are 
the irreducible polarizability:

</p>
<p align="center">

    <img src="./theory_mbt/tchi_def.png">
    
</p>


which describes the linear response of the density to changes in the total effective potential
(the superposition of the external potential plus the internal classical Hartree potential)
and the dynamically screened interaction, <i>W</i>, that is related to the bare Coulomb interaction, <i>v</i>, and
to the inverse dieletric function through:

<p align="center">

    <img src="./theory_mbt/W_def.png">
    

</p>


The dielectric matrix &#1108;(12) is related to the  irreducible polarizability <i>&#967;</i>(12) by the following relation:

<p align="center">

    <img src="./theory_mbt/e_from_tchi.png">
    
</p>


<p>
The pentagon sketched in the figure below shows how the various physical quantities are interrelated: 

</p>
<p align="center">


    <img src="./theory_mbt/hedin_pentagon.png">
     
</p>


The polarization function renormalises the bare interaction resulting in the screened interaction <i>W</i>(12). 
The screened interaction, <i>W</i>(12), the many-body propagator <i>G</i>(12), and the vertex function, &#915;(12;3), 
which describe the interactions between virtual hole and electron excitations [5],
are the essential ingredients for the determination of &#931;(12). 
<p>
The iteration starts by setting <i>G</i> = <i>G<sub>0</sub></i>. Then the set of equations should in principle be iterated until self-consistency in all terms is reached. 

<br><br><br></p>
<hr><br>


<h3><a name="gw_approximation"></a></h3>

<h3><b> 3. The GW approximation.</b></h3>


<p>The practical solution of Hedin's equations is extremely complicated
as
they are not just numerical relations but contain a functional
derivative in the equation for the vertex.
The direct evaluation of the vertex function is very challenging.
The set of equations can, however, be iterated assuming that only a few
iterations are actually needed to obtain physically meaningful results.
</p>
<p>
A widely used approach to the approximate solution of Hedin's equations is the so-called 
<i>GW</i> approximation[6], which consists in approximating the vertex function with a local and instantaneous function:

</p>
<p align="center">

    <img src="./theory_mbt/gamma_GW.png">

</p>


This approximated vertex, once inserted in the full set of Hedin's equations, leads to a considerable
simplification in the set of equations:

<p align="center">

    <img src="./theory_mbt/gw_pentagon.png">
     
</p>


<p>


Thanks to the neglect of vertex corrections, the irreducible polarizability <i>&#967;</i>(12) is now given by

</p>
<p align="center">

    <img src="./theory_mbt/RPA_chi.png">
     
</p>


which, once rewritten in terms of orbitals and energies, reduces to the RPA expression proposed by Adler [7] and Wiser [8].
<p>
In real space, the self-energy reduces to a simple direct product of the dressed electron
propagator, <i>G</i>(12), and the dynamically screened interaction, <i>W</i>(12):

</p>
<p align="center">

    <img src="./theory_mbt/gwsigma_rt_domain.png">
     
</p>


The self-energy, a simple product in the space-time domain, becomes a convolution when expressed in frequency-space:

<p align="center">

    <img src="./theory_mbt/self_def.png">
     

</p>


<p>
<!--dressed 
In practice, the electron propagator $G(12)$ is never explicitly known and has to be
approximated with the Green's function, $\Go(12)$, of an appropriate non-interacting system. 
-->
</p>
<p>
Ideally, the set of GW equations should still be iterated until self-consistency 
in all terms is reached; this is the fully self-consistent GW method (SCGW).
<!--
Hereafter this approach will be denoted with the name $SCGW$. 
$SCGW$ is the most satisfactory implementation from a theoretical point of view since it preserves energy, electron 
number and total momentum~\cite{Martin1959,Baym1961,Baym1962}.
-->
However SCGW calculations for real systems are still very challenging, and very few have been reported in the literature 
<!--
(see for example~\cite{Holm1998,Holm2000,Garcia-Gonzalez2001} for applications to 
the homogeneous electron gas and~\cite{Shirley1996,Schone1998,Kutepov2010} for calculations on simple metals and semiconductors).
-->
Moreover, the utility of fully SCGW results are still under debate within the scientific community.
</p>

<p>The problem is that self-consistency typically improves total
energies, but worsens spectral properties
(such as band gaps and optical spectra). Since obtaining the spectral
information is often the main reason for doing such difficult
calculations in the first place, many authors agree that a useful
self-consistent approach would need the inclusion of some kind of
vertex correction during the solution of the equations.
</p>

<p>
For this reason, the most common approach employed in the <i>ab initio</i> community
consists of using the best available approximation for <i>G </i>and <i>W</i> as a starting point 
and performing only a single-iteration of the parallelogram (the so-called one-shot <i>GW</i> method, or <i>G<sub>0</sub>W<sub>0</sub></i>).
In this case the self-energy is simply given by:

</p>
<p align="center">

    <img src="./theory_mbt/G0W0.png">
     
</p>


where <i>G<sub>0</sub><sup>KS</sup></i>(12) is the independent-particle propagator of the Kohn-Sham (KS) Hamiltonian,
and the screened interaction is approximated with the RPA calculated with KS energies and wave functions:

<p align="center">

    <img src="./theory_mbt/G0G0.png">
     
</p>


<br>
<br>
<br>
<hr><br>


<h3><a name="perturbative"></a></h3>

<h3><b>4. Perturbative approach.</b></h3>


Despite all the fundamental differences between many-body theory and DFT, 
the Kohn-Sham exchange-correlation potential can be seen as a static, local and 
Hermitian approximation to the self-energy. 
Indeed, in many cases the Kohn-Sham energies already provide a reasonable estimate of the band structure
and are usually in qualitative agreement with experiment. 
<p>
This observation suggests that a simple, albeit accurate, solution for the QP energies can be obtained using 
first-order perturbation theory, treating the exchange and correlation potential, <i>V<sub>xc</sub></i>, 
as a zeroth-order approximation to the non-local and energy dependent 
self-energy [11,12]

<!--
<p>
The GW approximation is usually formulated using the KS Green's function as the starting point 
and the KS exchange-correlation potential has to be subtracted later from the resulting self-energy.
are only a correction to the Kohn-Sham eigenvalues, 
-->
</p>
<p>
Under the assumption that the QP wavefunctions equal the KS orbitals,
we can expand the self-energy operator around <i>&#949;<sup>KS</sup></i> obtaining a closed expression for <i>&#949;<sup>QP</sup></i> :

</p>
<p align="center">

    <img src="./theory_mbt/qpene_perturbative.png">
    

</p>


where 

<p align="center">

    <img src="./theory_mbt/Z_def.png">
    
</p>
is the so-called renormalization factor. This corresponds to making a
Taylor expansion of the self-energy matrix element around the KS
energy, as depicted below.

<p align="center">

    <img src="./theory_mbt/self_energy_taylor.png">

    
</p>



<br>
<br>
<br>
<hr><br>


<h4><a name="RPA_Fourier_space"></a></h4>

<h4><b>5. The RPA polarizability in Fourier space.</b></h4>


In the reciprocal space and frequency domain (implying a Fourier transform (FT) of 
the real space coordinates and time variables), the independent-particle 
polarizability assumes the form:

<p align="center">

    <img src="./theory_mbt/chi0_sc_tr.png">


</p>


where only the transitions between valence (<i>v</i>) and conduction states (<i>c</i>) contribute
(for simplicity we have assumed a semiconductor with time-reversal invariance, 
the conventions used for the Fourier transform are discussed in the <a href="theory_mbt.html#notations">notation</a> section).
<p>
The number of bands used to compute the polarizability is specified by <a href="../input_variables/varbas.html#nband" target="kwimg">nband</a>,
while <a href="../input_variables/vargw.html#zcut" target="kwimg">zcut</a> gives the small complex
shift used to avoid the divergences in the denominators.
The frequency mesh is defined by the set of variables

<a href="../input_variables/vargw.html#nfreqre" target="kwimg">nfreqre</a>,
<a href="../input_variables/vargw.html#nfreqim" target="kwimg">nfreqim</a>,
<a href="../input_variables/vargw.html#freqremax" target="kwimg">freqremax</a>, and
<a href="../input_variables/vargw.html#freqremin" target="kwimg">freqremin</a>
(a number of more exotic grid choices are available through options beginning with 
<a href="../input_variables/vargw.html#gw_frqim_inzgrid" target="kwimg">gw_... or cd_...</a>).
</p>
<p>
<i>M</i> is a shorthand notation to denote the matrix element of a plane wave sandwiched 
between two wavefunctions (i.e. oscillator matrix elements).
The number of planewaves (PW) used to describe the wavefunctions is determined by 
<a href="../input_variables/vargw.html#ecutwfn" target="kwimg">ecutwfn</a>

while the number of <i>G</i>-vectors used to describe the polarizability 
(i.e. the number of <i>G</i> vectors in the oscillator matrix elements) is 
determined by <a href="../input_variables/vargw.html#ecuteps" target="kwimg">ecuteps</a>.
The oscillators are ubiquitous in the Many-Body part of ABINIT and their 
calculation represents one of the most CPU intensive part of the execution.
For this reason we devote section <a href="theory_mbt.html#oscillator_notes">6</a> 
to the discussion of some important technical details concerning their computation.
</p>
<p>
In principle, the set of <b>q</b>-points in the screening matrix is given by all 
the possible differences between two crystalline momenta of the wavefunctions 
stored in the KSS file, so it is controlled by the chosen <b>k</b>-point grid. 
The code, however, exploits the invariance of the two-point function under 
the action of any symmetry operation of the crystalline space group:

</p>
<p align="center"><img src="./theory_mbt/chi0_symmetry.png">

    <br>
</p>


so that only the <b>q</b>-points in the irreducible Brillouin zone (IBZ) have to be calculated explicitly. 
<!--
The use of symmetries is switched on/off with 
<a href="../input_variables/vargw.html#awtr" target="kwimg">awtr</a> (time-reversal), and 
<a href="../input_variables/vargw.html#symchi" target="kwimg">symchi</a> (crystal symmetries).
-->
<p>
In frequency and reciprocal space, the microscopic dieletric function is related to the
irreducible polarizability by the following relation

</p>
<p align="center">

    <img src="./theory_mbt/epsilon_vs_tchi.png">


</p>


from which the inverse dielectric function is obtained via matrix inversion.

Following Adler [7,8], the macroscopic dielectric function, &#1108;<sub>M</sub><sup>LF</sup>(&#969;), can be directly related to 
the inverse of the microscopic dielectric matrix by means of:

<p align="center">

    <img src="./theory_mbt/emacro_lfe.png">

</p>


The optical absorption spectrum -- the quantity one can compare with experiments --
is given by the imaginary part: Im[&#1108;<sub>M</sub><sup>LF</sup>(&#969;)].
<p>
Note that the equation above differs from 

</p>
<p align="center">

    <img src="./theory_mbt/emacro_nlfe.png">
    
</p>


due to the so called local-field effects introduced by the presence of the crystalline environment. 
These spectra, if calculated, are typically output as ...<b>_LF</b> and ...<b>_NLF</b> 

files during the course of a calculation.

<br>
<br>
<br>
<hr><br>


<h3><a name="oscillator_notes"></a></h3>

<h3><b>6. Notes on the calculation of the oscillator matrix elements.</b></h3>


Many body calculations require the evaluation of integrals involving the oscillator matrix elements

<p align="center">

    <img src="./theory_mbt/oscillator_longv.png">
    
</p>


where the <b>k</b>-point belongs to the full Brillouin zone.
<p>
These terms are evaluated by performing a Fast Fourier Transform (FFT) of the real space product of the two wavefunctions 
(the second expression in the equation above).
Thanks to the FFT algorithm, the CPU-time requirement scales almost linearly with the the 
number of points in the FFT box, moreover the code implements refined algorithms 
(for instance zero-padded FFTs, FFTW3 interface) to optimize the computation.
</p>
<p>
There can be a significant speed-up in this component depending on the numerical FFT library used. 
If possible, it should always be advantageous to link and use the FFTW3 library in GW calculations 
(controlled by setting 
<a href="../input_variables/vardev.html#fftalg" target="kwimg">fftalg</a> 312). The performance of
the various FFT libraries for a given type of calculation can be benchmarked with the <b>fftprof</b> utility.
</p>
<p>

For a given set of indeces (<i>b<sub>1</sub></i>, <i>b<sub>2</sub></i>, <b>k</b>, <b>q</b>), 
the calculation of the oscillator is done in four different steps:

</p>
<ol>
<li>
The two wavefunctions in the irreducible wedge are FFT transformed from the <b>G</b>-space to the real space representation, 
</li><li> 
The orbitals are rotated in real space on the FFT mesh to obtain the points <b>k</b> and <b>k-q</b> in the full Brillouin zone.
</li><li>
Computation of the wavefunction product.
</li><li>
FFT transform of the product to obtain <i>M</i>

</li>
</ol>

<p>
Each oscillator thus requires three different FFTs (two transforms to construct 
the product, one FFT to get M).
The number of FFTs can be significantly reduced by precomputing and storing in memory 
the real space representation of the orbitals at the price of a reasonable increase of 
the memory allocated. However, for very memory demanding calculations, the real space orbitals
can be calculated on the fly with an increase in computational time instead. This option 
is controlled by the second digit of the input variable 
<a href="../input_variables/vargw.html#gwmem" target="kwimg">gwmem</a>.
</p>
<p>
The third term in the equation defining the oscillators makes it clear that the product 
of the periodic part of the orbitals has non-zero Fourier components in a sphere 
whose radius is 2<i>R<sub>wfn</sub></i> where <i>R<sub>wfn</sub></i> is the radius of 
the <b>G</b>-sphere used for the wavefunctions (set by

<a href="../input_variables/vargw.html#ecutwfn" target="kwimg">ecutwfn</a>). To avoid 
aliasing errors in the FFT one should therefore use an FFT box that encloses the 
sphere with radius 2<i>R<sub>wfn</sub></i>, 
but this leads to a significant increase in the computing effort as well as in the 
memory requirements. The input variable 
<a href="../input_variables/vargw.html#fftgw" target="kwimg">fftgw</a> specifies
how to setup the FFT box for the oscillators and should be used to test how the aliasing 
errors affect the final results. The default setting of <b>fftgw 21</b> is safe, a setting 
of <b>fftgw 11</b> is fast but can be inaccurate, and a setting of <b>fftgw 31</b> gives 
the maximum possible accuracy at a significant computational cost.

<br><br><br></p>
<hr><br>


<h3><a name="hilbert_transform"></a></h3>

<h3><b>7. The Hilbert transform method.</b></h3>


<p>
The computational effort for the evaluation of the RPA polarizability 
with the Adler-Wiser expression scales linearly with the number of frequencies computed 
(<a href="../input_variables/vargw.html#nfreqre" target="kwimg">nfreqre</a> and
<a href="../input_variables/vargw.html#nfreqim" target="kwimg">nfreqim</a>),
albeit with a large prefactor which increases with the fourth power of the number of atoms. 
The main reason for the linear scaling is that the frequency dependence cannot be factorised 
out of the sum over transitions, hence a distinct and expensive summation over transitions 
has to be performed separately for each frequency.
</p>
<p>
This linear scaling represents a serious problem, especially when many frequencies are wanted, for
example when computing QP corrections within the contour deformation technique described in the <a href="lesson_gw2.html">second lesson</a> of the GW tutorial.

</p>
<p>
This computational bottleneck can be removed, under certain circumstances, by employing an efficient
algorithm proposed in [13] and subsequently revisited in [14], in which only the spectral function 

</p>
<p align="center">

    <img src="./theory_mbt/sf_chi0.png">
    
</p>


has to be evaluated in terms of electronic transitions between valence and conduction states.
The Dirac delta function can be approximated either by means of a triangular function centered 
at the energy transition following [14] or a gaussian approximant 
following [13] (see the related input variables 
<a href="../input_variables/vargw.html#spmeth" target="kwimg">spmeth</a>, and
<a href="../input_variables/vargw.html#spbroad" target="kwimg">spbroad</a>).
The spectral function is evaluated on a linear frequency mesh which covers the entire set
of trasition energies included in the calculation. 
The number of points in the mesh is given by
<a href="../input_variables/vargw.html#nomegasf" target="kwimg">nomegasf</a>.

<p>
The evaluation of the spectral function is rather efficient thanks to the presence of the 
delta-function in the expression above. For example, when 
<a href="../input_variables/vargw.html#spmeth" target="kwimg">spmeth</a>=1,
the CPU time required to compute the spectral function on an arbitrarily dense frequency 
mesh is just twice that required by a single static computation based on the standard 
Adler-Wiser expression.
</p>
<p>
The full polarizability is then efficiently retrieved by means of a less expensive frequency integration (a Hilbert transform):

</p>
<p align="center">

    <img src="./theory_mbt/hilbert_transform.png">
    
</p>


<p>
The price to be paid, however, is that a large table for the spectral function 
has to be stored in memory and a Hilbert transform has to be performed for each pair 
(<b>G</b><sub>1</sub>, <b>G</b><sub>2</sub>).
Since the computing time required for the transform scales quadratically with the number of 
vectors in the polarizability (governed by 

<a href="../input_variables/vargw.html#ecuteps" target="kwimg">ecuteps</a>),
the CPU time spent in this part will overcome the computing time of the standard Adler-Wiser 
formalism for large <a href="../input_variables/vargw.html#ecuteps" target="kwimg">ecuteps</a>.
A theoretical estimate of the crossover point is hard to give because it depends on many 
factors. However, if many frequencies are needed, such as for the evaluation of optical spectra, 
or accurate contour deformation integrations, or even mapping full grids in the complex plane, 
the Hilbert transform method can be significantly faster, and its use is well worth considering.

<br><br><br></p>
<hr><br>


<h3><a name="evaluation_gw_sigma"></a></h3>

<h3><b>8. Evaluation of the GW self-energy.</b></h3>


<p>
Following the standard approach, we separate the screened interaction into the static bare Coulomb term 
and a frequency-dependent contribution according to:

</p>
<p align="center">

    <img src="./theory_mbt/W_partition_simple.png">

    
</p>



where matrix notation is used.
<p>
This particular decomposition of <i>W</i>, once inserted in the convolution defining &#931;, 
leads to the split of the self-energy into two different contributions (exchange and correlation):

</p>
<p align="center">

    <img src="./theory_mbt/self_energy.png">
    

</p>



<p>
The exchange part is static and turns out to have the same mathematical structure as the Fock operator in Hartree-Fock theory, 
albeit constructed with quasiparticle amplitudes

</p>
<p align="center">

    <img src="./theory_mbt/self_x.png">
    
</p>


while the dynamic part &#931;<sub>c</sub>(&#969;) accounts for correlation effects beyond &#931;<sub>x</sub>.

<p>It is important to stress that, for computational efficiency, the
code does not compute the full self-energy operator by default.
Only its matrix elements for the states specified by <a href="../input_variables/vargw.html#kptgw" target="kwimg">kptgw</a> and 
<a href="../input_variables/vargw.html#bdgw" target="kwimg">bdgw</a>
are computed and used to obtain the QP corrections.
</p>
<p>
When expressed in reciprocal space, the diagonal matrix elements of the exchange part are given by:

</p>
<p align="center">

    <img src="./theory_mbt/self_x_mel.png">
    

</p>


The evaluation of these terms represents a minor fraction of the overall CPU time since only 
occupied states are involved. However we should always keep in mind that, due to the long 
range of the bare Coulomb interaction, the convergence with respect to the number of plane 
waves used in the oscillators <i>M</i> 
(<a href="../input_variables/vargw.html#ecutsigx" target="kwimg">ecutsigx</a>) is usually 
slow, much slower than the convergence of the correlation part, which is short-ranged. This plane
wave cutoff can be converged independently of others if nescessary, and given a much larger value 
in comparison to
<a href="../input_variables/varbas.html#ecut" target="kwimg">ecut</a>, 
<a href="../input_variables/vargw.html#ecutwfn" target="kwimg">ecutwfn</a> and 
<a href="../input_variables/vargw.html#ecuteps" target="kwimg">ecuteps</a>.
<p>
Another point worth noting is the presence in the expression of the Coulomb 
singularity for |<b>q</b>+<b>G</b>| &#8594; 0. From a mathematical point of view, 
the integral is well-defined since the singularity is integrable 
in three-dimensional space once the thermodynamical limit, <i>N</i><sub><b>q</b></sub> &#8594; &#8734;, is reached.
On the other hand, only a finite number of <b>q</b>-points can be used for practical applications,
and a careful numerical treatment is needed to avoid an exceedingly slow convergence with respect to the BZ sampling. 
To accelerate the convergence in the number of <b>q</b>-points, the code implements several 
techniques proposed in the literature. We refer to the documentation of 

<a href="../input_variables/vargw.html#icutcoul" target="kwimg">icutcoul</a> for 
a more extensive discussion.
</p>
<p>
The expression for the matrix elements of the correlation part is instead given by:

</p>
<p align="center">

    <img src="./theory_mbt/self_c_mel.png">
    
</p>


<!-- Safari does not recognize the code &#x1D4A5 for capital J -->
where all dynamical effects are now contained in the frequency convolution integral <i>J</i>.

<p>
The explicit expression for <i>J</i> depends on the method used to treat the screened interaction. 
The code implements four different plasmon-pole techniques to model the frequency dependence of
<b>W</b> in an efficient but approximate way, alternatively, it is possible to use the more sophisticated 
frequency integration of the contour deformation method [15]
for accurate QP calculations (see the related variables 
<a href="../input_variables/vargw.html#ppmodel" target="kwimg">ppmodel</a> and
<a href="../input_variables/vargw.html#gwcalctyp" target="kwimg">gwcalctyp</a>).
</p>
<p>
The double sum over <b>G</b>-vectors is performed for all the plane waves contained within
a sphere of energy

<a href="../input_variables/vargw.html#ecuteps" target="kwimg">ecuteps</a> 
(it cannot be larger than the value used to generate the SCR file).
For each state, the correlated matrix elements are evaluated on a linear frequency mesh centered 
around the initial KS energy and the derivative needed for the renormalization 
factor is obtained numerically (see
<a href="../input_variables/vargw.html#nomegasrd" target="kwimg">nomegasrd</a> and
<a href="../input_variables/vargw.html#omegasrdmax" target="kwimg">omegasrdmax</a>).

</p>
<p>
Note that here, in contrast to the exchange term, the sum over the band index <i>n</i> should extend up to 
infinity although in practice only a finite number of states 
can be used (specified by <a href="../input_variables/varbas.html#nband" target="kwimg">nband</a>).
</p>
<p>


<br><br><br></p>
<hr><br>


<h3><a name="plasmonpole_models"></a></h3>

<h3><b>9 Plasmon-pole models</b></h3>


One of the major computational efforts in self-energy calculations is represented by the calculation
of the frequency dependence of the screened interaction, which is needed for the evaluation of the convolution.
Due to the ragged behavior of <i>G</i>(&#969;) and <i>W</i>(&#969;)
along the real axis, numerous real frequencies are in principle
required to converge the results (note that the maximum frequencies
needed are now reported if <a href="../input_variables/varfil.html#prtvol" target="kwimg">prtvol</a> &gt; 9). 
On the other hand, since the fine details of <i>W</i>(&#969;) are integrated over, 
it is reasonable to expect that approximate models, able to capture the main physical features of 
the screened interaction, should give sufficiently accurate results with a considerably reduction 
of computational effort.
This is the basic idea behind the so-called plasmon-pole models in which the frequency dependence 
of <i>W</i>(&#969;) is modelled in terms of analytic expressions depending on coefficients 
that are derived from first principles, i.e. without any adjustable external parameters.

<p>
Four different plasmon-pole techniques are available in ABINIT and 
the input variable <a href="../input_variables/vargw.html#ppmodel" target="kwimg">ppmodel</a> selects the method to be used.
</p>
<p>
When <a href="../input_variables/vargw.html#ppmodel" target="kwimg">ppmodel</a>=1,2
the frequency dependence of the inverse dielectric function is modeled according to

</p>
<p align="center">

    <img src="./theory_mbt/ppmodel_imag.png">
    
</p>



<p align="center">

    <img src="./theory_mbt/ppmodel_real.png">
    
</p>


The two models differ in the approach used to compute the parameters. 
<a href="../input_variables/vargw.html#ppmodel" target="kwimg">ppmodel</a>=1
derives the parameters such that the inverse dielectric matrix is correctly reproduced at 
two different explicitly calculated frequencies: 
the static limit (&#969; = 0) and an additional imaginary point located at 
<a href="../input_variables/vargw.html#ppmfrq" target="kwimg">ppmfrq</a>. Unless the user overrides 
this, the default value is calculated from the average electronic density of the system.
The plasmon-pole parameters of 
<a href="../input_variables/vargw.html#ppmodel" target="kwimg">ppmodel</a>=2
are calculated so as to reproduce the static limit exactly and to 
fulfill a generalized frequency sum rule relating the imaginary part of the many-body inverse 
dielectric matrix to the plasma frequency and the charge density [12]

<p>
For a discussion of the models corresponding to <a href="../input_variables/vargw.html#ppmodel" target="kwimg">ppmodel</a>=3,4 we refer the reader to the original papers cited in the documentation
of the variable.

<br><br><br></p>
<hr><br>


<h3><a name="contour_deformation"></a></h3>

<h3><b>9. Contour deformation technique.</b></h3>


The contour deformation method was proposed in order to avoid having to deal with quantities close
to the real axis as much as possible [15].
The integral over the real frequency axis can be transformed into an integral over the contour 
depicted in red in the figure below. The integral over real frequency is traded with an integration 
along the imaginary axis plus contributions coming from the poles of the integrand lying inside 
the contour:

<p align="center">

    <img src="./theory_mbt/contour.png">

</p>


<p align="center">

    <img src="./theory_mbt/GW_CD.png">
    
</p>


In the above equation, the first sum is restricted to the poles lying inside the path 
<i>C</i><!-- Safari does not render it correctly.
&#x1D436;. 
-->.

<i>W</i><sub>c</sub>(<i>z</i>) represents the frequency dependent part of the screened interaction, whose expression in reciprocal space is given by:

<p align="center">

    <img src="./theory_mbt/Wc.png">

</p>


The integration along the imaginary axis is expected to converge quickly with respect to the 
number of  sampled frequencies since the integrand is typically very smooth.
Only the residues of the integrand have to be evaluated at the complex poles contributed 
by the Green's function whose frequency dependence is known.

<br>
<br>
<br>
<hr><br>


<h3><a name="notations"></a></h3>

<h3><b>10. Notation.</b></h3>


The following shorthand notations are employed:

<p align="center">

    <img src="./theory_mbt/notation_1.png">
    
</p>


where <i>v</i>(<b>r</b><sub>1</sub>,<b>r</b><sub>2</sub>) represents the bare Coulomb interaction, and &#951; is a positive infinitesimal. 

<p>

The Fourier tranforms for periodic lattice quantities are defined as

</p>
<p align="center">

    <img src="./theory_mbt/notation_2.png">
    
</p>


The volume of the unit cell is denoted with &#937;, while <i>V</i> is the
total volume of the crystal simulated employing Born-von Karman periodic boundary condition.
Unless otherwise specified, Hartree atomic units will be used throughout.

<hr>

<h3><p><b><a name="references">References</a>
</b></p></h3>


<ol>
<li><a name="Abrikosov1975">
  Abrikosov, Gorkov, and Dzyaloshinskii.
  Methods of quantum field theory in statistical physics,
  Dover, New York,
  (1975)
</a> </li><li><a name="Fetter1971">
  Fetter, and Walecka. 
  Quantum Theory of Many-Particle Systems,
  McGraw-Hill, New York,
  (1971)
</a> </li><li><a name="Mattuck1992">
  R.D Mattuck.
  A guide to Feynman diagrams in the many-body problem,
  Dover, New York,
  (1992)

</a> </li><li><a name="Aulbur2000">
  W. G. Aulbur, 
  L. Jnsson 
  and J. W. Wilkins, 
  Solid State Physics <b>54</b>, 1 (2000)
</a></li><li><a name="Onida1995">
  G. Onida et al.
  Phys. Rev. Lett.
  <b>75</b>, 818
  (1995)
</a> </li><li><a name="Hedin1965">

L. Hedin
Phys. Rev.
<b>139</b>, 
A796, (1965)
</a> </li><li><a name="adler">S. L. Adler, Phys. Rev. <b>126</b>, 413 (1962)</a> </li><li><a name="wiser">N. Wiser, Phys. Rev. <b>129</b>, 72 (1963)</a> </li><li><a name="2009">
  A. Kutepov, S. Y. Savrasov and G. Kotliar
  Phys. Rev. B
  <b>80</b>, 041103(R)
  (2009)
</a> </li><li><a name="Kutepov2009">
  A. Kutepov, S. Y. Savrasov and G. Kotliar
  Phys. Rev. B
  <b>80</b>, 041103(R)
  (2009)
</a> </li><li><a name="Hybertsen1985">
M.S. Hybertsen, and S.G Louie.
Phys. Rev. Lett.
<b>55</b>,
1418
(1985)

</a></li><li><a name="Hybertsen1986">
M.S. Hybertsen, and S.G Louie.
Phys. Rev. B
<b>34</b>,
5390
(1986)
</a></li><li><a name="Miyake2000">
T. Miyake, and F. Aryasetiawan.
Phys. Rev. B,
<b>61</b>,
7172
(2000)
</a></li><li><a name="Shishkin2006">
M. Shishkin, and G. Kresse.
Phys. Rev. B,
<b>74</b>,
035101
(2006)

</a></li><li><a name="Lebegue2003">
S. Lebegue, B. Arnaud, M. Alouani, and P.E. Bloechl.
Phys. Rev. B,
<b>67</b>,
155208
(2003)
</a></li>
</ol>

<script type="text/javascript" src="list_internal_links.js"> </script>

</body></html>