HiLAPWVersion 1.4— BASIC EDITION —Tamio OGUCHI
Osaka University, ISIR

1Introduction

In these days, band calculation is often called first-principles or ab initio calculation. Thismeans that a theoretical study on the basis of band calculation becomes calculating not
only band structure but also several physical quantities directly from first principles. There
have been three important issues in progress of the first-principles calculation. The first is
the establishment of the basic theory, typically density-functional and related theory. The
density-functional theory provides a practical approach of one-electron theory through localdensity approximation (LDA) and more recently generalized gradient approximation (GGA).
The second is the highly-developed band-structure calculation methods which make possiblea realistic calculation for a very complicated material system with reasonable computer resource. In the calculation, it is almost inevitable to assume numerical approximations andmodeling in addition to the basic ones such as the Born-Oppenheimer approximation and
LDA. Both kinds of the approximation are often used without making any distinction andmixed up in the everyday calculations. In order to discuss physics and validity of the basic
approximations, it is indispensable to minimize the numerical ones and at least to evaluate
their limitation and errors. Further development of the numerical techniques and algorithms
is highly desired for large-scale systems and complex phenomena to be investigated in recent
materials science. The third is rapid progress in high performance computing. A desktop
computer in every researcher’s oﬃce has computing power comparable to a computer called”supercomputer” ten years ago. (Concerning the first and second issues, beginner-course
students are recommended to study first fundamentals of the electronic structure theory [1].)Under these circumstances, a computer code for the first-principles calculation becomesa gigantic package including a lot of state-of-the-art techniques and related know-how. Usefulness and wide applicability of the first-principles calculation in materials science research
have been acknowledged. However, there has never been a general application code available
until quite recently as in quantum chemistry field. In this decade, a number of researchgroups have developed first-principles calculation codes and, in many cases, made them open
to public.
HiLAPW introduced in this article is a first-principles calculation package originally developed by a theory group in Hiroshima University. The current developments of the package are1

being made in collaboration between Osaka University and Hiroshima University. In section
2, the basic specifications of HiLAPW are explained. The package contains several executables,
which perform part of the first-principles calculation. The executables are linked via data
files, as described in section 3. In section 4, the way to install the package is given. Finally,some application examples are shown in the following sections.

2

Basic Specifications

The HiLAPW program package is designed to perform first-principles calculations based on
the density functional theory. Main features include the followings.

2.1

Linearized Augmented Plane Wave Basis Functions

The linearized augmented plane wave (LAPW) functions have high variational degrees of freedom and practically good eﬃciency without using pseudopotentials. Original non-linearized
APW was proposed by Slater [2] with muﬃn-tin potential approximation.
This program utilizes LAPW basis functions by Soler and Williams [5, 6], which areslightly diﬀerent in a practical sense from the original ones by Andersen [3] and Koelling andArbman [4].
To remove a possible appearance of the so-called ghost band, the radial energy-derivative
function Ṙl is forced to be exactly orthogonalized to the highest core state Rcl with keeping
the orthogonalization to the corresponding radial function Rl [7] as
˜ = Ṙ + αR + βR
Ṙ
l
l
cl
l

(1)

where the coeﬃcients α and β are determined by the orthogonal conditionsα=−⟨Rcl |Ṙl ⟩⟨Rcl |Rcl ⟩ − |⟨Rcl |Rl ⟩|2β = −α⟨Rcl |Rl ⟩

2.2

(2)
(3)

Scalar Relativistic Treatment

When the atomic number increases the relativistic eﬀects cannot be neglected. A proper
treatment of the relativity can be made with use of Dirac equation instead of Schrödinger
one. In this package, the scalar relativistic treatment is used including all the relativistic
eﬀects except the spin-orbit coupling (SOC) [8]. Additional inclusion of SOC in a secondvariation process provides an approximated but accurate description equivalent to solving
Dirac equation.

2.3

Self-Consistent Field Method

The Kohn-Sham equations should be determined in a self-consistent-field (SCF) manner andusually solved by an iterative scheme. The essential part of the iterative scheme is how toguess input electron density for next iteration from the input and output electron density
distribution functions in the previous iterations. Simple mixing method and extended Anderson method [9] can be employed depending on the SCF convergency behavior in diﬀerent
material systems.

