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.
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
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 ......
After running gap2_init, gwdir contains the following input files needed for running a GW calcuation:
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
%
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.
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
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 ToCWhen 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).
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).