HiLAPW— BASIC EDITION —Tamio OGUCHI
Osaka University, ISIR

1Introduction

In these days, band calculation is often called ﬁrst-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 ﬁrst principles. There
have been three important issues in progress of the ﬁrst-principles calculation. The ﬁrst 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 throughlocal-density approximation (LDA) and more recently generalized gradient approximation(GGA). The second is the highly-developed band-structure calculation methods whichmake possible a realistic calculation for a very complicated material system with reasonable computer resource. In the calculation, it is almost inevitable to assume numericalapproximations and modeling 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 and mixed up in the everyday calculations. In order to discuss physics andvalidity 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 powercomparable to a computer called ”supercomputer” ten years ago. (Concerning the ﬁrst
and second issues, beginner-course students are recommended to study ﬁrst fundamentalsof the electronic structure theory [1].)Under these circumstances, a computer code for the ﬁrst-principles calculation becomesa gigantic package including a lot of state-of-the-art techniques and related know-how.
Usefulness and wide applicability of the ﬁrst-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 ﬁeld. In this decade, a numberof research groups have developed ﬁrst-principles calculation codes and, in many cases,
made them open to public.
HiLAPW introduced in this article is a ﬁrst-principles calculation package originally developed by a theory group in Hiroshima University. The current developments of the package
are being made in collaboration between Osaka University and Hiroshima University. In
section 2, the basic speciﬁcations of HiLAPW are explained. The package contains several
executables, which perform part of the ﬁrst-principles calculation. The executables are1linked via data ﬁles, 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 Speciﬁcations

The HiLAPW program package is designed to perform ﬁrst-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 nonlinearized 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 Koellingand Arbman [4].

2.2Scalar 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) [7]. 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-ﬁeld (SCF) manner
and usually 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 extendedAnderson method [8] can be employed depending on the SCF convergency behavior indiﬀ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 [9] 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 ﬁnal
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

2.5

Atomic Forces

In pseudopotential methods with the plane wave basis, atomic forces, which are deﬁned
as a minus gradient of the total energy with respect to the atomic coordinates, are simply
given by 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
evaluated additionally to obtain the atomic forces accurately. In the package, Soler andWilliams’ formulation [5, 6] is used for the atomic forces. With the forces, structural
optimization and phonon 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 physicalquantities such as the electron density and the density of states. Blöchl’s tetrahedron
scheme [10] is utilized in the code. The other integration schemes with the occupation
given by 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,
the matrix element of the electric dipole with respect to the initial and ﬁnal states gives
the selection 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
attached with the package and FORTRAN90 compiler on user’s computer. List of the executables and related data ﬁles are summarized in Table 1.

4

Installation

The ﬁrst thing to install the package is to download the package hilapw.tar.gz to thehome directory on your computer and to extract it as# cd ~# gunzip hilapw.tar# tar xvf hilapw.tar
Hereinafter, the symbol # at the ﬁrst column denotes the prompt in the UNIX system andis not needed to be typed. The following directory structure under the home directory is
expanded with this.

3

Table 1: Executables and related data ﬁles in HiLAPW.
executable role
xsets
generation of the initial data

xlapw

xdoss

xnewa
xwbox
xpbox
xspin
xwcon

xsymm

xrept

SCF FLAPW calculation

input ﬁles
atomdata
atomdens
spgrdata
sets.inwavin
lapw.inoutput ﬁleswavoutsets.out

wavout
ekn
dis
foa
tau
ten
lapw.out
DOS calculation
wavin
pdos
doss.inpdos.index
doss.out
k-point update
wavin
wavoutnewa.out
electron density calculation on 3D mesh wavin
wbox.list
wbox.inwbox.out
potential calculation on 3D meshwavin
pbox.list
pbox.inpbox.out
addition of spin degree
wavin
wavoutspin.out
format conversion of wav ﬁlewavin
wavoutwavin.frm wavout.frm
wcon.inwcon.out
symmetry calculation
wavin
ekl
cml
symm.out
reorder of the eigenstates
ekl
eig1
cml
eig2
rept.out

