Gwex

From Wiki

Brief History

GWEX was named originally GWA and was written by Patrick Hahn as a part of his PhD thesis. Iniatially was intended to calculate the optical absorption spectra of molecules and clusters within the GW approximation (Green function, G, screened Coulomb potential W). Later, excitonic effects within the formalism of the Bethe-Salpeter Equation (BSE) for the polarization were taken into account. The latest version of the code (GWA 0.92) was not officially released or documented. It was intended to work both with a real-space code (North Carolina University, Prof. Bernholc) and with a modified version of the Vienna Ab Initio Simulation Package (VASP). The main routine of GWA 0.92 was written in C language with calls to FORTRAN routines. In the end of 2005, Patrick Hahn left the group and Luis Ramos started to use the code to calculate optical properties of Si nanocrystallites. Since the latest version was not documented and a bit experimental, changes had to be made and some parts of the code were re-written. Essentially variables which were not used were removed, variables and subroutines were renamed, formatting was changed, and the number of output files were extended. The current version of the code is completely written in FORTRAN and is no longer intended to work with the real-space code. As before, the code requires a modified version of the VASP.

Description

GWEX is written in fixed-format FORTRAN, though it needs a FORTRAN 90 compiler. The GWEX and the modified version of VASP do not work in parallel yet. However one will find many indications in the code that suggest a previous parallelization or an intention to do it in the future. To use the GWEX one needs to run the modified version of VASP previously. Since optical-matrix elements are required, PAW data sets are required. GWEX does not work with ultra-soft pseudopotentials. This code is intended to be used in localized systems, such as molecules and clusters. Therefore, it is assumed that there is only one k point in the KPOINTS file and this k point is the $Γ$ point. Only supercells with primitive vectors $(X,0,0)$ , $(0, Y,0)$ , and $(0,0, Z)$ are supported, but the code has only been used in the case where $X = Y = Z$ . In the present version of GWEX, the dimensions of the supercell are read directly from the POSCAR file present in the directory. The INCAR file of VASP should look like

PREC   = High
ENMAX  = 200
LREAL  = On
IALGO  = 48
ISMEAR = 0   ; SIGMA  = 0.001

IBRION =   1
NSW    =   0
POTIM  =  0.1
ISIF   =   0

LOPTICS=.TRUE.

NBANDS   = 50
NBVALOPT = 4    !  number of valence bands
NBCONOPT = 46  !  number of conduction bands
IQUASI   = 6

NGXGW = 48; NGYGW = 48; NGZGW = 48
NGXFGW= 48; NGYFGW= 48; NGZFGW= 48

The most important tag besides the tags for the generation of the file OPTIC are IQUASI=6 and NGXGW, NGYGW, NGZGW, and NGXFGW, NGYFGW, NGZFGW. If the values of the NG's are not large enough VASP usually stops. However, one should choose NGXGW, NGYGW, NGZGW, and NGXFGW, NGYFGW, NGZFGW as small as possible. The VASP run should be static, i.e., the atomic positions should be previously optimized. The number of empty bands should be at least twice the number of occupied bands, though it is convenient to have more than this number. After the VASP run, one finds additionally to the usual VASP outputs the following files

GVECTORS: contains the G vectors, their integer indexes, their components, their square of moduli, and their moduli.
Occ_num: contain the occupation number of each state. Not needed by GWEX.
Transmatrix: contains essentially the information that can be extracted from the binary file OPTIC, i.e., kpoint index, the indexes of the states of an optical transition, their eigenvalues and the real and imaginary parts of the operator $(d / d x, d / d y, d / d z)$ .
RHOMAT: EIGENVAL-like file. Not needed by GWEX.
VXCDATA_DIAG: contains the diagonal matrix elements of the exchange-correlation potential $V X C$ for all the bands, as it is calculated in VASP. In calculations where a diagonal approximation of the self-energy is needed, this file is used as input.
VXCDATA_NONDIAG: contains all the matrix elements of exchange-correlation potential $V X C$ (diagonal and off diagonal) for all the bands, as they are calculated in VASP. This file is read as input in the calculations where off-diagonal matrix elements of $V X C$ are required.

and the directory

TempBlochint: contains the binary files with the Bloch integrals. The files written under this directory are named TempBI_XXXX and TempBICOH_XXXX, the latter are necessary whenever a calculation takes advantage of the completeness identity of the Bloch integrals.

Input File for GWEX