2.4

Full Potential Method

The muﬃn-tin potential approximation assumed in the original APW and LAPW methods
often leads to a less accurate or unphysical solution in cases of sparsely packing or lowsymmetry systems. The full-potential method [10] assumes no shape approximation for theelectron density and potential functions and may be applicable to any crystal symmetryeven surface and interface systems. Note that the muﬃn-tin sphere is still assumed in the
full-potential method but only divides the regions where a plane-wave or spherical-wave
representation is adopted to express the electron density and potential functions. The final
results such as the total energy are unchanged by a diﬀerent choice of the muﬃn-tin sphereradius with well converged parameters for the representations.

2.5

Atomic Forces

In pseudopotential methods with the plane wave basis, atomic forces, which are defined as aminus gradient of the total energy with respect to the atomic coordinates, are simply givenby the Hellmann-Feynman theorem, leading to great success in the Car-Parrinello molecular dynamics simulations. In the LAPW method, where the augmentation functions have
explicit dependence of the atomic coordinates, the so-called Pulay term must be evaluatedadditionally to obtain the atomic forces accurately. In the package, Soler and Williams’ formulation [5, 6] is used for the atomic forces. With the forces, structural optimization andphonon mode calculation can be easily carried out.

2.6

Tetrahedron Integration Scheme with the Second-Order Corrections

Integration within the Brillouin zone is frequently performed to compute any of the physical quantities such as the electron density and the density of states. Blöchl’s tetrahedron
scheme [11] is utilized in the code. The other integration schemes with the occupation givenby the error function or the Fermi distribution function can be also used optionally.

2.7

Irreducible Representation of the Eigenstates

The group theory tells that each one-electron eigestate belongs to an irreducible representation of the k group. With the irreducible representation, one can understand how the
eigenstate is transformed by the symmetry operation (rotation or mirror). For example, thematrix element of the electric dipole with respect to the initial and final states gives theselection rule for photo-excitation spectra.

3

Executables

The HiLAPW program package is composed of kernel executables such as xsets and xlapw,
and some other optional ones. All the executables can be made with use of makefile attachedwith the package and FORTRAN90 compiler on user’s computer. List of the basic executables
and related data files are summarized in Table 1.
Table 1: Basic executables and related data files in HiLAPW.
executablerole(parallel ver.)xsets
generation of initial dataxlapw
(xplapw)

SCF FLAPW calculation

input files

output files

atomdata∗
atomdens∗
spgrdata alt∗
sets.inwavin
lapw.inlapwadv.in∗∗wavoutsets.out

wavout
ekn
dis
foa
tau
ten
lapw.out
xdoss
DOS calculation
wavin
pdos
(xpdoss)doss.inpdos.index
doss.out
xnewa
k-point update
wavin
wavoutnewa.out
xspin
addition of spin degree
wavin
wavoutspin.out
xsymm
symmetry calculation
wavin
ekl
cml
symm.out
xrept
reorder of the eigenstates ekl
eig1
cml
eig2
rept.out∗: Files on $HOME/hilapw/data are used unless provided in $CWD.
∗∗: Default values are assumed unless provided in $CWD.

4Installation

The first thing to install the package is to download the package hilapw.tar.gz on yourcomputer and to extract it as

# tar zxvf hilapw.tar.gz
Hereinafter, the symbol # at the first column denotes the prompt in the UNIX system andis not needed to be typed. The package is assumed to be installed on your home directory.
The following directory structure under the home directory is expanded with this.
hilapw/
1.4
bin@ -> current/binbz@ -> current/bz
configure.shcurrent@ -> 1.4
data@ -> current/data
Some of the tool commands and all of the executables after the proper compilation arelocated in the directory bin/. The directory bz/ contains the data for the typical Brillouinzone used for the band-structure plotting. In data/, the default data and some example data
are placed. To configure the information of the installed directory, type# cd hilapw# ./configure.sh
Before compilation, check the FORTRAN90 compiler name and its appropriate options available
on your computer and modify Makefile in hilapw/current, if necessary. (If you wish tocompile the parallel version of the package, compiler and option information in pMakefile arealso properly set.) In addition to the FORTRAN90 compiler, the libraries of linear algebra routines BLAS and LAPACK are required for installation of HiLAPW. After the proper specificationsin Makefile (and also pMakefile), one can compile by typing make to install.
# cd hilapw/current# make all
The same process is needed if you wish to compile the parallel version.
# make -f pMakefile all
For your calculations, let us create a new directory hilapw1 under your home directory.
# cd# mkdir hilapw1