4hilapw/
bin/
bz/
data/
lib/
ps/
sources/
Some of the tool commands and all of the executables after the proper installation 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 databases, atomdata,
atomdens and spgrdata, and some example archived data, Cu.tar, Si.tar and Fe.tar,
are placed.
Before compilation, check the FORTRAN90 compiler name and its appropriate options
available on your computer and modify makefile in three directories /hilapw/sources，/hilapw/lib and /hilapw/ps if necessary. In addition to the FORTRAN90 compiler, thelibraries of linear algebra routines BLAS and LAPACK are required for installation of HiLAPW.
After the proper speciﬁcations in makefile, one can compile by typing make to install.
# cd hilapw/sources# make clean# make all# make install
The same process is needed in the directories /hilapw/sources and /hilapw/lib.
# cd ~/hilapw/ps# make clean# make all# make install# cd ~/hilapw/lib# make clean# make all# make install
In order to make all the executables activated, the directory /hilapw/bin should be
added in your path. For example, when the csh or tcsh environment is used, the following
line after the last set path line in .cshrc can be inserted.
set path = ( . $home/hilapw/bin $path)
To activate it, just type.
# source ~/.cshrc
For example calculations, let us create a new directory hilapw1 under your home
directory.
# cd ~# mkdir hilapw1
5

5

Cu with fcc Structure

5.1

SCF Calculation

Under the directory hilapw1 just made above, another new directory Cu is created for
calculations of fcc Cu. By typing getdata, three database ﬁles atomdata，atomdens andspgrdata are imported from the default database directory ~/hilapw/data.
# cd hilapw1# mkdir Cu# cd Cu# getdata
To perform example calculations for fcc Cu, extract eleven ﬁles from the directory~/hilapw/data.
# tar xvf ~/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 particular
modiﬁer (A1 or A2) to the output ﬁles.

5.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) ﬁle tdos.ps as shown in Fig. 1. Here, the tool command getfermi
is for extracting the Fermi energy from a ﬁle. PSP is a general-purpose PS plotting toolwith the speciﬁcations of a standard input ﬁle, psp_tdos in this case.
For plotting partial DOS, since the required data are already generated and stored in
the ﬁle 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.

6

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

Figure 1: Total density of states in fcc Cu. The origin of energy is taken at the Fermi
energy.

5.3Band 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),
just type# PSP < psp_ek > ek.ps

6

Si with Diamond Structure

6.1

SCF Calculation

Following the previous example for fcc Cu, get started with some preparations.
# cd ~/hilapw1# mkdir Si# cd Si# getdata# tar xvf ~/hilapw/data/Si.tar
With the present data set, an example for calculating the total energy by changing thelattice constant for Si with diamond structure is provided.
# JOB-TEN

7

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

Figure 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 of s, p, d(eg ) and d(t2g ), respectively.
kz

ky
kx

Figure 3: Brillouin zone of the fcc structure.
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
ﬁtted to Murnaghan’s equation of state[11].
# GET-TEN > TEN# xefitm < TEN > fit_TEN
The executable xefitm is a tool program to ﬁt the data to Murnaghan’s equation of stateby the least-square method. Inside the output ﬁle fit_TEN, you can ﬁnd the results as
Coefficients for Murnaghan fitting
c1=0.3072618807D+07
c2=
0.1652098587D-02
E’=-1156.7411297535
V0=266.2218508318
E0=-1156.1556676775
8

fcc Cu

108
6ENERGY(eV)

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

Λ

Γ

∆

XSK

Σ

Γ

Figure 4: Band structure of fcc Cu. The origin of energy is taken at the Fermi energy.
B’
B0

=
=

4.0200000000
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 deﬁned as(

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

B0 = −

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 ﬁtted equation of state.

7

Ferromagnetic Fe with bcc Structure

7.1

SCF Calculation

Again, let’s create a new directory and take all the data ﬁles for the present example.
# cd ~/hilapw1# mkdir Fe# cd Fe# getdata# tar xvf ~/hilapw/data/Fe.tar
In this section, an SCF calculation is carried out and DOS and band structure areplotted in a similar manner to the calculations for fcc Cu. To start the SCF calculation,
type the shell-script9

-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.)

Figure 5: Volume dependence of calculated total energy in diamond Si. Dots represent
calculated data points and solid line is the ﬁtted equation of state.
# JOB-SCF
Detailed scripts in JOB-SCF should be checked carefully.

7.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
The output PS ﬁle tdos.ps is shown in Fig. 6.
Partial DOS can be obtained by# PSP < psp_pdos > pdos.ps
Figure 2 shows projected partial DOS on cubic harmonics, s，p，d(eg )，d(t2g ).

7.3

Spin Magnetic Moment

Results of xlapw are output in a ﬁle lapw.out. (With use of the prepared shell-scripts,
the ﬁle is rename as outA1 or outA2.)
One can ﬁnd the following output data in the ﬁle 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
10

D
3.86191
1.663135.52504
2.19878F

