GAP2: GW with Augmented Planewaves (version 2) (Last updated: June, 2021)

Table of Contents

Introduction

GAP2 is the second version of the all-electron GW implementation based on the full-potential linearized augmented planewave plus local orbital ((L)APW+lo) method, which handle core, semicore, and valence states on the same footing and therefore allows for a correct treatment of core-valence interaction. It is, therefore, able to handle a wide range of materials, irrespective of their composition. At present, GAP2 is interfaced to the WIEN2k code (currently mainly tested for version 14.2).

The main new features of GAP2 are the following:

An advantage of our code is the capability to explore d- and f-electron systems, materials traditionally categorized as strongly correlated. For such materials, a full-potential all-electron treatment is highly desirable. DFT with LDA or GGA exchange-correlation functionals fails dramatically for many such systems, and so may G0W0 carried out on top of LDA or GGA. Here, much of the problem may stem from the LDA/GGA starting point. As a first step towards establishing G0W0 for d and f-electron systems, we have implemented G0W0 based on LDA+U. This simple and effective approach has been applied to a series of prototype d and f-electron systems and shown to overcome the major shortcomings of LDA/GGA.

Return to ToC

Limitations

  1. The code aims at numerically accurate GW calculations, and as a price of that, the code is still not feasible for large systems. For moderate (high) accuracy, bulk systems with less than 20 (10) atoms/cell are probably affordable.
  2. Application of the code to low dimensional systems (surface, 2D systems, molecules, etc) is not encouraged, and can be VERY expensive to achive adequate numerical accuracy.
  3. The treatment of metallic systems is still in the experimental state, and has not been well tested.

Return to ToC

Installation

  1. Edit make.inc in terms of the system configuration. If using Linux-based systems, one only needs to change the LAPACK definition, and the use of INTEL MKL library is strongly recommended (other LAPACK packages is not tested). A sample of make.inc is also included for AIX and SunOS, but they have not been tested in this version.
  2. Compiling the code by "make": First generate the sequential version by "make seq" and then parallel version with "make para" (only subroutines related to parallelization are recompiled)
  3. Make sure the gap2c directory is included in PATH by  putting something like "export PATH=$HOME/prog/gw/gap2c:$PATH"  into ~/.bashrc, and then "source ~/.bashrc".

Return to ToC