5

Examples: Basic

5.1

Cu with FCC Structure

5.1.1

SCF Calculation

Under the directory hilapw1 just made above, another new directory Cu is created for calculations of fcc Cu.

# cd hilapw1# mkdir Cu# cd Cu
To perform example calculations for fcc Cu, extract files from the directory $HOME/hilapw/datalike# tar xvf $HOME/hilapw/data/Cu.tar
For starting the SCF calculation, just type the shell-script as# ./JOB-SCF
The students are recommended to look inside the script JOB-SCF and to understand how all
the processes are executed.
After the execution xlapw, the tool command LAcopy is used for attaching a particularmodifier (A1 or A2) to the output files.
5.1.2

Density of States

For the total density of states (DOS), type# xdoss# mv pdos pdosA2# getfermi outA2 > fermiA2# PSP < psp_tdos > tdos.ps
to get a PostScript (PS) file tdos.ps as shown in Fig. 1. Here, the tool command getfermi
is for extracting the Fermi energy from a file. PSP is a general-purpose PS plotting tool withthe specifications of a standard input file, psp_tdos in this case. For plotting partial DOS,
4.0

TOTAL DOS(/eV)

fcc Cu
3.0

2.0

1.0

0.0
-10 -8 -6 -4 -2 0 2 4
ENERGY(eV)

6

8 10

Fig. 1: Total density of states in fcc Cu. The origin of energy is taken at the Fermi energy.
since the required data are already generated and stored in the file pdosA2, type# PSP < psp_pdos > pdos.ps
to get projected DOS on cubic harmonics s, p, d(eg ) and d(t2g ) as shown in Fig. 2.

1.0

fcc CuDOS(/eV-atom-spin)

0.8
0.6
0.4
0.2
0.0
-10 -8 -6 -4 -2 0 2 4
ENERGY(eV)

6

8 10

Fig. 2: Partial density of states in fcc Cu. The origin of energy is taken at the Fermi energy.
Long broken, short broken, dotted and solid lines denote components of the partial DOS ofs, p, d(eg ) and d(t2g ), respectively.

5.1.3

Band Structure

Next, let us try to plot the band structure of fcc Cu.
# ./JOB-EK# ./JOB-SYM
The shell-script JOB-EK demands energy-eigenvalue calculations for k points along the highsymmetry lines of the fcc Brillouin zone in Fig. 3. The script JOB-SYM requests to extract
the irreducible representation for each eigenstate. To obtain the band structure (Fig. 4), justtype
# PSP < psp_ek > ek.ps
kz

ky
kx

Fig. 3: Brillouin zone of the fcc structure.

fcc Cu

10
8

ENERGY(eV)

6
4
2
0
-2
-4
-6
-8
-10
W Q L

Λ

Γ

∆

XSK

Σ

Γ

Fig. 4: Band structure of fcc Cu. The origin of energy is taken at the Fermi energy.

5.2

Si with Diamond Structure

5.2.1

SCF Total-Energy Calculations

Following the previous example for fcc Cu, get started with some preparations.
# cd $HOME/hilapw1# mkdir Si# cd Si# tar xvf $HOME/hilapw/data/Si.tar
With the present data set, an example for calculating the total energy by changing the lattice
constant for Si with diamond structure is provided.
# ./JOB-TEN
SCF calculations at eleven lattice constants between a = 5.20Å and a = 5.60Å are performed.
Using calculated results, total energies as a function of the unitcell volume are fitted to
Murnaghan’s equation of state[12].
# ./GET-TEN > TEN# xefitm < TEN > fit_TEN
The executable xefitm is a tool program to fit the data to Murnaghan’s equation of stateby the least-square method. Inside the output file fit_TEN, you can find the results as
Coefficients for Murnaghan fitting
c1=0.3072618807D+07
c2=
0.1652098587D-02
E’=-1156.7411297535
V0=266.2218508318
E0=-1156.1556676775
B’=
4.0200000000
B0=0.0066414363
97.6955282824 GPawhere V0 is the equilibrium unitcell volume in atomic units, B0 is the bulk modulus and B’is its volume derivative defined as(

d2 Etotal
B0 = B(V0 ) = V
dV 2(

B′ = −

1dB(V )
V
B0
dV

)

0

)
0