TOTAL DOS(/eV)

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

Figure 6: Total density of states in ferromagnetic bcc Fe. Upper and lower panels showdensity 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 deﬁned as an integrated
density of 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ﬃntin sphere (1Å) is 2.18µB , most of which comes from the d states. The s and p states are
weakly and negatively spin-polarized.
Note that the magnitude of the spin magnetic moment here is an integrated valueof the spin density, depending on the muﬃn-tin sphere radius assumed. In the case ofa ferromagnetic system, total spin magnetic moment per unitcell including contributions
outside the spheres can be evaluated by summing up the occupation for each spin bandand 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 ismostly localized within the muﬃn-tin sphere.
11

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

Figure 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.

7.4Band 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 ﬁle containing the band structure can be generated by# PSP < psp_ek > ek.ps
and shown in Fig. 9.
kz

ky

kx

Figure 8: Brillouin zone of bcc lattice.

12

bcc Fe

10ENERGY(eV)

5

0

-5

-10
Γ

∆ H G ND P Λ Γ Σ N

Figure 9: Band structure of ferromagnetic bcc Fe. The majority-spin and minority-spinbands are plotted with thick and thin lines, respectively. The origin of energy is taken at
the Fermi energy.

8

Input Data

In this section, three input data sets.in, lapw.in and doss.in are brieﬂy explained forthe fcc Cu example.

8.1

Input Data for xsets: sets.in

The input ﬁle 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
the line number of the ﬁle 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
06: 3.61 3.61 3.61
07: 90.0 90.0 90.0
08:-----atoms
09:1
10:Cu 1
11:0.0 0.0 0.0
12:-----k points
13:0
14:8 8 813The line 01 is a comment within 80 characters. This comment is also used as a headerline in some output ﬁles. The line 02 is a comment used for a data separator. The line 03
contains an integer number used as a ﬂag NSPIN of spin-polarization treatment. NSPIN=1
for non-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
the spin 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 ﬁle spgrdata 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 thea value. In Table 2, relations of the lattice constants to the lattice system and type arelisted. The line 08 is a comment used for a data separator. An integer number in line 09
Table 2: Lattice constants a, b, c, α, β, γ ，lattice system, lattice type. The c axis in monoclinic lattice is taken as the unique axis.
lattice systemcubic
tetragonal
orthorhombic
hexagonal
trigonal
monoclinic
triclinic

a
a
a
a
a
a
a
a

b
a
a
b
a
a
b
b

c
a
c
c
c
a
c
c

α
90
90
90
90
α
90
α

β
90
90
90
90
α
90
β

γ
90
90
90
120
α
γ
γ

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

indicates crystallographically independent atomic kinds NTTP. Two atoms which have thesame atomic number should be counted two when no symmetry operation exists between
both atomic sites. Line 10 contains name of the element and number of atomic sites of
one kind. In the following lines starting 11, atomic position vectors of all the atomic sitesof the kind are given in the conventional lattice units. The line 12 is a comment used fora 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 number of k-points to be calculated, and coordinates of k-pointin (2π/a, 2π/b, 2π/c) units with the sampling weight should be given in the following lines.

8.2

Input Data for xlapw: lapw.in

A ﬁle lapw.in is the input data for xlapw and contains parameters controlling the calculations.
01:HiLAPW 1.0
02:
20.0
03:
80.0
04:
-1
05:
0.005

emax
egmax
lbroaddeltae
14

06:
10
ne
07: c 0.0 c 0.0 c 0.0
ezr
08:
-1
atomic_loop_mode
09:
loop dta scf cmix smix field
10:
10 0.00
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 arethe cutoﬀ energy in Ry used in the wavefunction and the electron density, respectively. An
integer number in line 04 speciﬁes 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 methodwith 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, whichmust 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 numberof the valence states is proposed in the output ﬁle sets.out after the execution of xsets
as
NUMBER OF VALENCE STATESMINIMUM=
6
APPROPRIATE=
9
In line 07, the way to choose the energy parameter for each orbital-angular-momentum
wave of the atom kind is speciﬁed. Usually the example data set c 0.0 c 0.0 c 0.0 are
expected to work and all the energy parameters for s, p and d waves are taken at the centerof the occupied partial DOS. Sometimes, a higher energy parameter such as 1Ry above theFermi energy given as f 1.0 is needed for a particular partial wave in the case of a shallow
core state existing with the same orbital angular momentum. An integer number in line
08 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 arerelaxed according to a diﬀerent algorithm with some extra data speciﬁed (not describedhere because of advanced calculation). Line 09 is a comment. Input data in line 10 givecontrol parameters for SCF calculation. The example data denote 10 SCF iteration withthe electron density mixing of 0.2 by the 5th order extended Anderson scheme. In the caseof spin-polarized calculation, the spin-density mixing ratio is also speciﬁed.

