3.27 Periodic GW in FHI-aims

A periodic version of GW (more precisely the one-shot G0W0) has been implemented in FHI-aims. The current implementation is based on the k-space formalism, similar to the k-space implementation of periodic HF (can be invoked by setting use_hf_kspace to be true), and only works for insulating systems. So far, tests show that it is rather stable, but if you encounter any problem, consult the developers.

Complete descriptions of the key algorithm and implementation details can be found in Ref. [260]. Here we only present the central equation for computing the periodic G0W0 self-energy,

Σn,σG0W0(𝐤,iω) = 𝑑𝐫𝑑𝐫ψn,σ𝐤(𝐫)ΣσG0W0(𝐫,𝐫,iω)ψn,σ𝐤(𝐫)
= 12πm,𝐪μ,ν𝑑ωCn,m,σμ(𝐤,𝐤𝐪)W0,μν(𝐪,iω)Cm,n,σν(𝐤𝐪,𝐤)iωiω+μϵm,σ𝐤𝐪.

where the expansion coefficients Cn,m,σμ(𝐤,𝐤𝐪) are determined using the “RI-LVL” (or localized RI) approximation, similar to the periodic HF and hybrid function case discssued in Sec. 3.26.

At present, only the one-shot G0W0 scheme is implemented for periodic systems. The periodic G0W0 can be run in much the same way as the cluster (finite system) case. That is, one needs to specify a starting point by setting xc to a desired funcitonal. At the moment, LDA, GGAs and global hybrid functions can be employed as the starting points.

The periodic G0W0 is invoked by setting qpe_calc to gw_expt and setting k_grid to appropriate values. Similar to the cluster case, the periodic G0W0 self-energy is first computed on imaginary frequencies and then analytically continuted to the real frequency axis. As such, one needs to explicitly specify the type of analytical continuation procedure by setting the anacon_type to 0 (two-pole fitting) or 1 (Padé approximation.) The number of frquency points and the number of analytical continuation parameters can be specified by setting frequency_points and n_anacon_par to appropriate values. If not, the same default values as the clusters will be used for these two parameters. Further details can be found in Sec. 3.25.

As mentioned above, periodic GW calculations rely on the RI-LVL scheme. This implies one typically needs a larger set of auxiliary basis functions to achieve a satisfactory level of numerical accuracy. In practice, one usually needs to add so-called "for_aux" basis functions in species_default, that are used to generate additional auxliary basis functions. Tests show that the following “for_aux" basis setting usually works rather well, and is recommended as a first choice.

  for_aux     hydro 4 f  0.0
  for_aux     hydro 5 g  0.0