Figure 5 shows the volume dependence of the total energy obtained by typing# tail -103 fit_TEN > TEN2# PSP < psp_TEN > TEN.ps
where dots represent calculated data points and solid line is the fitted equation of state.
-1156.140

TOTAL ENERGY(Ry)-1156.142

Si-1156.144-1156.146-1156.148-1156.150-1156.152-1156.154-1156.156-1156.158-1156.160
220 230 240 250 260 270 280 290 300 310 320
VOLUME(a.u.)

Fig. 5: Volume dependence of calculated total energy in diamond Si. Dots represent calculated data points and solid line is the fitted equation of state.

5.3

Ferromagnetic Fe with BCC Structure

5.3.1

SCF Calculation

Again, let’s create a new directory and take all the data files for the present example.
# cd $HOME/hilapw1# mkdir Fe# cd Fe# tar xvf $HOME/hilapw/data/Fe.tar
In this section, an SCF calculation is carried out and DOS and band structure are plottedin a similar manner to the calculations for fcc Cu. To start the SCF calculation, type theshell-script# ./JOB-SCF
Detailed scripts in JOB-SCF should be checked carefully.

5.3.2

Density of States

全状態密度を計算するには In order to perform a DOS calculation and to get the total DOS
plot, type# xdoss# mv pdos pdosA2# getfermi outA2 > fermiA2# PSP < psp_tdos > tdos.ps

TOTAL DOS(/eV)

The output PS file tdos.ps is shown in Fig. 6. Partial DOS can be obtained by
4.0
bcc Fe
3.0
2.0
1.0
0.0
1.0
2.0
3.0
4.0
-10 -8 -6 -4 -2 0 2 4
ENERGY(eV)

6

8 10

Fig. 6: Total density of states in ferromagnetic bcc Fe. Upper and lower panels show densityof states for up and down spin bands, respectively. The origin of energy is taken at the Fermi
energy.
# PSP < psp_pdos > pdos.ps
Figure 2 shows projected partial DOS on cubic harmonics, s，p，d(eg )，d(t2g ).
5.3.3

Spin Magnetic Moment

Results of xlapw are output in a file lapw.out. (With use of the prepared shell-scripts, thefile is rename as outA1 or outA2.) One can find the following output data in the file as===== ELECPR
NUMBER OF ELECTRONS WITHIN MUFFIN-TIN
TYPE SPIN
TOTAL
S
P
1
UP
4.02803
0.091640.07447
1 DOWN
1.84418
0.093700.08735
U+D
5.872210.18534
0.16182
U-D
2.18385 -0.00206 -0.01287

D
3.86191
1.663135.52504
2.19878F

DOS(/eV-atom-spin)

1.0
0.8

bcc Fe

0.6
0.4
0.2
0.0
0.2
0.4
0.6
0.8
1.0
-10 -8 -6 -4 -2 0 2 4
ENERGY(eV)

6

8 10

Fig. 7: Partial density of states in ferromagnetic bcc Fe. Long broken, short broken, doted
and solid lines denote components of the density of states projected on cubic harmonics s，p，d(eg )，d(t2g ). Upper and lower panels show density of states for up and down spin bands,
respectively. The origin of energy is taken at the Fermi energy.
Here, the numerical data are the number of electrons, which is defined as an integrated densityof states at the Fermi energy for each partial (spin and orbital-angular-momentum projected)
state within the muﬃn-tin sphere. The spin magnetic moment inside the muﬃn-tin sphere(1Å) is 2.18µB , most of which comes from the d states. The s and p states are weakly andnegatively spin-polarized. Note that the magnitude of the spin magnetic moment here isan integrated value of the spin density, depending on the muﬃn-tin sphere radius assumed.
In the case of a ferromagnetic system, total spin magnetic moment per unitcell including
contributions outside the spheres can be evaluated by summing up the occupation for eachspin band and by subtracting them as----- BAND OCCUPATION
UP-SPIN
DOWN-SPIN
BAND
E-MIN
E-MAX
WEIGHT
E-MIN
E-MAX
WEIGHT1
0.065700.449101.00000
0.094160.55769
1.0000020.34063
0.52660
1.00000
0.463590.66445
1.00000
30.44597
0.68832
0.99998
0.556200.84957
0.8859040.517150.68832
0.999830.65867
0.84957
0.0078150.563640.708140.89049
0.71533
0.868920.00000
60.60697
1.38338
0.21600
0.782811.418460.00000
7
1.17953
2.48901
0.00000
1.26695
2.498210.00000
8
1.301822.48901
0.00000
1.381672.498210.00000
9
1.43833
2.50566
0.00000
1.49335
2.518260.00000
10
1.83505
3.374250.00000
1.92887
3.207810.00000------------------------------------------------------------------------------SUM
5.10629
2.89371
U-D
2.21257
to get the total spin moment of 2.21µB . Therefore, the spin magnetic moment in bcc Fe ismostly localized within the muﬃn-tin sphere.
5.3.4