8.3

Input Data for xdoss: doss.in

The ﬁle 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 ﬁle of xlapw. The range of each band can begiven in BAND OCCUPATION part of the output ﬁle. An integer in line 03 speciﬁes kind of15

the harmonics 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 bechanged by specifying the Euler angles (α, β, γ) in degrees in line 04. After the executionof xdoss, all the DOS data are output in the ﬁle pdos. An extra ﬁle pdos.index tells eachcolumn of the 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] D.D. Koelling and B.N. Harmon, J. Phys. C: Solid State Phys. 10, 3107 (1977).
[8] V. Eyert, J. Comp. Phys. 124, 271 (1996).
[9] M. Weinert, J. Math. Phys. 22, 2433 (1981).
[10] P.E. Blöchl, O. Jepsen and O.K. Andersen, Phys. Rev. B 49, 16223 (1994).
[11] F.D. Murnaghan, Proc. Natl. Acad. Sci. U.S.A. 30, 244 (1944).

16

HiLAPW— version 1.12 —Tamio OGUCHI
Osaka University, ISIRThe original edition of HiLAPW document describe practical procedures and examples
based on HiLAPW2002 (version 1.0) and additional issues introduced in version 1.12 aregiven in this article.

1

Additional Speciﬁcations and Tools in version 1.1

In HiLAPW version 1.1, several functionalities are included in addition to some patches toﬁx the bugs found in version 1.0.

1.1

Choice of Functional Form in the Local-Density or GeneralizedGradient Approximation

In version 1.0, it is not possible to choose the adopted functional form explicitly in thelocal-density approximation (LDA). On the occasion that the generalized-gradient approximation (GGA) is installed, the functional form used can be speciﬁed when xlapw is running.
Available functional form is written in the ﬁle excodata, in which the form in the ﬁrst line
is adopted. The default ﬁle excodata can be obtained by typing getdata in addition to
the previous database ﬁles.
1.1.1

Choice of Exchange-Correlation Functional: excodata

Before the execution of xsets, typing getdata gives three of the database ﬁles atomdata,
atomdens and spgrdata as described previously. In version 1.1, a new default ﬁle excodata
is also copied on your working directory. Exchange-correlation functional form used in
xlapw can be chosen by specifying the name of functional on the ﬁrst line of the ﬁleexcodata. The default excodata ﬁle is following:
LDA MJWGGA PBELDA BHLDA GLLDA VWNLDA PZ

1

In the present version 1.1, LDA functionals by Moruzzi-Janak-Williams (MJW) [1], BarthHedin (BH) [2], Gunnarsson-Lundqvist (GL) [3], Vosko-Wilk-Nusair (VWN) [4] and PerdewZunger [5] and GGA one by Perdew-Burke-Ernzerhof (PBE) [6] can be used.
The choice of the functional form is NOT a kind of adjustable parameters. Instead, itshould be used for checking the validity and reliability of the calculated results. Nevertheless, it is also of reality that either LDA or GGA works well for some particular systems
and physical phenomena and vice versa.

1.2

Modiﬁcations in Core Functions and Atomic-Force Calculation

Since GGA includes the second derivatives of the electron density, one has to pay more
attention to the representation of the electron density. Especially, as diﬀerent representations are assumed for the electron density inside and outside the muﬃn-tin spheres in
the full-potential method, a consistent numerical treatment must be paid to achieve thecontinuous representation of the electron density on the sphere boundary. Previously, thevalence electron density has been expressed in such a consistent manner but no such an
accurate treatment has been made on the core density, resulting to some serious problems
in the case of GGA.
Certain amplitude of a shallow core state may penetrate into the outer region of the
muﬃn-tin sphere. The amplitude outside the sphere was treated in the plane-wave representation and a constant density was assumed inside the sphere when the Fourier transformation is performed, leading to some spurious noise due to its discontinuity. A gaussian
type function is now used inside the sphere so as to fulﬁll the continuity on the sphereboundary. However, the previous treatment on the core density gives only very minor error in the atomic forces but the corrected representation requires an explicit considerationof the core density in the plane-wave part of the atomic forces. With the present corrections, the core density can be treated in GGA and atomic-force calculations in a consistentand accurate way.

