ecalj manualhttps://github.com/tkotani/ecalj
September 5, 2017
Abstract
ecalj is an all-electron full-potential electronic-structure calculation package, especially for
GW/QSGW calculations, in addition to its unique one-body problem solver which uses bothaugmented plane wave (APW) and the Muﬃn-tin oribital (MTO) simultaneously. Available
freely from github.
This is based on the formulation in Refs. [1] and [2] (See Kotani2114QSGWinPMT.pdf andKotaniKinoAkai2015FormulationPMT.pdf at ecalj/Document/Manual/) on top of developments [3]. By ecalj , we can do QSGW and related calculations based on the PMT=LAPW+LMTO
method. To get minimum usage for ecalj, read up to the Section.5 or Sec.8. (We can generatemodel Hamiltonian through Wannier functions(ecalj/MATERIALS/CuMLWF, and so on, but
not documented well yet...)Be careful. This document still contain many bugs... Let me know bugs in this manual (tryto ﬁx things step by step). We need your help on this manual.
————————
Requirement
To support of our ecalj activity, we need your acknowledgment to ecalj in your publicationsin the manner shown in ecalj/README.md;
https://github.com/tkotani/ecalj/#we-need-acknowledgment-as such as:[1] ecalj package at https://github.com/tkotani/ecalj/.
Its one-body part is developed based on Ref.[2].
[2] LMsuit package at http://www.lmsuite.org/.
Its GW part is adopted mainly from Ref.[1].
(as in the same manner of usual papers).
Contributions
Main contributors to the ecalj package are (only a list of people who directly related to thecode developments):
T.Kotani, Hiori Kino, Mark van Schilfgaarde, Sergey Faleev, Takashi Miyake,
Seungwoo Jang, Manabu Usuda. They also contribute to this documents.
Relation with LM suit at http://www.lmsuite.org/:
There is a full potential LMTO package (with a slightly diﬀerent implementation of PMT)
contained in lmsuit. In the package, there is the driver to use a modiﬁed version of the
GW code in ecalj (GW part in ecalj is occasionally adopted in lmsuit ). The FP-LMTO in
lmsuit and ecalj have some similarities (for input ﬁles and outputs). Thus some part of its
document in the web site is still applicable even to the ecalj package.
History
T.Kotani learned the LMTO-ASA-GW code by F.Aryasetiawan, which is on top of the Stuttgart’s
LMTO-ASA, mainly by van Schilfgaarde under O. Jepsen and O. K. Andersen1http://www2.fkf.mpg.de/andersen/LMTODOC/LMTODOC.html. A FP-LMTO package is originally given by M.Methfessel, Mark van Schilfgaarde and their collaborators. T.Kotani developed an all-electron full-potential GW method on top of the FP-LMTO during his subbaticalat 1999-2000 at the lab of M.Schilfgaarde. Then the QSGW method has developed with thehelp of S.Faleev and M.Schilfgaarde, ﬁrst published in 2004. Then T.Kotani was adding new
functionality such as spin ﬂuctuation, impact ionization. In year 2008, T.Kotani started the
PMT=LMTO+LAPW method to remove problem of the FP-LMTO. After T.Kotani moved
to Tottori-u at 2009, he spends time for the improvement of PMT method. T.Kotani and
H.Kino showed that the PMT method can describe even the diatomic molecules eﬃciently
and accurately.It turns out very limited number of APWs can remove such problems. Then
T.Kotani have developed PMT-QSGW method. It is now stable, easy, and reliable rather
than the previous version of QSGW on top of FP-LMTO.

Contents
1 Introduction
1.1 Uniqueness of the ecalj package. . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1.1 What do we expect for QSGW? . . . . . . . . . . . . . . . . . . . . . . .
1.2 Rule in this manual . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 What can we do with the ecalj package? . . . . . . . . . . . . . . . . . . . . . .

6
7
8
9
9

2 Install10
2.1 Binaries and Scripts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2 tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.3 Directory structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
3 Theory (note)

12

4 LDA/GGA calculations and Plots
14
4.1 Write crystal structure ﬁle, ctrls . . . . . . . . . . . . . . . . . . . . . . . . . . 14
4.2 Generate default ctrl from ctrls by ctrlgenM1.py . . . . . . . . . . . . . . . . . 17
4.3 crystal structure checker: lmchk . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
4.4 ctrl ﬁle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
4.5 Run LDA/GGA calculations, and get convergence . . . . . . . . . . . . . . . . 18
4.6 DOS, PDOS plot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
4.7 Band plot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
4.7.1 syml genearator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
4.8 Useful samples: ecalj/MATERIAL/ . . . . . . . . . . . . . . . . . . . . . . . . . 21
4.9 How to add spin-orbit coupling? . . . . . . . . . . . . . . . . . . . . . . . . . . 22
4.10 PROCASR (VASP format) generator . . . . . . . . . . . . . . . . . . . . . . . . 23
4.11 efermi.lmf . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
4.12 Eﬀective mass calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
4.13 density and |eigenfunction|2 plot . . . . . . . . . . . . . . . . . . . . . . . . . . 25
4.14 LDA+U, partially occupied core-hole . . . . . . . . . . . . . . . . . . . . . . . . 25
5 How to run QSGW calculation?
26
5.1 GWinput . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
5.2 Run gwsc script . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
5.3 Spectrum function: How to calculate hqn|Σ(ω)|qni . . . . . . . . . . . . . . . . 29
6 gwsc script to perform QSGW
33
6.1 outputs of gwsc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
6.2 Preparation stage of gwsc . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
6.3 Main stage of gwsc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
6.4 Other functions (or scripts) . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
7 Usage problems, QandA error messages.

236

8 Cautions for usage37
8.1 lmf –help . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
9 Wannier function
43
9.1 lwmatK1 and lwmatK2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
10 ctrl file details
44
10.1 How to set local orbitals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
11 GWinput details
50
11.1 generate a template of GWinput . . . . . . . . . . . . . . . . . . . . . . . . . 50
11.2 overview of GWinput . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
11.3 General section . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
11.4 <QPNT> section . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
11.5 set QPNT for eps mode (QforEPS section) . . . . . . . . . . . . . . . . . . . . . 58
11.6 <PRODUCT BASIS> section . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
11.7 ANFcond (we can skip here since we do not check this option now. Need fix this if necessary. 61
12 Main Output Files of GW part63
12.1 QPU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
12.2 XCU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
12.3 SEXU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
12.4 SEXcoreU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
12.5 SECU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
12.6 TOTE.UP (TOTE.DN) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
12.7 TOTE2.UP (TOTE2.DN) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
12.8 DOSACC.lda . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
12.9 DOSACC2.lda . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
12.10Core ibas* l*.chk . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
12.11VXCFP.chk . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
12.12The Fermi energies in this GW code. . . . . . . . . . . . . . . . . . . . . . . . . 64
13 mkGWIN lmf2 and its I/O Files66
13.1 echo 0|lmfgw . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
13.2 gwinit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
13.3 echo -100|qg4gw . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
14 gwsc script and its I/O Files
68
14.1 echo 0| lmfgw si . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
14.2 echo 1| qg4gw . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
14.3 echo 1|lmfgw si . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
14.4 lmf2gw . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
14.5 rdata4gw v2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
14.6 echo 1|heftet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
14.7 hchknw . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
14.8 echo 3|hbasfp0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
14.9 echo 0| hvccfp0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
14.10echo 3|hsfp0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
14.11echo 0|hsfp0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
14.12echo 1|hsfp0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
14.13echo 11|hx0fp0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
14.14echo 12|hsfp0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
14.15echo 0|hqpe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
15 Check list for convergence on GW calculations

374

16 Linear response calculations
76
16.1 eps lmfh, epsPP lmfh: the dielectric functions . . . . . . . . . . . . . . . . . . . 76
16.2 epsPP lmfh: the dielectric function(No LFC— faster) . . . . . . . . . . . . . . 77
16.3 How to calculate correct dielectric funciton? . . . . . . . . . . . . . . . . . . . . 77
17 Utility

83•Reference4Change LOG:
• Jan2017 Add minimum description for syml generator at Sec.4.7.1

5

1

Introduction

The ecalj is an ﬁrst-principles electronic-structure calculations package with some uniquefeatures. With ecalj , we can do not only standard calculations (LDA/GGA/LDA+U, relaxation of atomic positions), but also the quasiparticle self-consistent GW calculations(QSGW),
linear responses (charge and spin), Wannier functions (and U of them).
This is base on an unique one-body problem solver, the PMT method (=the Linearized
APW+MTO method) [1]. Thus we identify the QSGW method implemented in ecalj as the
PMT-QSGW method. Introduction to the PMT-QSGW is given in Sec.1.1. Today “QSGW”
is accepted as a standard procedure in the electronic structure calculations [4].
First, see README.md shown at https://github.com/tkotani/ecalj#ecalj- (or ecalj/README.md
in the package). Free to download ecalj package from it, and use it. The QSGW code isversion controlled by git.
The ecalj is related to a FP-LMTO package lmv7 seen at
http://titus.phy.qub.ac.uk/packages/LMTO/fp.html. The lmv7 and ecalj are branchedoﬀ at year 2009. After branched, main contributions are due to T.Kotani and Hiori Kino
(NIMS) until now. We added new features: simple install and test; all codes are in f90 (noC compiler); new methods, especially PMT-QSGW; MPI parallelization for QSGW; simple
usage with automatic setting of default ﬁles by python ver2; a small tool to convert VASPPOSCASR to ecalj, and so on. PMT-QSGW shows more stable convergence than the previous
version, FP-LMTO-GW 1 .

——————–
ecalj package mainly consists of two parts. One is one-body part(PMT part) (in ecalj/lm7K/),
the other is many-body part(GW part) (in ecalj/fpgw/).
One-body part (PMT part), based on [1]We can perform standard calculations such as LDA,GGA,LDA+U, atomic position relaxation,and so on. In addition, the PMT part has an interface to perform GW (and QSGW) calculation: the one-body part can include a given non-local exchange-correlation potential stored ina ﬁle sigm.*. 2 . The QSGW calculation is performed by a script gwsc, which has an iteration
loop calling the one-body program (lmf) and many-body part (GW part) alternatively. The
many-body part generate the ﬁle sigm.*. See Fig.6 and around.
Many-body part (GW part), based on [2]As the inputs for the GW calculation, we have to supply the eigenfunctions and the eigenvalues
from the one-body part to the GW part. The eigenfuncions re-expanded by the two typesof basis functions, the atomic-like argumentation functions in the muﬃn-tin(MT) spheres,and the plane-waves in the interstitial region, say, the interstitial plane-wave (IPW) hereafter.
IPW is deﬁned as the usual plane waves in the interstitial region, but zero within MTs’. The
IPWs + ”atomic like functions within MTs” make “a basis set to expand eigenfuncitons”. See
Eq.(17) of [2].
We need another basis set to expand “product of eigenfuncions”. That is, the mixed
product bases (MPB)[5, 6], which consists of the two kinds of bases (caution: do not mixed up with the basis set
for eigenfuncitons); (i)the local atom-centered functions conﬁned to MT spheres, so-called the
product basis; (ii) IPW. The product-basis are calculated from products of solutions to the
Schrödinger equation within the MT sphere. The Coulomb matrix v, the dynamically screened
Coulomb interaction W , and so on, are expanded in the MPB. It can virtually span all the
space made of product of eigenfunctions (but, in practical calculations, we need to use a small

1 In cases to treat magnetic systems which have intrinsic magnetic fluctuations, we may need to be careful about
initial condition or mixing procedure to get convergence. In cases, we need to start from LDA+U results as initialcondition from which we start QSGW. Let me know about such trouble.
2 In the case of using sigm.*, total energy shown in the console output (also in save.*) are just the indicator, notthe meaningful total energy

6

size of the bases to reduce computational time). We include full energy-dependence of W (ω).
See Sec.3 of [2].
Recently T.Kotani includes the Wannier function generator, which was originally developed
by T.Miyake and H.Kino on top of previous version of GW part. Thus the Wannier functions(including eﬀective interactions) can be generated in the PMT-QSGW. (U parameters in the
full-screening and cRPA).

1.1

Uniqueness of the ecalj package.

We will explain two unique points of ecalj .
PMT:Central part in an electronic structure packages is one-body problem solver. It means how
to calculate eigenvalues/eigenfunctions for a given one-body potential. Inversely, we have togenerate new one-body potential for given eigenfunctions/eigenvalues based on the density
functional theory (DFT) in the LDA or GGA (In the followings, LDA means both of LDA and
GGA). Then we can make the electron density self-consistent by iterations until converged,
and obtain total energy of ground states. Then we can calculate atomic forces by perturbation.
Based on such an one-body problem solver, we can implement kinds of methods; e.g, dielectric
function, magnetic susceptibility, transport and so on. Furthermore, we can implement higherlevel approximations such as the QSGW method explained below. An one-body problem solver(in linear methods) are characterized by(i) linear combinations of what basis set to represent eigenfunctions;(ii) how to represent electron density and one-body potential.
In ecalj , we use the PMT method [1, 7] as the one-body problem solver. The PMT method is
a new all-electron full potential method. It uses not only the augmented plane waves (APW)
but also the muﬃn-tin orbitals (MTO) together, in addition to the local orbital (lo’s), torepresent the eigenfunctions (no other methods use two kinds of augmented waves together).
Thus eigenfunctions are expanded in the linear combinations of the APWs, MTOs, and the
lo’s. The formulation is clariﬁed in Ref.[1]. Then the electron density and the one-bodypotential are given in the “3-components representation”. That is, the electron density (onebody potential) is divided into three components,“smooth part + onsite muﬃn-tin (MT) part − counter part”.
Here the counter part is in order to remove smooth part within MTs. This formalism (SolerWilliams formalism [8, 9]) is also used in the projected augmented wave (PAW) method such
as VASP.
We now usually use highly localized MTOs together with APWs of low energy cutoﬀ (3 ∼ 4Ry). 3 I think this is promising not only for eﬃcient DFT/QSGW scheme, but also for kindsof applications in future.
QSGW:
In ecalj , we can perform the GW calculation. The usual GW approximation is so-called“one-shot GW” starting from LDA. It usually only calculates diﬀerences between the quasiparticle energies (QPEs) and the LDA eigenvalues by a perturbation (only diagonal part ofself energy for the LDA eigenfunctions). Its ability is limited; it may fail when its starting
point (eigenfunctions and eigenvalues supplied by LDA) is problematic. This is the reason
why we originally develop the QSGW method. The QSGW now becomes popular and taken
as a possible candidate to go beyond current limitation of such GW and LDA/GGA [4]. In
principle, results given by QSGW do not depend on LDA anymore; the LDA are only used toprepare initial condition for self-consistency iteration cycle of the QSGW calculation 4 .
Usually the QPEs obtained by QSGW reproduce experiments better than LDA. For example, the band gap by GGA for GaAs is about 0.5 eV in contrast to the experimental value of
1.69 eV 5 . On the other hand, the QSGW predicts about 1.8∼ 1.9eV, a few tenth of eV larger
3 current implementation have not yet efficiently use this locality; this must allow us to speed up one-body problem
solver.
4 Exactly speaking, we use LDA idea for efficient implementation of QSGW; thus obtained results are slightly
dependent on the choice of LDA or GGA
5 We undo electron-phonon effect (0.06eV) and spin-orbit effect (0.11eV) from the true the experimental value

7

than experiment (for practical use, we sometimes use “hybrid functional between QSGW and
LDA” so as to obtain smaller band gap). Even in the case of NiO and so on, the QSGW gives
reasonable results (there is a tendency to give a little larger band gaps than experiments).
This is in contrast to the case of the one-shot GW applied to NiO, where we can not havegood agreement with experiments because the stating points in LDA is problematic.
The ecalj have other functions. LDA+U, atomic forces and relaxation (in GGA/LDA),
core level spectroscopy and so on. In addition, we can calculate dielectric functions andmagnetic responses from QPEs and the quasiparticle eigenfunctions given by LDA/QSGW.
But total energy in QSGW is still in research (shown total energies in QSGW calculations aredummy now).
The QSGW calculations are very time-consuming; roughly speaking, it takes 10 or moretimes expensive than usual one-shot GW (although we can reduce computational time bychoosing computational conditions). Thus the size of systems which we can treat is limited
to ten atom in a cell or something, say, with a node of 16 cores; computation may require aweek or so to have reasonable convergence. (heavy atoms require longer computational eﬀorts,
light atoms faster; non-magnetic systems are easier. We still have much room to accelerate
the method, but not have done yet so much. Minimum MPI parallelization is implemented).
The computational eﬀort is ∝ N 4 in the most time-consuming part of QSGW.

1.1.1

What do we expect for QSGW?

Let us recall hybrid functionals such as B3LYP, and LDA+U. In hybrid functional methods,
we use Vxc = (1-α)*LDA+α*(Fock exchange like term), where α is taken to be ∼ 0.25 usually
6. The α can be dependent on materials; for metals α should be almost zero. For larger band
gap insulator, α becomes larger. 7 Despite of success of the functional, its ability is limited.
For example, it is known that a hybrid functionals fail to describe metals such as bcc Fe. On
the other hand, we have LDA+U method which succeeded to describe materials including
localized electrons. However, it contains kinds of ambiguity and U is chosen by hand. The
important part of the hybrid functional methods and LDA+U is the non-local potential. It is
missing in the DFT. As we discussed above, they give some success but not satisfactory. We
somehow need to have a method to determine high-quality non-local potential (a substitutionof the exchange-correlation potential in LDA). It is the QSGW method.
Note two important aspects of non-local potential (missing character in the local potential
used in DFT). One is the onsite non-locality; it is also taken into account by LDA+U model.
However, note that relative shift of O(2p) band with respect to the center of 3d band is not in
LDA+U. The other is the oﬀ-site non-locality (mainly between nearest neighbors), which may
relate to LUMO-HOMO gap. A non-local potential can behave a projector which push down
only the HOMO states (valence band) to lower energy. This can be in the hybrid functional
but not in LDA+U.
In the QSGW, we determine such a non-local potential with the calculation of the GW
method, in a self-consistent manner (we repeat GW calculations until converged). We can
expect QSGW much more than hybrid methods/LDA+U. Roughly speaking, because the
QSGW automatically determine U of LDA+U, or alpha of the hybrid functionals. Moreaccurately speaking, we determine not only G0 but also W (the screened Coulomb interaction)
self-consistently. Here W corresponds to U and alpha. Thus QSGW gives reasonable results
even if it is applied to metals such as Fe. For systems with metallic screening, it gives small
non-locality (results are close to those of LDA). For systems with large band gap, QSGW gives
large enough non-locality (like 0.25*(Fock exchange)).
Since we now need to treat complex systems, e.g, metal on insulator, it is very essential totreat kinds of materials on a same footing.

8

The main purpose of QSGW is to determine an one-particle eﬀective Hamiltonian H 0, which describes the quasiparticle picture (or independent-particle picture) for the system

1.52 eV.
6 exactly speaking, we have range cutoff for Fock exchange term in the HSE functional in addition
7 If you use α = 1 (Hartree-Fock limit), the band gap of Si becomes 20eV or something.
8 people often pronounce thisH-naught

8we calculate. In other words, QSGW divides the full many-body Hamiltonian H into H =H 0 + (H − H 0 ). The screened Coulomb interaction W is determined self-consistently in the
QSGW iteration cycle.
In comparison with LDA, we see diﬀerences;• Band gap. QSGW tends to give slightly larger than experiments. It looks systematic.
• Band width. Usually, sp bands are enlarged (except very low density case such as Na).
This is the case for homogeneous electron gas. As for localized bands like 3d electrons,
they can be narrowed.
• Relative position of bands. e.g. O(2p) v.s. Ni(3d). More localized bands tends to get
more deeper. Exchange splitting between up and down (like LDA+U) get larger. In
cases such as NiO, magnetic moment become larger; closer to experimental values.
• Hybridization of 3d bands with others. QSGW tends to make eigenfunctions localized.
However, reality is complexed, and not so simple in cases.

1.2

Rule in this manual• This font is for executable ﬁle(program) or shell scripts.
• echo 3|hbasfp0 means doing hbasfp0 with the argument ’3’ supplied as the standard
input (read(*,*) in fortran).
• This font is for ﬁles, directories, contents of ﬁles, or variables used in codes.
• ctrl.si,rst.si and so on mean the case of Si. You may need to replace the extension si foryour case. (this extension is given by user. Lower case, number, and underscore [a-z0-9 ]are allowed.) In the followings, ctrl.* means a ﬁle wish such an extension.
• There are ﬁles named foobarU and foobarD, which are for up spin (isp=1) and for downspin(isp=2), respectively; for example, SEXU and SEXD. We sometimes use foobarU to
denote foobarU and foobarD together.
• k vector in the Brillouin zone is called as q or k.

1.3

What can we do with the ecalj package?

At Feb.2015, what we can do is as follow. We have limited parallelization. (e.g. k point
parallel).
• LDA/GGA LDA+U, calculations, atomic forces and relaxation. Spin-orbit is included
only for co-linear spin-density cases.
• Quasi-particle(QP) energy in the 1st-iteration from LDA. (one-shot GW )Make band plot for LDA and the QP energies.
• Spectrum function of the self-energy Σ. Life time (imaginary part) of QPs.
• Dielectric function, and its inverse. (including local-ﬁeld eﬀect or not).
• QP self-consistent GW (QSGW)• magnetic susceptibility.
• Wannier function. (not only one-body part, but also eﬀective interaction W and cRPA)

92

Install

Install and minimum tests are easy; even in a note PC, e.g., we can use gfortran in Ubuntu
14.04 on e.g., Thinkpad T420s for a test purpose. For productive runs, we may need multicores. Current implementation for parallelization by MPI is limited (not so much especiallyfor the dielectric function part yet). Thus, probably, it may be not so eﬃcient to use too many
cores.
Follow the instruction of ecalj/README.md.
or we can see the same one at https://github.com/tkotani/ecalj/#install-and-test. We
have a command ecalj/InsallAll.ifort and so on. This command installs ecalj and run aseries of install tests automatically.

2.1

Binaries and Scripts

Binary and Scripts contained in ecalj are• ctrlgenM1.py
Generate default input ﬁle ctrl.* from the structure ﬁle ctrls.* . The latter ﬁle only
contains information of crystal structure.
• lmfaSpherical atom calculation as initial condition, and core charge• lmchk
Check atomic positions, crystal symmetry, and computational conditions.
• lmf and lmf-MPIK LDA/GGA,LDA+U calculations. (or we can use Vxc in QSGW
instead). We mainly use lmf-MPI (k-parallel version) instead of lmf.
• PROCAR mode of lmf: Fat band mode.
mpirun -np 4 lmf-MPIK --mkprocar --band:fn=syml mgo gives PROCAR(Try an example /ecalj/MATERIALS/MgO PROCAR/. Run ./job at the directory).
• gwscQSGW calculation• job_band,job_fermisurface,job_tdos,job_pdos
band, fermi surface, tdos, pdos plot.
• epsPP lmfh
Dielectric funciton without local ﬁeld correction (LFC).
• eps lmfh
Dielectric funciton with LFC• epsPP lmfh chipm
Non-interacting transverse spin polarization.
• gw lmfh
One-shot GW calculation. This also show life-time of QPs (QPU lmf). (we need make itparallelized...)• genMLWF
Wannier functions and matrix elements of W on it. A implementation of cRPA included.
• dqpuA small python script to compare QPU.* ﬁles (eigenvalues are compared) numerically.
(Seungwoo say it cause “index error”; Need to ﬁx dqpu if something strange).

2.2

tests

Install.ifort run tests at ecalj/TestInstall.
In the following, si:gw_lmfh means ’>make si_gw_lmfh’ at ecalj/TestInstall/; this test is
performed with the Makeﬁle at the directory.

10

si:gw_lmfh/
si:gwsc/
gas:gwsc/
nio:gwsc/
fe_epsPP_lmfh_chipm:
gas:eps_lmfh/
gas:epsPP_lmfh/

Results: QPU: QPU,log.si: QPU,log.gas: QPU,log.nio: ChiPM*: EPS*: EPS*(These are just samples; not for practical calculations)2.3

Directory structureecalj├── InstallAll.ifort, Install.gfortran ! install and test├── Document/├── fpgw/! full potential GW code├── fpgw/Wannier ! Miyake’s Wannier code is reorganized here.
├── lm7K/! PMT method├── MATERIALS/ !job_materials.py contains examples.
├── StructureTool/ ! POSCAR converter, and a utility to├── TestInstall/! Install test; this is invoked from Install.*└── TOOLS/!Tools for developers

11

3

Theory (note)

Except the technical details, we need to know minimum for these theories.
• DFT in LDA/GGA.
• GW• QSGW
There are literatures as for GW. A recent one is Ref.[4]. In addition, it is better to know thebasics of the PMT method (LAPW+LMTO method) [1]. Here is a small note for GW .
Green function
The Green’s function G(r, r′ , ω) is the central quantity in the GW calculation. In the oneparticle theory (mean-ﬁeld theory/non-interacting case), the Schrödinger eq. is∂− H0 )ψ(r, t) = 0(1)∂t
Here H0 is the one-particle Hamiltonian which contains electrostatic potential plus exchangecorrelation potential Vxc (here we don’t care how it is given. This should be static but can be
non-local as Fock term). When we have (unkown) source term J(r, t) instead of 0 in the right∂−H0 ).
hand side of Eq. (1), we can calcuate ψ(r, t) by multipling inversion of the operator (i ∂t′This is the (non-interacting) Green funciton G(r, r , ω) deﬁned as(i∂− H0 )G = δ(t − t′ )δ(r − r′ ).
(2)∂t∂
As this shows, G is the inverse matrix of the operator (i ∂t− H0 ), (we pay attention to theboundary condition, especially for time-direction: retarted, advanced or time-ordered). Thus∂
− H0 ) = 1/(ω − H0 ). In other words, the Green function G is thewe can write G = 1/(i ∂tintegration kernel in order to convert a diﬀerential equation Eq. (2) to an integral equation.
This is the same as the case of Poisson equation (Laplacian) for the electrostatic problem;then
1/|r − r′ | is the Green function for the conversion (then we calculate electro static ﬁeld φ(r)for give source term ρ(r).
In the case of electro static problem, we have two ways to solve the problem. One is the
direct method to solve diﬀerential equation, the other is using such integration kernel. They
are equivalent, but the latter is easy to handle and suitable for numerical calculations (under
the assumption of “superposition low”, that is, linear response). Although G(r, r′ , t − t′ )contains time variables (to describe wave propagation), essentials are the same.
(i

Green function: many body caseG(r, r′ , ω) is the quantity which is well deﬁned even in the many-body perturbation theory. It
is deﬁned as the expectation value as h0|ψ̂(rt)ψ̂ † (r′ t′ )|0i (this is not accurate. See literatures..).
Here we use the second-quantized ﬁeld-operator ψ̂. This deﬁnition of G is a natural extensionof the one-particle case.
However, there are diﬀerences in the G of one-particle case and many-body case. In the
one-particle case, we have no interaction between electrons. Thus an electron added to theground state (which is by ﬁlling electrons up the Fermi energy) is the eigenstate. An electron
going on can not be an eigenstate, because of the correlation eﬀect and exchange eﬀect.
Let me explain the correlation eﬀect. In contrast to this, one-particle state (=the quasiparticle(QP)) in the many-body theory can be not the eigenstate; it moves in the sea of ohter
electrons and holes, which can be excited to other states, or an moving electron can move withpolalized cloud of other electrons and holes. Its motion is aﬀected by oher electrons. It willlose energy gradually; the QP can have life time. This may be identiﬁed as the correlationeﬀect.
Let me explain the exchange eﬀect. This is due to the Fermi statistics. The electron, which
you focus on, can not be distinguished from other electrons. (Thus mu-on can not feel thiseﬀect. G for mu-on in solids only contains pure correlation eﬀect). From a point of view, thisexchange eﬀect is interpreted as a hopping eﬀect, sudden jump from an electron to anotherelecron (actually the Hartree-Fock theory gives zero eﬀective mass at the Fermi energy formetallic systems).

12

We can include both of the eﬀects in the GW calculation at the lowest order. Exactly
speaking, these two eﬀects are really mixed up. However, in GW (since GW is at the lowest
order), these two are clearly separated.
The rigorous equation for G is


∂
(3)
i − H0 + (Σ(r, r′ , t − t′ ) − Vxc ) G = δ(t − t′ )δ(r − r′ ).
∂tOnly the diﬀerence is adding Σ(r, r′ , t − t′ ) − Vxc . In other words, all unknown eﬀects arepushed into this term. Many body eﬀect are pushed into (downfolded or projected into) the
self-energy (=dynamical one-particle potential) Σ(r, r′ , t − t′ ). We don’t treat many-bodyquantities such as X(r1 t1 , r2 t2 , r3 t3 ) directly in the GW calculations.
The GW method is how to give the Σ(r, r′ , t − t′ ) in the lowest order. Because of the longrange property of the Coulomb interaction, we need to have a special technique (not a simple
perturbation). Σ(r, r′ , t − t′ ) is calculated as the product of G × W , where G is (usually) justthe non-interacting Green function given in Eq. (2). W is the dynamically screened Coulombinteraction calculated in RPA.
QSGWNote that Eq. (3) completely cancels the eﬀect of Vxc (r, r′ ) since it is included in H0 . Thus
this division is not meaningful if we can solve a problem completely. However, what we can dois only the perturbation. Quantities for H0 is taken as the basic quantities. On it, we apply
the perturbation. Equivalently we can have a division H = H0 + (H + H0 ), where H is the
many-body Hamiltonian.
In priniple, we can start from any Vxc . However, as long as we use perturbation theory, we
have to minimize the eﬀect of Σ(r, r′ , t − t′ ) − Vxc . The degree of freedom of the choice of Vxc
is used to minimize it. This is the self-consistent pertubation theory idea for QSGW.
If Σ(r, r′ , t − t′ ) − Vxc give small eﬀect, we can use the concept of the QP (independent
particle picture) as the basis to evaluate physical quantities.
We have to determine best (or optimum) Vxc . How to do it? An idea is that we should
determine it so that the size of Σ(r, r′ , t − t′ ) − Vxc should be smallest (how to measure thesize?). Or, based on the Landau-Silin’s theory, the QPs contained in G should be reproducedby H0 . In anyway, this ends up with the self-consitent petrurbation theory. For trial Vxc (andelectron density), we have H0 ; from H0 we can calculate Σ(r, r′ , t − t′ ). Then we have todetermine new Vxc in a method. This is repeated until converged. No unique choice for the
method. Thus we testes some possible ways and now we usually use a standard procedure as
described in Ref. [2].

13

4LDA/GGA calculations and Plots

Calculations are performed by following steps. These steps are detailed in the following subsections.
To identify a set of ﬁles used for a material we calculate, we use an extension to ﬁles. For
example, ﬁles explained below are with extensions (only lower case allowed) of materials. For
example, ctrls.cu and ctrl.cu. In this case cu is the extension. Any extension works. Other
possible examples are ctrls.lagao3, ctrl.wgantest1, and so on.
1. Write crystal structure ﬁle ctrls.* , which contains crystal structure. It can be by hand,or convert it from POSCAR (in VASP). There is a tool to convert between POSCAR
and ctrls.* (See ecalj/StructureTool/README.txt for the tool. we have vasp2ctrl andctrl2vasp. Type them without arguments to see help.)
2. Generate ctrl.* from ctrls.* by a script ctrlgenM1.py, 9Here ctrl.* is the main control ﬁle which contains all required information to perform
calculations. ctrl.* contains not only the content of ctrls.* , but also other informationneeded for calculations. If necessary, we edit the generated ctrl.* ﬁle before next step.
3. Check crystal structure. lmchk is to conﬁrm the crystal structure (space-group symmetry and so on). lmchk is applied not to ctrls.* but to ctrl.* . ctrls.* never used in thefollowing steps.
CAUTION for a known bug! : If a crystal structure is only slightly diﬀerent froma structure with higher symmetry, the lmchk may give a wrong crystal symmetry. Insuch a case, you have to “standardize structure” by VESTA or some other tools. Thisoccurs e.g., when we use a structure numerically relaxed by VASP.
4. Run lmfa (calculations of spherical atoms (MT sites) in the cell). It also calculates coreeigenfunctions and valence electron charge to set up initial condition. Then we run main
calculation of LDA by lmf. It repeats iterations, and end up with converged results in
LDA. Main result (electron density satisfying self-consistency) is stored in restart ﬁlerst.* (binary ﬁle). It ﬁnished within a second.
5. Run LDA/GGA calculations. We can run the LDA/GGA calculation by lmf or lmfMPIK (-MPIK means kpoint parallel version).
6. Post processing.
Plot energy band, DOS, PDOS, by running scripts. We can use scripts job band forband plot (need syml.si ﬁle (symmetry line for band plot)). We also have job pdos,
job tdos, job fermi and so on for DOS, PDOS, fermi surfaces. Since we use gnuplot
to plot them, meanings of obtained data is apparently clear.

4.1

Write crystal structure file, ctrls

Let me show some samples of crystal structure ﬁles ctrls.* .
Cu: ~/ecalj/lm7K/TESTsamples/Cu/ctrls.cu------from here -----------------% const da=0 alat=6.798
STRUC ALAT={alat} DALAT={da}PLAT= 0.0 0.5 0.5
0.5 0.0 0.5
SITE
ATOM=Cu POS=0 0 0------to here ------------------

0.5 0.5 0.0

GaAs: ecalj/lm7K/TESTsamples/GaAs/ctrl.gaas------from here -----------------#id = GaAs
%const bohr=0.529177 a=5.65325/bohr
STRUC
9 ctrlgenM1.py exists originally at ecalj/lm7K/ (ctrlgenM1.py was already copied to your BINDIR= defined inecalj/Install.ifort in the installation).

14

ALAT={a}PLAT=0 0.5 0.5

0.5 0 0.5

0.5 0.5 0

SITE
ATOM=Ga POS=0.0 0.0 0.0
ATOM=As POS=0.25 0.25 0.25------to here -----------------SrTiO3: ecalj/lm7K/TESTsamples/SrTiO3/ctrls.srtio3------from here -----------------%const da=0 au=0.529177%const d0=1.95/au a0=2*d0 v=a0^3 a1=v^(1/3)
HEADER SrTiO3 cubic
STRUC
ALAT={a1} DALAT={da}PLAT=1 0 0 0 1 0 0 0 1
SITE
ATOM=Sr POS=1/2 1/2 1/2
ATOM=Ti POS= 000
ATOM=O POS=1/2 00ATOM=O POS= 0 1/2 0
ATOM=O POS= 0
0 1/2------to here -----------------Lines starting from ’#’ are neglected as comment lines. Lines starting from ’% const’deﬁne variables and set values (in these cases, da, alat, and bohr, and so on). Then thevariable alat is referred to as {alat}; in the cu case, {alat} means 6.798. Lines not startfrom ”#” nor ”%” are main content in ctrls.* . 10
Note that we have two tags of “categories” ”STRUC” and ”SITE”. (“HEADER” tag is
also; but it is just for user’s memo shown in console output). These tags should start from theﬁrst column. Thus ctrls.* is divided into multiple “categories”. In a category, we have “tokens”
such as ALAT, DLAT, PLAT. These under STRUC category. ALAT+DALAT specify unit oflength in this ctrl ﬁle. These are in a.u. (= bohr radius=0.529177Å).
The unit cell is given by PLAT (as noted, ALAT+DALAT as unit). In the above example of GaAs, three primitive cell vectors speciﬁed by nine numbers after PLAT=; they give
three primitive vectors; PLAT1=(0,0,0.5), PLAT2=(0.5, 0.0, 0.5), and PLAT3=(0.5, 0.5, 0).
DALAT is convenient to change lattice constant; but it is ﬁxed to be zero here; thus no eﬀect
in this example.
Note that SITE category can have multiple ATOM tokens. The number of ATOM token
under SITE should be the same as number of atoms in the primitive cell. In the case of GaAs;
SITE contain multiple ATOM tokens. POS= just next to ATOM is taken as subtokens under
ATOM token. 11 In cases, we specify such subtokens as SITE ATOM POS.
In the SITE category, we place atoms (MT names) in the primitive cell. In these cases we
use defaults atomic symbol (MT names) for ATOM. POS is in the Cartesian coordinate (in theunit of ALAT+DALAT).
To test ecalj, you may make a test directory and copy a ctrls.* to your directory. If you
have VESTA and ecalj/StructureTool/ installed, you can see its structure by$ viewvesta ctrls.cu(here $ means command prompt).
NOTE: As written in ecalj/README, you have to install VESTA and viewvesta.
Then set VESTA= at the top of ecalj/Structure/viewvesta, and make softlink
to it. The command viewvesta(~/ecalj/StructureTool/viewvesta.py) generate
POSCAR_cu.vasp ﬁrst, then send it to VESTA. viewvesta also accept POSCAR_cu.vasp
directly. Except names starting from ctrl and ctrls, viewvesta sends the nameto VESTA directly. We need extension ’.vasp’ to recognize it is written in VASP
10 For these variables, we can overlay values when we start programs. For example,’lmf -vdalat=0.1 si’ means thatalat is recorded in save.si file.
11 This may looks slightly uncomfortable since the end of range of ATOM is not clearly shown; it end just at thenext ATOM token or new category.

15

format. We have samples in ~/ecalj/StructureTool/sample.
A tool vasp2ctrl converts POSCAR .vasp to ctrls.. “–help” show a small help.
• ecalj/StructureTool/ is not tested well. Not believe it so much... We will ﬁx it
on your request. Another possible way is using cif2cell.
If you have a cif ﬁle, runcif2cell foobar.cif -p vasp --vasp-cartesian --vasp-format=5
And convert POSCAR to ctrls. cif2cell is available from github.
In ctrls.srtio3, we use an expression 1/2 to give POS. We can use mathematical expression instead of values. Mathematical expressions such as “+ − ∗/ sqrt(...)” are recognized.
(instead of 3**2, use 3^2. Use parenthesis, and no space for an expression). We can use default
atomic symbols (to check default atom name (MT name) type ctrlgenM1.py --showatomlist).
Instead of such default symbols, we can use your own symbol as
SITE
ATOM=M1 POS=1/2 1/2 1/2
ATOM=M2 POS= 000
ATOM=O POS=1/2 00ATOM=O POS= 0 1/2 0
ATOM=O POS= 0
0 1/2
SPEC
ATOM=M1 Z=38
ATOM=M2 Z=22
ATOM=O Z=8. Then we have to add extra category SPEC where we set Z number. (You can use Z=37.5
for virtual crystal approximation, however, you can not do it in ctrls now. Edit it in ctrl ﬁle.
Such a procedure will be explained in other place.xxx)
This is an example for Antiferro NiO:
#id = NiO
%const bohr=0.529177 a=7.88
STRUC
ALAT={a} PLAT= 0.5 0.5 1.0 0.5 1.0 0.5
SITE
ATOM=Niup POS= .0.0.0
ATOM=Nidn POS= 1.0 1.0 1.0
ATOM=O POS= .5.5.5
ATOM=O POS= 1.5 1.5 1.5
SPEC
ATOM=Niup Z=28 MMOM=0 0 1.2 0
ATOM=Nidn Z=28 MMOM=0 0 -1.2 0
ATOM=O Z=8 MMOM=0 0 0 0

1.0 0.5 0.5

In this case, we deﬁne Niup and Nidn sites. These are recognized as Ni atom because of given
Z number in SPEC. The subtoken MMOM=Ms,Mp,Md,Mf... re to specify number of magnetic
moments (µB ) for s,p,d,f channels (diﬀerence of up - down electrons within MT sites) as initialcondition. In this case, we set n(up)-n(down)=1.2 for Niup site for d channel. Even just one
ATOM name is given by yourself, all ATOM in SPEC should be given (in this case SPEC for
O should be given).
We can see other samples in ~/ecalj/lm7K/TESTsamples/*/ctrls.*. (we also have asample generator. See later section.) Note that ctrls ﬁle is jut in order to generate default ctrl
ﬁle in the followings. Not from ctrls but from ctrl, we can start calculations. (thus ctrls is notneeded if we prepare ctrl ﬁle directory).
It is possible to add RELAX= 0 0 1 after SITE_ATOM_POS; this means structure relaxation
along z-axis (also need to set DYN category
as seen at http://titus.phy.qub.ac.uk/packages/LMTO/tokens.html#DYNcat), but its defaults are given (but commented out) automatically in the ctrl ﬁle generated by the proceduredescribed in the following section). We detail it in other place xxx.

16

After ctrl.* is generated as shown below, we can run a command lmchk to check weather
crystal structure is correctly given or not. It ﬁnish in a second. It show symmetry information,and so on used in the calculation.
CAUTION!: Positions of atoms are not necessarily ﬁxed by ctrl.* when you restart calculation with rst.* ﬁle, because atomic positions are read from rst.*. We need to pay attention
when we use DYN option because lmf run may save relaxed atomic positions into the rst.*.
As “lmf –help” shows, lmf si --rs=1,1,1,0,0 can read atomic position from ctrl ﬁle.

4.2

Generate default ctrl from ctrls by ctrlgenM1.py

To run programs of lm7K (lmfa, lmf, lmchk) in ecalj, we need an input ﬁle ctrl.*, whichcontains not only structures but also other settings. To generate ctrl.* from ctrls.*, wehave a command ”ctrlgenM1.py” (written in python 2.x and call fortran programs(lmfa,lmchk)
internally). Two steps required to complete ctrl ﬁle: (i) we give reasonable options when werun ctrlgenM1.py. Then (ii) we may edit the ctrl ﬁle afterward. In anyway, ctrl ﬁle is thestarting point of calculations; ctrls is required just in order to generate ctrl.
At ﬁrst, try ctrlgenM1.py without arguments. It shows help. To generate ctrl from ctrl,
type$ ctrlgenM1.py cu --nk1=8Here cu specify ctrls.cu. The option –nk1=8 means the number of division of the Brillouinzone for integration. It means 8x8x8 division. If we like to use 8x8x4, we have to supply threearguments –nk1=8 –nk2=8 –nk3=4. The above command gives following console output.
$ ctrlgenM1.py cu --nk1=8=== INPUT arguments (--help gives default values) ===--help Not exist--showatomlist Not exist--nspin=1--nk=8--xcfun=vwn!(bh,vwn,pbe)--systype=bulk !(bulk,molecule)--insulator Not exist !(do not set for --systype=molecule)...
OK! A template of ctrl file, ctrlgen2.ctrl.cu, is generated.
As we see above, options which you speciﬁed are shown at the beginning of the console output(in this case –nk1=8). Others such as –nspin=1 are default settings. If we like to performspin-polarized calculations, we add other option ’–nspin=2’ asctrlgenM1.py nio --nspin=2 --nk1=6(NOTE: In the spin-polarized case, we need to set initial condition of size of magnetic momentat each atoms. Set it in ctrls.* as in the previous section, or edit MMOM of ctrl ﬁle(MMOM=s p d f ...) to be like MMOM=0 0 1.2.). The ctrlgenM1.py generates ctrl ﬁle named
as ctrlgenM1.ctrl.cu. To do calculations, copy it to ctrl.cu so that lmf can recognize it.
cp ctrlgenM1.ctrl.cu ctrl.cu

4.3

crystal structure checker: lmchk

Do lmchk to conﬁrm that we can let lmf know correct crystal structure. It also show crystal
structure informaiton, equivalent sites, site index and so on.
lmchk --pr60 cu(--pr# gives more informations if # is number)

Then it reads ctrl.cu. --pr60 is an option of verbose. Bigger number gives more information.
• Lattice info, Space group symmetry operations (in lmf format), and their generators(these operations can be generated from a few of them.)
See http://titus.phy.qub.ac.uk/packages/LMTO/tokens.html#SYMGRPcat+ about how
to represent the operations.

17• Show atomic positions in ctrl ﬁle.
• Tabulate MT radius and distance between atomic sites.
(lmchk –help shows help, but diﬃcult to see. Not need to read it ﬁrst.)
lmchk is also shows atom (MT site) id (position and class(equivalent positions). This isneeded to interpret PDOS.
CAUTION for a known bug! See 3 in Sec.4.:

4.4

ctrl file

It is not necessary to look into ctrl ﬁle ﬁrst, although some details are explained in thegenerated ctrl ﬁle. Please compare obtained results by lmf with those by other packages orliteratures; let me know if you ﬁnd something strange or your questions.
It is necessary to edit ctrl ﬁle to use full ability of lmf. For example, LDA+U, atomic
position relaxation, core level spectroscopy, Change setting of default MTO and lo, better
mixing procedure for stable convergence; higher accuracy, and so on.
But a few of ctrl ﬁle is easy to modify. Search these words and read explanations embedded
in ctrl ﬁle.
(1)XCFUN(choice of XC—it is not need to repeat ctrlgenM1.py). It is also possible to change number ofk points for sampling, to modify crystal structure slightly, and so on; all things needed are in
ctrl. It is not needed to repeat ctrlgenM1.py again.
(2)SOTo obtain correct dispersion around top of valence at Γ point for GaAs, we need to set SO=1
and NSPIN=2. QSGW calculation (by gwsc) do not allow this option now; Thus we run lmf(or lmf-MPIK) with such settings changed in ctrl ﬁle, after QSGW is converged.
lmf --input shows what can we write in ctrl ﬁle. But more than half are not for users,
but for developers (or irrelevant now).

4.5

Run LDA/GGA calculations, and get convergence

Here we show how to get converged results from a ctrl ﬁle.
At ﬁrst, we need initial guess of charge density. It can be given by a super position of
atomic charge density. To obtain the charge density, we solve atoms ﬁrst. It is by$ lmfa gaas | tee llmfa
It takes just a few seconds. Here tee is a command of Linux. It keeps console output (standard
output) to a ﬁle (llmfa in this case).
Then try$ grep conf llmfa. Then you see a key point that
conf:SPEC_ATOM= Ga : --- Table for atomic configuration --conf: isp l int(P) int(P)z
Qval
Qcore
CoreConf
conf:
1 0
4 0
2.000
6.000 => 1,2,3,
conf:
1 1
4 0
1.000
12.000 => 2,3,
conf:
1 2
4 3
10.0000.000 =>
conf:
1 3
4 0
0.000
0.000 =>
conf:
1 4
5 0
0.000
0.000 =>
conf:----------------------------------------------------conf:SPEC_ATOM= As : --- Table for atomic configuration --conf: isp l int(P) int(P)z
Qval
Qcore
CoreConf
conf:
1 0
4 0
2.000
6.000 => 1,2,3,
conf:
1 1
4 0
3.000
12.000 => 2,3,
conf:
1 2
4 3
10.0000.000 =>

18

conf:
1 3
4 0
0.000
0.000 =>
conf:
1 4
5 0
0.000
0.000 =>
conf:----------------------------------------------------This is an initial electron distribution, and how we divide core and valence. In this case core
charge Qcore are “6 electron for s channel=1s,2s,3s and 12 electron for 2s and 3p”. Qcore istreated by frozen core approximation. See Sec.2.5 in Ref.??. Qval means electrons for valences,p,d channels. The valence channels are 4s,4p,4d,4f (if we set EH=s,p,d,f) in this case The
int(P)z column is for local orbital. Thus we have 3d treated as local orbital. (ecalj allow add
one local orbital per l.)
The isp index means spin (1 or 2), since –nspin=1 (when we invoke ctrlgenM1.py) for
GaAs, no isp=2 exist. In summary we have 4s,4p,4d,3d,4f as valence. This means we use
corresponding number of MTOs and local orbitals.
After lmfa, let us start main calculation.
$ lmf cu
In unix, we ca save console output to llmf by $ lmf cu | tee llmf. As it starts iteration
calculations, it shows similar output again and again (it is a little too noisy now). Then you
end up with self-consistent result as......
it 8 of 30
ehf=-3304.895853
ehk=-3304.895853
From last iter
ehf=
-3304.895856
ehk=-3304.895855
diffe(q)= 0.000003 (0.000007)
tol= 0.000010 (0.000010)
c ehf=-3304.8958531 ehk=-3304.8958529
Exit 0 LMF
CPU time:
7.024sMon Aug 19 02:03:19 2013on

more=Fit 8 of 30 means it stop at 8th iteration, although we set maximum number of iteration
30. Note that this number is given by ITER_NIT=30 in ctrl.cu). ehf and ehk are the groundstate energy in Ry. They are calculated in a little diﬀerent procedure. Although they are
diﬀerent during iterations, it ﬁnally get to be the almost the same number. (But they can beslightly diﬀerent even converged for large systems. But you don’t need to care it so much).
NOTE: ehk:Hohenberg-Kohn energy, ehf: Harris-Faulkner energy.
“grep diﬀe lllmf” shows how the changed of total energy (and charges) during iteration.
diﬀe mean changes of energy with previous iteration, (q) is for electron density diﬀerence as
well. See also save.* ﬁle, which only show ehk and ehk obtained by each iteration.
“ grep gap llmf” shows how the band gap changes (in the usual setting), two same numbers
per iteration are shown now.
Thus we do have ground state energy. Although output of lmf is long, most of all are tomonitor convergence (we will shrink it). As long as it converged well, you don’t need to look
into it in detail. Eigenvalues are shown asbndfp: kpt 1 of 4, k= 0.00000 0.00000 0.00000
ndimh = 122-1.2755 -1.2008 -1.2008 -0.2052 -0.2052 -0.2052 -0.0766 -0.0766 -0.0766-0.0174 -0.0174 -0.0174 0.1094 0.1095 0.1095 0.2864 0.2864 0.4170
0.4170 0.4736 0.6445 0.6445 0.6445
This is at k= 0.00000 0.00000 0.00000 . (because of historical reason, two same bndfp: areshown in each iteration; two band path method). “lmf cu| grep -A6 BZWTS” shows the Fermi
energy (for insulator, we see band gap). Deep levels which gives little dependence on k arecore like levels. These are in Ry; zero level is not so meaningful (for convenience, it is simply
determined from the potential at MT boundaries).
rst.* contains is the main output which contains electron density. mix.* is a mixing ﬁle(which keeps iteration history). When you restart lmf again, it read rst.cu and mix.cu. If you
start from lmfa result, please remove them. We can do parallel calculation with lmf-MPIK,
we can invoke it with mpirun -np 8 lmf-MPIK cu. It should give the same answer.

19

4.6

DOS, PDOS plot

We already have script to plot dos, band, and pdos from the result of lmf self-consistentcalculations. We have scripts
job_tdos, job_band, job_pdos. Look into these scripts, and then you see how to plot them.
For total DOS plot, it is better to check ctrl ﬁle; BZ_TETRA=1(this is default; thus makesure that BZ_TETRA do not exist or BZ_TETRA=1). In addition, it might be better to enlarge
number of k point NKABC in ctrl ﬁle to have smooth curve. Then we do
job_tdos cuThis shows total DOS as The range of DOS and division of total DOS is given by DOSMAX

Figure 1: DOS(Cu)and NPTS in ctrl. Edit tdos.cu.glt for gnulot for your presentation. Please look into job tdos
ﬁle in your bin directory. It is a small script.
For PDOS plot,
job_pdos cu
It shows ﬁgures (number of ﬁgures are number of atoms in the cell) in gnuplot (they are
written in the same position on X-window; move top one a little). The command job_pdos isa little time-consuming because we use no symmetry to distinguish all lm channels. (PDOS
is not yet implemented for SO=1 case; spin-orbit coupling LṠ is added.) We can edit scriptof gnuplot (pdos.site*.*.glt) for your purpose. To plot again, run
gnuplot -persist pdos.site001.cu.glt
In principle, meanings of all data ﬁles are shown (see at the bottom of console output about
lm ordering in a line), thus not so diﬃcult to rewrite *.glt. For example, to plot eg and t2g
separately. (NOTE: site id is shown by lmchk).
WARNING: Usually lmf and so on recognize options such as -v option. For example, ’lmf gaas -vnspin=2’ or ’lmf gaas -vso=1’. This option changes values of variables deﬁned in% const section. This is recorded in save.* ﬁle, and also shown at the top of console output.
However, job tdos and so on, do not yet accept these options. Thus we need to modify ctrl
ﬁle without using -v option. Or you need to write these option to these command by hand
(we will ﬁx this problem in future.)

4.7

Band plot

For band plot, we have to set symmetry lines along which we plot eigenvalues. Collections
syml.* are in ecalj/MATERIALS/. (On Jan2017, we prepared automatic generator of syml.*
See Sec.4.7.1). Choose and modify one of them and rename it. I will gather other samples
soon. ’BZ wikipedia’ or something else will help you to interpret it.
To do band plot, we need syml.cu in your directory.
$ cp ~/ecalj/MATERIALS/Cu/syml.cu .
Then check syml.cu; it is

20

21
21
0

.5 .5 .5
0 0 0
L Gamma
0 0 0
1 0 0
Gamma X
(this is the terminator line)

We supply ten data for each lines (integer, three real, three real, two words). First line means,
we calculate eigenvalues for k points from k=(0.5,0.5,0.5) to k=(0,0,0). ”L Gamma” are namesof two end points (.5 .5 .5) and (0 0 0) in this case. These names are used in a gnuplot scriptfor band plot (bndplot.isp*.glt). Second line means, we calculate eigenvalues for k points fromk=(0,0,0) to k=(1,0,0). 3rd line means calculation just stop here. Units of k are in 2π/ALAT(or 2π/(ALAT+DALAT) if DALAT exist.). A line starting from ’#’ is neglected (comment line).
To do band plot, run$ job_band cu. This is for both nspin=1 and nspin=2 (These scripts try to determine the Fermi energy ﬁrst.
You may skip it in cases (but need to change the script)).

Figure 2: band plot(Cu)4.7.1

syml genearator

Jan2017, we added directory ecalj/GetSyml, which contains automatic generator of symmetryline ﬁle syml.* from ctrl.*). Simultaneously, we plot a picture of BZ showing the symmetry lines by matplotlib. See ecalj/GetSyml/README. We use open softwares spglib andseekpath (these should be installed separetely). We have to add acknoledgement to them (seeREADME) )when we use this generator.

4.8

Useful samples: ecalj/MATERIAL/

Not only ecalj/lm7K/TestSamples (some of them are by older version), We have a material
database in ecalj/MATERIALS/. Move to the directory, and type$ ./job_materials.py
Then it shows a help. You see...
=== Materials in Materials.ctrls.database are:===2hSiC 3cSiC 4hSiC AlAs AlN AlNzb AlP AlSb Bi2Te3
CdO CdS CdSe CdTe Ce Cu Fe GaAs GaAs_so GaN
GaNzb GaP GaSb Ge HfO2 HgO HgS HgSe HgTe InAs
InN InNzb InP InSb LaGaO3 Li MgO MgS MgSe MgTe
Ni NiO PbS PbTe Si SiO2c Sn SrTiO3 SrVO3 YMn2
ZnO ZnS ZnSe ZnTe ZrO2 wCdS wZnS
...

21

C

. For these simple materials (now 57 materials), input ﬁles can be generated, and run
them automatically by a command ./job_materials.py below. The ctrls are stored inecalj/MATERIALS/Materials.ctrls.database in a compact manner (in addition, options
passed to ctrlgenM1.py and options to lmf-MPIK are included). See ecalj/MATERIALS/README
about how to add new material to it; it is not diﬃcult. The command ./job_materials.py
gives ctrls.* for these materials from descriptions in the Materials.ctrls.database. And thenit generates ctrl ﬁle by calling ctrlgenM1.py internally, and run lmfa lmf-MPIK successively(when no –noexec).
Try ./job_materials.py Fe --noexec. (not fe but Fe as it shown above) at ecalj/MATERIALS/.
Then it makes a directory Fe/ and set ctrl.fe (also ctrls.fe) in the directory. Without ’–noexec’,
it does calculation for Fe successively. As for NiO and Fe, we see that ./job_materials.py
gives SPEC ATOM MMOM in generated ctrls and ctrl ﬁles. (Look into ctrls.fe; we needSPEC section when we add MMOM.)
Try job_materials.py GaAs Si.
Then directories GaAs/ and Si/ are generated. See save.* ﬁles containing total energies
iteration by iteration. Starting from ctrl.* in these directory, the command perform DFT
calculations (Console output is stored in llmf, save.* gives total energies. rst.* contains
self-consistent density, from which we can calculate energy bands and so on).
“./job_materials --all --noexec” generates ctrls and ctrl ﬁles of these materials. “./job_materials --all”
do self-consistent LDA calculations for materials (it takes an hour or more. To change thenumber of cores for lmf-MPIK, set option -np (number of core). See help of ./job_materials(type this without arguments).
To make band plot and so on for Fe, follow instructions already explained.
$ ./job_materials.py Fe (and need to type return)(If you like start over, remove Fe/ under it first).
$ cd Fe$ ./job_materials.py fe(but it might be better to do --noexec, and observe Fe/ctrls.fe andFe/ctrl.fe first. grep conf llmf shows the initial electron distribution).
$ cat save.fe (this shows total energies of each iteration. ’c ’ at
the first column gives converged result. ’h ’ is from atm file.)
If it does not ends with ’c ...’ line, something strange
occurs. see llmf (console out put of lmf is saved to llmf).
$ cp ../syml.fe .
$ job_band fe -np 4(As I said, this shell script do not yet accept
options to lmf. Look into the script).
(This calculate fermi energy first for safe; it takessome time)$ job_tdos fe$ job_pdos fe (as I said, this supress space-group symmetry, thus time consuming).
At the end of job_pdos, we show a help which pdos data is where (In pdos ﬁle, we have 26
numbers a line; ﬁrst is energy, 2-26 are pdos for s,p,d,f,g; which is which are shown in the help).
See joblmf ﬁle also (it contains options to invoke lmf. This is shown in save.*. In principle,
options in joblmf should be passed to band plot and so on. But not yet implemented (it is not
so diﬃcult; I have to do it).
After doing ./job_materials foobar, you may like move it back to original... In such acase, git works. At ecalj/, do$ mv MATERIALS MATERIALS.bk$ git checkout MATERIALSThen you can see MATERIALS/ is moved back to just downloaded one.

4.9

How to add spin-orbit coupling?

In the LDA level calculations, we can use

22

SO=0 (no SO), SO=1(LdotS), or SO=2(LzSz) schemes.
However, SO=1 is not the non-collinear method(z-direction is assumed to be spin direction; in cases
we may need to primitive cell so that the z-directionis the spin direction).
Within LDA, we have two possible ways.
(A) Do LDA and/or QSGW with SO=0 first.
Then apply the spin-orbit coupling by perturbation.
or (B) do self-consitent LDA calculations with SO=1 (or =2).
For semiconductors, we think you observe little differences.
Let me explain the case (A).
After converged with nspin=1, create new directory and copy
ctrl.gas, rst.gas, sigm.gas, QGpsi,
to it. Then we set
nspin=2
METAL=3 (usually this is default)
SO=1 (this is ldots calculation off-diagonal elements included).
Q=band (we do not change potential.)
in ctrl.gas.
Then run>lmf gas >& llmf_SO
You can see "band gap with SO" by> grep gap llmf_SO.
Then you can see two same lines.
VBmax = 0.101949 CBmin = 0.236351 gap = 0.134402 Ry = 1.82786 eV
VBmax = 0.101949 CBmin = 0.236351 gap = 0.134402 Ry = 1.82786 eV(two lines per iteration is shown in metel mode).
This is the band gap with SO as a first-order perturbation
on top of the "QSGW without SO". When you use ctrl fil e generated by
ctrlgenM1.py. You can do the above procedure with>lmf --rs=1,0 gas -vnit=1 -vso=1 -vnspin=2 -vmetal=3 --quit=band(--rs=1,0 read rst.gas but not write rst.gas. Run lmf --help.
The switch -v (-vso=1 in this case) replaces so=0 with so=1.
This is recorded in save.gas file).
For band plot, you can use the same procedure
for the case without SO. (Look into the shell script job_band.
You have to modify it so that’--rs=1,0 gas -vnit=1 -vso=1 -vnspin=2 -vmetal=3 --quit=band’is added to arguments for >lmf --band:syml ...).
(--quit=band is not necessary if we like to renew eigenfunctions
self-consistent. Anyway we expect little differences.)---\noindent QSGW with SO:\\
For given sigm file, it is possible to attain self-consisency with SO=1(or=2) with keeping sigm (then we
However, this imply that Vxc is fixed in QSGW,
it is not necessary better than the above procedure.

4.10

PROCASR (VASP format) generator

PROCASR mode for lmf (not yet in lmf-MPIK) Band weight decomposition. (Size of circles
show the size of components. Superposed on band plot). See ecalj/MATERIALS/MgO_PROCAR/README,

23

May19.2014

4.11

efermi.lmf

This is generated by lmf(lmf-MPI), which is used in the job band.

4.12

Effective mass calculation

New band plot and effective mass calculation (curvature):(see sample at ecalj/MATERIALS/mass_fit_test/ and its README).
-----------We now read syml.* in lm7K/fp/bndfp.F. Thus job_band is changed(nspin=1 or nspin=2 is automatically choosed by job_band command).
We do not use plbnd anymore. If necessary, you can modify"writeband" subrouitne in lm7K/fp/bndfp.F by yourself.
Follow steps to get effective mass...
1.
New syml read label of k point. It is shown by the gnuscript
file "bandplot.glt". Type "gnuplot -p bandplot.isp1.glt" and so on,when you like to remake band plot.
2.
New syml allow a special input suitable to determine effective mass
for semiconductor. An example of new syml is (this is a case of GaAs)------ syml example start -----------------#NOTE: ndiv2, ninit2 nend2 etolv etolc are for mass mode.
#ndiv qleft(1:3) qright(1:3) llabel rlabel ndiv2 ninit2 nend2 etolv(Ry) etolc(Ry)
5 0 0 0.5 .5 .5
Gamma L 257 1
32
0.1 0.01
5 0 0 0
1. 0 0Gamma X
5 0 0 0.75 .75 0
Gamma K 257 1
32
0.1 0.01------- end ---------------------New feature is start from the next columns of symmetry points labels.
In this example, 257 means a line connecting 0 0 0 (Gamma) and .5 .5 .5 (L) is divided
into 256. Then we only calculate from the 1st point to the 32th point
among 257 points.
Then we see not only band*.spin*, but also
we have files such as "Band001Syml001Spin1.dat"
files which contains eigenvalues for selected bands, that is,
we make a table of detailed plot of energy bands
whose energy E at Gamma (exactly speaking, at
left-end point), is evaltop-etolv (Ry) < E < econtop + etolc(Ry).
(As its head line shows, Band*Syml*.dat file contains data
isyml,iq, ib,isp, QPE-EF, QPE-QPE(start), |q|, mass=2*2*(QPE-QPE(start))/|q|**2).
here QPE(start) is QPE at the left-end of sym line.)
This sample is only for Gamma-L and Gamma-K.
Let us use this example for GaA and run job_band.
(e.g., job_band gaas -np 2 -vnspin=2 -vso=1;
note that -vfoobar=xxx replace value of foobar with xxx in ctrl.gaas)
Then you can see not only bandplot.isp1.glt.
Then you can see usual energy bands plot in addition to the dense energy bands
superposed on it. The dense energy bands are controlled by
the part "257 1
32
0.1 0.01". In the current version, set all lines24

of etolv and etolc are the same.
There are useful data, not only the eigenvalues in "Band001Syml*Spin*.mass"files. But, anyway, we have to determine effective mass from the eigenvalues
on the dense symmetry mesh.

3.
Least square fit by gnuplot.
We have an example is at ~/ecalj/MATERIALS/mass_fit_test
Look into ./job and run it.
See README in it.
---For your convenience, we have dE/dk in the bnd*.spin* files.
This is useful to determine the Fermi surface.
See the efermi.lmf to read the Fermi energy.

4.13

density and |eigenfunction|2 plot

We can generate density and>mpirun -np 24

P

n |Ψkn (r)|2by --density mode (after 18mar2016) Usages are

lmf-MPIK si --density

This is for density. We have smrho.xsf, and rho1mt*(1st components) and rho2mt*(2ndcomponents).
>mpirun -np 24

lmf-MPIK si --density,iq=1,ib=5

iq=1 means Gamma point. See index shown by the console output of lmf. The index ib is
the band index.
>mpirun -np 24

lmf-MPIK si --density,iq=1,ib=5,6,7

In this case we have sum for ib=5,6,7.
Use xcrysden as>xcrysden --xsf smrho.xsf
Then open menu ’Tools→DetaGrid’. Then plot isosurface and/or planes. For example, youneed to set isovalue, and/or turn on ’display color plane’ buttons.

4.14

LDA+U, partially occupied core-holeSeefile:///home/takao/ecalj/Document/BACKUP/MarksOriginalDoc/fp.html#ldaplusu
We need to add lines such as
IDU= 0 0 2 2 UH= 0 0 0.1 0.632 JH= 0 0 0 0.055
for each SPEC_ATOM, And initial occnum.foobar file.
An example is in ecalj/MATERIALS/GdNldau/
When you just like to generate initial condition for gwsc,
you have to remove (or comment out) IDU before the 1st iteration
with sigm file, because sigm may already can contain LDA+U kind of effect.
(Thus you may need to modify gwsc or stop it at the 0th iteration,and then remove IDU...)
Caution:We need the initial condition file such as occnum.gdn for LDA+U.
(you may need to set "% real" at the begninig of the file).

25

Note that definition of spherical harmonics is in ecaljmanual.pdf.
(real harmonics is usual ones used in jobpdos).
Look for the keyword ldau in fp/lmfp.F -> sudmtu.F which read occnum.gdn.

5

How to run QSGW calculation?

In the QSGW, we calculate a non-local exchange-correlation potential V xc (r, r′ ), by a procedure
LDA
of GW calculation (very time-consuming part). Then diﬀerence V xc (r, r′ )V −Vxc(r)δ(r−r′ )is stored into sigm.* ﬁle. The potential ﬁle sigm is a key to perform QSGW calculations as
LDAseen in Fig.6. The sigm contains static non-local potential ΣQSGW − Vxc.
Then, we again do one-body calculation by lmf (or lmf-MPIK) where we add this sigm
to one-body potential ; when we run lmf or lmf-MPIK(k-parallel mpi version), sigm is readand added to the one-body potential if we have HAM RDSIG=12 in ctrl.*. Thus this means
LDA
that we replace Vxc(r)δ(r − r′ ) with V xc (r, r′ ).
This iteration cycle is performed by a script “gwsc” as we explain later on. (In the default
setting of ctrl.* ﬁle, lmf try to read sigm.* ﬁle as long as it exists. If not, do lmf or lmf-MPIK
calculation.
To start QSGW calculation by gwsc, we need not only ctrl.si, but also another input ﬁleGWinput. Its template GWinput.tmp can be generated by the Step.2 as follows.
As a summary, you have to follow steps below in order to perform QSGW calculation.
Step 1. Perform from the Step.1 thru the Step.4 (up to lmfa) in Sec.4.5. (as same as the case of
LDA/GGA).
You don’t need to perform LDA calculation in advance, since gwsc perform LDA/GGA
calculation at its beginning. (It means that we start from the one-body Hamiltonian H0in LDA/GGA as initial condition.In cases, LDA/GGA give poor initial conditions for
QSGW; in such a case, we may need another trick to prepare starting point.).
[Caution: We have to use the same LMXA (l in the expansion of eigenfunctions in eachMT) for all the MT spheres. (This is due to historical reason; we may need to ﬁx this.)]xxx need to explain LMA xxx
Step 2. Run the script mkGWIN lmf2.
The purpose of this script is to get GWinput.tmp. Other ﬁles generated are not used in
the following stage.
Step 3. Edit GWinput.tmp and save it as GWinput.
GWinput is the input ﬁle describing the computational conditions for GW calculation.
Usually, the default setting gives reasonable results. To reduce computational time,
we may use lcutmx(atom)=2 for oxygen sites (this may be also other small atoms.).
grep lcut GWinput -A1 shows lcutmx for each atomic sites) These step 2. and step 3.
are just only to get GWinput.
Step 4. Run the script gwsc.

5.1

GWinput

In order to perform QSGW, one another input ﬁle GWinput (no extension) is necessary inaddition to ctrl.*. Thus all input ﬁles for QSGW is just two ﬁles, ctrl.* and GWinput. A
template GWinput can be generated by a script mkGWIN_lmf2. You may have to modify it incases for your purpose.
Let us start from ctrls.si;#id = Si%const bohr=0.529177 a= 5.43095/bohr
STRUC
ALAT={a}26PLAT=0 0.5 0.5

0.5 0 0.5

0.5 0.5 0

SITE
ATOM=Si POS=0.0 0.0 0.0
ATOM=Si POS=0.25 0.25 0.25
. Do ctrlgenM1.py si --tratio==1.0 --nk1=6 and copy ctrlgenM1.ctrl.si to ctrl.si. NOTE:
the option --tratio=1.0 means we use touching MT; this can be checked by lmchk si; since
defaults is almost unity (--tratio=0.97), this is irrelevant, just to explain options.
We have to write GWinput. The default is given automatically by a command mkGWIN_lmf2;$ lmfa si (lmfa is needed to do in advance).
$ mkGWIN_lmf2 si......
== Type three integers n1 n2 n3 for Brillowin Zone meshing for GW! ==n1=
Then it pause and ask numbers. You have to type three numbers as 2+ return + 2+return+2
return.
== Type three integers n1 n2 n3 for Brillowin Zone meshing for GW! ==
n1= 2
n2= 2
n3= 2
2 2 2
...(skip)...
OK! GWinput.tmp is generated!
Generated ﬁle is GWinput.tmp; you have to copy it to GWinput.
$ cp GWinput.tmp GWinput
These ’2 2 2’ you typed is reﬂected in a section ’n1n2n3 2 2 2 ’ in GWinput. This means 2x2x2(8 points in 1st BZ). You can edit it, and change it to e.g. ’n1n2n3 4 4 4’ if you like to calculate
self-energy on dense BZ mesh 8x8x8.
The template of GWinput is usually not so bad. But it may give a little expensive setting(or not very good enough in cases).

5.2

Run gwsc script

Let us perform QSGW calculation. For this purpose, we use a script gwsc. We need to do
lmfa in advance. Then do (not need to do lmf);
gwsc (number of iteration+1) -np (number of nodes) (id of ctrl)
If (number of iteration+1)=0, it gives one-shot calculation from LDA. But it is diﬀerent fromthe usual one-shot; since it calculates oﬀ-diagonal elements of self-energy also, we can plot
energy band dispersion. In cases (for usual semiconductors), it can give rather reasonable
results in comparison with experiments from practical point of view.
This is an example of one iteration of QSGW cycle. (now a little diﬀerent but essentiallysimilar)takao@TT4:~/ecalj/test1$ gwsc 0 -np 2 si
gwsc 0 -np 2 si### START gwsc: ITER= 0, MPI size= 2, TARGET= si--- No sigm nor sigm.$TARGET files for starting ------ goto sc calculation with given sigma-vxc --- ix=,0No sigm ---> LDA caculation for eigenfunctions
Start mpirun -np 2 /home/takao/ecalj/TestInstall/bin/lmf-MPIK si > llmf_lda
OK! --> Start echo 0| /home/takao/ecalj/TestInstall/bin/lmfgw si > llmfgw00
OK! --> Start echo 1|/home/takao/ecalj/TestInstall/bin/qg4gw > lqg4gw
OK! --> Start echo 1|mpirun -np 2 /home/takao/ecalj/TestInstall/bin/lmfgw-MPIK si> llmfgw01
OK! --> Start /home/takao/ecalj/TestInstall/bin/lmf2gw >llmf2gw
OK! --> Start echo 0|/home/takao/ecalj/TestInstall/bin/rdata4gw_v2 > lrdata4gw_v2

27

OK! --> Start echo 1| /home/takao/ecalj/TestInstall/bin/heftet > leftet
OK! --> Start echo 1| /home/takao/ecalj/TestInstall/bin/hchknw > lchknw
OK! --> Start echo 3| /home/takao/ecalj/TestInstall/bin/hbasfp0 > lbasC
OK! --> Start echo 3| mpirun -np 2 /home/takao/ecalj/TestInstall/bin/hvccfp0 > lvccC
OK! --> Start echo 3| mpirun -np 2 /home/takao/ecalj/TestInstall/bin/hsfp0_sc > lsxC
OK! --> Start echo 0|/home/takao/ecalj/TestInstall/bin/hbasfp0 > lbas
OK! --> Start echo 0| mpirun -np 2 /home/takao/ecalj/TestInstall/bin/hvccfp0 > lvcc
OK! --> Start echo 1| mpirun -np 2 /home/takao/ecalj/TestInstall/bin/hsfp0_sc > lsx
OK! --> Start echo 11| mpirun -np 2 /home/takao/ecalj/TestInstall/bin/hx0fp0_sc > lx0
OK! --> Start echo 2| mpirun -np 2 /home/takao/ecalj/TestInstall/bin/hsfp0_sc > lsc
OK! --> Start echo 0| /home/takao/ecalj/TestInstall/bin/hqpe_sc > lqpe
OK! --> == 0 iteration over ==OK! --> Start mpirun -np 2 /home/takao/ecalj/TestInstall/bin/lmf-MPIK si > llmf_gwscend.0
OK! ==== All calclation finished for gwsc 0 -np 2 si ====
Here echo (integer) is readin in at the beginning of the code. To see it, please look into gwsc
script (gwsc is at ecalj/fpgw/exec/ and copied to your bin/ by make install2). In anyway, thisconsole output shows calculations ﬁnished normally.
Now we get rst.si and sigm.si ﬁle which contains (static version of) self-energy minims
LDA
Vxc. What we did is the one-shot GW from LDA result; but note that we calculate not only
diagonal elements but also oﬀ-diagonal elements.
We can write energy dispersion (band plot) in the same manner in LDA. To do it, we need
rst.si, sigm.si, ctrl.si, QGpsi. (but QGpsi is quickly reproduced). After you have syml.si (e.g.
in ecalj/MATERIALS/), Do$ job_band si

Figure 3: Si, one-shot GW with off-diagonal elements
You can observe large band gap as shown in the Fig.5.2. (To see it again, gnuplot bnds.gnu.si -p.
All plots are in gnuplot, thus it is easy to replot it as you like).
We have QPU ﬁle (and also QPD for spin=2), which contains content of the diagonal partof self-energy. It will be explained elsewhere.
You can make total DOS and PDOS plot by$ job_tdos si$ job_pdos si
CAUTION:pdos plot is not allowed for so=1. (even tdos–¿ ask to t.kotani.)
To get ﬁnal QSGW results, we have to repeat iteration until eigenvalues are converged.
Note that total energy shown by console output llmf (and also shown in save ﬁle) is not so28meaningful in the QSGW; we just take it as an indicator to check convergence. Let us repeat
5 iteration more. ”-np 2” means one core to use.
$ gwsc 5 -np 2 si### START gwsc: ITER= 5, MPI size= 2, TARGET= si--- sigm is used. sigm.$TARGET is softlink to it ------ goto sc calculation with given sigma-vxc --- ix=,0
we have sigm already, skip iter=0---- goto sc calculation with given sigma-vxc --- ix=,1...(skeip here) ...
OK! --> == 5 iteration over ==OK! --> Start mpirun -np 2 /home/takao/ecalj/TestInstall/bin/lmf-MPIK
OK! ==== All calclation finished for gwsc 0 -np 2 si ====

si > llmf_gwscend.0

Note that we do skip 0th iteration (it is for one-shot from LDA) since we start from rst.si andsigm.si given by one-shot LDA. Thus we do just ﬁve iterations. Information of eigenvalues are
in QPU.(number)run ﬁles. (for magnetic systems with nspin=2), wee have QPD.(number)runalso). Check it by ls;$ ls QPU.*run
QPU0.run QPU.1run

QPU.2run

QPU.3run

QPU.4run

QPU.5run(These are overwritten when we again repeat gwsc; be careful.) Note that QPU0.run was old
one when you did 1-shot GW from LDA at the beginning. In anyway *.0run are confusingﬁles; remove them).
In order to check convergence calculations going well during iteration, do$ grep gap llmf*
This shows how band gap changes in llmf.*run ﬁles. In metal cases, we need to compare QPU
ﬁle, magnetic moment or grep ’[xc] save.*; this shows end of lmf iteration. Energy is not
so meaningful but can be indicator to convergence.
Let us check convergence of the QSGW calculations. For this purpose, it is convenient totake a diﬀerence of QPU(QPD) ﬁles by a script dqpu. These ﬁles are human readable. To
compare QPU4.run and QPU5.run, do$ dqpu QPU.3run QPU.4run
Then we see a list of numbers (these are the diﬀerences of values in QPU ﬁles). Then it showsat the bottom as
Error! Difference>2e-2 between:: sum(abs(QPU-QPD))= 0.05736

QPU.4runand

QPU.5runbut you don’t need to care it so much. You rather need to check the diﬀerence of values. I
can say most of all diﬀerence (especially around the Fermi energy are ) are almost 0.00eV or
0.01eV, we can judge QPEs are converged. If not converged well, you may need to repeat gwsc
again. (when the size of two QPU ﬁles are diﬀerent, dqpu stops.)

5.3

Spectrum function: How to calculate hqn|Σ(ω)|qni

We have an example at ecalj/MATERIALS/SiSigma/, where you can just type job. It calls ashell script gwsigma, which is just a modiﬁcation of gwsc for spectrum function plotting. If
you have sigm.*, it will automatically read it as in the case of gwsc.
By the script gwsigma, we calculate the diagonal elements iψ(q, n)|σc (ω)|ψ(q, n)h. Thus weneed to set q and band index n for which we calcualte. In addition, we need to set resolutionof ω.
1. Set <QPNT> section. This section is to set the q point, and band index for which we
calculate the self energy. In addition, energy mesh for plotting is set.

29

Figure 4: band plot(Si, QSGW one-shot test) and band(Si) (GGA)

Figure 5: band(GaAs), QSGW (test case)(A) Set q point set

If you set---------*** all q -->1, otherwise 0; up only -->1, otherwise 010---------You will have self-energy for all irreducible k points. This may be needed for A(omega).
or
You have to set all q points as---------*** q-points, which shoud be in qbz.,See KPNTin1BZ.
3<--- number of readin q point
1
0.0000000000000000
0.0000000000000000
0.0000000000000000 <--1st number is irrelevan
2-0.5000000000000000
0.5000000000000000
0.5000000000000000
30.0000000000000000
0.0000000000000000
1.0000000000000000---------To know allowed q points on regular mesh point, run the command "mkGWIN_lmf2", thensupply n1,n2,n3. The templete of GWinput.tmp contains all possible q points. Edit it.
NOTE:Anyq option can allow you to specify any q points by shifted mesh technique.
(if necessary, but only for some special purpose).
(B) Band index set

30

It is specified by the section---------*** no. states and band index for calculation.
2
4 5---------means the self-energy for band index 4 and 5. Just two bands.
If you like to plot self energy from 1 through 8, use*** no. states and band index for calculation.
8
1 2 3 4 5 6 7 8---------If you need 17 bands for example, it should be 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
in addition to 17 (number of bands at the first line).

(C) energy mesh set
At the bottom of <QPNT> section, we have---------------*****
0.01 2.0---------------Two real number should be supplied.
These are dwplot and omegamaxin, read in hsfp0.m.F by a line
read (ifqpnt,*,err=2038,end=2038) dwplot,omegamaxin
dwplot (=0.01 Ry) is mesh for self energy.
omegamaxin=(2.0 Ry) means the range "-2 Ry to 2 Ry" for self-energy plot.

Note that imaginary part of Sigma is given as the comvolution of ImW(omega) and the pole of Gree(esmr in GWinput gives energy smearing of the pole). Resolution for Im W (near omega=0) is by dw
I think that the reolution of self-energy is ~ 0.05 eV in the default setting.
This is because {\tt dw} \sim {\tt esmr} \sim 0.05 eV.

2

Run gwsigma. This will runecho 4| mpirun -np 24 hsfp0,
after dielectric funcition is calculated.
Then we have SEComg.UP (DN) files, Look for file handle, ifoutsec,
for the file in fpgw/main/hsfp0.m.F to see format for the file.
(not hsfp0.sc.m.F but hsfp0.m.F). Search a line
open(ifoutsec,file=’SEComg’//sss) (around hsfp0.m.F L1052)
You can find that we use folloing lines to plot SEComg.*.
---------------write(ifoutsec,"(4i5,3f10.6,3x,f10.6,2x,f16.8,x,3f16.8)")&
iw,itq(i),ip,is, q(1:3,ip), eqx(i,ip,is),&(omega(i,iw)-ef)*rydberg(), hartree*zsec(iw,i,ip) !,sumimg
This means we use energy in eV.
iw:omega index
itq(iq): band index specified by <QPNT>
ip:
k point index specified by <QPNT>
is:spin indexq:q vector (cartesian in 2pi/alat)eqx:eigenvalue in eV. (I think relative to the Fermi energy)(omega(i,iw)-ef)*rydberg(): omega relative to the Fermi energy
hartree*zsec(iw,i,ip):
Self energy. real and imaginary part.(complex, two values)
You can only repeat echo 4| mpirun -np 24 hsfp0

31

when you change setting in <QPNT> section.
3 Example. There is an example MATERIALS/SiSigma/
plot ’SEComg.UP’ u ($9):($10) w l,’’ u ($9):($11) w l
can give a plot for Re (Sigma_c(omega)) and Im(Sigma_c).
9th: energy in eV
10th: real part
11th: imag part

(omega(i,iw)-ef)*rydberg()
Re hartree*zsec(iw,i,ip)
Im hartree*zsec(iw,i,ip)

4 To get integrated spectrum function (DOS), we need to superpose all the spectrum
function (All q points and all band index). Be careful about the degeneracy (multiplicity)
for each q points. You have to build it from SEComg ﬁle. To know the multiplicity, searchfollowing lines ofkeyword Multiplicity in the console output of qg4gw (lqg4gw).
Anyway, consider about “is it worth to do?” To conﬁrm your result, use sum rule (sumof spectrum weight). And pay attention to the relation between real and imag parts(Hilbert transformation).

32

6

gwsc script to perform QSGW

6.1

outputs of gwsc

Figure 6: Shell script gwsc to perform QSGWcontains an iteration loop to make sigm (and eigenvalues, eigenfunctions) self-consistent. The iteration loop is written in ecalj/fpgw/exec/gwsc (abash script). Exactly speaking, we have to pass all the required data (not only eigenvalues andeigenfunctions, but also crystal structure, q+G vectors, symmetry information and so on) to GW
LDA
part. The purpose of the GW part is to calculate ΣQSGW − Vxc.
When gwsc runs normally, it gives console output as follows. This is a case for ctrl.gaas
for>gwsc 10 -np 24 gaas
. Without arguments, typing gwsc shows usage as>An example of usage: gwsc 5 -np 4 si, where 5 means 5+1 iterations. We recommend you do look into the script gwsc. It uses run arg, which is a special subroutineof bash; but not so diﬃcult to understand it. (In the followings, we assume /home/binx/ isyour bin directory at which we have all binaries for ecalj.)### START gwsc: ITER= 10, MPI size= 24, TARGET= gaas--- No sigm nor sigm.$TARGET files for starting ------ goto sc calculation for given sigma-vxc --- ix=,0No sigm ---> LDA caculation for eigenfunctions
OK! --> Start echo --- | mpirun -np 24 /home/binx/lmf-MPIK gaas > llmf_lda
OK! --> Start echo 0 | /home/binx/lmfgw gaas > llmfgw00
OK! --> Start echo 1 | /home/binx/qg4gw > lqg4gw
OK! --> Start echo 1 | mpirun -np 24 /home/binx/lmfgw-MPIK gaas > llmfgw01
OK! --> Start echo --- | /home/binx/lmf2gw > llmf2gw... (preparation stage ends here; start main stage) ...
OK! --> Start echo 0 | /home/binx/rdata4gw_v2 > lrdata4gw_v2
OK! --> Start echo 1 | /home/binx/heftet > leftet
OK! --> Start echo 1 | /home/binx/hchknw > lchknw
OK! --> Start echo 3 | /home/binx/hbasfp0 > lbasC
OK! --> Start echo 3 | mpirun -np 24 /home/binx/hvccfp0 > lvccC
OK! --> Start echo 3 | mpirun -np 24 /home/binx/hsfp0_sc > lsxC
OK! --> Start echo 0 | /home/binx/hbasfp0 > lbas
OK! --> Start echo 0 | mpirun -np 24 /home/binx/hvccfp0 > lvcc
OK! --> Start echo 1 | mpirun -np 24 /home/binx/hsfp0_sc > lsx
OK! --> Start echo 11 | mpirun -np 24 /home/binx/hx0fp0_sc > lx0
OK! --> Start echo 2 | mpirun -np 24 /home/binx/hsfp0_sc > lsc
OK! --> Start echo 0 | /home/binx/hqpe_sc > lqpe
OK! --> Start echo --- | mpirun -np 24 /home/binx/lmf-MPIK gaas > llmf_oneshot... (this is the end of main stage) ...
== 0 iteration over ==---- goto sc calculation for given sigma-vxc --- ix=,1
OK! --> Start echo --- | mpirun -np 24 /home/binx/lmf-MPIK gaas > llmf
OK! --> Start echo 0 | /home/binx/lmfgw gaas > llmfgw00... (lines here omitted ) ...
OK! --> Start echo 0 | /home/binx/hqpe_sc > lqpe== 1 iteration over ==... (lines here omitted ). ..
== 2 iteration over ==... (this repeat until ITER= 10(the first argument to gwsc) ) ...

33

This shows that gwsc invoke lmf-MPIK,lmfgw,qg4gw,... successively. echo 3|hbasfp0
means running a fortran program hbasfp0 with the argument ’3’ from the standard input(read(*,*) in fortran code). We can divide these successive calls to “preparation stage” and“main stage”. Preparation stage is just to prepare eigenfunctions and so on which are required
for the “main stage” of GW calculation. At the end of “main stage”, we have the potential
ﬁle sigm.
As it shows, console output are going to l* ﬁles.

6.2

Preparation stage of gwsc

At the end of this stage, we get required eigenfunctions, BZmesh data, and so on, which arerequired for the main stage. Note that echo 0 | lmfgw means supply an integer to the fortranprogram lmfgw from standard input by read(*,*).
• lmf-MPIK (k-parallel version of lmf)
This is the one-body calculation for given sigm.gaas. At the beginning, we do not havesigm.gaas. In this case lmf-MPIK just perform LDA/GGA calculation.
• echo 0 | lmfgw:
Get some small information ﬁles to start qg4gw. If you type just lmfgw, and observe
what occurs. It shows a menu and pauses (asking you to supply an integer); then we
supply 0 in this case. (If we do echo 0 | lmfgw, no pause occurs.)• echo 1 |qg4gw : Get k points used in the GW calculations and the corresponding G
vectors. And what is the irreducible k point (See console output of qg4gw. gwsc keeps
it in lqg4gw).
Since we use the oﬀset-Gamma method for BZ integration for G × W , we need shiftedmesh points to calculate W at oﬀset-Gamma points. The q vectors of oﬀset-Gamma
method is in Q0P ﬁle. (If you have two points in Q0P, we see two shifted mesh points
in addition to regular mesh points.) Remember that cutoﬀ of G is given by QpGcut psi
and QpGcut cou in GWinput. (Based on the experiences, we use smaller QpGcut cou to
reduce computational time. Explained in other section xxx).
LDA• echo 1 | lmfgw-MPIK : Calculate the LDA eigenfunctions, eigenvalues, and hψ|Vxc|ψiat the irreducible k points (shown at the bottom of output lqg4gw of qg4gw.)• lmf2gw: store these data into DATA4GW V2 and CphiGeig, whose I/O is controlled bya key subroutine gwinput v2.f.

6.3

Main stage of gwsc

We can start the main stage of GW calculation from these ﬁles;GWinput DATA4GW_V2 CphiGeig QGpsi QGcou Q0P QIBZ SYMOPS BZDATA HAMindex CLASS;
these ﬁles contains eigenfunctions and so on in the manner of Eq.(17) of [2], eigenvalues andother required information. This is the starting point of the GW calculation.
• GWinput : computational conditions.
• DATA4GW V2 : Crystal structures and so.
• CphiGeig : Eigenvalues and Eigenfunctions• QGpsi : q and G vector for the eigenfunctions(q means k in the previous section),• QGcou : q and G vector for the Coulomb matrix• Q0P : q points near q=0 instead of q=0 (oﬀset Gamma points)• QIBZ : irreducible q points (This is also contained in BZDATA).
• SYMOPS : point group operation• BZDATA : q points date (and tetrahedron weights if necessary) for BZ integrals.
• HAMindex :Hamiltonian index, all required complex index for Hamiltonian of PMT

method. (See the top of subroutine write hamindex in lm7K/subs/m hamindex.F. This is also in

34

fpgw/gwsrc/m hamindex.F. Identical ﬁles are in two diﬀerent directory— it should be avoided
in future.)• CLASS : class information or atomic sites (equivalent sites).
With these ﬁles, we do the main stage as• rdata4gw V2: Read DATA4GW V2 and so on, and decompose it into ﬁles required in
the followings. And calculate PPOVL* ﬁles (overlap matrix of interstitial plane waves.
Because of technical reasons some diﬀerent types of PPOVL* with q-point index).
• heftet : Get the Fermi energy EFERMI by tetrahedron method. It is used in hx0fp0.
• hchknw : stores the number of required ω points along real-axis into NW.
(NW is not essentially used, but is supposed to exist in the followings.)• echo 3|hbasfp0: gives the product basis for Core exchange.
• echo 0|hvccfp0: gives the Coulomb matrix for the Core exchange.
• echo 3|hsfp0 sc : gives the Core exchange part of the self-energy. 12• echo 0|hbasfp0: gives the product basis.
• echo 0|hvccfp0: gives the Coulomb matrix v.
• echo 1|hsfp0 sc : gives the exchange part of the self-energy.
• echo 11|hx0fp0 sc : gives the correlated part of the screened Coulomb interaction
W − v.
• echo 2|hsfp0 sc : gives the correlated part of the self-energy.
• echo 0|hqpe sc: gather data and write down ﬁnal results into sigm, QPU, TOTE.UP
ﬁles.

6.4

Other functions (or scripts)

In addition to gw lmfh, there are some other additional scripts and functions.
• gw lmfh : The one-shot GW calculation. Lifetime(impact ionization rate) of QPs.
• gwsc : QSGW calculation explained here
• epsPP lmfh, eps lmfh : Dielectric function without or with local-ﬁeld eﬀects.
• run-mode 4 of hsfp0: to plot the spectrum function Σ(ω).(need to be ﬁxed again probably).
• epsPP lmfh chipm : non-interacting spin susceptibility. One-degree of freedom like
Rigid moment approx. After it ends, you need to do calj_nlfc_metal and/or calj_summary_mat
to get the full spin susceptibility.
• genMLWF : Wannier funciton and its matrix elements of the Screened Coulomb interaction.

12 Correlation part due to cores is neglected. In future, we will switch to a version without PB for core part to
reduce computational time and for numerical accuracy.

35

7

Usage problems, QandA error messages.

1.Bandplot for FSMOMMETHOD/=0
Even when you use FSMMOMMETHOD/=0 in GWinput for gwsc,
you need to set FSMOMMETHOD=0 (or comment it out) when you run job_band_nspin2.
[If you run job_band_nspin2 with FSMOMMETHOD/=0, it make a shift(adding bias magnetic field).]
2.Note that ctrlgenM1.py automatically set this for --systype=molecule.
Then we have
TETRA=0
N=-1#Negative is the Fermi distribution function W= gives temperature.

W=0.001 #W=0.001 corresponds to T=157K as shown in console
In addiiton, FSMOM (n_up-n_down) is needed (FSMOMMETHOD=1)if wehave magnetic moment.
3. core>evalence message.
Ecore is grater than Evalence.
For save, we do not allow this.
Complare ECORE file and valence levels, shown in log file orconsole output.
4. If you see a error message from lmf (e.g., internally called in the gwsc script).;
Exit -1 rdsigm: Bloch sum deviates more than allowed tolerance (tol=5e-6)
You have to enlarge RSRNG so that lmf finsh normally.
See FAQ at the bottom of http://titus.phy.qub.ac.uk/packages/LMTO/gw.html
E.g. try RSRNGE=8 for Ni...
5. Back ground charge and fractional Z.
You can use fractional numbers for ATOM_Z, and also can set
valence charge by BZ_ZBAK (I removed BZ_VAL).
You see console out put, e.g,"Charges:valence19.80000

hom background.20000nucleii-28.00000

deviation from neutrality:

cores

8.00000

0.00000. This is a case with BZ_ZBAK=.2.
NOTE: at the first iteration, Charges: shows such as
Charges:valencehom background8.00000
0.12300cores

20.00000

nucleii-28.00000

deviation from neutrality: 0.12300

because of the initial condition by superposition of atoms. It show
deviation seems nonzero. But charge should be conserved from thenext iteration.
6. Not converged in metal. --->mixing may help
For example, if you try metal such as Bi2Sr2CuO6, it may fail at LDA/GGA level.
Then use ITER MIX=A2,b=.2. or something (.2 means it only mix 20% of output to givenew input for next iteration). Then I see convergence. (b is the
mixing parameter.

36

7. Use PZ or not.
If spillout of core is not so small (more then 0.05 or something.),
it is better to use PZ(lo). Treat the core as valecne.
Bi4d is such a case. Maybe use PZ=0,0,4.9
8. Core treatment
See 10.1103/PhysRevB.76.165106 (Eq.35 and after).
Now I usually not use CORE2 (CORE1 only).
9. ERROR EXIT! rgwina: 2nd wrong l valenceThis may be because you use wrong GWinput.
Back it up. And run mkGWIN_lmf2 (any n1 n2 n3 is fine).
10. Known bug.
Error occurs when system is anisotroic such as CuAlTe2.
Temporary fix is "Add token NPWPAD=100 in HAM category".
(guess of used APW fails (more than expected)).
CuAlTe2,CuGaTe2 cases.
11. Known bug
a little unstable when metal GGA, especially when we have large empty regions.

8

Cautions for usage
1. == not meaningful total energy in QSGW===
Total energy shown in QSGW mode in current version is not meaningful. (just treat as
an indicator to convergence).
2. == Do we use VWN or GGA for QSGW? ===In principle, QSGW results should not depends on VWN or GGA (XCFUN=1 or 103 inctrl). But there is minor dependence, because
1. frozen core density.
2. core eigenfunctions.
3. radial basis functions
4. Slight numerical reason(This is probably because Sigma-interpolation procedure But not exactly ﬁgured out yet→ aﬀect about 0.02eV as for band gap for GaAs. ). In anyway, use VWN (HAM XCFUN=1)
as standard. And such technical things aﬀects, 0.05 eV level of error for band gap.
3. EH and EH2 : For si, if EH and EH2 are the same, the following error occurred.
fexit,fexit2,fexit3 error retval=-1

Exit -1 zhev_tk2: nev /=nevx something wrong.
The large EH, EH2 get to be meaningless. We usually use up to ∼ 2. (If you use verylarge EH such as E∼ 10, I am not so sure weather it is )4. The options about the product basis within MT. (SeungWoo’s memo)<PRODUCT_BASIS>
tolerance to remove products due to poor linear-independency

37

0.100000D-02 ! =tolopt
When the product basis are made, we may have poorly linear independent basis. For
example, one in the set {f1 , f2 , ..., fn } would be almost give by a linear-combination ofothers. We need to make the linear-independent set. Therefore, after calculating theoverlap matrix hfi |fj i. We do diagonalization, then we remove eigenvectors corresponding to small eigenvalues than tolopt. See the hbasfp0 command in gwsc
lcutmx(atom) = maximum l-cutoff for the product basis.
4

4

4

2

2

4

4

For φ1 × φ2 case, |l1 − l2 | ≤ ltot ≤ |l1 + l2 |. So ‘lcutmx’ changes the maximum cutoﬀ forthe ltot . The order is the same as the order of atoms in the ctrl ﬁle.
atom

l

nnvv

nnc !

1

0

3

3

1

1

3

2

1

2

2

1

‘atom’ means the atom number identiﬁed in the ctrl ﬁle.
‘l’ is the angular momentum quantum number.
‘nnvv’ is the number of radial functions (valence) on the augmentation-waves.
‘nnc’ is the number of radial functions for core.
The latter two ones, ‘nnvv’ and ‘nnc’, will be understood more clearly if you see thefollowing ones.
atom

l

n

occ unocc! Valence(1=yes,0=no)

1

0

1

1

1

! 5S_p

1

0

2

0

0

! 5S_d

1

0

3

1

1

! 4S_l

1

1

1

1

1

! 5p_p

1

1

2

0

0

! 5p_d

1

1

3

1

1

! 4p_l-----Above options are about the product basis set within MT (Valence).
‘atom’ and ‘l’ are explained above. ‘nnvv’ for ‘atom = 1 and l = 0’ was ‘3’ so this case
we have 3 basis (‘n = 1, 2, 3’)‘n’ is the degree of freedom of the radial function, φ. ‘n = 1’ means φ, ‘n = 2’ means φ̇,and ‘n = 3’ means kind of φ̈, which the dot above the letter represents the diﬀerentiation
with respect to the energy. So ‘n = 1 and 2’ is related to the linearization of the radialfunction and ‘n = 3’ is the local orbital which is restricted in MT. The local orbital canbe modiﬁed changing ‘PZ’ in the ctrl ﬁle. Finally, the number of the basis set which isneeded for expanding eigenfunctions is (l + 1)2 × n.
‘occ’ and ‘unocc‘ mean that we use only ones that checked as ‘1’, in other words weneglects ‘0’ cases for making product basis. Be careful for confusion with name ‘occ’
and ‘unocc’. These don’t mean that occupied or unocc. When making product basis,
M = φ1 × φ2 , ‘occ’ corresponds to φ1 and ‘unocc’ to φ2 . For example,
atom

l

n

1

0

1

occ unocc
1

1! 5S_p! Valence(1=yes,0=no)

2

3

1

0

1

! 4f_p-----38

If the options are like the above, the product basis will be consists of (φ1 = φatom=1,l=0 )×(φ2 = φatom=1,l=0 ), (φatom=1,l=0 × φatom=2,l=3 ). As you can see, (φ1 = φatom=2,l=3 ) isskipped.
In the ctrl ﬁle, ‘EH’ controls the l part. As for ‘EH’, (s, p, d, f) are used but GWinput
ﬁle uses (s, p, d, f, g). ‘EH’ : HEAD part. ‘GWinput’ : contains TAIL part... need more
explanation.
atom

l

n

1

0

1

occ unocc
00

ForX0 ForSxc ! Core (1=yes, 0=no)0

0! 1S -----

1

0

2

0

0

0

0

! 2S

1

0

3

0

0

0

0

! 3S

Above options are about the product basis set within MT (Core).
‘nnc’ for ‘atom = 1 and l = 0’ was ‘3’ so this case we have 3 basis (‘n = 1, 2, 3’)Finally, for the convergence check, we can modify the following three things, (i) tolerance,(ii) lcutmx, and (iii) occ and unoccu.

== one show QSGW (not one-shot GW) ==one-shot QSGW can be useful in cases.
As it contains off-diagonal part, we can resolve band tanglement
problem in Ge (no band gap).
== Restart calculation in lda ==lmf(lmf-MPIK) read rst.* in defaluts.
rst contains electron density.
If rst is already converged, it stops after two iteration.
rst contains atomic positions.
So, in order to read atomic positions change in ctrl,
Use options shown in lmf --help.
== Restart calculation in qsgw ==To remove mixsigm* (mixing for sigm), maybe required.
== iteration check ===First, watch console output of gwsc (do redirect to output file)
Need to check OK! signs arrayed on 1st columns.
gwsc iteration is time cosuming,
So we need to check calculations are normally going on or not.
Memory inefficiency.
Set ’KeepEigen off’ an ’KeepPPOVL’ off.
In fact, out code is still inefficient for memory usage.
grep gap llmf ---> minimum gap at mesh point.
see save.* ,or grep ’[xc] ’ save.*the end of iteration of lmf is shown as x or c.
(if failed, QPU file.

39

dqpu QPU.4run QPU.3run
As for usual semi-conductor, accuracy abou t0.1 eV is limit of current implementation.
Set vwn (xcfun=1) looks better (stable) for GW.
$grep rms lqpe*
shows
...

rmsdel=2.44D-04...

rmsdel=4.91D-03...

rmsdel=2.44D-04...

rmsdel=3.37D-04

If rsmdel is getteing to be smaller, it is on convergence path.
(but in magnetic cases, it may give be too good even not yet going to
be converged..., beccause magnetic energy is so small)
grep diffe llmf---> difference of energies of each iteration.

ehf (harris energy)
ehk (Hohenberg kohn energy)== emax cutoff for APWs. ==We can not use so many APWs in current version,
because of overcompleteness (this is because null vector within MTs),
In anyway, use pwemax=3 as standard (test it with 4 or 5).
To avoid failure of calculation, we may use smaller MT radius foralkali, and alkali-earth elements.
In feature, I think we can introduce pseudopotentials for these atoms only.
== Check Used MTO
Near beginig of console output, what MTO you use is shown as: (GaAs case).
sugcut:make orbital-dependent reciprocal vector cutoffs for tol= 1.00E-06

spec

l

rsm

eh

gmax

last term

cutoff

Ga

0*

1.13

-1.00

6.579

1.19E-06

1459

Ga

1*

1.13

-1.00

7.028

1.26E-06

1807

Ga

2*

1.13

-1.00

7.475

1.09E-06

2109

Ga

3

1.13

-1.00

7.920

1.06E-06

2637

Ga

0*

1.13

-2.00

6.579

1.19E-06

1459

Ga

1*

1.13

-2.00

7.028

1.26E-06

1807

Ga

2

1.13

-2.00

7.475

1.09E-06

2109

As

0*

1.18

-1.00

6.300

2.13E-06

1243

As

1*

1.18

-1.00

6.720

1.26E-06

1471

As

2*

1.18

-1.00

7.140

1.37E-06

1837

As

3

1.18

-1.00

7.558

1.05E-06

2229

As

0*

1.18

-2.00

6.300

2.13E-06

1243

As

1*

1.18

-2.00

6.720

1.26E-06

1471

As

2

1.18

-2.00

7.140

1.37E-06

1837== gwsc cause error stop.
Have you ever changed MTO setting? Consistent with GWinput?

40

== QSGW for Fe.
It is better to use 3p as core. Furthermore, 3d+4d as valence is better.
Thus we need to set PZ=0,3.9,4.5
I also got aware that emax_sigm should be large enough (4$\sim$5 Ry)
to have smooth band dispersion. n1n2n3 can be 10x10x10.
== RSRNGE: enlarge RSRNGE ===Use RSRNGE=10 or so (in cases, RARNGE=20 or more is required),
for large number of k points. Try and enlarge it if it fails with amessage "Exit -1 rdsigm: Bloch sum deviates more than allowed tolerance (tol=5e-6)".
We will have to make it automatic in future.
== Q0P check
In cases, it is better to use Q0Pchoice=2 instead of default Q0Pchoice=1.
(For slabs, Q0Pchoice=2 may be better; need check more. In anyway,
it is problematic to use unbalanced k points for anisotropic cell).
See Copmuter Physics Comm. 176(2007)1-13).
=== When calculation in LDA level fails ===when calculation fails in LDA level.
(1) smaller MT(2) fewer PW. smaller pwemax.
(3) core as semicore.

=== LDA+U ===not yet written...
=== MAE by rotating crystal ===(we have a sample at lm7K/TESTsmaples/MAEtest/, but only in GGA/LDA).
=== spin wave ===J calculation.

====
If not stable convergence in gwsc, try to set
mixbeta 0.5(and/or mixpriorit 3 or something)at the begining of sigma.

=======
cleargw (directory):
This command clean up up intermediate files under (directory).
This recursively into deeper level. Be careful, or edit it.
I use it as ’>cleargw .’41

-----------------------------------Magnetic moment within MTs are shown as-----------------------------------charges:

old

new

smooth

17.240314

17.240740

mmom

0.000024-0.000010

site

1

mmom
site

2

mmom
site

3

mmom
site

4

mmom

6.207135

6.206590

1.062276

1.062991

6.207115

6.206834-1.062323-1.062958

1.172718

1.172918

0.000011-0.000011

1.172718

1.172918

0.000011-0.000011...

<--- here<--- here

In this case, MTsite1 has 1.062991 and MTsite2 has -1.062958.
>grep ’lin mix’ -A30 llmf
can take out this message (if console output is in llmf).

----------------------------ORBITAL MOMENT in pertubation:----------------------------Try>lmf nio --rs=1,0 -vso=1 --quit=band >llmf
After converged, try>grep IORBTM -A20 llmf
Then llmf shows shows orbital moment in first order perturbation.
(Here --rs=1,0 read rst.* file but not change it. See >lmf --help.
--quit=band means quit just after band calculation.)== EPS mode,
Check Im part of chi0 is smoothly damping at high energy (typically
1Ry or larger enengy range). If there is some large Im part remains,something strange (usually due to orthogonality problem of
eigenfunctions when you set low q).
Related source codes are in ecalj/lm7K/ .
A command ecalj/lm7K/ctrlgenM1.py can generate ’standard input file (ctrl file)’
just from a given crystal structure file called as ctrls file.
Binaries are lmf and lmf-MPIK (MPI k-parelell verion).

42

8.1

lmf –help

lmf –help show option of –rs=(ﬁve numbers); this let lmf know how to read atm.* ﬁle whichis the initial atom ﬁle by lmfa.

9

Wannier function

We can generate Wannier functions (maximally localized Wannier Functions or similar) bya script genMLWF. It automatically perform cRPA calculation sucessively. (If it is notnecessary, insert ’exit’ in genMLWF, after it performs lmaxloc2).
Try to run examples in ecalj/MATERIALS/Sample_MLWF/. Read README in it. To run thescript genMLWF, we need to get GWinput by editing GWinput.tmp. (mkGWin_lmf2 contains
default Wannier section). In addition, we have some settings (energy windows and so on).
This is the example of the initial conditions for Cu case. 5 is the number of Wannier function.
The most left one means φ index and the right one of it is φ̇ index. They are written in the@MNLA CPHI ﬁle.
Then we can run genMLWF. After it ﬁnished, we can analyze it results. (if you don’t
need Wannier funciton plot, You can skip a line of wanplot in genMLWF. Then we don’t needto set vis_wan_* options.)

9.1

lwmatK1 and lwmatK2

If you input the following command>grep Wan lwmatK*You will get the following results. (This case : Cu cases)
lwmatK1:

Wannier1124.644475

0.000000 eV

lwmatK1:

Wannier12

24.644576

0.000000 eV

lwmatK1:

Wannier13

25.471361

0.000000 eV

lwmatK1:

Wannier1

4

24.644575

0.000000 eV

lwmatK1:

Wannier1525.470946

0.000000 eV

lwmatK2:

Wannier1

10.000000 eV-21.263759-0.000000 eV

lwmatK2:

Wannier120.000000 eV-21.263839

0.000000 eV

lwmatK2:

Wannier130.000000 eV-21.931033-0.000000 eV

lwmatK2:

Wannier1

40.000000 eV-21.263839-0.000000 eV

lwmatK2:

Wannier150.000000 eV-21.930702-0.000000 eV### Wanneir Branch now under developing (imported from T.Miyake’s Wannier and H.Kino’s).
A. make at ecalj/fpgw/Wannier/ directory, and do make, and make install.
(need to check Makefile first). You first have to install fpgw/exec/ in advance.
B. Samples are at these directories.
MATERIALS/CuMLWFs (small samples),
MATERIALS/CuMLWF/
MATERIALS/CuMLWFs/
MATERIALS/FeMLWF/
MATERIALS/NiOMLWF/
MATERIALS/SrVO3MLWF/
C. With GWinput and ctrl.*, run>genMLWF

43

at these directories.
In GWinput, we supply settings to generate Wannier funcitons. (Sorry,not documentet yet..)D. After genMLWF, do>grep Wan lwmatK*then compare these with Result.grepWanlwmatK
These are onsite effective interactions (diagonal part only shown).
*.xsf are for plotting the Maximally localized Wannier funcitons.
Anyway, documentaion on Wannier is on the way half.
Time consuming part (and also the advantage) is for effective interaction in RPA.
Look into the shell script genMLWF; you can skip last part if you don’t need the effective interaction.

10

ctrl file details

A ctrl ﬁle is usually generated from a ctrls ﬁle by the ctrlgenM1.py (a crystal structure ﬁle is
not “ctrl” but “ctrls”.). It contains self explanation. Here we give complementary explanations
to it. Let us Look into a ctrl ﬁle. This is a head part of ctrl.cu generated by ctrlgenM1.py:### This is generated by ctrlgenM1.py from ctrls### For tokens, See http://titus.phy.qub.ac.uk/packages/LMTO/tokens.html.
### Do lmf --input to see all effective category and token ###### It will be not so difficult to edit ctrlge.py for your purpose ###VERS

LM=7 FP=7

IO

SHOW=T VERBOS=35 TIM=2,2# version check. Fixed.

# SHOW=T shows readin data (and default setting at the begining of
console output)# It is useful to check ctrl is read in correctly or not(equivalent with --show option).
# larger VERBOSE gives more detailed console output.
SYMGRP find# ’find’ evaluate space-group symmetry automatically.
# Usually ’find is OK’, but lmf may use lower symmetry...
Note that # means comment lines. We can also use lines starting from % const ... to deﬁnevariables and set constant.
We see “categories” such as VERS, IO, and so on. The beginning of categories are startingfrom the ﬁrst column. Under categories, we have ”tokens” such as VERBOSE. Thus we specify
full name of token VERBOSE under category IO as IO_VERBOSE.
• IO_TIM is for debugging. It shows which subroutines are called and so on. Bigger number
shows deeper subroutines.
• SYMGRP is a category without token under it; we set generators of space group (See explanation in previous paragraph). When we set find, it automatically calculate symmetryof crystal lattice. If we like to enforce symmetry, set some of generators which are shownby lmchk.
• We see ctrls is embedded in the ctrl by ctrlgenM1.py.
... (skip) ...
% constda=0 alat=6.798

STRUC

ALAT={alat} DALAT={da}44PLAT=

0.0 0.5 0.5

NL=4 NBAS= 1
SITE

0.5 0.0 0.5

0.5 0.5 0.0

NSPEC=1

ATOM=Cu POS=0 0 0... (skip) ...
NL, NBAS(number of SITE) and NSPEC(number of SPEC) are automatically added by
ctrlgenM1.py. It is possible to deform unit cell by adding optional tokens (it is possible torotate PLAT for magnetic anisotropy calculation). See http://titus.phy.qub.ac.uk/packages/LMTO/tokens.html#
For new calculations, it is better to ﬁnd some examples ﬁrst.
• SITE category: As for MT sites, we have two categories. (1)SPEC(species) and (2)SITE(specify
centers of atoms(species) in primitive cell). As for SPEC, we specify MTs(radius, Z,
MTOs on it) appeared in the cell. These are deﬁned subtokens under SPEC ATOM=foobar(we have multiple SPEC ATOM=foobar).
Then we place these MTs at SITE sections in the cell. At SITE, we specify atomic
sites (What SPEC ATOM is placed to positions by POS) in a primitive cell. We set
POS= by direct form (Cartesian) but with the unit of ALAT+DLAT. Total number of SITE
(number of tokens SITE ATOM) is the number of atoms in the primitive cell. Setting
POS= under SITE ATOM=foobar means that we place MT named as foobar deﬁned inSPEC ATOM=foobar. In addition, we can set SITE ATOM RELAX, if you like to ﬁnd
relaxed structure (we simultaneously set DYN category) in LDA. As for relaxation, see
LaGaO_relax/ctrl.lagao example, and read
http://titus.phy.qub.ac.uk/packages/LMTO/tokens.html#DYNcat.
The SITE ATOM=foobar (with same foobar with diﬀerent POS) are not necessarilyequivalent with respect to the space group operation of a system. Thus SITE_ATOM=foobar
are divided into “classes” which are connected by the operation. The lmf automaticallyjudge “classes” (see also info by lmchk). Thus not need to specify it, but it may be betterto check it. A sample is lmchk lagao at ~/ecalj/lm7K/TESTsamples/LaGaO_relax• SPEC category: In ctrls, we have not yet speciﬁed contents of SPEC; we have just given
default symbols or only Z= when we use non-default names (shown by ctrlgenM1.py –showatomlist). The command ctrlgenM1.py adds default SPEC sections.
We have some SPEC_ATOM, under which we give subtokens such as SPEC_ATOM_R(MTradius), SPEC_ATOM_Z(nucleus charge), cutoﬀ parameters of angular momentum, and soon. These SPEC_ATOM is refereed to in SITE.
An example of SPEC category is
SPEC
ATOM=Fe Z=26 R=1.70
KMXA={kmxa}

LMX=3 LMXA=4 NMCORE=1

PZ=0,3.9,4.5
EH=-1 -1 -1 -1

RSMH=0.85 0.85 0.85 0.85

EH2=-2 -2 -2

RSMH2=0.85 0.85 0.85

MMOM=0 0 2 0
ATOM=... (then the similar block of ATOM= are repeated.)...
Under the token ATOM=Fe, we have subtokens SPEC_ATOM_Z,SPEC_ATOM_R, and so on.

45Subtokens Z= is the nucleus charge and R= MT radius. Note that Fe is just a name todistinguish MT sphere in the cell. If you set SPEC ATOM Z=27, it is recognized as Co(since Z=27). LMX=3 is the maximum l of MTOs. Thus maximum l of MTO is l=3. The
maximum of l to expand electron density and potential within MT is LMXA (in contrast to usual LAPW), we can use quite small LMXA such as LMXA=4. NMCORE=1
means we calculate core density without non magnetically-polarization. This can reduce
computational confusion.
PZ is to set local orbital (if not, no local orbitals). EH and RSMH are to specify ﬁrst setof MTOs.(We can check how local orbitals are set by lmfa explained in the next section).
EH2 and RSMH2 are to specify second set of MTOs.
After PZ=, we have three numbers. These are numbers for s,p,d,f,g,... channels. Zeromeans not exist. You can use space or comma(,) as delimiter. Here not only the integer part of principle quantum number, but also the fractional part should be supplied (If PZ=0,3,4, it does not work.) Now PZ=3.9 for p and PZ=4.5 for d. Thismeans we use local orbital for 3p, and local orbital for 4d (fractional parts (continuous principle quantum number) are large ∼ 0.9 for core like orbital, and smaller
for extended orbital ∼ 0.3 or something. See Logarithmic Derivative Parameters at
http://titus.phy.qub.ac.uk/packages/LMTO/lmto.html). This is a little confusing,thus we will explain this in appendix. See Sec.??.
EH(damping factor), and RSMH (where the smooth Hankel function bent) determines
MTOs (or its envelope function as a smooth Hankel function). We now set four numbers
for them. Thus we set MTOs s,p,d,f with EH=-1 and RSMH=0.85. Our current testshows that RSMH is one half of R (that is, 0.85=1.70/2, but minimum RSMH is 0.5)
and not need to be dependent on s,p,d,f. (If LMX=2, s,p,d are allowed and no f MTOs.)
EH is -1; not need to change except test purpose. In a similar manner, EH2 and RSMH2
for second set of MTOs are given. Just three numbers means these for s,p,d.
MMOM=s,p,d,f... gives initial condition of magnetic moment in µB (number of up-down
electron).
In cases such as As, the local orbital given by default ctrl is responsible of rather deepcore, and it is not need to be treated as valence electrons. In such a case, we don’t need
local orbital.
In the case of AntiFerro-II NiO, it contains two NiO in a primitive cell. Thus it is reasonable to have two SPEC ATOM as Ni1 and Ni2, although subtokens under ATOM=Ni1
and ATOM=Ni2 (e.g. SPEC_ATOM_EH for them) are the same except initial condition ofmagnetic moment of MMOM=s,p,d,f... See example of NiO.
The minimum help of call Category token subtoken are listed with minimum explanation
with$ lmf --input
It gives a long output. But many of them are experimental and not need to manage them. A
part of it is46

Token

Input

cast

(size,min)-----------------------------------------... ...
STRUC_ALAT

reqd

r8

1,

1

Scaling of lattice vectors, in a.u.
... ..
This is an minimum explanation of it. ”reqd” means ”required” (no default). r8 means it readwith real number, 1,1 means that ALAT=xxx should contain one number minimum (max is
also one) (See also STRUC PLAT, and so on).
There are kinds of examples in ecalj packages. Please look into their ctrls.* and ctrl.*
These are in lm7K/TESTsample/* and ecalj/CMDsampls. In addition, ecalj/MATERIALS
contain many samples (need a command); see a later subsection.
As for what is shown in $ lmf --input, most of important tokens are already describedin the ctrl ﬁle generated by ctrlgenM1.py. So, we don’t need to care many options shown by
it.
But we have not yet explained some useful features; STRUC category to deform crystal; DYN category for dynamics; LDA+U treatment; Adding background charge; Core-Holetreatments. We will prepare examples for it if requested.
http://titus.phy.qub.ac.uk/packages/LMTO/tokens.html#STRUCcat
http://titus.phy.qub.ac.uk/packages/LMTO/tokens.html#DYNcat
For QSGW calculation:
We need a setting in ctrl ﬁle to read sigm ﬁle (HAM SIGP). It is simpliﬁed now,and not need to care it so much. As we set RDSIG=12 in defaults, lmf read sigm
ﬁle and add it to one-body potential as long as sigm.* exist.
NOTE for old users: We now set SIGP[MODE=3 EMAX=9999.] in ctrl ﬁle to read
self-energy in lmf (or lmf-MPIK). This is because we use very localized MTOs (similar with the Maxloc Wannier). Our test shows reasonable results and this simplify
algorithms. In my previous version, we asked you to use SIGP[MODE=3 EMAX=2.0]
where EMAX is a little (0.5Ry) less than emax_sigm. If something strange occurs,
try this setting).
• In principle, QSGW result should not depended on the choice of XCFUN. However, it can aﬀect slightly. In our tests, it seems slightly better to use VWN (XCFUN=1) for QSGW calculations. (BUT need to check more...)10.1

How to set local orbitals

As we stated, do "lmfa |grep conf" to check used MTO basis.
We have to set SPEC_ATOM_PZ=?,?,?
(they ordered as PZ=s,p,d,f,... ) to set local orbitals.
lmv7 (originally due to ASA in Stuttgart) uses a special terminology"continious principle quantum number for each l", which is justrelatated to the logalismic derivative of radial funcitons at MT

47

boundary. It is defined as
P= principleQuantumNumber + 0.5-1/pi*atan(r* 1/phi dphi/dr),where phi is the radial function for each l.

For example,

P= n.5 for l=0 of free electron (flat potential) because phi=r^0,
P= n.25 for l=1 because phi=r^1;
P= n.147584 for l=2 because phi=r^2; P=, n.102416 n.077979 for l=3,4.
(Integer part can be changed). See Logarithmic Derivative Parameters in
http://titus.phy.qub.ac.uk/packages/LMTO/lmto.html#section2
Its fractioanl part 0.5-atan(1/phi dphi/dr) is closer to unity forcore like orbital, but closer to zero for extended orbitals.
Examples of choice:Ga p: in this case, choice 0 or choice 2 is recommended.
We usually use lo for semi-core, or virtually unoccupied level.
(0)no lo (4p as valence is default treatment without lo.)
3p core, 4p valence, no lo: default.
Then we have choice that lo is set to be for 3p,4p,5p.
(1)3p lo ---> 4p val (when 3p is treated as valence)
3d semi core, 4d valence
Set PZ=0,3.9(P is not requied to set. *.9 for core like state. It is just an initial condition.)(2)5p lo ---> 4p val (PZ>P)
Set PZ=0,5.5
5.5 is just simply given by a guess (no method have yet
implemented forIf 5.2 or something, it may fails
because of poorness in linear-dependency. We may need to observe
results should not change so much on the value of PZ.
(3xxx)4p lo ---> 5p val (we don’t use this usually. this is for test purpose)
4p lo, 5p valence
Set PZ=0,4.5 P=0,5.5 (In this case, set P= simultaneously).
(NOTE: zero for s channel is to use defalut numbers for s)
Ga d: (in this case, choice 0 or choice 1 is recommended).
(0)no lo (3d core, 4d valecne, no lo: default.)
Then we have choice that lo is set to be for 3d,4d,5d.
(1) 3d lo ---> 4d val
Set PZ=0,0,3.9(when 3d is treated as valence)(P is not required to set)(2) 5d lo ---> 4d val(PZ>P)

Set PZ=0,0,5.5(3xxx) 4d lo ---> 5d val(this is for test purpose)

Set PZ=0,0,5.5 P=0,0,4.5(NOTE: zerofor s,p are to use defalut numbers )If you like to read from atm.ga file instead of rst file(if exist).

48

You have to do lmf --rs=1,1,0,0,1, for example. See lmf --help
Becase rst file keeps the setting of MTO, thus change in ctrl is notreflected without the above option to lmf.
=============================================================

49

11

GWinput details

11.1

generate a template of GWinput

As in the previous section, we need two input ﬁles ctrl.si and GWinput. In principle, these twodetermines ﬁnal results uniquely. A template GWinput.tmp is generated by mkGWIN lmf2.
Required ﬁles are
Input ﬁles
• ctrl.si : The control ﬁle for PMT method.
(Recently modiﬁed mkGWIN lmf2 runs lmfa internally. If you use older version, do lmfa
in advance).
Output ﬁles
• GWinput.tmp : A ﬁle including computational conditions for the GW calculation. In addition, it speciﬁes the k points for which you calculate the QP energy.
When you run mkGWIN lmf, it asks you to supply three numbers for BZ integration as== Type three integers n1 n2 n3 for Brillouin Zone meshing for GW! ==n1=
Then you need to type a number e.g. as ”2 Return ” for n1. Then you need to repeat itfor n2 and n3 as
n1= 2 Returnn2= 2 Returnn3= 2 Return. These numbers speciﬁes what k points in BZ is used for BZ integration (In this case,2 × 2 × 2 = 8 k point in the 1st BZ is used. Roughly speaking, we need 4 × 4 × 4 (or 6 × 6 × 6)
to get band gap for Si and so on, with ≈ 0.1 eV accuracy.)
Then you have to edit GWinput.tmp and copy it to GWinput. We details the GWinput inlater chapter.
We need to repeat mkGWIN lmf2 when you change MTO sections in ctrl ﬁle (adding PZcase, and so on).

11.2

overview of GWinput(Because of historical reason, input ﬁle is diﬀerent from ctrl.*).
The main input ﬁles is GWinput. This controls the setting of GW calculation. The ﬁleGWinput consists of structures as
keyword1 data1
keyword2 data2...
In each lines, it consists of keyword and data. Data can be single or plural. As for keywords,
upper case or Lowercase is not distinguished. All keywords should start from 1st column (nospace at head). Order of lines are irrelevant. As for logical variable, you can use anything“true, yes, on, 1, T” for .true., and anything “false, no, oﬀ, 0, F” for .false.
Or we have “tag sections” in GWinput speciﬁed by <PRODUCT_BASIS>, <QPNT>, <PBASMAX>,<QforEPS>, and <QforEPSL>.

(<PRODUCT_BASIS> is requires for all kinds of calculations.

<PBASMAX> is optional. <QforEPS> and/or <QforEPSL> are required for epsilon mode). It
is like<PRODUCT_BASIS>
tolerance to remove products500.100000D-02 ! =tolopt
lcutmx(atom)
3 3
atom

l

...
</PRODUCT_BASIS>. In these tag sections, you have to keep format for its own (usually numbers are read by freeformat read(*,*)).
The fundamental readin routine for GWinput is a subroutine getkeyvalue deﬁned in
gwsrc/keyvalue.F. getkeyvalue is a general and convenient readin routine in full use of
the f90 features. Read a head part of the ﬁle and try to do ”grep getkeyvalue *.F” in gwsrc/
or main/ so as to see how to use it (test routine is main/kino_input_test.F.)
Soche the GWinput consists of three sections

1. General section
2. <PRODUCT_BASIS> section
3. <QforEPS>,<QforEPSL> section (only eﬀective for dielectric function mode)
4. <QPNT> section (only eﬀective for one-shot mode)
5. <PBASMAX> section (optional)
We will explain each by each in the followings.

51

11.3

General section

In general section, it looks like! ##### From GWIN0 ################
n1n2n3

1

QpGcut_psi

4.000 !(See unit_2pioa for unit) |q+G| cutoff for eigenfunction.

1

1 ! for BZ meshing in GW

QpGcut_cou

3.000 !(See unit_2pioa for unit) |q+G| cutoff for Coulomb and W.

unit_2pioa off ! off --> a.u.; on--> unit of QpGcut_* are in 2*pi/alatalpha_OffG

1.000 !(a.u.) Used in auxially function in the offset-Gamma method.

!emax_chi0

99999.000 !(Ry) emax cutoff for chi0

emax_sigm

3.000 !(Ry)(Optional)

emax cutoff for Sigma! ##### FREQUENCIES from GWIN_V2 ################
dw

0.005000 !(a.u.) energy-mesh (bin width size) along real axis.

omg_c
!

0.040 !(a.u.) energy-mesh is twiced at omg_ccoaser mesh for higher energy. Width get to be doubled at omg_c.

iSigMode

3 ! QSGW mode switch for gwsc. use =3.

niw

10 ! Number of frequencies along Im axis. Used for integration to get Sigma_c! E.g. try niw=6 and niw=12
delta-0.10D-05 !(a.u.)

Broadening of x0. negative means tetrahedron method.

! used by hx0fp0. You get smeard x0 witth abs(delta).
deltaw0.020000 !(a.u.) Mesh for numerical derivative to get the Z factor

esmr

0.003000 !(Ry) used by hsfp0. Keep esmr smaller than band gap for insulators! Poles of G^LDA are treated as if they have width esmr in hsfp0.
! Change esmr for metals.
GaussSmear on

See DOSACC*---especailly around Ef.

! Gaussian or Rectangular smearing for Pole of G^LDA with esmr for hsfp0.

—————————————————

1. BZ integration.
n1n2n3 3 integers as N1 , N2 , N 3 (no default); They are ≥ 0.
Brillouin Zone mesh for integration is determined by keywordsBZmesh and n1n2n3. Current version only allow regular mesh point including Gamma point for G × W . But
not that Chi RegQbz below allow you to use oﬀ-Gamma mesh for W (ω) (and dielectric
function mode).
We usually take ’4 4 4’, ’6 6 6’ or ’8 8 8’ for GaAs. For metal such as Fe, ’10 10 10’ ormore is better.
Chi RegQbz (on or oﬀ)
Chi RegQbz = on (default): Use regular mesh (including gamma) for eps calculation.
Chi RegQbz = oﬀ : Use oﬀ-Gamma mesh (Not including gamma) for eps calculation.
(In cases, Chi RegQbz off gives faster convergence as for n1n2n3; not only for GW, butalso for dielectric functions eps lmfh.)
2. Plane wave (q + G) cutoﬀ
QpGcut psi 1 real (no default)
QpGcut Cou 1 real (no default)
unit 2pioa 1 logical (no defalt)
We have two cutoﬀ for q + G. QpGcut psi is the cutoﬀ of |q + G| for the IPW in the
expansion of the eigenfunctions. QpGcut Cou is for the IPW of the interactions v, D, W .

52

Its unit is speciﬁed by unit 2pioa ; ”oﬀ” means unit in a.u. and ”on” means in unit of
2π

alat

. (alat is length scale unit in ctrl.*).

Rule of thumb: QGcut_psi is a little (usually 0.5 or so) larger than QpGcut_cou. It
becomes accurate if we use large QpGcut_cou. But it enlarge size of IPW(interstitialplane wave) part of Mixed product basis. For test, try 2.7, 3.2, 3.7 for QGcut_cou (andadd 0.5 or 1 for QGcut_psi). Larger one is expensive.
We expand eigenfunctions in the Muﬃn-tin division of the space. See Eq.[?] in Ref.[?];
in the current GW implementation, we use very simple form of eigenfuncitons (not by
the 3-component formalism in the [?]).
Thus the form of expansion is just related to the division of space; not directly related
to the diﬀerence among LAPW, LMTO, and PMT.
3. Eigenfunctions within MTs (no parameters setting for them).
The radial functions (phi and phidot for each l), corresponding to the true parts, (= 2ndcomponent in the 3-component formalism [?]), are automatically determined already inthe one-body part of program lmf-MPIK.
4. Cutoﬀ for used bands.
emax chi0 : 1 real (optional,default=∞), in Ry
emax sigm : 1 real (optional,default=∞; We usually use 3 Ry).
emax_sigm is the maximum limit of the self-energy (measured from the Fermi energy).
See the paper [?] which shows how the results are aﬀected by emax sigm. But in cases,
small emax_sigm can give poor dispersion curve (slightly unnatural behavior) because of
sudden cutoﬀ by emax sigm. However, we like to use smaller value to reduce computational time.
That is, larger is better, but expensive (And note that we simultaneously need to use
empty spheres when we use large emax sigm, as shown in [?], ).
Generally speaking, accuracy less than ∼0.1eV (for bandgap) is allowance of current
technique. Probably, it may be possible to have better accuracy, but it may ask us torepeat many calculations with changing conditions to conﬁrm stability.
nband chi0 : 1 integer (optional,default=∞)
nband sigm : 1 integer (optional,default=∞)
These specify how many bands you use in hx0fp0 (for chi0) and in hsfp0 (for sigma).
Higher bands above them are neglected.
5. Energy mesh related parameters.
— new after march2016 —HistBin dw : 1 real (a.u.). Mesh width along real axis for W (ω).
HistBin ratio : 1 real (a.u.).
HistBin_dw and HistBin_ratio specify real space bins which we accumulate imaginary
part weight of polarization functions. The bins are (see frhis in ecalj/fpgw/gwsc/m_freq.F)
ωi = b ∗ (exp(a ∗ (i − 1)) − 1),(4)

where [ωi , ω( i+1)] (i = 1, 2, ...,nwhis+1) is the i-th bin. HistBin_dw is bin width at ω = 0.
The ratio ω( i + 1)/ω( i) for large ω is exp(a)=HistBin_ratio. This choice of gettingcoarser at high energy is because we think W (ω) around ω ∼ 0 gives most important53

contribution to the GW approximation. If histogram bins are too wide, dielectric function
can be less accurate, but results may be not so much aﬀected.
In GW calculation, the Plasmon pole is important. It is determined not by the Drude
weight; it usually gives small contribution to the Plasmon pole (for example, Si is described well by a Plasmon pole model, but Si has no Fermi surface). We expect that the
GW results are not so sensitive to the choice of HistBin_dw, HistBin_ratio usually.
We may use ﬁne mesh when we plot quantities such as W (ω) near ω = 0.
The ecalj gives W̄ (ω = 0) ∼ 0 for metal; where W̄ (ω) is the eﬀective interaction averaged
in the Γ-cell [10]. And W¯(ω get closer to v for larger ω.
Default value of them are given in fpgw/gwsrc/m_freq.F. They arecall getkeyvalue("GWinput","HistBin_ratio",oratio, default=1.03d0)
call getkeyvalue("GWinput","HistBin_dw",dw, default=1d-5) !a.u.
These are safer settings. Thus new mkGWIN_lmf2 (GWinput.tmp generator) gives default
lines in GWinput as
HistBin_dw

2d-3 ! 1d-5 for metal!(a.u.) bin width along real axis at omega=0.

HistBin_ratio 1.08 ! 1.03 for check. frhis(iw)= b*(exp(a*(iw-1))-1),
where a=ratio-1.0 and dw=b*a. This is given in the conv2gwinput.F.
— old version. before 14march2016 –dw : 1 real (a.u.). Mesh width along real axis for W (ω).
omg c : 1 real (a.u.).
dw and omg c determines ω mesh along real axis to calculate W (ω)...
————————(WARNING! Some of my examples may show as if they are in ”(Ry)”. But they are Wrong!)

delta : 1 real (a.u.). We usually use very small number as -1d-8 for gw lmfh, eps modeand so.
xxx does this really make the stabilization? xxx This is the size of δ in denominatorof Π (EQ.xxx). But (I think) we can (or can not) use it so as to make broadening fortheoretical test (maybe not exactly corresponding to δ). or when you make calculation
stabilized (Takao need to check this point again xxx!)[Old note. Need check xxx: In gw lmf, it is used for broadening of x0 when it call
hx0fp0. Then delta is δ is EQ.32. The sign of delta is just used as a ﬂag whether you
use the tetrahedron method of dielectric constant [11] or not; minus sign means ”Use thetetrahedron method for D”; plus sign means you do it by simple sum. You can usually
use this default setting. But it might be possible to use a larger value to smear the ﬁne
structures on the energy-dependence of W in cases. This might be necessary if W is soenergy-dependent and dw is not so small to resolve the structure —but I don’t know.]
niw : 1 integer.
Number of integration points along the imaginary axis(FIG.1) to get Σc .

See rou-tines wint* called from sxcf*.F, which is called from the main routine hsfp0.m.F (orhsfp0.sc.m.F in the QSGW case). The integration points are iω ′ (n) = i(1/x(n) − 1),
where x(n) is the usual Gaussian-integration points for the interval [0,1]. In addition,54

we give the special analytical treatment for the peaky part at ω ′ = 0. Out tests shows
niw=6 for Si is good enough for 0.01 eV accuracy. The convergence as for niw is quite
good. This integration scheme has been developed by Ferdi Aryasetiawan. The number of points should be the one of 6,10,12,16,20,24,32,40,or 48. It is because we use asubroutine gauss in /gwsrc/mate.F prepared by Ferdi. We will replace better one in
future. See II-F in Ref.I.
GaussSmear : 1 logicalesmr : 1 real (Ry). Used by hsfp0 (and hsfp0.sc for QSGW).
Poles of the Green function GLDA are treated as if they have width esmr in hsfp0. If
GaussSmear is on, each pole of GLDA is smeared by a Gaussian function with σ =esmr
in the calculation of hsfp0. If GaussSmear is oﬀ, we assume rectangular smearing for thepoles. Usually it is necessary to take rather smaller value than band gap for insulators.
Try to use 0.003 or so in the case of Si and GaussSmear=on.
For metal, this esmr is somehow related to how we capture the Fermi surface; In principle,
we have to take the limit n1n2n3 → ∞ and esmr → 0). However, we may inevitably use
some ﬁnite esmr to make calculations converged.)deltaw : 1 real (a.u.) only for one-shot case.
in EQ.8. We calculate hψ kn |Σ(ǫkn +deltaw is the interval for the numerical derivative ∂Σ(ω)∂ω
deltaw)|ψ kn i and hψ kn |Σ(ψ kn − deltaw)|ψ kn i in addition to hψ kn |Σ(ǫkn )|ψ kn i. Fromthese values, we can calculate two Z (or second-derivative of Σ(ω)), as shown in SECU.
It will help to see whether the used deltaw is O.K. or not.
6. Oﬀset-gamma point.
Q0PChoice 0 :1 integer
Q0P Choice gives how to determine the oﬀset gamma points. Initially we take them as
1: —q— is ten times smaller than regular mesh.(default)
2: —q— is average in the Gamma cell (cell of BZ including Gamma point).
Then we choose only inequivalent q points based on the point group symmetry. Obtainedoﬀset gamma points is given in Q0P ﬁle.
alpha offG : 1 real (a.u.)alpha offG corresponds to α in EQ.48. alpha offG=1d0 is usually good in the sensethat it seems to be almost a limit at α → 0. So you can usually ﬁx it as alpha offG=1d0,
and check the convergence as for n1n2n3.
7. core orthogonalization (default=oﬀ)
CoreOrth 1 logical — recently, this option is not maintained — Better to use local
orbital instead, so that core charge do not spill out.
If this is on, we enforce cores orthogonalized to valence φ and φ̇ (these appear in II-C
in Ref.I). This procedure enforce the correct orthogonal condition, thus we have correct
behavior for the dielectric function at q → 0. However, it may deform core functions
too much, especially in the case of shallow 3d (or maybe 4d) cores. So we don’t recommend use this option, even though then the orthogonality condition is somehow broken.
Anyway you can check weather it aﬀects to results or not by this switch.
8. QP self-consistent GW.
iSigMode 1 integer (no default).
This is required for QSGW calculation by the script gwsc. We have some possible ways
to make GW self-consistent (how we determine Vxc from calculated Σ(ω)). We now
mostly use iSigMode=3.

55

Σnn′ (ǫn ) + Σnn′ (ǫn′ )(mode-A in [?]).
2
1: Use Σnn′ (EF ) + δnn′ (Σnn′ (ǫn ) − Σnn′ (EF )) (mode-B in [?]).
3: Use Re

5: Use δnn′ Σnn (ǫn ) (Eigenvalue-only self-consistency, keeping the eigenfuncions as given)).
See /gwsrc/sxcf_fal2.sc.F, which is called from the main routine hsfp0 sc (this is theroutine to calculate self-energy)).
9. Others
KeepEigen 1 logical (default=on)
These are for memory usage. When KeepEigen is on, eigenfunctions (Eigen) are kept inmemory during calculation. If you have not enough memory in your machine, use themoﬀ. Then you can save memory usage. However, then we may have too frequent accessto ﬁles. So %CPU might get lower. Be careful to use these options.
Verbose 1 integer (default=0) If 0, it gives minimum standard output. If 40 or higher,
it shows too much output. (these verbosity control is not well-organized yet).
10. LFC@Gamma , EIBZmode , multitet are for test purpose.xxx

11.4<QPNT> section(only for one-shot GW. Not suitable to make band plot in BZ.) This section is to specify theq points and bands index for which you calculate the QP energies (QPE). An example is<QPNT>--- Specify the q and band indeces, for which we evaluate the self-energy --*** all q -->1, otherwise 0; up only -->1, otherwise 000*** no. states and band index for calculation.
3
15 16 17*** q-points, which shoud be in qbz. See KPNTin1BZ.
2
10.0000000000000000
0.0000000000000000
0.0000000000000000
20.2886751358562925-0.5000000000000000
0.0000000000000000
30.0000000000000000
0.0000000000000000
0.3073140749846343</QPNT>
Numbers are read by free format read(5,*), thus the numbers should be separated by space.
At the next line to the ﬁrst ***, you have to give two numbers used as ﬂags. Both of themtakes 0 or 1. 1st one is whether you calculate QPE for all q points (in IBZ) or not. If it is 1,
you calculate QPE for all q. If it is 0, you calculate them only for q points speciﬁed within
this ﬁle. In the case of metal where you want to calculate the Fermi energy for QPE, you needto calculate all the eigenvalues somehow above the Fermi energy (If you put 1, it is safer but
too time-consuming). The second number is whether you calculate QPE for both spins or not.
It is usually 0. In the case of antiferro material, it should be 1.
From the next line to the second ***, you have to specify the states for which you calculate the QPE. In this example, you calculate the 3 bands of QPE for 15th, 16th, and 17th
eigenfunctions (they are ordered from the bottom).
From the next line to the third ***, you have to specify the q points. The ﬁrst numbersof each line are dummy. In this case, you calculate QPE for two q points. The third q pointis neglected because 2 is given at ﬁrst.
When you generate GWinput.tmp, you see all the possible q points are listed (these q
points should be a part of the regular mesh points).
In the QSGW mode (gwsc), this section is neglected (then we calculate all QPE on regularmesh points); so its hsfp0 sc part is quite expensive (usually it takes time more than hx0fp0).
Additional Note ————–

56

QPNT nbandrange num1 num2 (two integers).
This override setting in <QPNT>. (I think this switch may still work, but not checked recently).
AnyQ on (default is oﬀ)
If this is on, you can specify any Q point which is not on the mesh point. For the purpose, weneed to prepare eigenfunctions at extra k points. But it is automatic. In order to make thecomputation eﬃcient. Even in this case, from the computational view, it is better to choose
q on the two times ﬁner divided mesh (or three times ﬁner divided k mesh). This is used forFig.6 in Phys. Rev. B 74, 245125 (2006).

57

11.5

set QPNT for eps mode (QforEPS section)

For eps modes (scripts eps_*, which are for linear responses. See Sec.16), you have to specify
q point in the following ways.
1. QforEPSIBZ on
Then all Q point in IBZ are used.
2. Use section as<QforEPS>
0d0 0d0 0.01d0
0d0 0d0 0.02d0
0d0 0d0 0.04d0
0d0 0d0 0.08d0</QforEPS>
In addition, you can specify Q points as<QforEPSL>
0d0 0d0 0d0

1d0

0d0

0d0 8

0d0 0d0 0d0

.5d0

.5d0

0d0 8

</QforEPSL>This is along the line— 8 point along the line (not left-end q; so omitting 0 0 0). The ﬁrst
line means line (0d0 0d0 0d0)—(1d0 0d0 0d0) is divided to 8. So we have 7 points, (0.125 00), (0.25 0 0),... (1 0 0).

58

11.6

<PRODUCT BASIS> section

This section is to deﬁne product basis to expand W and so. Numbers are read by free formatread(5,*),thus the numbers should be separated by space. The line number in this section ismeaningful (you can not add comment lines).

<PRODUCT_BASIS>
tolerance to remove products due to poor linear-independency
0.100000D-04 ! =tolopt; larger gives smaller num. of product basis. See lbas and lbasC, which are outp
lcutmx(atom) = maximum l-cutoff for the product basis. =4 is required for atoms with valence d, like N
4 3
atom
l nnvv nnc ! nnvv: num. of radial functions (valence) on the augmentation-waves, nnc: num. fo
1
0
2
3
1
1
2
2
1
2
3
0
1
3
2
0
1
4
2
0
2
0
2
1
2
1
2
0
2
2
2
0
2
3
2
0
2
4
2
0
atom
l
n occ unocc ! Valence(1=yes,0=no)
1
0
1
1
1
! 4S_p ----102
10! 4S_d1
1
1
1
1! 4P_p1
1
2
00! 4P_d12
1
1
1! 4D_p12
2
00! 4D_d1
2
3
1
1! 3D_l1
3
1
0
1! 4f_p13200! 4f_d141
00! 5g_p14200! 5g_d201
11! 2S_p ----20200! 2S_d2
1
1
1
1! 2P_p21
2
00! 2P_d221
1
1! 3d_p22200! 3d_d2
3
1
0
1! 4f_p23200! 4f_d241
00! 5g_p24200! 5g_datomln occ unocc ForX0 ForSxc ! Core (1=yes, 0=no)
1
0
1
0
0
0
0
! 1S ----1
0
2
0
0
0
0
! 2S
1
0
3
0
0
0
0
! 3S
1
1
1
0
0
0
0
! 2P
1
1
2
0
0
0
0
! 3P
2
0
1
0
0
0
0
! 1S ----</PRODUCT_BASIS>
This section is read in the free format in fortran. So, e.g., 0.01 works as same as 0.10000D-01.
The line order is important (you have to keep the order given by GWinput.tmp). Be carefulatom atom id—lmf may re-order it and pass it to gw code. Look into LMTO ﬁle (generated
by mkGWIN lmf2); which contains crystal structure information after such re-ordering bylmf. I used ! to make clear that things after ! are comments. But ! is not meaningful – justthe expected numbers of data separated by blank(s) are read for each line from the beginningof lines.
• 0.100000D-02 ! =tolopt controls a number of Product basis to expand the Coulombinteraction within MTs. tolopt is a criterion to remove the poorly linear-independent
product basis. Note that the product basis, which is to expand the Coulomb interaction,is diﬀerent from the basis to expand eigenfunctions. In our experience, 0.100000D-02(=0.001) is not so bad. If you like to reduce computational time use 0.01 or so, but alittle dangerous in cases. With 0.0001, we can check stability on it.
(note: By supplying multiple numbers, we can specify tolopt atom by atom. Remember
lmchk gives atom ID.)59

• lcutmx(atom) is the l cutoﬀ of product basis for atoms in the primitive cell (do lmchk
for atom id). In the case of Oxygen, we can usually use lcutmx=2 (need check by
the diﬀerence when you use lcutmx=2 or lcutmx=4). Then the computational time is
reduced well.
• (dec2014:<PBASMAX> is not checked recently; see fpgw/main/hbasfp0.m.F and fpgw/gwsrc/basnfp.F).)
You can use <PBASMAX> section to override this setting. It is given as<PBASMAX>1

5 5 5 3 32

5 5 3 2 33

3 3 2 2 2</PBASMAX>
The ﬁrst number is for atom index (ﬁxed), and other are product basis for each l channel.
• The integer numbers in 4th line lcutmx gives the maximum angular momentum l for the
product basis for each atomic site. In our experience, lcutmx=4 is required when thesemi-core (or valence ) 3d electrons exist and we want to calculate the QP energies ofthem.
• Keep a block starting from ” atom l nnvv nnc ...” as it originally generated in GWinput.tmp. It just shows that how many kinds of radial functions for cores and valence
electrons for each atom and l. nnvv=2 in the case of φ and φ̇; nnvv=3 in the case to add
the local orbital in addition.
• There are two blocks after the line ”atom l n occ unocc :Valence(1=yes, 0=no)’ andafter ”atom l n occ unocc ForX0 ForSxc !

Core (1=yes, 0=no)’. These are used tochoose atomic basis to construct the product basis. The product basis are generated fromthe products of two atomic basis.
GWinput.tmp generated by mkGWIN lmf2 contains labels on each orbitals as 4S_p,
4S_d, 4P_p... Here 4S_p is for φ4s ; 4S_d for φ̇4s ; 3D_l for φlocal
3d . Capital letter just after
the principle-quantum number means the orbital is used as ‘Head of MTO’; lowercase
means just used only as the ‘tail of MTO’.
The switches for columns labeled as occ and unocc. take 0 (not included) or 1 (included).
With the switch, we can construct two groups of orbitals, occ and unocc. In this samcorecore
ple GWIN V2 as for atom 1, {φ4s , φ̇4s , φ4p , φ4d , φlocal
3d , φ3s , φ3p } consist the group occ,

and {φ4s , φ4p , φ4d , φlocal
3d , φ4f } consists the group unocc. So the any product of combicore
core
local
nations {φ4s , φ̇4s , φ4p , φ4d , φlocal
3d , φ3s , φ3p } × {φ4s , φ4p , φ4d , φ3d , φ4f } are included as

for the basis of the product basis. As for atom 2, {φ2s , φ2p , φ3d } × {φ2s , φ2p , φ3d , φ4f }are included.
• Core section: (not worth to read, since we currently use no CORE2, A=B=C=0.)
Each line of the last section of Product BASIS forms
atomln occ unocc
ForX0 ForSxc :CoreState(1=yes, 0=no)
1
2
1
A
x
B
C
At ﬁrst you have to understand the concept of CORE1 and CORE2 in EQ.35 Ref.I.
However, in our recent calculations, we do not use “CORE2” generally. So, in such acase, set A=B=C=0. And treat shallow cores (above Efermi−2Ry or so ) as valence electronby “local orbital method” in lmf.
• Be careful. Current version is inconvenient... Need to repeat mkGWIN lmf2 to generate
GWinput template when you add PZ (local orbital).

60

[( Note: you can skip here if you don’t use CORE2.)
Each of A,x,B,C takes 0 or 1.

There are some possible combination of theseswitches;
1. If you take ( A x B C )= (1 0 1 1), then the core is included in core2. Inother words, this core is treated in the same manner of the valence electron.
2. If you take ( A x B C )= (0 0 0 0), then the core is included in core1. The(exchange only) self-energy related to this core is included in SEXcore.
C is the key switch which determine whether it is included in core1 or core2.
There could be another option.
3. If you take ( A x B C )= (1 0 0 1). This core is in core2. But it is notincluded in the calculation of D and W . This core is only included for SEX
and SEC calculations.
These three kinds of choices are reasonable ones but we can consider some anotherchoice. In the following, we show how these switches (A,B,C) aﬀect executions
called from gw lmfh (essentially as same as gw lmf).
• hbasfp0(mode 3) :Product basis for exchange due to core.
We include the C=0 cores as a part of the product basis as if A=1 x=0.
• hsfp0(mode 3): exchange mode for core.
Σx only due to the C=0 cores are calculated.
• hbasfp0 (mode1): Product basis.
Only see the switch A and x. The product basis is generated from (occupied
× unoccupied), where A=1 core is included as one of the occupied basis.
• hsfp0 (mode 1): exchange mode.
Only see the switch C. Σx due to valence and due to C=1 cores are calculated.
• hx0fp (mode 1): W − v calculation.
Only see the switch B. W is calculates using all the valence and B=1 cores.
• hsfp0 (mode 2): correlation mode.
Only see the switch C. Σc due to valence and due to C=1 are calculated.
• After you perform gw_lmfh or anything, you ﬁnd output ﬁles lbas by hbasfp0 (mode1),
and/or lbasc by hbasfp0 (mode3) for core. These contains important information about how
many and how product basis are chosen. E.g. ’grep nbloch lbas’ shows how many productbasis are used in the calculations.

11.7 ANFcond (we can skip here since we do not check this option now.
Need fix this if necessary.
This ﬁle is used in hx0fp0 in the calculation of W − v (or rather Π in the program) to specify
the antiferro condition.
Note : Now only for the case that (a translation vector + spin ﬂip) is a symmetry operation.
This should be given by hand. For the cases of not antiferro, this ﬁle should not exist.
Even if ANFcond does not exists for antiferro case, hx0fp0 works but it requires about two
time computational eﬀorts.

61

The existence of this file means the Antiferro condition is used for x0k
Product basis B({\bf r}-{\bf a}) is translated to B({\bf r}-{\bf a}-Af})= B({\bf r}-{\bf a}’-T_0})
1d0 1d0 1d0! Af=Antiferro translation vector in Cartesian.
1 2
2 1
3 4
4 3
The ﬁrst line speciﬁes the Antiferro translation vector. From the second line, we specify
that atom i in the primitive cell is mapped to what atom j(i) in the cell with the opposite
spin by the translation. In this case, j(1) = 2, j(2) = 1, j(3) = 4, j(4) = 3. You have to becareful as for the true atomic position used in the GW calculations can be diﬀerent from thegiven atomic positions in ctrl.MnO. The true atomic positions is written in LMTO.
In the case of one-shot GW (gw lmf and gw lmfh), it may be better to set ”up only” QPE,
so that you only calculate QPE of up spins at the same time.
In the case of gwsc, we just calculate QPE for up spins automatically (QPNT section isneglected).

62

12

Main Output Files of GW part

12.1

QPU

This is the main output13 in human readable format.
An example of one-shot GW by gw lmfh) is (In the case of QSGW, Z (Z = 1) is notshown):===============================================================
quasiparticle energies MAJORITY===============================================================
E_shift= 0.4263273221017709D+00 0.6075150850568627D+00 0.7046628446164018D+00 eVqstate SEx
0.0 0.0 0.0 1 -29.56
0.0 0.0 0.0 2 -30.52
0.0 0.0 0.0 3 -20.67...

SExcore SEc
vxc
dSE dSEnoZ eLDAeQP eQPnoZeHF Z
2Z*Simg ReS(elda-1.97 10.40 -20.22 -0.52 -0.90 -19.08 -19.42 -19.71 -30.81 0.58 0.95
-21.12-2.24 10.09 -21.53 -0.70 -1.14 -18.06 -18.58 -18.93 -29.72 0.61 0.96
-22.66-1.87
5.97 -16.85 0.19
0.28 -7.20 -6.83 -6.65 -13.32 0.67 0.66
-16.57

From the 6h line, we have the eigenvalue data. All of the unit of energy is in eV. We shouldnote that the zero-level of these values eLDA eQP eQPnoZ can be changed by hqpe. This eLDA
- E shift are the eigenvalues relative to a Fermi energy determined by the smearing method.
Detailed value of eLDA is in TOTE2.UP. Detailed value of eLDA- E shift is in TOTE.UP.
q : k vector
state: Band index n, which is from the lowest eigenvalue (not include cores).
SEx: = = hΨkn |Σcore2+valence(r, r′ )|Ψkn i
x
SExcore: = hΨkn |Σcore1(r, r′ )|Ψkn ix
SEc: = hΨkn |Σcore2+valence(r, r′ , ǫn (k))|Ψkn ic
LDA
vxc: LDA exchange correlation energy. hΨkn |Vxc([ntotal ], r)|Ψkn i

dSE: Znk × dSEnoZ
LDA
dSEnoZ: hΨkn |Σcore1(r, r′ ) + Σcore2+valence(r, r′ , ǫn (k))|Ψkn i − hΨkn |Vxc([ntotal ], r)|Ψkn ixxc= SEx + SExcore + SEc - vxc
eLDA: LDA eigenvalues. ǫn (k)eQP: QP energy. ǫn (k)+dSE
eQPnoZ: QP energy without Z. ǫn (k)+dSEnoZ
eHF: HF energy of 1st iteration. ǫn (k)+SEx + SExcore -vxc
Z: Z factor. Znk
2Z*Simg: Quasi-particle life time. 2Znk × ImhΨkn |Σcore2+valence(r, r′ , ǫn (k))|Ψkn ic(Is this really the usual deﬁnition of the life time?—don’t believe me)
ReS(elda): RehΨkn |Σcore1(r, r′ ) + Σcore2+valence(r, r′ , ǫn (k))|Ψkn i
x
xc

12.2

XCU

LDA exchange-correlation. Detailed data of above vxc.

12.3

SEXU

Exchange part of the self-energy due to valence electrons. Detailed data of above SEx.

12.4

SEXcoreU

Exchange part of the self-energy due to core. Detailed data of above SExcore.
13 Note that QPU also implies QPD and so on. U is for up D is for down spins.

63

12.5

SECU

Correlation part of the self-energy. Detailed data of above SEc.

12.6

TOTE.UP (TOTE.DN)

This is a central output. It contains LDA and QP energies. These values are relative to a
Fermi energy determined by the smearing method. It contains two kind of QP energies QP
QPnoZ. The ﬁrst line contains the Fermi energy in Ry determined by the smearing method. It
is also shown in the end of DOSACC.lda.

12.7

TOTE2.UP (TOTE2.DN)

This is a central output. It contains zero-level shifts from TOTE.UP. The ﬁrst line contains
the Fermi energy in eV (= the Fermi energy in TOTE.UP but it is in Ry) and three energy
shifts E shift, which are the same values in the 4th line of QPU.

Note that all *.chk ﬁles are just to check calculations (not read in by successive executions).

12.8

DOSACC.lda

This lists all the eigenvalues in ascendant order. States with almost the same eigenvalues aredegenerated states. The 4th column contains number of electrons up to the eigenvalue.

12.9

DOSACC2.lda

This is similar with DOSACC.lda. But we remove the degeneracy.

12.10

Core ibas* l*.chk

Used core eigenfunctions.

12.11

VXCFP.chk

This contains eigenvalues and hψkn |Vxc |ψkn i in both units, Ry and eV. See below.

12.12

The Fermi energies in this GW code.

smear
tetra
We mainly have two kinds of Fermi energy EFEERMIEFEERMI.

1. At ﬁrst eigenvalues given by lmfgw is in VXCFP.chk. You can see%head VXCFP.chk### LDA exchange correlation ####qvec
ikp iband
eigen
0.0000 0.0000 0.0000 1 1-0.96932423...

VXC(ntotal)-1.00727912

VXC(nvalence)
0.00000000

These are raw values. TOTE contains the eigenvalues but relative to a Fermi energy
smear
EFEERMI
which is determined by the smearing method. It is also shown at the top partof output ﬁles lsx sf and lsc sf. And you also see the value at the end of DOSACC.lda.
This is the head of TOTE.UP;

64

eigen(eV)-13.18843159

VXC-

%head TOTE.UP
43
0.0000000
0.0000000...

80.0000000
0.0000000

8.520283353474250E-003
0.0000000
1
1 -0.1330435686590073D+02 -0.1322984339282777D+02 -0.1
0.0000000
2
1 -0.7555264915356062D+00 -0.6267595395613325D+00 -0.6

smearHere EFEERMI=8.520283353474250E-003. From the second lines, they are LDA eigensmear
values and QP energies (Z included and Z=1); they are relative to the EFEERMI.
smear-13.18843159 eV - EFEERMI(which should be translated into in eV) = -0.1330435686590073D+02

eV. Here -13.18843159 is the value in VXCFP.chk shown above.
tetra
2. There is the another Fermi energy EFEERMI, which is used by mode 11 (or mode 1) of

hx0fp0 in gw lmfh. It is determined by heftet and stored in EFERMI.
3. hqpe gives TOTE2.UP and QPU. They contains the same values. You can see eLDA
not only in QPU but also in TOTE2.UP. At top lines of TOTE2.UP, you see%head TOTE2.UP
43
8 0.1159252712507000D+00
0.0000000
0.0000000
0.0000000
10.0000000
0.0000000
0.0000000
20.0000000
0.0000000
0.0000000
30.0000000
0.0000000
0.0000000
4
..

eQP

eQPnoZ

0.7555207081466229D+00 0.6267572296579150D+00 0.60094
1 -0.1254883615775411D+02 -0.1260308616316985D+02 -0.1
1 -0.5783388983382487D-05 -0.2309903417430093D-05 -0.1
1 -0.1369098933889923D-05 -0.6195200397129952D-07 0.2
10.0000000000000000D+00 0.0000000000000000D+00 0.0

smear,where a number in ﬁrst line EFEERMI= 0.1159252712507000D+00 eV =8.520283353474250E-

003 Ry, the same as the previous one. This is a case when you did hqpe with augment
4 (it means we set the 4th-band eigenvalue zero). Another 3 values in the ﬁrst line areshifts from TOTE. Shown eshift(eLDA) = 0.7555207081466229D+00 eV. E.g., the secondline shows-0.1254883615775411D+02 eV= -0.1330435686590073D+02(in TOTE) + eshift(eLDA)
eV.
When you do hqpemetal, three shifts at the ﬁrst line in TOTE2.UP is determined so as
to give the eigenvalues relative to the Fermi energies shown in EFERMI, EFERMI.QP1,
and EFERMI.QPz=1. These are Fermi energies by tetrahedron method.

As for gwband lmf, it recalculates eigenvalues for all q along SYML. Thensmear
the default ”zerolevel” = EFEERMI- eshift(lda). Because the eigenvalues given by this band-mode are presumably the same, we haveShown LDA eigenvalue= -13.18843159(raw data by band mode—same as that in VXCFP.chk) - zerolevel= (-13.18843159 - EFERMIsmear) + eshift(lda).
= -0.1330435686590073D+02( this is in TOTE.UP) + eshift(lda)= -0.1254883615775411D+02( this is in TOTE2.UP).
It means that values in TOTE2.UP recovers. But if raw data by band mode is diﬀerent
from it, these is a trouble. It does not recover the values in TOTE2.UP(=QPU).
As for the QPE, we calculate the diﬀerence from LDA values in TOTE2.UP at ﬁrst, andadd the diﬀerence to the Shown LDA eigenvalue.

65

Z

13

mkGWIN lmf2 and its I/O Files(QPNT.chk contains irreducible k point for given n1 n2 n3; KPTin1BZ.gwinit.chk contain all
k points in Brillouin Zone).
The purpose of the script mkGWIN lmf2 is to give a template GWinput.tmp. The script
is complicated because of historical reasons. However, its essential is simple; we calls three
executions in this script asecho 0 | lmfgw si
echo 1 | gwinit
echo −100 | qg4gw. We explain each by each.

13.1

echo 0|lmfgw— makes SYMOPS LATTC CLASS NLAindx.
Input ﬁles
• GWIN0 : This is a ﬁle, which contains your supplied n1 n2 n3 when you invoke the script.
This ﬁle is given within the script of mkGWIN lmf2 (as ”here document”).
cat <<EOF >GWIN0n1 n2 n3
$n1 $n2 $n3
cut
4.0 3.0
alpha1
Number of bands
999 99999.0
999 3.0
EOF• ctrl.si : Master input ﬁle of lmf calculation.
Output ﬁles
• LATTC : contains the information of primitive translation vectors, lmxa and konf. See ??
• SYMOPS : The point group operations. See ??
• CLASS : Equivalent atomic positions are called as ’class’. This small ﬁle contains a map
between atomic site and ’class’.
• NLAindx : This ﬁle contains indexes (pvalence , l, a ) for orbitals in the MT. (pvalence is radialfunction index, a is atomic site index). Eigenfunctions are expanded in this order.
• ldima : Number of MTOs for each atomic site. (this is used only from hqpe sc—QSGWmode).
• ves* : not meaningful at this stage• rhoMT* : not meaningful at this stage

13.2

gwinit— Get GWIN V2.tmp and QPNT.tmp
Input ﬁles
• GWIN0 :• LATTC :• SYMOPS :• NLAindx :Output ﬁles
• GWIN V2.tmp : A part of GWinput.tmp• QPNT.tmp : A part of GWinput.tmp• (KPNTin1BZ.gwinit.chk) : check KPNT in the 1BZ.
If SYML exist, gwinit gives also a template QPNTforSYML.tmp suitable for such SYML.
Here SYML specify how to plot the energy band. See explanation for bandplot script.

66

Note that LATTC SYMOPS CLASS NLAindx are overwritten when you execute gw lmfh
because we repeat echo 0|lmf at the head of gw lmfh.

13.3

echo -100|qg4gw— Generate GWinput.tmp
Input ﬁles
• GWIN0 : (copy of GWIN0.tmp by gwinit)• GWIN V2 : (copy of GWIN V2.tmp by gwinit)• QPNT : (copy of QPNT.tmp by gwinit)
Output ﬁles
• GWinput : (this is copied to GWinput.tmp)
This command ”echo -100|qg4gw” is a ﬁle converter from these two ﬁles into GWinput. And
it is copied to GWinput.tmp. (mkGIN lmf keeps GWinput if it exist before you invoke it.).

67

14

gwsc script and its I/O Files

In gwsc, we have a loop of QSGW self-consistency. Look into the gwsc script. In each
iteration, we perform these fortran programs;
NO_MPI=0 #this is used for non-mpi versions of fortran program.
### self-consistent calculation for given Sigma(self-energy) ###run_arg ’---’ $MPI_SIZE $nfpgw/lmf-MPIK

llmf $TARGET### Preparation stage ##############################################
argin=0;
argin=1;
argin=1;

run_arg $argin $NO_MPI $nfpgw
run_arg $argin $NO_MPI $nfpgw
run_arg $argin $MPI_SIZE $nfpgw
run_arg ’---’ $NO_MPI$nfpgw/lmfgw
llmfgw00 $TARGET/qg4gw
lqg4gw#Generate requied q+G v/lmfgw-MPIK llmfgw01 $TARGET/lmf2gw
llmf2gw #reform data for gw### Main stage of gw ################################################
argin=0;
argin=1;
argin=1;

run_arg $argin $NO_MPI$nfpgw
run_arg $argin $NO_MPI $nfpgw
run_arg $argin $NO_MPI $nfpgw/rdata4gw_v2 lrdata4gw_v2 #prepare files/heftetleftet # A file EFERMI for hx0fp0/hchknw
lchknw # A file NW, containing nw## Core part of the self-energy (exchange only) ##argin=3;
argin=3;
argin=3;

run_arg $argin $NO_MPI $nfpgw
run_arg $argin $MPI_SIZE $nfpgw
run_arg $argin $MPI_SIZE $nfpgw/hbasfp0/hvccfp0/hsfp0_sc

lbasC # Product basis generation
lvccC # Coulomb matrix for lbasC
lsxC # Sigma from core1## Valence part of the self-energy Sigma ##argin=0; run_arg $argin $NO_MPI $nfpgw/hbasfp0
argin=0; run_arg $argin $MPI_SIZE $nfpgw /hvccfp0
argin=1; run_arg $argin $MPI_SIZE $nfpgw /hsfp0_sc
argin=11; run_arg $argin $MPI_SIZE $nfpgw /hx0fp0_sc
argin=2; run_arg $argin $MPI_SIZE $nfpgw /hsfp0_sc
argin=0; run_arg $argin $NO_MPI $nfpgw/hqpe_sclbas # Product basis generationlvcc # Coulomb matrix for lbas
lsx # Exchange Sigma
lx0 $lx0_para_option #x0 partlsc #correlation Sigmalqpe #all Sigma are combined.

run arg:
Here a subroutine of bash run arg was used, which is given in ecalj/lm7K; it just invoke
a command with the argument $argin (this is read by read(*,*) in fortran). In cases with
MPI_SIZE/ = 0, mpirun is invoked. Console out put go to l* ﬁles. For example,
argin=2; run_arg $argin $MPI_SIZE $nfpgw /hsfp0_sc

lsc #correlation Sigma

invokes hsfp0 sc with argument ’2’ by mpirun with the -np $MPI SIZE. Console outputs are
written into logﬁles such as lqpe. $nfpgw contains path to the execution binaries.

In the following, We explain input/output ﬁles for each fortran program. Note that “echo0|lmfgw” means invoking lmfgw with argin=0.

68

14.1

echo 0| lmfgw si

See Sec.13.1.

14.2

echo 1| qg4gw

This makes q points, and G vectors for these q. (q was k in previous sections.) Main routineof qg4gw is fpgw/main/qg4gw.m.F and calls fpgw/gwsrc/mkqg.F
Input ﬁles
• GWinput :• LATTC :• SYMOPS :Output ﬁles
• QGpsi : (bin) q and G vector for the eigenfunctions.
• QGcou : (bin) q and G vector for the Coulomb matrix• Q0P : oﬀset-Γ points which are the replacement of the q=0 points. See section??.
• QIBZ : q points in the Irreducible BZ.
• BZDATA : (bin) BZ data for integration (include tetrahedrons if necessary). See e.g.
main/hx0fp0.sc.F and search ”call read BZDATA”, which is a readin routine of this ﬁle deﬁned in rwbzdata.F.
• KPTin1BZ.mkqg.chk : list of q in the 1st BZ for check.
• QBZ : q point in the 1st BZ.
• EPSwklm : Required information for the BZ integration (mainly in order to evaluate theweight in the Γ cell). See Eq.xxx in [?].

14.3

echo 1|lmfgw si

Calculate eigenfunctions, eigenvalues and hψ|HKS |ψi
Input ﬁles
• ctrl.si :• rst.si : (bin) Restart ﬁle of the lmf calculation. It contains all information• sigm.si : (bin) If this exist and• QGpsi,QGcou,Q0P : :• NLAindx : :Output ﬁles
• gwa.si : (bin) atomic data
• gwb.si : (bin) band data
• gw1.si : (bin) hψ|HKS |ψi• gw2.si : (bin) hψ|HKS − Vxc (ntotal )|ψi.
• vxc.si,evec.si : (bin) used in hqpe.sc.m.f as ”v xc” and ”evec”).
vxc.si contains hψ|Vxc (ntotal )|ψi including oﬀ-diagonal part. evec.si contains eigenfunctions.
• normchk.si : norm check (only for check) This is like this> head -20 normchk.si#
IPW
IPW(diag)
Onsite(tot)
Onsite(phi)Total
0.436015
0.805123
0.563972
0.562573
0.999988
0.339134
0.620353
0.660515
0.656881
0.999649
0.339133
0.620353
0.660516
0.656882
0.999649
0.339133
0.620353
0.660516
0.656882
0.999649
0.507738
0.648515
0.492040
0.487673
0.999778...
This check is sometimes important for debugging and to determine the cutoﬀ parameter
QGcut psi . The ﬁrst line (corresponding to 1st band of 1st q point) means that total normalization almost unity = 0.999988 = 0.436015 + 0.563972. Because we expand the MTO
by IPW, the normalization is a bit diﬀerent from unity, especially for higher bands. You can
see that it get closer to unity for larger QGcut psi, though it does not reach to unity becauseof some contribution of the higher angular momentum contribution within MT. [Values of
Onsite(phi) are not correctly shown in the case when you use local orbital.]
Due to historical reason, data in vxc.si and exec.si and others contains duplicated data.

69

14.4

lmf2gw

All the required information are stored into DATA4GW V2 and CphiGeig.
Input ﬁles
• gwa.si :• gwb.si :• gw1.si :• gw2.si :• Q0P :• CLASS :• NLAindx :Output ﬁles
• DATA4GW V2 : (bin) Main data for GW calculations.
I/O of DATA4GW V2 is controlled by gwinput.f, which contains detailed information.
• CphiGeig : (bin) Eigenfunctions for GW calculations.
• VXCFP.chk : Eigenvalue and Vxc check (only used for check)
It is like this;### LDA exchange correlation ####qvec
ikp iband
eigen
0.0000 0.0000 0.0000 1 1-0.68505346
0.0000 0.0000 0.0000 1 2
0.19292662
0.0000 0.0000 0.0000 1 3
0.19292763
0.0000 0.0000 0.0000 1 4
0.19292777...

VXC(ntotal)-0.91850436-0.99853478-0.99853469-0.99853461

VXC(nvalence)
0.00000000
0.00000000
0.00000000
0.00000000

eigen(eV)-9.32070032
2.62492096
2.62493477
2.62493664

VXC(ntotal)(eV)-12.49698668-13.58586453-13.58586334-13.58586222

Here VXC(nvalence) is not used now. The eigenvalue in eigen is in Ry.

— This is the end of the preparation stage. —From here, the main stage.

14.5

rdata4gw v2— Read DATA4GW V2 and some ﬁles, and decompose it into ﬁles required in the following
GW steps. (checked! dec2014)
Input ﬁles
• GWinput :• DATA4GW V2 :• CphiGeig :• QGpsi :• QGcou :• Q0P :• QIBZ :• SYMOPS : points group operations.
Output ﬁles
• hbe.d : data size
• Core ibas* l*.chk : core eigenfunctions just for check.
• LMTO : basic date for the crystal.
• EValue : (bin) valence eigen value• ECORE : core data and core eigenvalues• CPHI : (bin) Coeﬃcients of eigenfunctions as for the atomic-like argumentation waves in
MTs’.
• GEIG : (bin) Coeﬃcients of eigenfunctions as for IPW.
• PHIVC : (bin) All the radial functions.
• @MNLA CPHI : index set for CPHI. This is not refereed just a check write.
• @MNLA core : index set for core. This is not refereed just a check write.
LDA• VXCFP : (bin) this is for diagonal elements of Vxc(ntotal ).
• PPOVLI.* : (bin) Overlap matrix of IPW. xxxxxxxxxx

70

VXC(nvalence)
0.00000000
0.00000000
0.00000000
0.00000000

• PPOVLG.* : (bin) PPOVLG Overlap matrix xxxxxxxxxxoxf IPW. not exactly the theoverlap matrix. see around line 500 in rdata4gw.m.f• PPOVL0 : (bin) xxxxxxxxxxxxx• HVCCIN : (bin) Required inputs for hvccfp0. Information in this ﬁles.
• NQIBZ : q point info. Only used for parallel test mode.
• normchk.dia : Norm check. These numbers should be almost the same as those in normchk.si
These ﬁles are input for the following steps. The name of ﬁle fooU means that it relates to
up-spin. We have fooD ﬁles in the case of spin-polarized calculation with nspin=2.

14.6

echo 1|heftet— Get the Fermi energy EFERMI by tetrahedron method. It is used in hx0fp0.
Input ﬁles
• EVU :
• BZDATA :• GWinput :• ECORE : (dummy)• SYMOPS : (dummy)• LMTO :• hbe.d :Output ﬁles
• EFERMI : contains Fermi energy given by the tetrahedron method. It is used in hx0fp0
but not in hsfp0.
• DOSACC.lda,DOSACC2.lda : They are lists of the all the eigenvalues from the bottom.
DOSACC2.lda is a list to show only the un-degenerated eigenvalues. They are just check write.
But it is an indicator for you to determine esmr in GWinput.

14.7

hchknw— Calculate the required number of ω points along real axis.
This NW is not essentially used in gw lmfh (but required as a dummy ﬁle). Only used in
gw lmf.
Input ﬁles
• BZDATA :• GWinput :• ECORE : (dummy)• SYMOPS : (dummy)
Output ﬁles
• NW : contains number of ω points.

14.8

echo 3|hbasfp0— Make product basis.
Mode 3 is for the core states. It generate a product basis on each MT suitable to expand
to calculate the exchange part due to core. See explanations for the input ﬁle of GWinput.
Input ﬁles
• LMTO :• PHIVC :• GWinput :Output ﬁles
• BASFP* : (bin) Product basis functionsR• PPBRD V2 * : (bin) Radial integrals on each MT, symbolically written as φ(r)φ(r)B(r)dr• PHIV.chk : Valence radial functions (for check).

71

14.9

echo 0| hvccfp0— Calculate the Coulomb matrix in the Mixed basis
Input ﬁles
• HVCCIN :• Q0P :
• BASFP* :Output ﬁles
• VCCFP : The Coulomb matrix expanded in the mixed basis• Mix0vec : This is used only for dielectric-constant calculation (mode 2 or 3 of hx0fp0). This
contains the expansion of the plane wave exp(iqr) in the mixed basis. See Usuda’s note.

14.10

echo 3|hsfp0— Exchange part of the self-energy for the core
Input ﬁles
• GWIN V2,LMTO,ECORE :• CLASS :• hbe.d :• Q0P :• PPBRD V2 * :• CPHI :• GEIG :• VCCFP :• PPOVL :Output ﬁles
• SEXcoreU : The core part of the exchange self-energy for q and band index speciﬁed in<QPNT>. See ??.

14.11

echo 0|hsfp0— Make product basis for the valence part.

14.12

echo 1|hsfp0— Exchange part of the self-energy for the valence part.
Output ﬁles
• XCU : The LDA exchange self-energy for q and band index speciﬁed in <QPNT>. See ??.
• SEXU : The valence part of the exchange self-energy for q and band index speciﬁed in<QPNT>. See ??.

14.13

echo 11|hx0fp0— Screened Coulomb interaction W (Sergey mode)
Input ﬁles
• GWinput, LMTO, ECORE, EVU :• NW : dummy• hbe.d :• Q0P :• PPBRD V2 * :• CPHI,GEIG :• PPOVL :• VCCFP :• ANFcond : (optional) This ﬁle is to specify antiferro condition. This should not exist forother cases. This ﬁle should be given by hand.
Output ﬁles
• WV.d : size of the dielectric function• WVR : (bin) W − v in the expansion of mixed basis along the real axis

72

• WVI : (bin) W − v in the expansion of mixed basis along the imaginary axis

14.14

echo 12|hsfp0— Correlation part of the self-energy(Sergey mode)
Input ﬁles
• GWinput, LMTO, ECORE, SYMOPS : These are readin by genallcf v3.
• CLASS, hbe.d, EVU, Q0P :• PPBRD V2 * :R
Radial integrals on each MT, symbolically written as φ(r)φ(r)B(r)dr. These are generated by hbasfp0.
• CPHI,GEIG :• PPOVL :• WV.d, WVR, WVI :Output ﬁles
• SECU : The correlation part of the self-energy for q and band index speciﬁed in <QPNT>.
See ??.

14.15

echo 0|hqpe— Summarize the output
Input ﬁles
• SEXcoreU,XCU,SEXU,SECU : See ??.
Output ﬁles
• QPU : The QP energies and related value summary in human interface. See ??.
• TOTE : The detailed values of the QP energies. See ??.
• TOTE2 : The detailed values of the QP energy. See ??. This is used for bndplot.
NOTE: For example, if you do {echo 4$|$hqpe}{hqpe}, it just shift zero level of QPE, sothat 4th line (counted from top) eigenvalue (in QPU) is to be zero.

73

15

Check list for convergence on GW calculations

Results could be dependent on cutoﬀ parameters in GWinput. In my opinion, generally speaking, it is so easy to have convergence more than ∼ 0.1eV for band gap...
• Number of k points n1n2n3 .
Probably 4 × 4 × 4 (or 6 × 6 × 6) are reasonable choice for insulator in the case of two
atoms such as GaAs. In other words, ‘”used periodic cell volume” = 4 × 4 × 4× (# ofatoms) × (Volume per atom)”. For example, for 2 × 23 = 16 atoms per cell, we can use
2 × 2 × 2 instead of 4 × 4 × 4 (this is the same as in the case of LDA).
For metals, 12 × 12 × 12 or more per atom may be required. But it depends on case bycase.
• MTO’s and APWs. Basis for eigenfunctions.
If we like to get “best converged results, we may need to use large enough MTO’s. (Ithink the default setting is reasonable; but there is a room to change MTO settings in
ctrl.*.). But, in cases, we can not use large APW cutoﬀ (pwemax) more than 3.0 ∼ 4.0
Ry because of poor linear dependency of basis set (calculation in LDA level fails). In
such a case, we need to use “smaller MT radius R=”. Then we may need semi-core as
local orbital when “its spillout in the outside of MTs is too large”.
In GWinput, we can set the number of unoccupied states which you take into account byemax chi0 , emax sigm , nband chi0 , and nband sigm . But we usually unset them
except emax sigm for gwsc.
(We may need a kind of completeness of the basis set; the completeness could be important from the view of ‘Coulomb hole’ picture.)
NOTE: Current PMT-QSGW method [?] expand the static version of self-energy of
QSGW just in the basis of MTO’s (no APWs). Thus the expansion can be unsatisfactory(it depends on case by case, and required convergence). In such a case, we inevitably
have to use empty spheres (MTO’s) which is for empty region.
• Cores.
We usually use all cores as core1 (exchange only core). If necessary, it is better to treat
shallow cores by the local orbitals (such cores are treated as valence).
—followings are obsolate —If we treat cores by core2 (not only for exchange, but also in dielectric functions),
we have to be careful about the core wave orthogonalization with respect to valence
eigenfuncitons; this is a little complicated; probably it is better not to use core2).
(This is related to CoreOrth (only for core2). If it aﬀects so much, D function might
be too poor due to the poor orthogonality condition between core and valence.) (NOTE:CoreOrth is not maintained recently)————• QpGcut psi IPW cutoﬀ to expand eigenfuncitons in the interstitial region. We usuallyuse QpGcut\_psi 4.0. Usually not so bad. Larger is better but expensive. You may test
calculations with QpGcut\_psi 3.0, and how much diﬀerence of results.
• QpGcut cou IPW cutoﬀ to expand the Coulomb interaction in the interstitial region(Mixed product basis(MPB) consist of this IPW and PB). interstitial region. We usuallyuse QpGcut\_cou 3.0. Usually not so bad. Larger is better but expensive. You may test
calculations with QpGcut\_cou 2.5, and how much diﬀerence of results.
• Product basis section.
At least, lcutmx=4 will be necessary for atoms with d electrons. (but lcutmx=2 lookreasonable for oxygen, lcutmx=1 for hydrogen. Need check).
• esmr .
In our experience, esmr=0.003000 (default) is reasonable. But there is a room to checkstability on it for metals. In principle, for larger n1n2n3, we can use smaller esmr.

74• dw, omg cniw

It will be worth to try to check how much the results changed due to them. But usuallydw=0.005, omg_c=0.04 is not so bad a choice. As for niw=10 seems to be not so bad
usually, but it is safer to check the convergence on it (test cases with niw=6,10,12,16).
• deltaw . ∼ 0.01 a.u. will be not so bad. See two Z values shown in SXCU. It is betterto try to check how about the dependence on this.
• chi regqbz off . (on is default). If oﬀ, we use oﬀ-Gamma mesh (Gamma point is between mesh points) for dielectric functions when we perform GW or QSGW. “chi regqbzoﬀ” may accelerate convergence on number of k points.
• PWMODE=1 in ctrl This used ﬁxed number of G vector determined at q=0

75

16

Linear response calculations

With these scripts for linear response calculations, eps*, we can calculate q-dependent dielectric funciton ǫ(ω, q) (and v, W ) (and χ for spin ﬂuctuation). But (because of numericalreason), we can not use q = 0 limit. (if |q| is too small, we have numerical problem, zero
divided by zero, because we have not implemented the version to use q = 0.)• eps lmfh
Dielectric function epsilon with local ﬁeld correction. Expensive calculation (we mayneed to reduce number of wing parts in future...).
• epsPP lmfh
epsilon without local ﬁeld correction. 1 − heiqr |v|eiqr iheiqr | χ0 |eiqr i• epsPP lmfh chipm

For spin susceptibility. This essentially calculate non-interacting spin susceptibility.
Then it is used for the calculation of full spin susceptibility with util/calj_*.F programs(small quick programs). See spin wave paper. See spin susceptibility section Sec.??.
• (not maintained now; we will recover this)eps lmfh chipm
This gives full non-interacting spin susceptibility. Testing. We have to determine U(Stoner I) for the determination of full spin susceptibility. TDLDA? or so?
• (This is old mode --- not maintained) epsPP lmfh chipm q
For spin susceptibility, spin susceptibility heiqr |χ(q, ω)|eiqr i In this script, You have toassign that isp=1 is majority, isp=2 is minority. This is with long wave approximation.
———————• We use the histogram method (the Hilbert transformation method); we ﬁrst calculate its
imaginary parts with the tetrahedron technique for dielectric functions. Then we get its realpart by the Hilbert transformation.
You need to choose HisBin dw,HisBin ratio . The width of histogram bins are getting larger
when omega gets larger. dw is the size of histogram-bin width at omega=0. At omega=omg c,
its width gets twiced.
To plot dielectric function with reasonable resolution, it might be better to set dw 0.001
and omg_c 0.1 for example. You may have to choose small enough omega for spin wave mode
as 0.001 Ry (Or smaller). omg c is given like 0.05 Ry or so. But sometimes it can be like 1Ry.
• epsPP lmfh only calculation an matrix element of dielectric funciton for exp(iqr). Thus very
faster than eps lmfh mode. It uses a a special product basis set for cases without inversion(problem is in how to expand exp(iqr) in the MPB; the product basis is not from phi andphidot, but from spherical Bessel functions).
In *_lmfh_* modes( I now use little for *_lmf_* modes), you can use small enough delta.
Use small enough delta (=-1e-8 a.u.) for spin wave modes (also you can use it for
dielectric function and GW). This is necessary because pole is too smeared
if you use larger delta.

16.1

eps lmfh, epsPP lmfh: the dielectric functions

You can invoke the script, e.g. as ”eps lmfh si”.
—————Specify q point in <QforEPS> or so. Mesh for ω is speciﬁed by dw, omg c .
The obtained data are in EPS*.dat and EPS*.nlfc.dat. EPS*.nlfc.dat contains the result
without local-ﬁeld correction EPS*.dat contains the result with local-ﬁeld correction (this is
generated only for eps_lmfh. Both of them contains
q(1:3), ω, Re(ǫ) Im(ǫ), Re(1/ǫ), In(1/ǫ)
in each line.
• This code works OK only for q is near 0. Be careful for q → 0 limit. Too small q cangive strange spectrum at high energy (real part is aﬀected by it)76

Because q → 0 gives too large cancellation eﬀects (the denominator and numerator
go to zero—it means we need very accurate orthogonalization between occupied andunoccupied states). This is a kind of disadvantage of our method (though there is anadvantage— our code can calculate dielectric function even for metal as far as you use
large enough number of k point.)• The calculate of dielectric functions usually requires so many k point. For example, forSi, n1 n2 n3 = 4 4 4 is too small. It gives too large dielectric constants ∼ 19.4 though
the converged value should be ∼ 13. (we need 10x10x10 or more like 20x20x20 for somereasonable results). For GaAs, we observed that reasonable ǫ(ω) requires rather large
number of q points like 15x15x15 or 20x20x20 for n1n2n3 . This is too time-consuming
to get result (but you can use “very small product basis”(just sp polarization for thispurpose; it makes speed up so much). Or, you can calculate ”ǫ(ω) without LFC”. See
section for eps PP lmfh.
• Core orthogonalization problem (only when core2 is used)—- CoreOrth is not maintained recently — CoreOrth gives so serious eﬀect for ǫ(ω),
if you include some cores as ”core2” in the product basis setting. (This means that you
includes transitions from ”core2 to valence” in the calculation of ǫ(ω)).
Then you have to use ”CoreOrth on”. Without it, you will have rather large imaginary
part at rather high energy Such transitions from core to higher valence bands is artiﬁcial
due to the incomplete orthogonality between core and the higher bands. However, shallower d semi-core might be deformed too much by this option. Try to plot Core *.chk
ﬁles, which contains core radial functions. Anyway, it is better to treat shallow core as
valence by “local orbital”.

16.2

epsPP lmfh: the dielectric function(No LFC— faster)

You can calculate ǫ without LFC by epsPP lmfh. It is very faster than eps lmfh.
To calculate ǫ(q, ω) without LFC accurately, the best basis set for the expansion of the
Coulomb matrix within MT is apparently not the product basis, but the Bessel functions
corresponding to the plane waves exp(iqr). We use such a basis in this mode. However,
our experience shows that the changes are little even with the usual product basis (we don’t
describe this here).

16.3

How to calculate correct dielectric funciton?

(this subsection is essentially OK... but need to clean it up. dec2014)
There are prolems to calculate correct epsilon.
At first, we talk about epsPP_lmfh, which is No LFC. Main problem are----------------1.Convergence for number of k point(specified by n1n2n3).
Roughly speaking, 20x20x20 is required for not-so-bad results for Fe and Ni.
It is better to do 30x30x30 to see convergence check.
However, in the case of ZB-MnAs (maybe because of simple structure around Ef),
it requires less q points.
figs are for GaAs.
fig001: n1n2n3 convergence for Chi_RegQbz = on case.
fig002: n1n2n3 convergence for Chi_RegQbz = off case.
(Chi_RegQbz in explained in General section in this manual).
As you see, k points convergence looks a little better in Chi_RegQbz=off(mesh not including gamma). However a little ploblem is that its thereshold around
0.5eV is too high and slowly changing.

77

fig003: Alouanis’(from Arnaud) vs. ‘‘Chi_RegQbz = on’’ vs. ‘‘Chi_RegQbz = off’’
As you see, the threshold of the Red line (20x20x20 Chi_RegQbz=on) and Alouani’s
are almost the same, but the red line is too oscilating at the low energy part.
On the other hand, ‘‘Chi_RegQbz = off’’ in Green broken line is not so satisfactoryat the low energy part.
fig.gas_eps_kconf.pdf shows the convergence behevior of epsilon for

2.$q \to 0$ convergence (this is related to whether Chi_RegQbz=on or off).
If you use very small q like q=0.001 is GaAs, it can cause a problem.
Use q=0.01 or larger (maybe q=0.02 or more is safer).
Very small q can give numerical error for high-energy region.
In fig004, we show the high energy tail part of Im $\epsilon(\omega)$ for GaAs case.
At q=0.01 (this means q= 2*pi/alat * (0 0 0.01)), the imaginary partis a little too large . Less than 80eV, q=0.02 gives good results when compared withother high q results, though it still has noise above 80eV.
In fig005, I showed the same results compared with Alouani’s (his is up to 40eV).
Both gives rather good agreements. As you see, q=0.06 or above might be necessary
to get reasonable convergence for high energy part abouve 40eV.
We have to be careful for this poorness in high energy part--- it may effectlow-energy Re[$\epsion$] through KK relation. However this can be very smallehough.
In fig.gas_eps_qconv.jpg, we checked the convergence of eps (\omega=0,q) for q \to 0.
As you see, it gives convergence, however, q=0.01 is a little out ofcurve---this should be because of the poorness in the high energy part.
so q=0.02 or q=0.03 is safer, and you can get eps within 1 percenr accuracy.
3. Including Core for dielectric constant is dangerous.
It can cause very poor results if you include core part in GWinput.
You need to include core just as valence (with local orbital).
In fig008, we showed core effects. It starts from \approx 16eV
(this is core to conduction transition).
fig007 showd the check about the q point dependence---even with large q,
it would not change.
These shows that the core excitation can have larger energy range.
This is in contrast to the valence case(then the most of excitaion is limited to less than 10eV).
We have to be careful for such high-energy exciation... The LMTO basis might
be not so good for high energy.
4. basis set.
Use QpGcut_psi \approx 3.0 a.u. or so (as same as GW calculation).
In the case of epsPP* mode,
QpGcut_cou can be very small--- In our codes now,ngc>=1 should be for all q vector shown in lqg4gw02 (output of echo 2|qg4gw).
[In principle, it should be only for the q vector for which we calculate epsilon.
But there is a technical poorness in our code--(maybe) a problem here; the plane-wave part of the eigenfunction generated
in lmfgw is not correctly passed to lmf2gw when ngc=0].

78

-- eps_lmfh: including LFC ---------------------------------To include eps with LFC, do eps_lmfh.
But lcutmx=2 seems to be good enough to get 0.5 percent error (maybe better than this).
Test it 10x10x10 or so. (I need to repeat if necessary).
Further you can use smaller QpGcut_cou like 2.2 or so,with rather smaller product basis (up to p timed d, not including f).
Note: epsPP_lmfh is designed to use good basis to calculate eps
without LFC. This is usually in agreement with what you obtained by eps_lmfh;
however it can give slight difference when you use small product basis.

---Summary -------------------So in conclusion, I think a best way to do is
1. set q=0.02 [q=2pi/alat(0 0 0.02)] or so for GaAs case.
If you want to check, do q=0.03 and q=0.06 also.
‘‘Chi_RegQbz = off’’ is better for matrials like GaAs with direct gap.
2. You can use small QpGcut_cou but all ngc should be one or more.
3. As for the Product basis setting in epsPP* scripts, only
lcutmx and tolerance (this can be like 0.001 or so) are relevant.
E.g. set lcutmx=4 or so.
4. Do nk=20 18 16 and take interpolarion to determine eps(omega=0, q=0).
5. To get eps with LFC, set QpGcut_cut as xxx, and set lcutmx=2 where(occupied sp) \timex (unoccupied spd) are included.
But correct EPS*.nolfc.d is rather from epsPP_lmfh script.

gaseps0 ChiRegQbz=on (default)
50"gasepstest0202020/EPS0001.nolfc.d" using ($4*rydb):($6)"gasepstest0181818/EPS0001.nolfc.d" using ($4*rydb):($6)"gasepstest0161616/EPS0001.nolfc.d" using ($4*rydb):($6)"gasepstest0101010/EPS0001.nolfc.d" using ($4*rydb):($6)

45
40
35

---plot------

30
25
20
15
10
5
0
0

2

4

6

79

8

10

12

14

ﬁg001
gaseps0 ChiRegQbz=off
50"gaseps0bz2161616/EPS0001.nolfc.d" using ($4*rydb):($6)"gaseps0bz2181818/EPS0001.nolfc.d" using ($4*rydb):($6)"gaseps0bz2202020/EPS0001.nolfc.d" using ($4*rydb):($6)

45
40
35

---plot------

30
25
20
15
10
5
0
0

2

4

6

8

10

12

14

ﬁg002
gaseps0 ChiRegQbz=on(bz2) vs. off vs. Alouani’s
50"gasepstest0202020/EPS0001.nolfc.d" using ($4*rydb):($6)"gaseps0bz2202020/EPS0001.nolfc.d" using ($4*rydb):($6)"e2ofGaAsRPAwithoutLFPAW.txt" using ($1):($2)

45
40
35

---plot------

30
25
20
15
10
5
0
0

2

4

6

ﬁg003

80

8

10

12

14

gaseps0bz21616164 High energy part (numerical error)
0.4

gaseps0bz2161616 q=0.01
q=0.02q=0.03q=0.06q=0.12q=0.180.35

0.3---plot------

0.25

0.2

0.15

0.1

0.05

0
0

20

40

60

80

100

120

140

ﬁg004
ChinoRBZ=T 161616 noGd3d compared with Alouai ---Tail part0.4
q=0.02q=0.03q=0.06q=0.12
Alouani’s

0.35

0.3

---plot------

0.25

0.2

0.15

0.1

0.05

0
0

20

40

60

ﬁg005

81

80

100

120

140

/home/takao/DATA/gaseps0bz2161616 4 6
0.4
alouanis161616high q=0.02
161616 q=0.12
0.35

0.3

---plot------

0.25

0.2

0.15

0.1

0.05

0
20

40

60

80

100

120

ﬁg007
Effects of Gd3d
5gaseps0bz21616162noGa3d q=0.02
gaseps0bz21616162
q=0.02

4---plot------

3

2

1

0
0

5

10

15

20

ﬁg008

82

25

30

35

40

45

50

17

Utility

1. structure tool viewvesta,ctrl2vasp,vasp2ctrl
2. Other calculation examples A recent paper by Deguch et al [12] contains a link to asample package.
3. cif to PROCAR./cif2cell --vasp-cartesian --vasp-format=5 cifs/BaTiO3_cubic.cif83Index
Mix0vec, 61
gwinit, 55
hbasfp0, 60, 61
NLAindx, 55, 58, 59
hchknw, 60
NQIBZ, 60
heftet, 60
NW, 60, 61
hqpe, 62
PHIV.chk, 60
hsfp0, 61, 62
PHIVC, 59, 60
hvccfp0, 61
PPBRD V2 *, 60–62
PPOVL, 61, 62
hx0fp0, 61
PPOVL0, 60
lmf2gw, 59
PPOVLG.*, 60
lmfgw00, 55, 58
PPOVLI.*, 59
lmfgw01, 58
Q0P, 26, 58, 59, 61
qg4gw, 56
QBZ, 58
qg4gw00, 58
QGcou, 26, 58, 59
rdata4gw v2, 59
MNLA core, 59
QGpsi, 26, 58, 59
MNLA CPHI, 59
QGpsi,QGcou,Q0P, 58<PRODUCT BASIS> section, 48
QIBZ, 26, 58, 59<QPNT> section, 45
QPNT, 56(KPNTin1BZ.gwinit.chk), 55
QPNT.tmp, 55
ANFcond, 61
QPU, 52, 62
ANFcond (we can skip here since we do not check this option
Need fix this if necessary., 50
SECU,now.
53, 62
BASFP*, 60, 61
SEXU, 52, 61
BZDATA, 26, 58, 60
SEXcoreU, 52, 61
CLASS, 27, 55, 59, 61
SEXcoreU,XCU,SEXU,SECU, 62
CLASS, hbe.d, EVU, Q0P, 62
SYMOPS, 26, 55, 58–60
CPHI, 59, 61
TOTE, 62
CPHI,GEIG, 61, 62
TOTE.UP (TOTE.DN), 53
Core ibas* l*.chk, 53, 59
TOTE2, 62
CphiGeig, 26, 59
TOTE2.UP (TOTE2.DN), 53
DATA4GW V2, 26, 59
VCCFP, 61
DOSACC.lda, 53
VXCFP, 59
DOSACC.lda,DOSACC2.lda, 60
VXCFP.chk, 53, 59
DOSACC2.lda, 53
WV.d, 61ECORE, 59, 60
WV.d, WVR, WVI, 62
EFERMI, 60
WVI, 62
EPSwklm, 58
WVR, 61
EVU, 60
XCU, 52, 61
EValue, 59
ctrl.si, 40, 55, 58
GEIG, 59, 61
gw1.si, 58, 59
GWIN0, 55, 56
gw2.si, 58, 59
GWIN V2, 56
gwa.si, 58, 59
GWIN V2,LMTO,ECORE, 61
gwb.si, 58, 59
GWIN V2.tmp, 55
hbe.d, 59–61
GWinput, 26, 56, 58–60
ldima, 55
GWinput, LMTO, ECORE, EVU, 61
normchk.dia, 60
GWinput, LMTO, ECORE, SYMOPS, 62
normchk.si, 58
GWinput.tmp, 40
rhoMT*, 55
General section, 42
rst.si, 58
HAMindex, 26
set QPNT for eps mode (QforEPS section), 47
HVCCIN, 60, 61
sigm.si, 58
KPTin1BZ.mkqg.chk, 58
ves*, 55LATTC, 55, 58
vxc.si,evec.si, 58
LMTO, 59, 60

84

References[1] Takao Kotani, Hiori Kino, and Hisazumu Akai. Formulation of the augmented plane-wave and muﬃntin orbital method. Journal of the Physical Society of Japan, 84(3):034702, March 2015.
[2] Takao Kotani. Quasiparticle self-consistent GW method based on the augmented plane-wave
and muﬃn-tin orbital method. J. Phys. Soc. Jpn., 83(9):094711 [11 Pages], September 2014.
WOS:000340822100029.
[3] Takao Kotani, Mark van Schilfgaarde, and Sergey V. Faleev. Quasiparticle self-consistent gw method:
A basis for the independent-particle approximation. Physical Review B, 76(16):165106, 2007.
[4] Fabien Bruneval and Matteo Gatti. Quasiparticle self-consistent GW method for the spectral properties of complex materials. In Cristiana Di Valentin, Silvana Botti, and Matteo Cococcioni, editors,
First Principles Approaches to Spectroscopic Properties of Complex Materials, volume 347, pages
99–135. Springer Berlin Heidelberg, Berlin, Heidelberg, 2014.
[5] T. Kotani and M. van Schilfgaarde. All-electron GW approximation with the mixed basis expansion
based on the full-potential LMTO method. Solid State Communications, 121(9-10):461–465, 2002.
[6] Christoph Friedrich, Markus Betzinger, Martin Schlipf, Stefan Blügel, and Arno Schindlmayr. Hybrid
functionals and GW approximation in the FLAPW method. Journal of physics. Condensed matter :an Institute of Physics journal, 24(29):293201, July 2012.
[7] Takao Kotani and Mark van Schilfgaarde. Fusion of the LAPW and LMTO methods: The augmented
plane wave plus muﬃn-tin orbital method. Physical Review B, 81(12):125117 (2010) [5 pages], March
2010.
[8] José M. Soler and Arthur R. Williams. Simple formula for the atomic forces in the augmented-planewave method. Phys. Rev. B, 40(3):1560–1564, Jul 1989.
[9] José M Soler and Arthur R Williams. No Title. Phys. Rev. B, 42:9728, 1990.
[10] Takao Kotani, Hiori Kino, and Hisazumu Akai. Formulation of the Augmented Plane-Wave andMuﬃn-Tin Orbital Method. Journal of the Physical Society of Japan, 84(3):034702, March 2015.
[11] J. Rath and A.J. Freeman. Generalized magnetic susceptibilities in metals: Application of the analytic
tetrahedron linear energy method to Sc. Phys. Rev. B, 11:2109, 1975.
[12] D Deguch, K Sato, H Kino, and T Kotani. Accurate energy bands calculated by the hybrid quasiparticle self-consistent GW method implemented in the ecalj package. Journal of the Physical Society ofJapan, page accepted, 2016.

85