Band StructureBand structure of ferromagnetic bcc Fe is calculated and plotted along the high-symmetrylines of the bcc Brillouin zone (Fig. 8).
# ./JOB-EK# ./JOB-SYM
The PS file containing the band structure can be generated by# PSP < psp_ek > ek.ps
and shown in Fig. 9.
kz

ky

kx

Fig. 8: Brillouin zone of bcc lattice.

6

Input Data

In this section, three input data sets.in, lapw.in and doss.in are briefly explained for the
fcc Cu example.

6.1

Input Data for xsets: sets.in

The input file sets.in includes data about the target material system such as crystal structure, atomic information and k-point grid. (In the following data, 01: and so represent theline number of the file for sake of explanation and is not included in the actual data.)
01:fcc Cu
02:-----nspin
03:1
04:-----space group
05:Fm-3m

10

bcc FeENERGY(eV)

5

0

-5

-10
Γ

∆ H G ND P Λ Γ Σ N

Fig. 9: Band structure of ferromagnetic bcc Fe. The majority-spin and minority-spin bands
are plotted with thick and thin lines, respectively. The origin of energy is taken at the Fermi
energy.
06: 3.61 3.61 3.61
07: 90.0 90.0 90.0
08:-----atoms
09:1
10:Cu 011:0.0 0.0 0.0
12:-----k points
13:0
14:8 8 8
The line 01 is a comment within 80 characters. This comment is also used as a header line
in some output files. The line 02 is a comment used for a data separator. The line 03
contains an integer number used as a flag NSPIN of spin-polarization treatment. NSPIN=1 fornon-spin-polarized and NSPIN=2 for spin-polarized calculation. In non-spin-polarized case,up and down spin states are treated equally and all the occupations are multiplied by thespin degeneracy of two when summing up the states. The line 04 is a comment used fora data separator. In line 05, name of the space group of the system is declared. Only the
space group names included in the file spgrdata alt are allowed to use. In line 06 the lattice
constants a, b, c are given in Å units and α, β, γ are in degrees. If the lattice constants
are incompatible with the space group given above, the lattice system of the space group
has the priority. For example, if a diﬀerent b value from a is given, b is replaced with the avalue. In Table 2, relations of the lattice constants to the lattice system and type are listed.
The line 08 is a comment used for a data separator. An integer number in line 09 indicates
crystallographically independent atomic kinds NTTP. Two atoms which have the same atomic
number should be counted two when no symmetry operation exists between both atomic

Table 2: Crystal system, lattice constants a, b, c, α, β, γ, lattice type. The b axis inmonoclinic lattice is taken as the unique axis here and the other choices may be allowed ina wider definition of space group.
crystal system (lattice)
cubic
tetragonal
orthorhombic
hexagonal
trigonal (hexagonal)
trigonal (rhombohedral)monoclinic
triclinic

a
a
a
a
a
a
a
a
a

b
a
a
b
a
a
a
b
b

c
a
c
c
c
c
a
c
c

α
90
90
90
90
90
α
90
α

β
90
90
90
90
90
α
β
β

γ
90
90
90
120
120
α
90
γ

lattice types
P, I, F
P, I
P, I, F, C
P
P
R
P, B
P