2

Additional Speciﬁcations and Tools in version 1.11

In HiLAPW version 1.11, new functionality for input method of atomic positions is includedin addition to some patches to ﬁx the bugs found in version 1.1.

2.1

Input Method of Atomic Positions

When the atomic positions are given in the previous version, the position vectors of all
atoms in unitcell should be input and their consistency to the space group is checked.
It is highly desired to simplify the input method in cases of many-atom systems and/or
displacing the atoms for a given space group. Therefore, it is now possible in this version
to input only one position vector for one kind of atom.
For example, the input data sets.in for Y2 C3 in the previous version are following.
Y2C3 #220 I-43d a=8.18976 x(Y)=0.05017 x(C)=0.29481-----nspin
12

-----space group
I-43d
8.18976 8.18976 8.18976
90.0 90.0 90.0-----atoms
2Y
8
0.05017 0.05017 0.05017 x,x,x
0.30017 0.30017 0.30017 1/4+x,1/4+x,1/4+x
0.55017 0.44983 -0.05017 1/2+x,1/2-x,-x
0.80017 0.19983 0.69983 3/4+x,1/4-x,3/4-x-0.05017 0.55017 0.44983 -x,1/2+x,1/2-x
0.69983 0.80017 0.19983 3/4-x,3/4+x,1/4-x
0.44983 -0.05017 0.55017 1/2-x,-x,1/2+x
0.19983 0.69983 0.80017 1/4-x,3/4-x,3/4+x
C 12
0.29481 0.00000 0.25000 x,0,1/4-0.29481 0.50000 0.25000 -x,1/2,1/4
0.54481 0.50000 0.25000 1/4+x,1/2,1/4
0.45519 0.00000 0.25000 3/4-x,0,1/4
0.25000 0.29481 0.00000 1/4,x,00.25000 -0.29481 0.50000 1/4,-x,1/2
0.25000 0.54481 0.50000 1/4,1/4+x,1/2
0.25000 0.45519 0.00000 1/4,3/4-x,00.00000 0.25000 0.29481 0,1/4,x0.50000 0.25000 -0.29481 1/2,1/4,-x0.50000 0.25000 0.54481 1/2,1/4,1/4+x
0.00000 0.25000 0.45519 0,1/4,3/4-x-----k-points
04
4
4
The present version 1.11 may accept the previous data and also the following simpliﬁed
data format as
Y2C3 #220 I-43d a=8.18976 x(Y)=0.05017 x(C)=0.29481-----nspin
1-----space group
I-43d
8.18976 8.18976 8.18976
90.0 90.0 90.0-----atoms
2Y
00.05017 0.05017 0.05017 x,x,xC
00.29481 0.00000 0.25000 x,0,1/4
3

-----k-points
04
44
where the value 0 following the element name denotes only one position vector for one kindof atom.
Note that since all the atomic position vectors are generated according to the input
atomic position and the space group, the generated position vectors should be checked.

3

Additional Speciﬁcations and Tools in version 1.12

In HiLAPW version 1.12, a modiﬁcation in linealization for preventing the ghost band ismade in addition to some patches to ﬁx the bugs found in version 1.11.

3.1

Modiﬁcation in Linear Method

In linear method [7, 8], the energy dependence of the radial function Rl (r, E) used as abasis set can be represented eﬃciently. This is well understood from its Taylor expansion
Rl (r, E) = Rl (r, El ) + (E − El )Ṙl (r, El ) + ·(1)

with use of the radial and energy-derivative functions Ṙl (r, El ) = (dRl (r, E)/dE)El solved
at a particular energy El . Alternatively, the energy-derivative function is orthogonalizedwith the corresponding radial function normalized inside the MT sphere∫ S
0

r2 drRl (r, El )Ṙl (r, El ) = 0(2)

and the radial and energy-derivative functions provide a minimal basis for a given orbital
angular momentum component l.
In linear augmented plane wave method, the plane waves are augmented with the radialand energy-derivative functions and the radial function with any logarithmic derivative D
can be constructed with their linear combination such as
Rl (r, D) = Rl (r, El ) + ω(D)Ṙl (r, El )(3)This is quite convenient for the connectivity to the plane wave on the sphere. Especially, theexpectation value inside the sphere with respect to the radial function with the logarithmic
derivative D is shown to be correct up to the third order.
hRl (D)|H − E|Rl (D)iS = O(E − El )4(4)