How to use the code for GW calculations?

  1. Run a full SCF calculation using WIEN2k in the directory denoated as w2kdir henceforth
  2. Prepare gw input files by running the shell script in the w2kdir: "gap2_init -f casename -d gwdir -nkp number_of_kpoints_for_gw ......". gwdir is the working diretory for the gw calculation to be carried out. It is strongly recommended that gwdir be different from w2kdir. More detailed information about the options for gap2_init can be obtained by "gap2_init -h"
  3. Modify gwdir/gw.inp if necessary, and then run "gap2.x" or "gap2-mpi.x" (using mpirun) in the gwdir.
  4. GW calculations are usually very cost-demanding, and therefore parallelization is often necessary to obtain useful results within reasonable time. Three-level parallelization is implemented, including the parallelization over the loop of q, k and the summation of upccupied states. In most cases, the default setup based on the number of processes (nproc) and the number of q/k-points should be most effective. It is possible to change the default parallelization setup by setting the input block in gw.inp.
  5. Some useful information is given in the main output file, case.outgw. Further analysis can be done by using the C-shell script gap2_analy (see post-processing).
  6. To add high-energy local orbitals on top of standard (L)APW+lo basis (the SCF calculation still uses the default (L)APW+lo), use the option "-newin1 <nLO>" (see Phys. Rev. B, 93,115203(2016)) for more details) when running gap2_init. You are prompted to edit case.in1, but in general you just check whether there is anything abnormal and do not need to do anything.
  7. To add high-energy local orbitals to larger l channels than l=3, additional care has to be taken.
    • a) Make some changes in WIEN2k as explained as followings:
      • Enter the directory when WIEN2k is stored, and copy SRC_lapw1 and SRC_nmr to SRC_lapw1_lom10 and SRC_nmr_lom10, respectively.
      • Enter SRC_lapw1_lom10, set LOMAX = 10 (and NSLMAX = 11 in old versions, e.g. 14.2)  in param.inc, recompile the lapw1 executables ("make all"), and then copy them to the upper directory by adding a flag "_lom10" to each of them, i.e. cp lapw1 ../lapw1_lom10, etc.
      • Enter SRC_nmr_lom10, set LOMAX = 10 in param.inc (modules.F in more recent versions of wien2k, e.g. 18.2) , and make similar operations as above.
      • Some changes need to be made in x_lapw and x_nmr_lapw. Here are modified x_lapw and x_nmr_lapw for WIEN2k 19.2 (grep 'JH' to find out what changes have been made that need to be made for other versions of WIEN2k).  
    • b) Use “-lom10” when running gap2_init. In this case, you are prompted to edit case.in1 twice.

    You are promted to edit case.in1 twice. For the first time you need to modify it based on up to which l  you would use HLOs. Taking Si as an example, if you want to add HLOs up to l=5, the original Si.in1 reads

         WFFIL  EF=.353692248950   (WFFIL, WFPRI, ENFIL, SUPWF)
         9.00       10    4 (R-MT*K-MAX; MAX L IN WF, V-NMT
         0.30    2  0      (GLOBAL E-PARAMETER WITH n OTHER CHOICES, global APW/LAPW)
         0    0.30      0.000 CONT 1
         1    0.30      0.000 CONT 1
         K-VECTORS FROM UNIT:4   -9.0       2.0    7   red   emin/emax/nband      

    and it has to be modified to

         WFFIL  EF=.353692248950   (WFFIL, WFPRI, ENFIL, SUPWF)
         9.00       10   10 (R-MT*K-MAX; MAX L IN WF, V-NMT
         0.30    6  0      (GLOBAL E-PARAMETER WITH n OTHER CHOICES, global APW/LAPW)
         0    0.30      0.000 CONT 1
         1    0.30      0.000 CONT 1
         2    0.30      0.000 CONT 1
         3    0.30      0.000 CONT 1
         4    0.30      0.000 CONT 1
         5    0.30      0.000 CONT 1
         K-VECTORS FROM UNIT:4      -9.0    1000.0              7 red          
    For the second time to edit case.in1, again you just check whether there is anything abnormal and in most cases you don't need to do anything. Return to ToC

    How to do CRPA calculations to estimate the Hubbard U

    To run constrained random phase approximation (cRPA) to obtain the Hubbard U (L. Vaugier, H. Jiang and S. Biermann, Phys. Rev. B 86, 165105(2012)), the following procedures should be taken:

      The code supports doing  cRPA calculations with either the projector-based Wannier functions (dmftproj, a module of in the TRIQS package)  or the maximally localized Wannier functions by using wien2wannier and wannier90. 
      To use the later, you must have first installed wannier90. The main options related to CRPA are followings:
       -w <wf_method>   # CPRA: control how to obtain Wannier projectors
                         1 -- use dmftproj
                         2 -- use MLWF obtained by wannier90+wien2wannier
    
       -win  <proj bot top>  # the basic input parameters to generate the WF projectors
                proj -- a string defining the projections in the form "atom1:orb1;atom2:orb2", where orb can be the following
                     "s/p/d/f"   -- d/f orbital in the cubic spherical harmonic
                     "ds/fs"     -- d/f orbitals in the sphericalharmonics representatioin
                     "dt2g/deg"  -- the t2g/eg manifold of the d orbitals
               bot top -- indicate the range of bands considered for WF.
                      for wf_method=1 the lower and upper energy bound (in eV) w.r.t. the Fermi energy for wf_method=2 band indices for bottom and top bands
                      if '-win ...' is not present, WF projectors are generated in the default way.
                         for wf_method==1, it is assumed that case.indmftpr is already prepared
                         for wf_method==2, gap2_w2w is run interactively
    
    Once you run gap2_init successfully, you need to edit gw.inp manually to set the transitions that are going to be excluded.
    
    For example, for NiO, to do a d-d model CRPA calculation,
    1) run
            gap2c_init -emax 6.0 -nkp 216 -d ../NiO-crpa-e6-k6c -w 2 -win Ni:d 8 12
        to generates projector-based Wannier functions by calling wien2wannier+wannier90.
    
    2) edit gw.inp to set the window for the transitions that are going to be excluded
    
      %eps_mask
          1  | 0.0    # iop_mask_eps |  wt_mask
          5  | 8      # noc_excl | iocc_excl
          5  | 8      # nunoc_excl | iunoc_excl
      %
    
    If you want to study the full frequency dependence of U, you need also change the block FreqGrid, 
    
    %FreqGrid                         # Frequency grid parameters
     1   |  2   |  4.0   |  1.E-8  |  0    # iopfreq | nomeg | omegmax | omegmin | nomeg_blk
    %
    
    The default as shown above is to only calculate U at zero frequency (1.E-8) and a large frequency (4 Hartree). You may want to modify the block as
    
    %FreqGrid                         # Frequency grid parameters
     1   |  100   |  2.0   |  1.E-8  |  0    # iopfreq | nomeg | omegmax | omegmin | nomeg_blk
    %
    
    3) run the gap2
    After you finish the calculations, you can get the CRPA results mainly from case.Umat, which reads like
     # bare Coulomb
    #Averaged U(full) U(diag) J (eV)=     24.8024   26.2366    0.8984
    # U_{m1,m2} (il1=1, il2=1)
       26.946936   24.632813   24.083271   25.001160   25.001160
       24.632813   26.946934   25.307122   24.389233   24.389233
       24.083271   25.307122   25.763042   23.877974   23.877974
       25.001160   24.389233   23.877974   25.763042   23.877974
       25.001160   24.389233   23.877974   23.877974   25.763042
    # J_{m1,m2} (il1=1, il2=1)
       26.946936    1.157061    1.131355    0.675925    0.675925
        1.157061   26.946934    0.524114    0.979545    0.979545
        1.131355    0.524114   25.763042    0.953544    0.953544
        0.675925    0.979545    0.953544   25.763042    0.953544
        0.675925    0.979545    0.953544    0.953544   25.763042
    
    # W at omega (eV)=    0.000000
    #Averaged U(full) U(diag) J (eV)=      6.1850    7.4628    0.8278
    # U_{m1,m2} (il1=1, il2=1)
        7.509596    5.529655    5.498480    6.226036    6.226305
        5.529655    7.516499    6.473899    5.743091    5.743442
        5.498480    6.473899    7.431092    5.738599    5.738903
        6.226036    5.743091    5.738599    7.428170    5.737279
        6.226305    5.743442    5.738903    5.737279    7.428650
    # J_{m1,m2} (il1=1, il2=1)
        7.509596    0.991821    1.037988    0.643204    0.643189
        0.991821    7.516499    0.511595    0.906458    0.906434
        1.037988    0.511595    7.431092    0.879218    0.879222
        0.643204    0.906458    0.879218    7.428170    0.879201
        0.643189    0.906434    0.879222    0.879201    7.428650
    
    # W at omega (eV)=  108.845584
    ......
    

    Return to ToC

    Input files

    After running gap2_init, gwdir contains the following input files needed for running a GW calcuation:

    • gw.inp -- the master input file (see below).
    • case.struct - the WIEN2k struct file (directly copied from w2kdir);
    • case.energy - WIEN2k energy file (directly copied from w2kdir);
    • case.in1 - WIEN2k in1 file (directly copied from w2kdir/case.in1[c], but may have different parameters as that used in SCF calculations.
    • case.vector[/up/dn] - WIEN2k vector file
    • case.vsp[/up/dn] - WIEN2k vsp
    • case.vxc - copied from case.r2v[/dn] file)
    • case.core - core energies and wavefunctions (case.scfc and case.corewf combined).

    The master input file is named as gw.inp, which is generated when running the script gap2_init. In most cases, one does not need to make any changes to gw.inp. A sample of gw.inp reads as followings.

    #Initialization options: gap2c_init -emax 1000.0 -nkp 0 -d ../bro-k122-nlo1 -newin1 1
    Task    = 'gw'
    %gw
      0 | 0                 # iop_sxc | iop_vxc 
    %
    
    %BZConv                 # BZ convolution options
     "tetra" | "imfreq"
    %
    
    %FreqGrid                         # Frequency grid parameters
     3  | 16  | 0.42  |  0.0  |  0     # iopfreq | nomeg | omegmax | omegmin | nomeg_blk
    %                                 # iopfreq= 1 (equally spaced), 2 (Gauss-Laguerre) or 3 (double Gauss-Legendre) 
    CaseName = "bro"             # case name, used as prefix for input/output files
    dftpkg  = 0                  # indicate which package is used for DFT calculations ( 0 -> wien2k ) 
    Restart = F
    Task    = "gw"               # Option for task
    SavDir ="./tmp"              # directory to store temporary data that might be reused 
                                 # once you are sure you are not going to reuse them, you should remove this directory 
                                 # since it may occupy large space especially with large N_k
    
    emingw = -2.0                # emingw and emaxgw ( unit --  Ry. ) are used to control the range of bands 
    emaxgw =  2.0                # for which GW correction are going to be calculated. Only states whose LDA energies 
                                 # falls between E_Fermi+emingw and E_Fermi+emaxgw are calculated 
    
    nspin = 1                    # 1 for spin-unpolarized and 2 for spin-polarized calculations 
    ComplexVector = F            # T for for sytems without inversal symmetry 
    UseScratch = 2               # 0 - if not using scratch space in which zzk are all kept in core 
                                 #     and minm is always calculated from scratch. 
                                 # 1 - use the scratch and different vectord files for different processes
                                 # 2 - use the scratch and different processes share the single vectord 
    SymVector = T                # whether use symmetrized output from wien2k, i.e. indicating whether KS energies 
                                 # and vectors in the irreducible (T) or full (F) BZ are written in case.energy/vector
    
    Minm_mblksiz = 48            # block size to split m-index for the calculations of minm 
    
    %SelfEnergy             # option for correlation self-energy 
      2  |  0   |   1       #  |    |   
    %                       # Number of poles ( previous maxexp + 1, valid range: 2.. nomeg/2 ) 
                            # iopes: 0/1/2/3 - without or with itereration > 
                            # iopsac:0/1     - Pade's approximation / multipole fitting 
    
    emaxpol = -1.0E10       # the upper bound for unoccupied states considered in polarizatioin 
    eminpol = -1.0E10       # the lower bound for deep core states considered in polarization 
    emaxsc  = -1.0E10       # the upper bound for unoccupie states considered in the self-energies
    barcevtol = 0.5         # tolerance used to reduce the bare Coulomb matrix eigenvectors as the basis set 
    MB_ludot = F            # whether considering the contribution of udot when setting up the mixed basis
    MB_emin = -1.E10        # only core states with energy above MB_emin are considered when setting up MB basis 
    MB_emax = 20.0          # LO orbitals with energy higher MB_emax are excluded when setting up MB basis 
    
    LOmax = 3               # the LOMAX used in WIEN2k(lapw1,nmr) 
    
    %MixBasis               # Mixed basis parameters
    0.7                     # Q_mb  
    2    |    1.E-4  |  0   # lmbmax  | wftol  | lblmax 
    %
    
    ºreCoul               # Bare Coulomb interaction
    2.0  |  1.E-15
    %
    
    nvel = 192.0            # the number of valence electrons
    
    %kmesh
     1 | 2 | 2               # the number of k-points along each axis
     0 | 0 | 0 | 2           # shift from k=0 (in integer coordidate) and their common division
    %
    

    Return to ToC

    Output files

    The temporary files generated during the calculation are stored in gwdir/tmp, and can be removed once all calculations have been done. In gwdir, the following new files will be created after the calculation.

    • case.outgw - the main output file
    • case.eqpeV_[GW/GW0] - G0W0 and GW0 quasi-particle energies in units of eV.
    • case.eqpeH_[GW/GW0] - G0W0 and GW0 quasi-particle energies in units of Hartree (mainly used for post-processing purpose.
    • case.sxc_nn - diagonal elements of the exchange and correlation self-energies (mainly used for post-processing purpose).
    • case.vxc_nn - diagonal elements of the DFT exchange-correlation potential
    • case.emac - the macroscopic dielectric function, which is essentially useless for a GW calculation since it is along the imaginary frequecy and it is usually not converged with respect to the number of k-points

    Return to ToC

    Post-processing

    The main results of a G0W0/GW0 calculation by running the GAP2 code are a set of quasi-particle energies on an equally spaced k-mesh. To calculate densities of states or plot the band structure diagram along high symmetry directions in the BZ, one usually needs QP energies at k-points that are not considered in GW calculations. Calculating QP energies at an arbitrary k is possible, but computationally expensive. Currently, the GAP2 code uses the Fourier interpolation approach (Pickett et al. Phys. Rev. B 38, 2721 (1988)) to obtain QP energies on k-points that are not included in the original k-mesh. For that purpose, the C-shell scripts gap2_analy and gap2_gwnvf can be used.

    The script gap2_gwnvf generates two files, namely case.energy_gw and case.vector_gw, which contain interpolated quasiparticle energies, using data from case.energy and case.vector in the current directory. These QP interpolated energies and vectors can then be used for further analysis or other calculations.

    The script gap2_analy first call "x lapw1" to calculate KS energies and vectors on a new k-mesh, and then use gap2_gwnvf to obtain QP corrected energies and vectors on the same k-mesh, from which either DOS ("x tetra") or band structure diagrams ("x spaghetti") can be obtained.

    Note that both gap2_analy and gap2_gwnvf should be run in w2kdir. More information about these two scripts can be obtained by running them with a command-line argument "-h".

    Return to ToC

    Consideration of spin-orbit coupling

    As mentioned previously, the SOC effects can be taken into account only at the perturbative level, which is essentially just using GW corrected Kohn-Sham vectors without SOC as the input to do a one-shot SOC calculation. This is realized by using the post-processing script gap2_analy with the option "-so0". But before running gap2_analy with "-so0", you first need to prepare necessary input files by calling "initso_lapw".

    Return to ToC

    Citing GAP2

    When publishing the work that has used the GAP2 code, please cite the following papers

    Hong Jiang and Peter Blaha, Phys. Rev. B, 93,115203(2016).

    Hong Jiang, Ricardo I. Gomez-Abal, Xinzheng Li, Christian Meisenbichler, Claudia Ambrosch-Draxl, and Matthias Scheffler,Computer Phys. Commun.,184,348 (2013).

    Return to ToC

    Obtaining GAP2

    The GAP2 code is available to the licensed WIEN2k users free of charge. The latest stable version can be downloaded here (last update: Aug 10, 2021). 

    Return to ToC

    Publications using FHI-gap/GAP2

    Publications related to external use of the code

    (Please send the information on your new publication when the GAP/GAP2 code has been used in your work).
    1. Ferenc Karsai, Paul Tiwald, Robert Laskowski, Fabien Tran, David Koller, Stefanie Gr?fe, Joachim Burgd?rfer, Ludger Wirtz, and Peter Blaha, F center in lithium fluoride revisited: Comparison of solid-state physics and quantum-chemistry approaches, Phys. Rev. B 89, 125429 (2014).

    Publications by the code developers and collaborators

    1. Min-Ye Zhang and Hong Jiang*,Electronic Band Structure of Cuprous and Silver Halides: a Numerically Accurate All-Electron GW StudyPhys. Rev. B. 100, 205123(2019) arXiv:1906.02472v1. (2019).
    2. Hong Jiang,Revisiting the GW Approach to d- and f-electron Oxides, Phys. Rev. B, 97, 245132(2018).
    3. P. Delange, S. Backes, A. van Roekeghem, L. Pourovskii, H. Jiang, S. Biermann, Novel approaches to spectral properties of correlated electron materials: From generalized Kohn-Sham theory to screened exchange dynamical mean field theory, J. Phys. Soc. Jp. 87, 041003(2018).
    4. S. K. Panda, H. Jiang, S. Biermann, Pressure dependence of dynamically screened Coulomb interactions in NiO: Effective Hubbard, Hund, intershell and intersite components, Phys. Rev. B 96, 045137(2017) arXiv:1612.07571.
    5. Feng Wu, Huihui Wang, Yu-Chen Shen, and Hong Jiang*,Electronic properties of ionic surfaces: a systematic theoretical investigation of alkali halides, J. Chem. Phys. 146, 014703 (2017).
    6. Hong Jiang and Peter Blaha, GW with linearized augmented planewaves extended by high-energy local orbitals, Phys. Rev. B, 93,115203(2016).
    7. Ambroise van Roekeghem,Thomas Ayral,Jan M. Tomczak, Michele Casula,Nan Xu,Hong Ding,Michel Ferrero,Olivier Parcollet,Hong Jiang,and Silke Biermann,Dynamical correlations and screened exchange on the experimental bench: spectral properties of the cobalt pnictide BaCo2As2, Phys. Rev. Lett. 113, 266403 (2014) arXiv:1408.3136.(2014).
    8. Wei An, Feng Wu, Hong Jiang, Guang-Shan Tian and Xin-Zheng Li*, Systematic investigation on topological properties of layered GaS and GaSe under strain, J. Chem. Phys. 084701 (2014).
    9. Philipp Hansmann, Loig Vaugier, Hong Jiang and Silke Biermann, What about U on surfaces? Extended Hubbard models for adatom systems from first principles, J. Phys.: Condens. Matter 25 094005 (2013).
    10. Hong Jiang, Ricardo I. Gomez-Abal, Xinzheng Li, Christian Meisenbichler, Claudia Ambrosch-Draxl, and Matthias Scheffler, FHI-gap: a GW code based on the All-electron augmented plane wave method, Computer Phys. Commun.,184, 348 (2013).
    11. Loig Vaugier, Hong Jiang and Silke Biermann, Hubbard U and Hund's exchange J in transition metal oxides: screening vs localization tends from constrained random phase approximation, Phys. Rev. B 86, 165105(2012).
    12. Hong Jiang, Patrick Rinke and Matthias Scheffler, Electronic properties of lanthanide oxides from the GW perspective, Phys. Rev. B, 86, 125115 (2012).
    13. Hong Jiang, Electronic Band Structures of Molybdenum and Tungsten Dichalcogenides by the GW Approach J. Phys. Chem. C 116,7664-7671 (2012).
    14. Xin-Zheng Li, Ricardo Gomez-Abal, Hong Jiang, Claudia Ambrosch-Draxl, and Matthias Scheffler, Impact of widely used approximations to the G0W0 method: an all-electron perspective, New J. Phys. 14, 023006 (2012).
    15. Huihui Wang, Feng Wu and Hong Jiang, Electronic band structures of ATaO3 (A=Li, Na and K) from first-principles many-body theory, J. Phys. Chem. C 115, 16180 (2011).
    16. Hong Jiang, Structural and electronic properties of ZrX2 and HfX2(X=S, and Se) from first principles calculations J. Chem. Phys. 134, 204705 (2011)
    17. Hong Jiang, Ricardo I. Gomez-Abal, Patrick Rinke, and Matthias Scheffler, First-principles modeling of localized d states with the GW@LDA+U approach, Phys. Rev. B 82, 045108 (2010).
    18. Hong Jiang, Ricardo I. Gomez-Abal, Patrick Rinke, and Matthias Scheffler, Quasiparticle Electronic Band Structure of Zirconia and Hafnia Polymorphs, Phys. Rev. B 81, 085119(2010).
    19. Hong Jiang, The GW Method: Basic Principles, Latest Developments and Its Applications for d-and f-Electron Systems (Invited Review), Acta Phys.-Chim. Sin.,26, 1017(2010).
    20. Hong Jiang, Ricardo I. Gomez-Abal, Patrick Rinke and Matthias Scheffler, Localized and itinerant states in lanthanide oxides united by GW@LDA+U , Phys. Rev. Lett. 102, 126403 (2009). pdf
    21. R. Gómez-Abal, X. Li, M. Scheffler, and C. Ambrosch-Draxl,Influence of the core-valence interaction and of the pseudopotential approximation on the electron self-energy in semiconductors. Phys. Rev. Lett. 101, 106404 (2008).

    Return to ToC