sites. Line 10 contains name of the element and number of atomic sites of one kind. In thefollowing lines starting 11, atomic position vectors of all the atomic sites of the kind are givenin the conventional lattice units. The line 12 is a comment used for a data separator. An
integer number in line 13 denotes how to set up the k-point grid. A Uniform grid including(excluding) Γ point is constructed for 0 (-1). If it’s a positive integer, it represents the numberof k-points to be calculated, and coordinates of k-point in (2π/a, 2π/b, 2π/c) units with thesampling weight should be given in the following lines.

6.2

Input Data for xlapw: lapw.in

A file lapw.in is the input data for xlapw and contains parameters controlling the calculations.
01:HiLAPW 1.0
02:
20.0
03:
160.0
04:
-1
05:
0.005
06:
0
07:
-1
08:
loop dta
09:
10 0.00emax
egmax
lbroaddeltaeneatomic_loop_mode
scf cmix smix field
5 0.20 0.00 0.00 0.00

The line 01 is a comment within 80 characters. Real numbers in lines 02 and 03 are thecutoﬀ energy in Ry used in the wavefunction and the electron density, respectively. An
integer number in line 04 specifies the method of k-point integration. The number is −1 forthe tetrahedron method with corrections, 0 for the standard linear tetrahedron method, 1
for the broadening method with the error function and 2 for the broadening method with
the Fermi function. In the case of the broadening method, a broadening width in Ry should
be given in line 05. In line 06, number of the valence states should be given, which mustbe equal to or greater than half of the number of the valence electrons in insulating systems
and must be suﬃciently larger than that in metallic cases. An appropriate number of thevalence states is proposed with 0 for the input. An integer number in line 07 controls atomic
force calculation. If it’s −1, no atomic force is calculated. Atomic forces are calculated
for 0 and greater. When it’s a positive integer, atomic positions are relaxed according to adiﬀerent algorithm with some extra data specified (not described here because of advanced
calculation). Line 08 is a comment. Input data in line 09 give control parameters for SCF
calculation. The example data denote 10 SCF iteration with the electron density mixing of
0.2 by the 5th order extended Anderson scheme. In the case of spin-polarized calculation,
the spin-density mixing ratio is also specified.

6.3

Input Data for xdoss: doss.in

The file doss.in contains the control parameters when calculation the density of states.
01:doss.in
02:-0.2 1.8 2000
03:1
04:0 0 0
The line 01 is a comment within 80 characters. In line 02, the energy range in Ry andenergy mesh number are given for DOS calculation. You may get some information about
the position of the bands from the output file of xlapw. The range of each band can begiven in BAND OCCUPATION part of the output file. An integer in line 03 specifies kind of theharmonics on which the partial DOS is projected. Spherical harmonics for 0, cubic harmonics
for 1 and hexagonal harmonics for 2. The principal axes of the harmonics can be changed
by specifying the Euler angles (α, β, γ) in degrees in line 04. After the execution of xdoss,
all the DOS data are output in the file pdos. An extra file pdos.index tells each column ofthe data corresponds to which partial DOS.

References[1] For example, see: R.M. Martin, ”Electronic Structure” (Cambridge University Press,
2004).
[2] J.C. Slater, Phys. Rev. 51, 392 (1937).
[3] O.K. Andersen, Phys. Rev. B 12, 3060 (1975).
[4] D.D. Koelling and G. Arbman, J. Phys. F: Metal Phys. 5, 2041 (1975).
[5] J.M. Soler and A.R. Williams, Phys. Rev. B 40, 1560 (1989).
[6] J.M. Soler and A.R. Williams, Phys. Rev. B 42, 9728 (1990).
[7] M. Weinert, G. Schneider, R. Podloucky, and J. Redinger, J. Phys.: Cond. Mater 21,
084201 (2009).

[8] D.D. Koelling and B.N. Harmon, J. Phys. C: Solid State Phys. 10, 3107 (1977).
[9] V. Eyert, J. Comp. Phys. 124, 271 (1996).
[10] M. Weinert, J. Math. Phys. 22, 2433 (1981).
[11] P.E. Blöchl, O. Jepsen and O.K. Andersen, Phys. Rev. B 49, 16223 (1994).
[12] F.D. Murnaghan, Proc. Natl. Acad. Sci. U.S.A. 30, 244 (1944).