However, this does not mean the above setting is optimal, or always sufficient. If in doubt, one is encouraged to vary the effective charge parameter Z (the value “0.0" above) to check how the obtained results (e.g., the band gap) change.

The outputs from periodic G0W0 calcualtions are the G0W0 quasiparticle energies within a window of energy levels at pre-set 𝐤 points. The specification of the 𝐤 points can be done in two ways. In the first way, one can directly compute the quasiparticle energy band structures along a set 𝐤-point paths. The way to set up the 𝐤-point paths is exactly the same as in KS-DFT calculations, through the keyword output band. In the second way, one can choose to print out the quasiparticle energy level on a sef of 𝐤 points belonging to the uniform 𝐤 grid, that is used in the preceding SCF calculations. This is set through the keyword output gw_regular_kgrid. At each 𝐤 point, the set of energy levels for which one asks for quasiparticle energy calculations can be specified through the keyword state_lower_limit and state_upper_limit, in exactly the same way as in the cluster case. If state_lower_limit and state_upper_limit are not explicitly specified, G0W0 quasiparticle calculations will be done for the states close to the Fermi level. You can also use state_limits_ewin to choose the states to calculate by setting an energy window.

Tags for general section of control.in:

 

Tag: qpe_calc(control.in)

Usage: qpe_calc gw_expt
Purpose: If set, and if a finite 𝐤 grid is also set via the keyword k_grid, periodic G0W0 calculation will be started.

Here we choose to use a different keyword gw_expt instead of gw as in the cluster case, in order to emphasize that this is still an experimental version. In future we may change the setting gw_expt back to gw.

 

Tag: state_limits_ewin(control.in)

Usage: state_limits_ewin ref e_lower e_upper
Purpose: This keyword offers one way to control the states to calculate quasi-particle energies besides the keywords state_lower_limit and state_upper_limit. When enabled, FHI-aims will only calculate the quasi-particle energies for the bands with at least one state located within the energy window determined by ref, e_lower and e_upper.

ref: edge, fermi or abs. It sets the energy reference of the later two energy bound options.

  • edge: use the VBM (HOMO) and CBM (LUMO) for the reference of lower and upper bound, respectively.

  • fermi: use the internal Fermi energy for the reference of both bounds.

  • abs: no reference is used, i.e. the energy bounds are absolute energy values.

e_lower and e_upper: The lower and upper energy bound with respect to the reference (see the ref option above). They are given in units of eV.

Example: edge -10.0 10.0 will set the energy window between 10.0 eV below the VBM and 10.0 eV above the CBM.

By default, energy window is set to 20 eV below VBM and 10 eV above CBM, i. e. edge -20.0 10.0. This default is overwritten when either state_lower_limit and state_upper_limit is explicitly set. When both state_limits_ewin and state_lower_limit/state_upper_limit are specified, state_limits_ewin will be respected.

After the energy window is set, quasi-particle energies are calculated for bands with at least one state located within this energy window.

Note that, at the moment, this keyword only works in periodic GW, i.e. when qpe_calc is set to gw_expt.

 

output sub-tag: gw_regular_kgrid(control.in)

Usage: outputgw_regular_kgrid
Purpose: If set, the quasiparticle energy levels at the regular 𝐤 grid will be printed out. However, one does not automatically print out the information on all 𝐤 points (which would be often too much), but rather for a number specified through output k_eigenvalue number. That is, the 𝐤 points where quasiparticle energy levels are printed out are precicely the same as those on which the preceding KS/HF eigenvalues are printed.

 

output sub-tag: gw_nscf_exx_only_band(control.in)

Usage: outputgw_nscf_exx_only_band
Purpose: If set, the non-self-consistent exact-exchange (EXX) band structure will be printed in Non-SCF_EXX_only_bandXYYY.out on the same band path and with the same format as the regular band output bandXYYY.out. There is no computational overhead by switching this option on as the EXX contribution has already been calculated to solve the quasi-particle equation.

 

Tag: periodic_gw_optimize(control.in)

Usage: periodic_gw_optimizewo_inverse avoid_hang freqbatch:16 blacs_block_size:128 lvlkq_block_size:16 pol_s_block_size:512
Purpose: Keywords to change the technical details of the periodic GW calculations. wo_inverse is to switch off the Cholesky matrix inverse algorithm (the current default) and use the spectrum decomposition algorithm. avoid_hand is to use a more stable version of message passing. freqbatch:xxx is to use an extra layer of parallization over frequency integration where xxx is the number of sub mpi groups. blacs_block_size:xxx is the BLACS block size of the parallel linear algebra. lvlkq_block_size:xxx and pol_s_block_size:xxx are used to control the batch sizes for g_times_w and polarisability calculations. The default is hard coded to make sure of using  800 MB RAM per MPI.

 

Tag: periodic_gw_optimize_single_precision(control.in)

Usage: periodic_gw_optimize_single_precisiong_times_w polarisability self_energy
Purpose: Keywords to switch on the single precision calculations in certain region of the code.

 

Tag: periodic_gw_optimize_init(control.in)

Usage: periodic_gw_optimize_initinverse scale fft
Purpose: Keywords to switch on the optimized version of the init steps for periodic GW. fft requires the code to be linked with an optimized FFT library.

 

Tag: periodic_gw_optimize_lvl_init(control.in)

Usage: periodic_gw_optimize_lvl_init boolean
Purpose: Keywords to switch on the optimized version of the initialization of periodic LVL triple coefficients.
Note: This switch typically reduce the initialization time by more than half, and works well for systems with large number of atoms. It is introduced to accelerate the periodic GW calculation, but should also help for other post-SCF calculations involving RI-LVL. One caveat is that peak memory usage at this stage can increase in calculation using a number of k-points when it is switchted on. If out-of-memory break is observed during initialization, we suggest the user turn off this switch. The users are welcome to report the case to the developers to improve this part of code. boolean: .true. or .false. Default is .true.

 

Tag: periodic_gw_optimize_use_gpu(control.in)

Usage: periodic_gw_optimize_use_gpug_times_w polarisability inverse pxgemm
Purpose: Keywords to switch on the GPU accelerated on certain steps in periodic GW calculations.

 

Tag: periodic_gw_head_correction(control.in)

Usage: periodic_gw_head_correction switch q0x q0y q0z
Purpose: Control whether the head of polarization matrix in Coulomb-diagnoal basis at the Γ point should be corrected by the macroscopic dielectric function. Users usually should not change it unless knowing exactly what the purpose is.

boolean: boolean, default .true. to switch on this correction. This is shown to improve the k-point convergence for quasi-particle energies in the current implementation.

q0x , q0y , q0z: float numbers, altogether defining the direction in which the macroscopic dielectric function is calculated from dielectric tensor. By default, all of them are set 1, i.e. we are taking the (1,1,1) direction, which is equivalent to taking the isotropic average in systems with a diagonal dielectric tensor.

 

Tag: periodic_gw_discard_wc_gamma(control.in)

Usage: periodic_gw_discard_wc_gamma boolean
Purpose: Discard the correlation self-energy contribution from Wc(𝐪=0). This will lead to slow convergence with respect k-grids of the quasi-particle band structures, but avoid the use of different Coulomb operators when building dielectric matrix and self-energy from W at the Gamma point. This option was added for a detailed analysis of the effect of RI functions and for experimental use. Default is .false., i.e. Gamma point contribution to self-energy is always included. In general, users should not change this, otherwise band structures well below convergence would be expected.
boolean: .true. or .false.
Note: currently a hack is used for this purpose: zero out Wc(𝐪=0,iω) when adding its contribution to Σc. So this option has negligible impact on the efficiency.

 

Tag: periodic_gw_use_average_inverse_dm_gamma(control.in)

Usage: periodic_gw_use_average_inverse_dm_gamma boolean
Purpose: compute the inverse dielectric matrix at the Γ point as an average over the Brillouin zone of the Born-von-Karman cell.
Option boolean: .true. or .false.; default .false.
By default, the dielectric matrix at the Γ point, i.e. at the long-wavelength limit, is actually computed using response function 𝝌(𝐪=0) and singularity-lifted bare Coulomb matrix 𝐕(𝐪=0). The resulted matrix is then corrected by the macroscopic dielectric function evaluated at certain 𝐪 direction based on 𝐤𝐩 theory. When this keyword is switched on, the program will account for correction from all directions through a spherical average.
It has been tested in cubic systems that this flag significantly speed up the convergence of quasiparticle band gap with respect to k-grid and extra for_aux functions (see the beginning of this section for discussion).
Note: this keyword is experimental. It does not work with the optimize keywords in periodic_gw_optimize except for inverse.

The inverse dielectric matrix at the Γ point under Coulomb eigenvectors representation is replaced by an average around the Γ point (frequency argument is left out for clarity):

𝜺v1¯= 1Vγγd𝐪𝜺v1(𝐪) (3.97)
1Vγγd𝐪^𝜺v1(𝐪^)0qmax(𝐪^)dqq2
= 1Vγγd𝐪^𝜺v1(𝐪^)qmax3(𝐪^)

where the subscript v indicates the Coulomb eigenvector representation and γ stands for a small region enclosing the Γ point, currently chosen as the Brillouin zone of the Born-von-Karman cell with volume Vγ. qmax(𝐪^) is the distance to the Γ point from the border of γ along solid angle 𝐪^. 𝜺v1(𝐪^) is calculated by block matrix inversion of

𝜺v(𝐪^)=[H(𝐪^)𝐰(𝐪^)𝐰(𝐪^)𝐁] (3.98)

where the head term H(𝐪^) is a scalar, the wing term 𝐰 a column vector of size Nv1 and the body term 𝐁 a (Nv1)×(Nv1) matrix, where Nv is the number of Coulomb eigenvectors (see prodbas_threshold).

The expressions for the head and wing terms are

H(𝐪^)= α=x,y,zβ=x,y,zq^αHαβq^β (3.99a)
wλ(𝐪^)= α=x,y,zq^αwλα (3.99b)

where

Hαβ(iω)=δαβ4πN𝐤Vcell 𝐤σ{nf(ϵn𝐤σ)pα,nn,σ𝐤pβ,nn,σ𝐤ω2 (3.100)
+m<n2(fm𝐤σfn𝐤σ)pα,mn,σ𝐤pβ,mn,σ𝐤[(ϵm𝐤σϵn𝐤σ)2+ω2](ϵm𝐤σϵn𝐤σ)}

and

wλα(iω)= vλμXμλwμα(iω) (3.101)
= 4πVcell2vλN𝐤μXμλ𝐤σmnCmnμ(𝐤,𝐤)pα,nm,σ𝐤[ϵn𝐤σϵm𝐤σ]2+ω2fn𝐤σ(1fm𝐤σ)

where vλ and 𝐗λ is the Coulomb eigenvalue and eigenvector, respectively. pα,mn,σ𝐤 is the momentum (or velocity) matrix in Kohn-Sham basis, briefly discussed in Sec. 3.42.

The head tensor Hαβ in Eq. (3.100) and wing vector wμα in Eq. (3.101) will be written to head_tensor.dat and wing_vector_vertical.dat, respectively. Both files are unformatted and can be processed using the script read_head_wing.py under the utilities folder.

The idea basically follows Freysoldt et al,[95] namely Equations (27), (33), and Appendix B therein. Differences in our implementaion include:

  • In [95], the basis corresponding to the head element is uniquely defined in the plane-wave basis framework, i.e. the G=0 plane wave. Such plane-wave basis is not directly available in FHI-aims. Instead, we approximate by the eigenvector with the largest eigenvalue of the singularity-lifted Coulomb matrix in auxiliary basis.

  • The conventional reciprocal-space algorithm for periodic GW is adopted in our implementation, thus we do not transform to W(r,r) in context of space-time approach in [95].

  • The q-average is taken during the calculation of screened Coulomb calculation, instead of the self-energy integral. This is reasonable in the current FHI-aims periodic GW implementation using cut Coulomb to compute screened Coulomb from the dielectric matrix, which is regular and smoothly varying near the Γ point.

  • Currently we do not rely on expanding the geometrical terms using spherical harmonics. To ensure numerical precision, the Lebedev grid with highest order (131, 5810 grid points) is used. This can be improved in the future.

 

Tag: periodic_gw_output_dielectric_matrix(control.in)

Usage: periodic_gw_output_dielectric_matrix ikpoints ifreqs
Purpose: Write dielectric matrix ϵ(ki,ωj) at specific k-point and frequency-point to file.

ikpoints: a list of integers (ikpoint, starting from 1) separated by comma. This option specifies indices of k-points at which ϵ should be dumped. Note that the indices are counted in the full BZ. When symmetry is used, ϵ is only calculated within the irreducible zone, thus the matrices at k-points outside the irreducible zone will not be written.
ifreqs: a list of integers (ifreq, starting from 1) separated by comma. This option specifies indices of frequency points at which ϵ should be dumped.

The dielectric matrices will be written in files named after dielectric_matrix_kpt_{ikpoint}_ifreq_{ifreq}.csc. The matrices are stoed in ELSI CSC format, and can be read using the scripts under the utilities/elsi_matrix folder.

Example: periodic_gw_output_dielectric_matrix 1,2 1,2,3 asks the code to write at most 6 files.

Note: At the time of writing, this keyword does not work with restart_periodic_gw and optimize keywords in periodic_gw_optimize except for inverse.

 

Tag: periodic_gw_force_finite_size_correction(control.in)

Usage: periodic_gw_force_finite_size_correction bool
Purpose: flag to enforce applying finite-size correction to periodic GW calculations. By default, finite-size correction (see keywords periodic_gw_head_correction and periodic_gw_wing_correction) is switched off for low-dimensional systems, identified by k-grid setup. When periodic_gw_force_finite_size_correction is set to .true., finite-size correction will be applied anyway.
bool: .true. or .false., default to .false.

 

Tag: periodic_gw_self_energy_restart(control.in)

Usage: periodic_gw_self_energy_restart key
Purpose: Control how to restart the periodic GW calculation by reading/writing the correlation self-energy on imaginary axis from files.

Note that this keyword is different from restart_periodic_gw which is used to restart from an interupted calculation.

This keyword is used only when you want to redo the analytic continuation and solving quasi-particle equation with different analytic continuation parameters. This requires the same system, including the number of spin, basis, QP states, k-points and frequencies points.

At present the file names to read/write are hard coded:

  • For self-energy on kgrid, the relavant file is self_energy_grid.dat.

  • For self-energy on band paths, the relavant file is self_energy_band.dat.

key: read, write, or read_and_write. By default this keyword is switched off.

When set to read, the program will look for the above files (depending on whether you want to calculate QP energies on k-grid or band-paths). If they are found, the evaluation of correlation self-energy on imaginary axis will be skiped, and the program will proceed to analytic continuation and calculation QP energies. If they are not found, the program will raise error and stop.

When set to write, the program will perform the normal GW calculation, and write Σc(iω) to above files.

When set to read_and_write, the program will first check the above files as using read. If they are found, the program will skip Σc(iω) part, similar to read. Otherwise, the program will turn off the reading and do the normal GW calculation.

By default, the output files are read/written in binary (unformatted) format by direct access. Plain-text format is available only for reading by adding a third argument plain, but is not recommend due to the numerical sensitivity of AC on input.

 

Tag: periodic_gw_skip_completed_bands(control.in)

Usage: periodic_gw_skip_completed_bands bool
Purpose: flag to check if GW band output files exist at persent working directory, and skip the re-evaluation of corresponding bands when they are found.
bool: .true. or .false., default to .false.

When set to .true., the program will check if GW_band1<iband>.out exists. If it does, it will not re-compute the iband-th band and proceed to the next. If spin polarization is switched on, GW_band2<iband>.out will be checked as well, and the evaluation will be skipped only when both files are found. We suggest to use this keyword together with periodic_gw_self_energy_restart and/or restart_periodic_gw.

 

Tag: restart_periodic_gw(control.in)

Usage: restart_periodic_gw <dirname> binary
Purpose: keyword to enable restart of periodic GW calculation from a previous abrupted job.
Options:
<dirname>: the name of directory to save/load the checkpoint files.
binary: (optional) flag to control whether the checkpoint files are unformatted binary.

When swichted on, FHI-aims will search for directory <dirname> and checkpoint files for periodic GW calculations therein. If they are found, the calculation will load data from checkpoint files and continue. Otherwise the calculation will create the directory and write data files under it. For example, restart_periodic_gw restart_pgw will save/load checkpoint files under the directory "restart_pgw". By default, the checkpoint files are in human-readable format. If the optional flag binary is also appended, the checkpoint files will be written as unformatted binary to reduce disk space usage.

A few notes:

  • Each process will write/load its own checkpoint file. The restart functionality implemented in restart_periodic_gw requires the same number of MPI tasks used in the previous job and the restart run to ensure the saved data are correctly loaded. Otherwise, the program will stop and warn the user.

  • Only the band data are saved in the checkpoint files in current implementation. Therefore, restart_periodic_gw does not work with output gw_regular_kgrid

The following section summarizes the different optimizations available for the periodic GW calculations.

Single Precision Optimizations (periodic_gw_optimize_single_precision)
Keyword Description
periodic_gw_optimize_single_precision g_times_w Use single precision for g_times_w
periodic_gw_optimize_single_precision self_energy Use single precision for self energy
periodic_gw_optimize_single_precision polarisability Use single precision for polarisability
Initialization(periodic_gw_optimize_init)
Keyword Description
periodic_gw_optimize_init inverse Use Cholesky decomposition for matrix inversion
periodic_gw_optimize_init fft Use optimized FFT library for Spherical Bessel Transformation. Requires FFT library.
periodic_gw_optimize_init scale Use optimized logsbt_scale subroutine
periodic_gw_optimize_lvl_init <.true./.false.> Use optimized LVL triple coefficient initialization
K-grid Symmetry(periodic_gw_optimize_kgrid_symmetry)
Keyword Description
periodic_gw_optimize_kgrid_symmetry all Use all symmetry for reduced kgrid
periodic_gw_optimize_kgrid_symmetry inverse Use only inverse symmetry for reduced kgrid
periodic_gw_optimize_kgrid_symmetry none no symmetry used
CPU Optimizations(periodic_gw_optimize)
Keyword Description
periodic_gw_optimize inverse Use Cholesky decomp for dielectric matrix inversion. Default uses eigen decomposition
periodic_gw_optimize avoid_hang Broadcast one by one to avoid hang behaviors for lvl_kq
periodic_gw_optimize lvlkq_block_size:<value> Use lvlkq block size of <value> for self-energy calculation
periodic_gw_optimize pol_s_block_size:<value> Use polarisability s_block size of <value> for polarisability_kpoints calculation
periodic_gw_optimize blacs_block_size:<value> Use BLACS block size of <value> for pzgemm
periodic_gw_optimize freqbatch:<value> Use freq batch with the number of batch of <value>
GPU Optimizations(periodic_gw_optimize_use_gpu)
Keyword Description
periodic_gw_optimize_use_gpu g_times_w Use GPU for g_times_w calculation
periodic_gw_optimize_use_gpu polarisability Use GPU for polarisability calculation
periodic_gw_optimize_use_gpu inverse Use GPU for parallel inverse (requires slate library)
periodic_gw_optimize_use_gpu pxgemm Use GPU for parallel matrix multiplication (requires cosma library)