An example of the the input file INPUTGWA is given below

7
4 8 104
80 80 80
400 0.1 20.0

# Calculation_mode (see below)
# valence_bands  gw_bands  total_bands
# NGX NGY NGZ
# Number_of_time_steps energy_step(eV)  energy_limit(for non-linear energy_step)
#
#
#  1:  time integration for self-energy with using completeness where appropriate
#      ( Eq.VII-8 of P. Hahn's thesis )
#  2:  time integration for self-energy WITHOUT using completeness
#      ( Eq. III-89 of P. Hahn's thesis )
#  3:  static part of self-energy with completeness where appropriate
#      ( Eq. VII-6 of P. Hahn's thesis ) COH is diagonal
#  4:  static part of self-energy WITHOUT using completeness
#      ( Eq. III-100 of P. Hahn's thesis ) COH is diagonal
#  5:  diagonal of static part of self-energy with completeness
#  6:  diagonal of static part of self-energy WITHOUT using completeness
#  7:  BSE  (needs one of options from 1 to 6 first)

This input file is as simple as possible. There are essentially two kinds of GW calculations or in other words, two kinds of approaches to calculate the self-energy operator $Σ$ . The ones which take advantage of a completeness identity of the Bloch integrals and the ones which does not. Details can be found in Patrick Hanh's PhD thesis (in German) or in Molecular electronic excitations calculated from a solid-state approach: Methodology and numerics, P.H. Hahn, W.G. Schmidt, and F. Bechstedt, Phys. Rev. B 72, 245425 (2005). In the different calculation modes above the self-energy $Σ$ is either approximated within the COHSEX approximation (Coulomb-hole and Screened-exchange) $Σ = Σ C O H + Σ S E X$ without dynamical corrections (options "3" to "6") or approximated by $Σ = Σ C + Σ X$ with dynamical corrections for the correlation term $Σ C$ included. The exchange part $Σ X$ is frequency indepedent and therefore calculated exactly. The correlation part can be separated further in a static and a dynamical part $Σ C = Σ C - s t a t + Σ C - d y n$ .

The calculation mode is selected in the first line of INPUTGWA. It may be important to keep the output screen of the each run. Options "1" and "2" include the dynamical part of the correlation in the self-energy. In these calculation modes the self-energy is always a diagonal operator. The off-diagonal elements of the self-energy operator are not taken into account in the calculation. Options "1" and "2" are the best estimate of the self-energy operator within the GWEX code. The integration over the frequency is made in such a way that it is possible to perform an iteration with the eigenvalues of the quasiparticle eigenvalue problem. Such iterations are formally needed, since the self-energy operator depends on the energies one has to find in the quasiparticle eigenvalue problem. In previous versions of the code the part corresponding to the Bloch integrals in the integral was stored in the disk for each frequency. In the new version we keep these data in the memory.

Options "3" to "6" do not include any dynamical part of the self-energy. Although they are significant faster than options "1" and "2", they are not complete. One important point to be observed is that whenever the quasiparticle eigensystem is solved in the non-diagonal case (options "3" and "4"), there is the risk of an unphysical lift of the states degeneracies. This can be either due to the small size of the supercell or due to the number of unoccupied states, which have an influence on the calculation. Option "3" and "4" compute the SEX (screened exchange) part of the self-energy completely, i.e. with diagonal and off diagonal elements, together with a diagonal COH (Coulomb hole) part. This is done in this way because we noticed that the off diagonal terms of the COH part leads to unphysical breaking of degeneracies. Options "5" and "6" correspond to a diagonal approach for both SEX and COH. These are very fast calculation modes and can be used in tests or when one has a good estimate for the dynamical part of the correlation, e.g. a good estimate of the derivative of the self-energy versus energy.

After running options "1" and "2", the following files are written:

