- Introduction
- Limitations
- Installation
- How to use the code for GW calculations?
- How to do CRPA calculations to estimate the Hubbard U
- Input files
- Postprocessing
- Consideration of spin-orbit coupling (SOC)
- Citing GAP2
- Obtaining GAP2
- Publications using FHI-gap/GAP2

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:

- MPI-based three-level parallelization.
- GW using LAPW enhanced with high-energy local orbitals (HLOs) (see, H. Jiang and P. Blaha, Phys. Rev. B, 93,115203(2016))
- Perturbative consideration of spin-orbit coupling (SOC)
- Constrained random phase approximation for first-principles determination of the Hubbard U using the projected Wannier functions (with the interface with dmftproj) and maximally localized wannier functions (MLWF) with the interface to wannier90 by wien2wannier.

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 *G*_{0}*W*_{0} 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 G_{0}W_{0} for d and f-electron systems, we have implemented *G*_{0}*W*_{0} 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.

- 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.
- 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.
- The treatment of metallic systems is still in the experimental state, and has not been well tested.

- 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.
- 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) - 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".

- Run a full SCF calculation using WIEN2k in the directory denoated as
henceforth**w2kdir** - 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*" - Modify
*gwdir/gw.inp*if necessary, and then run "*gap*2.x" or "*gap2-mpi*.x" (using mpirun) in the*gwdir*. - 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*. - Some useful information is given in the main output file, case.outgw. Further analysis can be done by using the C-shell script
(see post-processing).*gap2_analy* - 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. - 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 readsWFFIL 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/nbandand it has to be modified to

WFFIL EF=.353692248950 (WFFIL, WFPRI, ENFIL, SUPWF) 9.00 10

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**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### 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 ......### 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

, which is generated when running the script*gw.inp*. In most cases, one does not need to make any changes to*gap2_init**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 % 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

### Post-processing

The main results of a

*G*_{0}*W*_{0}/*G**W*_{0}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 scriptsand**gap2_analy**can be used.**gap2_gwnvf**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

and**gap2_analy**should be run in**gap2_gwnvf***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 papersHong 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).### 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).

### 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).- 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

- Min-Ye Zhang and
**Hong Jiang***,*Electronic Band Structure of Cuprous and Silver Halides: a Numerically Accurate All-Electron GW Study*, Phys. Rev. B.**100**, 205123(2019) arXiv:1906.02472v1. (2019). - Hong Jiang,
*Revisiting the GW Approach to d- and f-electron Oxides*, Phys. Rev. B,**97**, 245132(2018). - 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).
- 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. - 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).
- Hong Jiang and Peter Blaha,
*GW with linearized augmented planewaves extended by high-energy local orbitals*, Phys. Rev. B, 93,115203(2016). - 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). - 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). - 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). - 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). - 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). - Hong Jiang, Patrick Rinke and Matthias Scheffler,
*Electronic properties of lanthanide oxides from the GW perspective*, Phys. Rev. B, 86, 125115 (2012). - Hong Jiang,
*Electronic Band Structures of Molybdenum and Tungsten Dichalcogenides by the GW Approach**J. Phys. Chem. C*116,7664-7671 (2012). - 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). - 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). - Hong Jiang,
*Structural and electronic properties of ZrX2 and HfX2(X=S, and Se) from first principles calculations**J. Chem. Phys.*134, 204705 (2011) - 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). - 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). - 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). - 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 - 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).

- a) Make some changes in WIEN2k as explained as followings: