This guide documents the structure of the GWinput file used by the GW package. The GW suite does not use the ctrl file, but reads tags from GWinput. You can, however, create a template GWinput from the ctrl file, as explained for example, in this tutorial. This is particularly convenient because GWinput is not easy to read, or to construct by hand.
Note: in future GWinput will be subsumed into the ctrl file.
This documentation denotes Ref. I as a reference to PRB76 165106 (2007).
Table of Contents
- How input is organized
- The general category
- The PRODUCT_BASIS category
- The QPNT category
- The PBASMAX category
- The QforEPS and QforEPSL categories
How input is organized
File GWinput is the main put file for the GW package. Its structure has the tag-value style roughly similar to the ctrl file, but it is much less sophisticated. All but the first category is delimited html-style, e.g.
<PRODUCT_BASIS> ... </PRODUCT_BASIS>
Categories must be ordered as follows:
(general, no name or delimiter) <PRODUCT_BASIS> ... </PRODUCT_BASIS> <QPNT> ... </QPNT> <PBASMAX> ... </PBASMAX> <QforEPS> ... </QforEPS>
Some lines within a category are comments; some require a fixed order though the first (general) category do not. The user is strongly advised to make a GWinput template using
lmfgwd --job=-1, and edit it.
<PBASMAX> is optional; <QforEPS> is needed only to compute macroscopic the dielectric function or magnetic susceptibility at specified k points.
The general category
The general category has the structure
Some tags are optional.
data may be a string, a real number, an integer or a logical, and sometimes a vector of numbers. Lines beginning with, or text following, a ! are comments, and ignored by the reader. Below is a sample general category. The first column was added to label lines; they are described following the box.
!Verbose 41 ! 0-->default; 100-->debug GWversion 12 ! 0-->Sep12 compatible. !Q0P_Choice 1 ! 1(default):Near q->0 limit 2:1/q^2 average in Gamma region !CoreOrth off ! off --> Do not orthogonalize core to valence (default) !EXonly .15 ! for exchange-only calculations !EIBZmode off ! turn off to suppress use of symmetrization in calculating polarization !TimeReversal off ! when time-reversal symmetry is broken KeepEigen on ! keep eigenfunctions in memory KeepPPOVL on ! keep PPOVL in memory !Chi_RegQbz on ! Offset Gamma point mesh for chi. on => no offset BZmesh 1 ! Offset Gamma point mesh for Sigma=iGW WgtQ0P 0.01 ! Weight used when BZmesh is 2 NormChk 0 ! 1,2 writes norm check files (diagonal or full) n1n2n3 4 4 4 ! for GW BZ mesh QpGcut_psi 2.7 ! |q+G| cutoff for eigenfunction QpGcut_cou 2.2 ! |q+G| cutoff for coulomb int. unit_2pioa off ! off --> units of 2 preceding Gcut are a.u.; on--> units are 2*pi/alat alpha_OffG 1 ! (a.u.) parameter (auxiliary function, offset-Gamma method) !nband_chi0 999 ! nband cutoff for chi0 (Optional) !emax_chi0 999. ! emax cutoff for chi0, Ry (Optional) !nband_sigm 999 ! nband cutoff for Sigma (Optional) emax_sigm 2.5 ! Energy cutoff for Sigma, Ry (Optional) dw 0.01 ! mesh spacing along Real axis (Hartree) omg_c 0.1 ! mesh spacing increases linearly with w: becomes 2*dw at omg_c iSigMode 3 ! QSGW mode switch (QSGW only) niw 6 ! # freq. on Im axis; used for integration to make Sigma_c delta -1e-4 ! delta-function broadening for calc. x0, a.u.. delta<0 => tetrahedron deltaw 0.02 ! width in finite diff for sigma energy derivative, a.u. esmr 3e-3 ! Broadening in the poles of G(LDA) (hsfp0) GaussSmear on ! on => broadening of poles in G(LDA) by Gaussian !mixbeta .25 ! mixing of input, output sigma for self-consistency !mixpriorit 2 UseBSEW 0 ! include BSE (ladder diagram) contributions to W NumCondStates 8 ! (BSE only) number of conduction states to include in 2-particle hamiltonian NumValStates 3 ! (BSE only) number of valence states to include in 2-particle hamiltonian MinOmega_BSE 0.0! (BSE only) minimum frequency for BSE diectric function for optics MaxOmega_BSE 2.0! (BSE only) maximum frequency for BSE diectric function for optics Qvec_bse 1 1 1 ! (BSE optics only) q-vector for single-shot BSE dielectric function ImagEner1 -0.08 ! (BSE only, Ry, default=0.02) broadening at the minimum frequency for optics ImagEner2 0.107 ! (BSE only, Ry) broadening at the maximum frequency for optics; it varies linearly between minimum and maximum frequency x0_small_BSE off! (BSE only) x0_small_BSE on is a special purpose mode; See documentation below NumValBSE_no 9 ! (BSE only) Used when x0_small_BSE is on. See documentation below numOmega_BSE 999! (BSE only) NLF_corr_BSE on ! (BSE only)
Tags by referenced by line:
1 Specifying the number of divisions for a uniform mesh in Brillouin zone. Requires 3 integers: it is the analog of BZ_NKABC in the ctrl file.
2 Controls changes in functionality as the GW code has evolved. GWversion consists of a compound of digits.
- 1s digit
- include local field correction when calculating W(q=0); same as LFC@Gamma, Sep12 version of aux function
- Same as 1, but Aug 14 auxiliary function
- 10s digit
- number of intermediate states are computed relative to , not (this is a bug fix)
- alternative analytic auxiliary functions to regularize singularity at q-> (q0irre.F)
- 100s digit
- alternative estimate for average SE high lying states (hqpe.sc.m.F)
- rotate sigm in basis of LDA eigenfunctions (hqpe.sc.m.F)
- 1000s digit
- Use high accuracy but slow Bessel functions (not usually necessary)
GWversion=12 is the recommended choice. Use GWversion 0 for best compatibility with the Sep 2012 version of code.
3 Controls how the offset gamma points are determined. Q0P_Choice should be one of:
- 6 points along the reciprocal lattice vectors ( for each vector). This the default.
- 6 points along Cartesian axes (± for each vector)
- 4 points in Qlat(1) and Qlat(2) (± for each vector)
- 2 points in Qlat(3) (± for each vector)
These points may be reduced by point group symmetry. Offset gamma points are generated by qg4gw and stored in file Q0P.
For slabs, Q0Pchoice 2 may be better; but take care has this has not been fully analysed. Regardless, it is problematic to use unbalanced k points for anisotropic unit cells. see Computer Physics Comm. 176, 1 (2007).
One-dimensional chains seem to be well represented by Q0P_Choice 3 together with BZmesh 2.
4 If on cores are orthogonalized to valence partial waves and . Enforcing orthogonalization condition is needed to ensure the correct behavior in ε as q→0. However, orthogonalization may excessively deform shallow core orbitals, so usually it is safer not to use this switch. It is useful, however, to check how orthogonalization affects results. If it does, you can alternatively include these shallow cores in the secular matrix using local orbitals, which is very reliable and has no problems with orthogonalization.
5 Use EXonly for exchange-only calculations. The value following EXonly tells the code how much Hartree Fock character to admix with LDA. Thus use EXonly 1 for Hartree Fock.
6 By default, symmetry operations are used in calculating the polarizability. We follow a procedure similar to that of Friedrich et al as described in Phys. Rev. B81, 125102. Use EIBZmode off to suppress use of symmetrization.
7 Usually the polarizability assumes time-reversal symmetry Use TimeReversal off if it is broken.
8,9 To reduce disk access, eignvectors and the overlap matrix of the product basis are read and retained in memory by default. You can cause the GW code to read them from disk every time one is needed. This conserves memory but increases the (slow) disk access.
Use KeepEigen off to always read eigenfunctions from disk;
Use KeepPPOVL off to always read PPOVL from disk.
10 Use this switch to stagger the k mesh in calculating the dielectric function. It should be on (the default) or off. off offsets the mesh to avoid the point, analogous to what BZ_BZJOB does for lmf and related codes; see ctrl file input (but note that BZJOB can accept up to three switches, Chi_RegQbz only one).
At times the static limit of macroscopic dielectric function εM(ω→0) can converge very slowly with the number of k divisions nk, so that getting the converged result can be painful. It may happen (and often does) that εM(ω→0) converges with nk from below with the staggered mesh and from above with the usual one. In such a case the two meshes bound εM(0) for finite nk, and should converge in the 1/nk→0 limit. By plotting εM(0) as a function of 1/nk for both meshes at several values of nk, you can extrapolate each of them to 1/nk→0. If the two meshes extrapolate to the same point, it can be taken as a reliable value for εM .
An instance of this is NH4PbI3. The Figure shows εM(0) plotted for the two meshes. Even while εM(0) is not well converged for either mesh (especially the regular one), extrapolating to 1/nk=0 gives much confidence in the converged result.
See also the Brillouin integration page.
11 affects the k mesh in calculating the self-energy, analogous to BZ_BZJOB in the ctrl file. It should be 1 (the default) or 2: 2 offsets the mesh to avoid the point. 1 is recommended.
12 The total weight for the offset gamma method when BZmesh 2 (and is effective only for that case). More precisely it is fraction of the weight given to a regular mesh points (n1n2n3)−1. In principle, the final result should not depend on WgtQ0p. It is usually OK to take 0.01 (default), but there is a little dependence and you can take smaller values. When WgtQ0p is very small, e.g. 10−5, the QP levels become sensitive to the errors in normalization of the eigenfunction. You can improve the normalization with a larger QpGcut_psi. Normalization is printed out in normchk.dia.
13 Set to 1 or 2 causes the code to write norm check files (diagonal or full).
14 Number of k divisions along each reciprocal lattice vectors.
15 G cutoff (actually q+G cutoff) in the plane wave expansion of the eigenfunctions.
16 G cutoff (actually q+G cutoff) in the plane wave expansion of two particle objects such as and .
17 specifies whether QpGcut_psi and are in QpGcut_cou and are in atomic units (off), or in units of 2π/a (on).
18 corresponds to in Eq. 48 of Ref. I. alpha_offG 1d0 corresponds to offset points close to zero, and is usually a good choice. You can check convergence in *n1n2n3, varying this number.
19 if present, limits the number of bands used in computing ε(ω).
20 if present, limits the energy in computing ε.
21 if present, limits the number of bands used in computing Σ.
22 if present, limits the energy in computing Σ.
23,24 dw and omg_c fix the frequency mesh for real-axis integration of the polarizability, Eq. 32 of Ref. I, in atomic units. For small the mesh points are linearly spaced, separated by dw. If omg_c is not present, mesh points are linearly spaced on the entire grid. Otherwise the spacing gradually increases with , becoming twice dw at . Points on the mesh are interpolated to determine . Point is given by
is approximately dw/2.
25 Defines kind of quasiparticle self-consistency (how static Σ is extracted from the dynamical one).
iSigMode=3: Use (Eq. 10, mode-A in Ref. I).
iSigMode=1: Use (Eq. 11, mode-B in Ref. I).
iSigMode=5: Use (Eigenvalue-only self-consistency).
There is no “right” choice but 3 has the best formal justification.
26 Number of integration points along the imaginary axis to get , as described in II-F in Ref.I. The integration points are , where is the set of points for a Legendre quadrature in the interval [0,1]. (There is also special treatment of the sharply peaked contribution near .) Convergence is rapid in niw. Tests show niw=6 is sufficient to reach 0.01 eV accuracy in QP levels in Si.
27 Small imaginary part in Eq. 32 of Ref. I, in Hartree. turns on tetrahedron integration.
28 When is required, e.g. for the Z factor, it is computed numerically with a finite difference. deltaw is the energy spread used in the finite difference.
29,30 Broadening of the poles in the self-energy makers hsfp0 and hsfp0_sc. Poles of the Green function are smeared out by a width esmr. If GaussSmear is on, the smearing has the shape of a Gaussian function; otherwise a rectangular smearing is used. For insulators the value should be smaller than the bandgap (you can use small numbers like 0.001 or even smaller). For metals, the number should be a little larger. Some experimenting is needed, as it depends on the case. For Fe, esmr=1e-2 seems to be a good value.
31,32 In QSGW, self-energies are mixed with prior iterations to stabilize convergence. Provided a file sigm is available on disk (taken to be , i.e. the input ), hqpe_sc assembles results from the self-energy maker hsfp0_sc to make a new , and writes to sigm a new estimate for . It takes the linear combination
If mixpriorit is nonzero, hqpe_sc looks for pairs from prior iterations in file mixsigma. If they are available, Anderson mixing is used to estimate a new from a linear combination of and pairs from mixpriorit prior iterations (or as many as available, if fewer than mixpriorit). If mixpriorit is zero, is written. Anderson mixing used is the same as lmf uses to mix input and output charge densities, and some description can be found in the input file manual. Note in particular the discussion of parameters tj found the hqpe_sc output, file lqpe.
33-44 These are switches connected with the BSE solver that adds ladder diagrams to the polarizability. To add this correction, set UseBSEW to a nonzero value. You also need NumCondStates and NumValStates. Ladders are calculated after W is calculated in the RPA. They can be added in two contexts: when calculating the optics, i.e. the macroscopic dielectric function, and to modify the screened coulomb interaction W, when computing the self-energy Σ.
For optics, ε(ω) is computed i the window MinOmega_BSE < ω < MaxOmega_BSE, for q→0 along the vector given by Qvec_bse. The spectrum is broadened by a broadening which increases linearly in ω from the initial value ImagEner1 at MinOmega_BSE, to a final broadening at ImagEner2.
x0_small_BSE on turns on a special purpose mode that calculating ε, it subtracts the contribution to χ0 from the subset of transitions from the full χ0, and then adds on the BSE contribution. NumValBSE_no tells the BSE solver how many valence states not to include when calculating W.
The PRODUCT_BASIS category
An instance of this category, taken from NiO, is shown below. (The first column was added to label lines; they are described following the box.) NiO has two Ni atoms and two O atoms. Thus the information for sites 1 and 2 should be identical, as and also for sites 3 and 4.
Note: when generating a template for GWinput, lmfgwd automatically generates equivalent product basis information for equivalent sites. Whenever you edit GWinput by hand it is a good idea to check for consistency. lmfgwd can perform a number of checks for you, including checking equivalency of equivalent sites. To run a check, do
<PRODUCT_BASIS> ! Product basis block tolerance = minimum eigenvalue in PB overlap to reduce linear dependency 3e-4 3e-4 1e-3 lcutmx(atom) = l-cutoff for the product basis 4 4 3 3 atom l nnvv nnc ! nnvv: num. radial functions (valence) for augmentation-waves. nnc = num. for core. 1 0 2 3 1 1 3 1 1 2 3 0 1 3 2 0 1 4 2 0 ... 3 0 2 1 3 1 2 0 3 2 2 0 3 3 2 0 3 4 2 0 ... atom l n occ unocc :Valence(1=yes, 0=no) 1 0 1 1 1 ! 4S_p * 1 0 2 0 0 ! 4S_d 1 1 1 1 1 ! 4P_p 1 1 2 1 1 ! 4P_d 1 1 3 1 1 ! 3P_l 1 2 1 1 1 ! 3D_p 1 2 2 0 1 ! 3D_d 1 2 3 0 1 ! 4D_l 1 3 1 0 1 ! 4F_p 1 3 2 0 0 ! 4F_d 1 4 1 0 0 ! 5g_p 1 4 2 0 0 ! 5g_d ... 3 0 1 1 1 ! 2S_p * 3 0 2 0 0 ! 2S_d 3 1 1 1 1 ! 2P_p 3 1 2 0 0 ! 2P_d 3 2 1 1 1 ! 3D_p 3 2 2 0 0 ! 3D_d 3 3 1 0 1 ! 4F_p 3 3 2 0 0 ! 4F_d 3 4 1 0 0 ! 5g_p 3 4 2 0 0 ! 5g_d ... atom l n occ unocc ForX0 ForSxc :CoreState(1=yes, 0=no) 1 0 1 0 0 0 0 ! 1S * 1 0 2 0 0 0 0 ! 2S 1 0 3 1 0 1 1 ! 3S 1 1 1 0 0 0 0 ! 2P ... 3 0 1 0 0 0 0 ! 1S * ...
2 This line is used as a marker to identify the next line (keyword is tolerance), so it should be present.
3 A list of floating-point numbers. The size of list should range between 1 and 2×lmax+1 numbers, where lmax is the largest integer in line 5. The code will populate any numbers omitted with the last number you supply.
These numbers refer to tolerances in the overlap matrix of product basis functions. Product basis functions of different l are automatically orthogonal, so each l can be treated independently of the others. The overlap matrix constructed from product basis for a particular site and l is diagonalized to form an orthonormal set. Orthogonalized functions with small eigenvalues add little to the hilbert space of functions. Those with eigenvalues below tolerance are discarded, both reducing the size of the basis and avoiding numerical issues stemming from linear dependencies.
For small l it is safer to use small tolerances; but for higher l tolerances can be coarser for the same accuracy. Computational time is more sensitive to the higher l tolerances because they have more functions of different m per radial function.
4 This line is used as a marker to identify the next line (keyword is lcutmx(atom)), so it should be present.
5 specifies the l cutoff lmax for each site. Product functions are expanded to l = 2×lmax. You must supply one lmax for each site. lmax=3 is usually good enough for light sp elements, lmax=4 for transition metals. But as the nucleus becomes heavier and larger, lmax needs to be made somewhat larger. We have not tested the convergence for f shell elements, but lmax=6 is almost certainly safe. Smaller lmax reduces the computational cost, at the expense of accuracy.
By choosing rough tolerances considerable savings in time may be realized for systems with few atoms. But for larger systems, e.g. with 8 atoms or more, the time is limited by the plane wave part of the mixed basis, and there is little advantage to playing with this tolerance. The values autogenerated by lmfgwd –job=-1 are usually pretty safe.
6 This line is used as a marker to identify the next section (keyword is atom), so it should be present.
7-18 This table specifies the number of partial waves in the valence nnvv, and number of core levels nnc, one line for every site and l. nnvv will customarily be two for partial waves and , or three for channels that also have a local orbital . For the first (Ni) atom, there are three s core levels, but only a core because the is included in the valence as a local orbital (note nnvv=3 for the p). O (site 3 or 4) has only one core level, the .
This table is made automatically by lmfgwd –job=-1, and it is tied to the basis. If you change the basis, be sure to check that the table is still valid (lmfgwd –job=-2) or regenerate the table.
19 This line is used as a marker to identify the next section (product basis for valence states, keyword atom), so it should be present.
12-18 This table specifies which partial waves should be used to construct the product basis. There is one line for every l, valence partial wave, and every site. The lines in this table must synchronize with the preceding one. Recall that two partial waves combine to make a single product basis function. From the table in this section, two lists of functions are assembled: one for the first function which enters into the bra (|…>) part of the matrix element and another for the second which enters into the ket (<…|) part of the matrix element.
Column l refers to l quantum number of the partial wave
Column n refers to the type of partial wave (1 for , 2 for , 3 for ).
Column occ if nonzero, is added to the “bra” list
Column unocc if nonzero, is added to the “ket” list.
The remainder of the line is not used, but the template maker supplies a label so you can readily identify the kind of orbital it is. _p refers to principal wave (), _d refers to energy derivative () _l refers to local orbital ().
We have found that including the adds little to the effectiveness of product basis, so the template maker leaves them out. You can always include them by putting 1 in the occ and/or unocc columns.
In the sample section above, the default was altered. For Ni partial waves entering into the “bra” set is : , and the “ket” list is .
This basis is likely larger than needed, e.g. the inclusion of orbitals, and in the unocc column are probably not necessary.
44 This line is used as a marker to identify the next section (product basis for core states, keyword atom), so it should be present.
45-51 This section is the analog for core states of the preceding table. Columns l, n, occ, unocc serve the same function as for the preceding section. Column ForX0 tells the GW code to include this state in the calculation of the polarizability. There is an orthogonalization issue, however, see line 4 in the general section for a discussion.
Column ForSxc tells the GW code to include this state in the calculation of the self-energy.
Return to Table of Contents
The QPNT category
This category is designed for one-shot GW, to select which k points and QP levels for which to make the self-energy.
Note: if your system is metallic there are restrictions: you must calculate the self-energy for all k points and occupied states, in order to obtain the Fermi level.
A typical instance of this category is shown below. (The first column was added to label lines; they are described following the box.)
<QPNT> ! Specify particular k-points and bands (one-shot mode) --- Specify qp and band indices at which to evaluate Sigma *** Sigma at all q -->1; to specify q -->0. Second arg : up only -->1, otherwise 0 0 0 *** no. states and list of band indices to make Sigma and QP energies 4 4 5 6 7 *** q-points (must belong to mesh of points in BZ). 3 1 0.0000000000000000 0.0000000000000000 0.0000000000000000 3 -0.5000000000000000 0.5000000000000000 0.5000000000000000 7 0.0000000000000000 0.0000000000000000 1.0000000000000000 ... </QPNT>
The category as written above applies to 1-shot GW or QSGW calculations. (Extra information must be given when generating the dynamical self-energy; this will be described shortly.)
Lines preceding ‘***’ or ‘###’ are not used (only line 2 in this case). Lines beginning with one of these are markers that flag the start of a subsection following the line.
3-4. This is the marker for the first subsection and its contents. The two integers on line 4 should be 0 or 1. If the first is 1, QP levels are calculated for all q in the irreducible BZ. If it is zero, QP levels are computed only for a list following the third instance of *** (lines 8-13 in this case). For a metal, you should always pick 1 so the Fermi level can be found.
The second integer is usually 0. Choosing 1 causes the code to calculate QP levels for the first spin only. (Useful for an antiferromagnetic system such as NiO).
5-7. This is the marker for the second subsection and its contents. Its purpose is to select particular QP levels which calculate. (You can of course calculate them all but it is faster to make only a subset.) Line 6 tells the reader how many levels to calculate; line 7 lists the levels. It is an error for the list in line 7 to be smaller than number of levels you select.
8-13. This is the marker for the third subsection and its contents. Its purpose is to select the q points. (If the first integer in line 4 is 1, this subsection is not used.) Line 9 tells the reader how many points q you want QP levels for; call this nq. Then follow nq lines specifying the q points.
- The first number is not used
- Each k point must be a point on the mesh you specify; see n1n2n3, line 14 of the general category. The template maker generates a list of all the irreducible points (it was 8 points in this case). You can use a text editor to shuffle the rows (as was done here, to select the Γ, L and X points of a cubic material).
Extra input for the dynamical self-energy
hsfp0 –job=4 is a special mode that makes the dynamical self-energy. In this mode the reader parses lines in the QPNTS category until one is found begining with ‘****’ (the line must be present). Two numbers must appear in the next line. The two lines should be inserted before the normal part of this category described above.
An example is:
***** ---Specify the q and band indices for which we evaluate the omega dependence of self-energy --- 0.01 2 (Ry) ! dwplot omegamaxin(optional) : dwplot is mesh for plotting.
0.01 and 2 are the energy spacing and range over which is calculated. In this mode the usual dw and omg_c (lines 23 and 24 in the general section). The entire range is twice the second number: it extends both above and below the Fermi level.
The PBASMAX category
Not documented yet.
The QforEPS and QforEPSL categories
These categories apply to modes to calculate dielectric or magnetic response functions. Use one or the other.
Response functions are calculated at selected k points. If you add a line
then the list becomes all the q in the irreducible BZ.
Otherwise, use a QforEPS or QforEPSL category. In the former, you specify a list of q points, e.g.
<QforEPS> 0d0 0d0 0.01d0 0d0 0d0 0.02d0 0d0 0d0 0.04d0 </QforEPS>
Alternatively specify q points along symmetry lines (useful for plotting a response function along particular symmetry lines). For example, for a cubic system, the following yields seven q points from X to Γ and another seven from Γ to L (the left end point is omitted).
<QforEPSL> 1d0 0d0 0d0 0d0 0d0 0d0 8 0d0 0d0 0d0 .5d0 .5d0 .5d0 8 </QforEPSL>