However, in linear method, the so-called ghost band may happen to appear and become
very serious problem, particularly in the case of systems with a shallow core state. To
eliminate the ghost band, the energy parameter El must be set to a higher energy value
and the energy representation of the radial function may be worse in the valence band
region.
One of the reasons of the ghost band originates in the fact that the energy-derivative
function in the valence band region may be a possible bound state similar to the shallow
core state. Figure 1 shows the logarithmic derivatives of the s and p channels in the4

LOGARITHMIC DERIVATIVES(Ry)

20
Ti-s

Ti-p

10
0
-10
-20
-5 -4 -3 -2 -1 0 1 2 3 4 5ENERGY(Ry)-5 -4 -3 -2 -1 0 1 2 3 4 5ENERGY(Ry)

Figure 1: The logarithmic derivative of the s and p channels in Ti. Broken and dottedlines denote that for the radial and energy-derivative functions. Solid lines indicate thatfor the orthogonalized energy-derivative function.
valence band region in Ti, which has shallow core 3s and 3p states. A possible boundstate may appear in the energy region where the logarithmic derivative is negative. Theenergy-derivative functions with s and p waves become negative in the valence band region.
Let us consider to correct the energy-derivative function to prevent the appearanceof the ghost band. Well localized core states with no amplitude on the MT sphere are
rigorously orthogonalized to the radial and energy-derivative functions. However, the orthogonality is broken when a core state has ﬁnite amplitude penetrating outside the MT
sphere, which actually happens for shallow core systems. If this is the case, the energyderivative function has a negative logarithmic derivative in the valence band region andmight give rise to a bound state (ghost band) similar to the shallow core state. So, theorthogonality to the core state should be recovered to eliminate the ghost band. Suppose
˜ orthogonalized not only to the shallow core
the following new energy-derivative function Ṙlstate Rcl but also to the corresponding radial function.
˜ = Ṙ + αR + βR
Ṙ
l
l
cl
l

(5)

where coeﬃcients α and β are determined by the orthogonal conditionsα=−hRcl |Ṙl ihRcl |Rcl i − |hRcl |Rl i|2(6)

β = −αhRcl |Rl i(7)

˜ has a non-negative value ofAs expected, this orthogonalized energy-derivative function Ṙl
the logarithmic derivative in the valence band region and shows very smooth behavior as
a function of energy (see solid lines in Fig. 1). In an actual band structure calculation
with the new energy-derivative function, the energy parameters are always chosen within
the valence band region without any appearance of the ghost band.

5

4Installation

Downloaded the complete set of HiLAPW version xx.xxx, hilapw_xx.xxx_tar.gz, is placedon the home directory, unzipped and extracted by typing# cd ~# gunzip hilapw_xx.xxx_tar# tar xvf hilapw_xx.xxx_tar
to get the following directory structure under the home directory:
hilapw/xx.xxx/
bin/
bz/
data/
lib/
ps/
sources/sources_fermi/bin/bz/
current/
data/
The directory current is symbolically linked to that with the latest version name (xx.xxx
in the present case). The directories bin, bz and data are also symbolically linked to thecorresponding directories under current. Compilation by typing make is just in a similarmanner to the previous version.
# cd hilapw/current/sources/# make clean# make all# make install
The codes for band ﬁtting described in the advanced edition of HiLAPW document arelocated in the directory hilapw/current/sources_fermi/. After completion of the compilation above, the band-ﬁtting codes can be compiled as follows:# cd hilapw/current/sources_fermi/# make clean# make all# make install
The setting for path is just the same as before. For example, if you use csh new settingin .cshrc is given as
set path = ( . $home/hilapw/bin $path)6References[1] V.L. Moruzzi, J.F. Janak and A.R. Williams, Phys. Rev.B 12, 1257 (1975).
[2] U. von Barth and L. Hedin, J.Phys.C:Solid State Phys. 5, 1629 (1972).
[3] O. Gunnarsson and B.I. Lundqvist, Phys. Rev. B 13, 4274 (1976).
[4] S.H. Vosko, L. Wilk and M. Nusair, Can. J. Phys. 58, 1200 (1980); S.H. Vosko and
L. Wilk, Phys. Rev. B 22, 3812 (1980); G.S. Painter, Phys. Rev. B 24, 4264 (1981).
[5] J.P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).
[6] J.P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 77, 3865 (1996); 78, 1396
(1997).
[7] O.K. Andersen, Phys. Rev. B 12, 3060 (1975).
[8] D.D. Koelling and G. Arbman, J. Phys. F: Metal Phys. 5, 2041 (1975).

7