SIGMA_RE: contains $ε λ + R e < λ | Σ(ε) - V X C | λ >$ ( $ε λ$ are the KS eigenvalues) versus energy for each band. The file format is suitable for the 2D-plot packages xmgrace (or xmgr) and probably gnuplot (default: 5000 energies per band).
SIGMA_IM: contains $I m < λ | Σ(ε) - V X C | λ >$ versus energy for each band. The file format is suitable for xmgrace (or xmgr) and probably gnuplot (default: 5000 energies per band).
SIGMA_ENER: contains the spectral function (see P.H. Hahn's thesis or the reference mentioned above) for each band in the same format as for SIGMA_RE and SIGMA_IM.
SIGMAC_DIAG: contains the static contribution of the correlation $Σ C - s t a t$ .
SIGVXC_DIAG: contains $Σ C - s t a t + Σ X - V X C$ .

After running options "1" to "6", the following files are written:

HF_DIAG or HF_NONDIAG: the so-called Hartree-Fock-like eigenvalues. Currently these files are empty after the run due to some experiments with parallelization. Anyway, the diagonal terms are always displayed in a table in screen output together with KS eigenvalues and other quantities.
QUASIVECXXXXYYYY: list the number of states, eigenvalue and eigenvector of the quasiparticle eigenvalue problem. XXXX is the total number of states (bands) and YYYY is the number of states considered in the self-energy. This file is needed to run option "7".
SEXcomp or SEXwocomp: diagonal elements of the SEX part of the self-energy for the cases with (comp) and without completeness (wocomp).
COHcomp or COHwocomp: diagonal elements of the COH part of the self-energy for the cases with (comp) and without completeness (wocomp).

and the KS spectrum and QP (GW) spectrum are in the files:

EPS_KSYYYY: contains the macroscopic dielectric function within the Kohn-Sham (KS) indepedent-particle approximation in a 3-column format: Photon energy, $R e ε(ω)$ , $I m ε(ω)$ (default: 5000 energies, Lorentzian broadening factor 50 meV). YYYY is the number of bands considered in the self-energy.
EPS_QPYYYY: contains the macroscopic dielectric function within the quasiparticle approximation (GW) in a 3-column format: Photon energy, $R e ε(ω)$ , $I m ε(ω)$ (default: 5000 energies, Lorentzian broadening factor 50 meV). YYYY is the number of bands considered in the self-energy.

The second line of INPUTGWA contains three integer numbers that control respectively the number of occupied bands, the number of bands considered in the calculation of the self-energy, and the total number of bands. The number of occupied bands is fixed and should correspond one half of the electrons in the system. The number of bands considered in the calculation of the self-energy is usually chosen as twice the number of occupied bands or more. The third number can be used to control the number of bands considered in the calculation of the dielectric matrix. Since the diagonalization routines of VASP do not calculate the highest energy states accurately, it is not prudent to include all the bands anyway. By reducing the number of total bands one can save computer time, since the calculation of the dielectric matrix is the most time demanding part of the code. However, one should be careful when reducing the number of bands used to calculate the dielectric matrix, because this may have direct influence on the screening and quasiparticle energies.

After running options "1" to "6", one can select option "7" in INPUTGWA and calculate the spectrum including excitonic effects within the Bethe-Salpeter-Equation (BSE) formalism. The only change in this calculation mode with respect to the first version of the code is that not all the states considered in the self-energy are considered in the solution of the BSE. Due to problems of representation of unoccupied extended states by plane waves, not all the unoccupied states should be taken into account. The determination of which states should be included in the BSE calculation could be done by inspecting the constant value of the average electrostatic at the vacuum region and excluding all the unoccupied states above it. Although this is a more sensible procedure, the number of states found would have to be introduced as a new parameter in INPUTGWA. For that reason we followed a more straightforward procedure, which automatically sets the number of states in the BSE. We look up the KS eigenvalues and determine the number of states below the zero energy. Then we take into account in the solution of the BSE as many states of the quasiparticle eigensystem as this number. After running option "7" the following file is written

EPS_BSEXXXXYYYY: contains the macroscopic dielectric function including quasiparticle (GW) and excitonic effects (BSE) in a 3-column format: Photon energy, $R e ε(ω)$ , $I m ε(ω)$ (default: 5000 energies, Lorentzian broadening factor 50 meV). XXXX is the total number of bands and YYYY is the number of bands considered in the BSE calculation.

In the third line of INPUTGWA the NGX, NGY, and NGZ used in the VASP calculation should be given. In the fourth line of INPUTGWA are the parameters which control the integration over the frequency in the options "1" and "2": The number of steps or number of frequencies considered $(N s t e p s)$ in the integration, the energy step in eV $(E s t e p)$ , and the energy limit in eV for a linear increase of the energy step $(E l i n e a r)$ . The latter parameter was introduced because at high energies the terms contribute less to the integral. After this limit the energy step starts to vary exponentially like $E s t e p (1.01) n$ where $i n t (E l i n e a r / E s t e p) < n < N s t e p s$ .