ESopt Manual
ver 1.2.1.E

August 2011 (English version)

Division of Frontier Materials Science,
Department of Materials Engineering Science,
Graduate School of Engineering Science,
Osaka University

1ESopt is a program package for the electronic structure calculation based on thedensity functional theory (DFT). The program adopting the plane-wave expansion
method is originated from the package developed in the Institute for Solid State
Physics (ISSP), University of Tokyo, which is called “opt”. This open software isprovided as “ISSP FPSOPT” by ISSP. Our ESopt, which is a revised version ofopt, was produced by choosing functions selectively and by applying the Fortran 95
scheme. We also added another diagonalizing routine to the original version.
The original “opt” has two characteristic features.
• High readability of the source code
• Realization of a CG (conjugate-gradient method) to optimize all degrees of
freedom of the wavefunction.
The former advantage is succeeded in ESopt. To solve the memory allocation problem by Fortan 77, some routines became lengthy and hard for easy reading in theoriginal program. This point is avoided using Fortran 90. The latter specialty wasdeveloped by Dr. Tadashi Ogitsu, who is now in the Lawrence Livermore National
Laboratory (LLNL), US. Compared to the ordinal CG method using the band-byband minimization in the band structure calculation, the method by Ogitsu perform the simultaneous optimization of all degrees of freedom for the Kohn-Sham
wavefunction. Thus, “opt” is a special program, whose implementation realizes theoriginal idea of the CG method as it is.
This manual is provided in order to achieve the next purpose. We give the basic
equations and the real implementation of the steps in this manual. By comparing
the mathematical representation and the program code, we can give beginners achance to perform the computation, understanding the calculation steps in theirmind. For this purpose, we do not intend to describe the theory completely, but
we select the minimum contents. However, we introduce some speciﬁc examples
showing real implementation adopted in “opt” aiming at deep understanding of thereaders.
ver. 1.0.0 The ﬁrst Fortran 90 version adopting only Sopt functions optimizing
atomic positions.
ver. 1.1.0 The revised version realizing Allopt functions performing stress calculations and cell optimization.
ver. 1.2.0 The version combining the band-structure calculation.

1

Expression of the total energy

In this section, we give the expression of the total energy in the plane-wave expansion method. We describe the formula when the norm-conserving pseudo-potential2

is adopted. This allows readers to understand practical calculation methods in the
DFT simulation. We deﬁne RI , ZI , V̂Ipseudo (r − RI ) as the atomic coordinate, theionic charge, and the pseudo-potential of the I-th atom.
We utilize the next expression for the total energy of a crystal with the volumeΩ.
Etotal = Ekin + Eel−el + Exc + Eel−ion + Eion−ion .
(1)Here, Ekin is the kinetic energy of electrons, Eel−el represents a part of the electronelectron interaction called the Hartree energy, which is the classical Coulomb energy for the electron charge density. We also have the exchange-correlation energy
Exc , the static electron-ion Coulomb interaction Eel−ion , and the ion-ion repulsive
Coulomb energy Eion−ion . We assume that the ion core gives the spherical potential.
Each contribution has the next expression in the real space.
Ekin =
Eel−el∑ ∫k,n,σ Ω

dr ϕ∗k,n,σ (r)()

1− ∇2 ϕk,n,σ (r) ,2

1 ∫∫n(r)n(r′ )=dr dr′,2 Ω|r − r′ |

(2)
(3)

∫

Exc =
Eel−ion =Ω

dr εxc (n(r))n(r) ,

∑ ∑∫

Ω

k,n,σ I

Eion−ion =(4)

dr ϕ∗k,n,σ (r)V̂Ipseudo (r − RI )ϕk,n,σ (r) ,1 ∑ ZI ZJ.
2 I,J |RI − RJ |

(5)
(6)

Here, n(r) is the single-particle density given as,n(r) =∑nσ (r) ,

(7)

σ=↑,↓

nσ (r) =∑|ϕk,n,σ (r)|2 ,(8)

εk,n,σ ≤EF

The Kohn-Sham orbital represented by a wave function, ϕk,n,σ (r), is given for then-th band with the spin σ at the sampled k point, k. The energy of the orbital isgiven by εk,n,σ . The Fermi energy EF is for the imaginative independent Fermion
system.

2

Dual space formalism

In this package, each wavefunction is expanded in the plane wave basis by adopting
the periodic boundary condition. The Bloch theorem and the Fourier transformation
allows us to represent the wavefunction in a series expansion of plane waves as,
ϕk,n,σ (r) =∑ϕk,n,σ (G) exp (i(k + G) · r) .

G

3

(9)

Here, the summation with respect to G is taken on the reciprocal lattice.
In a real calculation of atoms and molecules, we often consider a periodic structure of atoms and/or molecules, allowing the periodic boundary condition. By thisway, we can accelerate the simulation using the speed of the plane wave expansion
method. We will explein some extra limitations of the method in the latter part.
We call the unit of the periodic structure “the unit cell”. Introducing a large ﬁnite
integer N , we assume that the unit cells of N 3 are periodically aligned in a lattice.
The repeating unit cells are obaying the total periodic boundary condition. Thus
we have double periodicity. We introduce a volume Ωcell as the volume of the unit
cell. The total system has the volume Ω = N 3 Ωcell .
We consider the Fourier transformation. We will introduce a cutoﬀ so that a ﬁnite
descrete sum of G will appear. To derive formulae of the tansilation, however, it isbetter to derive the equations for continuum system with coordinates r as continuum
variables. Therefore, in eq. (??), the summation of G is an inﬁnite summation. The
electron charge density should have the same periodicity as the crystal. Thus wehave the Fourier transform of the density.
1 ∫nσ (G) =drnσ (r) exp (−iG · r) ,(10)
Ωcell Ωcellnσ (r) =∑nσ (G) exp (iG · r) ,

(11)

G

In the calculation step, both of the real space and the reciprocal space (or the
momentum space) are utilized. The integration for the exchange-correlation energy
is usually performed in the real space. the kinetic energy is much easily evaluatedin the momentum space. To perform both of numerical integration, we adopt two
representations in the real and reciprocal spaces. The method to perform these
integration by transforming information from each other is called “the dual spaceformalism”. This idea with the fast Fourier transformation (FFT) allows us toreduce the calculation steps.
The solution of the Kohn-Sham equation is represented as the wavefunction
ϕk,n,σ (G) in the reciprocal space. The single-particle density given by the wavefunctions determines each energy term in the density functional. To realize thecalculation, we should have representations of the density both in the real and reciplocal spaces. The procedure is given as follows.
To have a numerics, G is restricted in a ﬁnite region in the whole reciplocal space.
More precisely, the cutoﬀ energy Ecut is introduced for the plane-wave expansion.
The Fourier coeﬃcient of the wavefunction ϕk,n,σ (G) is assumed to be zero, when
|G|2 > Ecut . The energy cutoﬀ determines a subspace of the reciprocal space. We
deﬁne a ﬁnite reciprocal mesh, in which |G|2 ≤ Ecut holds. （Fig. ??）Owing to thediscrete Fourier transformation, a corresponding lattice structure is given in the real
space. This latter mesh is called the real mesh. Points on the real mesh is denotedas rj . Now we have the next procedure.
4

((6

'EWV

+PXGTUG((6/QOGPVWO5RCEG

4GCN5RCEG

図 1: A schematic viewgraph explaining the real mesh and the reciprocal meshtransformed by the Fourier transformation. In the reciprocal space (the momentum
space), the wavefunction has a ﬁnite coeﬃcient for the plane wave with G satisfying√|G| ≤ Ecut .

1. The wavefunction ϕk,n,σ (G) in the reciprocal space is inverse-Fourier-transformed
to have another expression ϕk,n,σ (rj ) given by Eq. (??) in the real mesh.
2. Using ϕk,n,σ (rj ), the electron density nσ (ri ) is calculated using Eq. (??).
3. The Fourier representation of the density nσ (ri ) is given by FFT as a series
expansion.
In the last step, we utilize,nσ (G) =∑nσ (rj ) exp (−iG · rj ) .

(12)

j

3

Expression of the energy

Following the procedure given in the last section, we have {ϕk,n,σ (G)}, {ϕk,n,σ (rj )},{nσ (rj )}, {nσ (G)}. The actual calculation is done using the formulae given in thissection. Each contribution for the energy is given by an expression per a unit cell.

3.1

Kinetic energy

This term is obtained in the reciprocal space. The kinetic energy Ēkin per a unit
cell is,Ēkin =∑ ∑ ∫k,n,σ G,G′

dr
Ωcell

|k + G|2
exp (i(G − G′ ) · r) ϕ∗k,n,σ (G′ )ϕk,n,σ (G)25=Ωcell ∑ ∑|k + G|2 |ϕk,n,σ (G)|2 .
2 k,n,σ G(13)This expression may be used as it is for the discrete representation by rewriting
the sumation with a ﬁnite sumation with respect to G. The calculation of the
kinetic energy is possible by evaluating an approximated expression using the ﬁnite
diﬀerence method, since the total calculation cost will not increase so much. The
expression and the treatment is much easy in the momentum space .

3.2

Hartree energy

This term is also evaluated in the reciprocal space. We deﬁne Ēel−el as the Hartree
energy per a unit cell.
1
Eel−elN3∫∫1n(r)n(r′ )=drdr′
2 Ωcell|r − r′ |
Ω
∫
∫
1
1 ∑
=
dr
ds
nσ (G)nσ (G′ ) exp (iG · r + iG′ · (r − s))2 Ωcell|s| G,G′ ,σ,σ′ΩĒel−el ==∫Ωcell ∑1nσ (G)nσ (−G) ds exp (iG · s) .
2 G,σ,σ′|s|ΩHere, the divergence appearing when G → 0 also appears in Eion−ion , Eel−ion andtotally cancelation of these three diverging terms happens. Thus, the diverging
term is omitted in this stage. For Ēel−el , the term for G = 0 corresponds to thisdivergence. Thus, we remove this term at G = 0. Furthermore, when N → ∞, ifG ̸= 0∫1
4π
ds exp (iG · s) = 2 .
|s|GΩWe can derive nσ (−G) = n∗σ (G), since the density n(r) is real. If we write n(G) =∑σ nσ (G), then we have the next expression.
Ēel−el =Ωcell ∑′ 4π|n(G)|2
.
2 G
G2

∑

(14)

Here, ′G speciﬁes a summation except for G = 0.
When the evaluation is done in the real space, the calculation requires a six
2) numerical steps with Nmesh as the number ofdimensional integral and O(Nmeshreal-mesh points. When it is in the reciprocal space, the deterministic process is the
FFT steps of O(Nmesh log Nmesh ) numerical steps. Thus, the calculation is easily
done in the reciprocal space.

6

3.3

Exchange-correlation energy

There are several known expressions for the exchange-correlation energy density
εxc . In this package, we adopted the Perdew-Wang expression (PW92). Comparedto the Perdew-Zunger expression (PZ81), PW92 applies the spin-polarization dependence by Vosko et al. and it is given by weighted parametrization considering
the error in each calculation point of the Monte-Carlo calculation Basically, theseparameters are determined refering the Diﬀusion Monte-Carlo calculation of theuniform electron gas system.
∫Ēxc =Ωcell

drεxc (n(r)) n(r) .

(15)

The exchange-correlation energy density εxc may be separated into the exchange
contribution εx and hte correlation energy density εc .
εxc = εx + εc .

(16)

In ESopt, only the spin-paramagnetic states are treated. The exchange energy
density is given as
(

)1

3
9π 3
εx (rs ) = −4πrs 4()1
3 3
rs = −,
4πn(17)while expression for the correlation energy density is omitted.

3.4Electron-ion interaction

The charge of ions are treated as the classical charge distribution. For each chargeof an ion is described by the pseudo-potential. In this package, Eel−ion is calculatedby the potential V̂Ipseudo . The eﬀect of core electrons is included in the pseudopotential, so that the state of the valence electrons is explicitly determined in thepseudo-potential method. In other word, the charge distribution of core electrons is
ﬁxed. In this package, the energy of the pseudo-potential is given as a summationof the local potential and the non-local potential.
non−local
local
Eel−ion = Eel−ion+ Eel−ion

=

∑∫
I

+

Ωc elldr n(r)VIlocal (|r − RI |)∑ ∑ ∫k,n,σ I,l̸=l0Ωc elldr ϕ∗k,n,σ (r)δ V̂Il (r − RI )ϕk,n,σ (r) .

(18)Here, l denotes the angular momentum. The potential adopted in the package,
one of components with l = 0, 1 or 2 is included in the local potential, where the7speciﬁed local component is denoted as l0 . The other components are treated by
the non-local potential.
VIlocal (r) = VIcore + δVIl0 .
(19)
Here, VIcore is the potential by the core electrons.
In the reciprocal space, Eel−ion is given as,
Eel−ion = Ωcell

′ ∑
∑
G

+ Ωcell

SI (G)n(G)VIlocal (G) + Ωcell n(0)

I

∑ ∑ ∑
k,n,σ G,G′ I,l̸=l0

∑

αI

I
∗
l
SI (G − G )ϕk,n,σ (G)δ V̂I (k + G, k + G′ )ϕk,n,σ (G′ ) .
′(20)
SI (G) = exp (−iG · RI ) .

(21)

Here, VIlocal (G) is the Fourier component of the local potential of the I-th atom,
and SI (G) is the structure factor.
We noted that divergence in Eel−el , Eion−ion , and Eel−ion are perfectly cancelled,when they are summed. Introduction of the pseudo-potential method, however,causes a residual value and the expression of the total energy has a shift. This shift
given as the limit of G → 0 for the expression of the energy is obtained as the
local
diﬀerence between the local component Eel−ion
in the limit of G → 0 and the limit
2of the Coulomb potential −4πZI /Ωcell G for the I-th atom. The value is given as
αI ,
)
(
4πZI
.
(22)
αI = lim SI (G)VIlocal (G) −G→0Ωcell G2
The energy of the non-local component is calculated by O(N 2 ) operation steps for
N plane waves. The value is reduced to O(N ) by adopting an approximation given
by Kleinman and Bylander, which is called the K-B approximation. We adopt theformula by the K-B approximation given as,δ V̂ non−local =∑ |δ V̂ l ϕl,m ⟩⟨ϕl,m δ V̂ l |⟨ϕl,m |δ V̂ l |ϕl,m ⟩

l,m

.

(23)

Here ϕl,m denotes a pseudo wave function with the angular momentum l with amagnetic quantum number m. The eq. (??) is exact, if the wavefunctions ϕl,m form
the complete set. In the K-B approximation, since the functions are given by thepseudo wavefunction, the expression becomes an approximated formula in general.
In a large number of systems, this approximation is known not to cause any badeﬀect.
However, in some systems, a side eﬀect by the K-B approximation is known toappear. in a general expression, the non-local potential is given as a local expression
in the radial direction, which is called the semi-local potential. When the K-B
approximation is adopted, non-locality arises also for the radial direction. Therefore,8the Wronsky theorem becomes inapplicable. Thus, for the K-B approximation, thepotential may have a false bound state, which may be energetically lower than atrue bound state. This abnormal bound state is called the ghost state.
In the K-B approximation, the non-local potential δ V̂ l (k + G, k + G′ ) is given by
the next expression.
∫1
dr exp (−i(k + G) · r) δ V̂ l (r)P̂l exp (i(k + G) · r)Ω Ω
4π=
(2l + 1)Pl (cos θk+G,k+G′ )Ω∫δ V̂ l (k + G, k + G′ ) =×cos θk+G,k+G′
χlI (k + G)∞0

dr r2 jl (|k + G|r)jl (|k + G′ |r)δ V̂ l (r)4πχlI (k + G)χlI (k + G′ )′=(2l + 1)Pl (cos θk+G,k+G ), (24)ΩXIl(k + G) · (k + G′ )=,(25)|k + G||k + G′ |
=

XIl =

∫ ∞
0

∫ ∞
0

dr r2 jl (|k + G|r)δ V̂ l (r)Rl (r) ,(26)

dr r2 δ V̂ l (r)|Rl (r)|2 .

(27)

Here, Rl (r) is the radial part of the pseudo wavefunction ϕl,m = Yl,m (θ, ϕ)Rl (r) forthe pseudo-potential δ V̂ l (r).

3.5

Ion-ion interaction

This contribution is evaluated by the Ewald summation.
Eion−ion =1 1 ∑ ZI ZJ
2 Ncell I,J |RI − RJ |

=

Zn Zm
1 ∑
2 n,m,L |dn − dm − L|(erfc(κ|dn − dm − L|) erf(κ|dn − dm − L|)
1 ∑
Zn Zm
=
+
2 n,m,L|dn − dm − L|
|dn − dm − L|

)

.(28)

Here, L is the translation vector for the periodic system. The coordinate RI for
an atom is rewritten as RI = dn + L with the position of the atom in a unit cell
dn . The ﬁrst term in eq. (??) is obtained in the real space adn the second term is
calculated in the reciprocal space.
1 ∫dr f (r) exp (−iG · r − µr)f (G) = lim √
µ→0
Ωcell
2
4π e−(G/2κ)=.
Ωcell G2
9(29)

Here, the Fourier coeﬃcient f (G) behaves in the limit of G → 0 as,
[

(

]

)

4π
G 2
1
f (G) ≃1−
+ O(G4 )
Ωcell
2κ
G2
[
]
4π
1
1
=
−
+ O(G2 )Ωcell G2 2κ

(30)

The ﬁrst diverging term is cancelled with Eel−el and Eel−ion so that the second term
remains. Therefore, the expression of Eion−ion without the divergence is given using
the Ewald parameter γewald as,
γewald =

∑

[

Zn Zm

n,m

+

(

2
4π ∑′ e−(G/2κ) iG·(dn −dm )1e−2Ωcell GG
4κ∑′ erf(κ|dn − dm − L|)L|dn − dm − L|

)

]

2κ
− √ δn,m .
π(31)

When the number of atoms in a unit cell becomes large, the calculation time
for the Ewald summation becomes non-negligible. To save the calculation cost, we
may optimize the numerical summation required to achieve convergence. In eq.
(??), we have cutoﬀs for summations in the real and reciprocal space. They are
estimated as 1/κ and 2κ. When a side of a unit cell is L, the length of a reciprocallattice vector becomes G0 = 2π/L. In both spaces, the maximum values of lattice
vecrors required for the convergence in γewald are given as cutoﬀ radii. We haverelations for these values as Lmax = LNRmax ∝ 1/κ and Gmax ∝ G0 NGmax ∝ 2κ.
Therefore, since the total number of( the) lattice points in each space is estimated as
3
κL
1
3
3
Ntotal ∝ NRmax
we can minimize the calculation steps by+ NGmax∝ (κL)
3 +ππtaking κ ≃ L . This value becomes κ ≃ (Lx LyπLz )1/3 , when the cell is a rectangularparallelepiped.

4

Hellmann-Feynman ForceIn this section, we give the expression of the inter-atomic force. The force FI
acting on the I-th atom is given by taking the derivative of the total energy withrespect to RI .
FI = −
= −

dE
dRI
∑
i

[

(

)

(

)

d
d
⟨ϕi |
H |ϕi ⟩ +⟨ϕi | H|ϕi ⟩ + ⟨ϕi |H
dRI
dRI

(

)]

d
|ϕi ⟩]
dRI(32)

Here, i(= (k, n, σ)) is a set of quantum numbers specifying a Kohn-Sham wavefunction. If the Kohn-Sham Hamiltonian is written as H, we have the total energy∑
in a simple expression of E = i ⟨ϕi |H|ϕi ⟩. In general, we have contribution from10the wave function for the derivative, the calculation cost becomes very big. When
wavefunctions ϕi give the orthonormalized complete set of the eigen states of H,
however, the second term of eq. (??) becomes zero.
∑

[(

)

d
⟨ϕi | H|ϕi ⟩ + ⟨ϕi |HdRI
i∑ d⟨ϕi |ϕi ⟩
=
εi
=0
dRI
i

(

d
|ϕi ⟩
dRI

)]

(33)

Thus, the expression for the inter-atomic force becomes,
FI = −

∑

⟨ϕi |

i

dH
∂E
|ϕi ⟩ = −
dRI
∂RI

(34)

The equation (??) has a meaning that the physical force on the atom is given by thepartial derivative of the energy with respect to the atomic coordinate, when eachwavefunction is the eigen state of the Hamiltonian. This conclusion is sometimes
called as the Hellmann-Feynman theorem. Using the force theorem, we can obtain
the physical force acting on the atom rather easily. However, even wehn the theorem
is adopted, we need to have enough accurate convergence of the total energy to thetrue ground state energy. In addition, obtaining the force with enough accuracy ismuch harder than calculating the total energy with the same accuracy. The reason is
because the error in the force becomes the ﬁrst order in the error in the wavefunction,
while the error in the total energy is the second order.
When the wavefunction is expanded in the plane waves, we have a special featurethat the atomic coordinates do not explicitely appear in the expression. In this condition, when the derivative of the energy in our method using the norm-conserving
pseudo-potential is obtained, only two terms from Eel−ion and Eion−ion have ﬁnite
values in the derivative.
FI = −∂(Eel−ion + Eion−ion ) .
∂RI

(35)

∂Eel−ionis given by two contributions from the local component and the∂RI
non-local component.
Here,

5

Quantum Stress

The expression of the internal stress is given by the variational method similarly
to the force. Let’s consider an electronic state with a unit cell with the volumeΩcell . The total energy is given as E({ψk,i (r)}, {RI }). The variational principle
tells that the ground state energy is given as the minimum of E({ψk,i (r)}, {RI })by optimizing the variables, {ψk,i (r)}, {RI }. To obtain the internal stress, we need11to consider the symmetrized strain tensor ε by eliminating the lotational degrees
of freedom as the second order tensor. By the scaling relation with respect to theinﬁnitesimal change in ε, we have the next expression.
RI → R′I = (1 + ε)RI ,ψk,i (r) → ψk,i (r′ ) = det(1 + ε) 2 ψk,i (r).
1

The ground state energy shift as,
Etot → E(ε) = Etot + ∆E(ε).
This shift ∆E(ε) is expanded in a series with respect to ε. The stress tensor σ isgiven as the expansion coeﬃcient in the ﬁrst-order term.
∆E(ε) = −Tr(σε)Ωcell + O(ε2 ),1 ∂E(ε).
σαβ = −
Ωcell ∂εαβ ε→0(36)

When the strain and the stress become isotropic, we have the well-known expressiondEof P = −.
dΩcellThe strain tensor εαβ is a symmetric tensor and G behaves as (1 − ε)G in theﬁrst order. Moreover, Ωcell changes as det(1 + ε)Ωcell . Since Ωcell n(G) and SI (G)are invariant, we have the next expression.
el−el
non−local
kin
xc
local
α
ewald
σαβ = σαβ
+ σαβ
+ σαβ
+ σαβ
+ σαβ
+ σαβ
+ σαβ
.

kin
σαβ
=2

∑

|ψi (k + G)|2 (k + G)α (k + G)β ,

k,G,i

1
el−el
σαβ
=−

′
∑
4π|n(G)|2

2 G

xc
σαβ
= δαβ

(

|G|2

)

2Gα Gβ− δαβ ,
|G|2

1 ∫drn(r)(µxc (n) − εxc (n)),Ωcell Ωcell

local
σαβ
=−

′
∑

n∗ (G)SI (G)G,I
non−local
σαβ= −2∂VIlocal (G)
,
∂εαβ

′
′
l
∗
′ ∂δ V̂I (k + G, k + G )SI (G − G )ψk,i (G)ψk,i (G )∂εαβk,i G,G′ I,l

∑ ∑ ∑

′

α
= δαβ n∗ (0)
σαβ

∑

αI ,

I
ewald
σαβ
−

1 ∂γewald.
Ωcell ∂εαβHere, we introduced a relative coordinate qI for the atomic coordinate, which behaves RI (t) = (1 + ε)qI (t).
126Structure of the package

The structure of ESopt has a form explained in this section at present, although anessential change might happen in future. When ESopt.tar.gz is unpacked, you havea directory ESopt. In this directory, you have subdirectories named INPUT_DATA,
mknon, commons, params. We explain contents of these directories.

6.1

ESopt

In this main directory, program ﬁles (*.f90) and makefile are contained as well
as some script ﬁles to create working directories. In ESopt, input ﬁles should be
prepared in the directory INPUT_DATA. Then, you have to do make in the main
directory Excecution of the calculation is done in the working directory, whose nameis speciﬁed by a line of makefile.

6.2

INPUT DATAIn this directory, all input ﬁles are stored. Users have to create ﬁles the main input
ﬁle CORD and some controlling parameters in *.CNTL. In CORD, the structure of thematerial and the conditions for the calculation are given by parameters, which arerequired also in the compiling step. The control ﬁles *.CNTL determine a parameter
set, which are read from these ﬁles at run time. These controlling parameters maychange the action of the excecutable ﬁle called “opt”. In this directory we have also
the pseudo-potential parameter ﬁles and data ﬁles specifying the sampling k-point
and k vectors (or lines in the Brillouin zone) used in the band structure calculation.

6.3

mknon

In this directory, there are sub-programs to determine a parameter set required in
the actual computation. Because the programs detemine internal parameters, theuser is just requested to deﬁne essential input parameters. The ﬁrst sub-programkbexe obtains the mesh for FFT, the number of Fourier components, and the cutoﬀ
using the size of the unit cell. The number of states as well as the total charge is
given by a summation of the valence electrons in each atom.
The subprogram outputs some of parameters speciﬁc to a system in a data ﬁleparameters.data, which is written in a Fortran unformatted form, Others are stored
in input.dat, which is also an unformatted ﬁle. The second subprogram mknon
obtaines the FFT mesh in a rectangular distribution in the momentum space and alist of the Fourier components in the cutoﬀ radius by deﬁning correspondence to themesh points, and outputs the arrays in gpt_list, whic is also unformatted. These
sub-programs run automatically, when the user executes make.
13

6.4

commons and params

In this package, variables and arrays are transfered among subroutines by thefollowing rule. The essential arrays are deﬁned as the global arrays in a module by
the deﬁnition of the public array, which are refered in each subroutine. Some ofvariables are written in a common statement. This is because our revision of theoriginal Fortran 77 package is not complete. Using parameters determined by mknon,
the size of each array is obtained to use in the allocation of the global array. The
data mknon/parameters.data is read by parameters.f90, while arrays are given
by globalarray.f90. Some variables commonly used in subroutines are written asa header ﬁle stored in commons. In params, some parameters including the LDAparameters, which are independent of material, are stored.

7

Overview of the action
A typical procedure to perform computation using ESopt is the following.
1. By creating INPUT_DATA/CORD, the user set the system parameters including
atomic coordinates, atomic species, and the unit cell, as well as the essentialconditions for the cutoﬀ radius, and the k-sampling method.
2. The user should check and write destination, which determines the name of
the working directory and is written in the line DESTDIR=destination at thehead of makefile. In this instruction, we omit a description how to handle
the machne dependence.
3. Do make.
4. Next, cd destination. Set the control ﬁles, SOPT.CNTL, ALLOPT.CNTL, FILE_IO.CNTL,
for the execution properly.
5. Run the prgram by ./opt, or by submitting job by qsub qsub using the job
scheduling system, SGE. In the last action, the ﬁrst qsub is the command,
while the second qsub is a script prepared when you run make.

8

Inputs and outputs

As explained in the last section, the main steps are 1) prepare INPUT_DATA/CORD,
makefile, 2) make, 3) change your directory to the working directory, 4) prepare*.CNTL. Then the preparation for the run is completed. Now instruction of inputs
and explanation of outputs are given.

14

8.1

CORD

The format of INPUT_DATA/CORD are the following. This example describes Si in
the diamond structure.
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ <<<<<<<<<<<<<<<<<<< Basic Inputs >>>>>>>>>>>>>>>>>>>>>>>>+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#Atom coordinates relative to unit vectors,
Atom,mode.
0.250000000D+00 0.250000000D+00 0.250000000D+00 Si o
0.500000000D+00 0.500000000D+00 0.500000000D+00 Si o
#^
^ ^
^ ^
^ ^^ ^
This part determines the atomic coordinates in the internal coordinate, which is
the relative coordinate with respect to the unit vectors. The input ﬁeld starts
from # and ends at the second #. The number of lines gives the number of atoms.
Each component of the relative coordinate has to be written by numbers ﬁlling an
area from ^ to ^. This is a ﬁxed format. The atomic species are given by two
characters identifying the lower-case letter and the chapital letter. The user specify
with the value as Si. A small letter in mode speciﬁes the motion of the atom, where
the atom with o is optimized and that with f is ﬁxed. The number of valence
electrons and the pseudo potential parameters like l0 for each element is written in
INPUT_DATA/TABLE_OF_ELEMENTS.
================<< Other basic inputs >>=========================________________________________________________________________| unit_vec: Primitive vectors of unit cell in real space [Ang]|| kp_file:
Sampling type of primitive F.B.Z.
||[file name with 6 char.]|| E_cut:
Cutoff energy for plane-wave [Ry]|| band_calc: Choice [yes|no|metal|bands|mband]|| f_crt:Critical force [Ry/a.u.] for judging convergency. || f_tol:Parameter for judging convergency of electronic||
deg. freedom.(Relative )|| calc_mode: Calculation mode.(bas[ic] or ext[entended])|| ideb:ideb = 1; Running verbose mode.
|~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~&basic_input
unit_vec = 2.715, 0.000, 2.715,
15

2.715, 2.715, 0.000,
0.000, 2.715, 2.715,
kp_file= ’DP10GM’,
E_cut= 8.0,
band_calc = ’metal’,
f_crt= 5.0D-03,
f_tol= 1.0D-8,
calc_mode = ’bas’,
ideb
= 0
&end
This section is given by the namelist format of Fortran.
• unit_vec : The unit vectors for the unit cell are given in
a1x a1y a1z
a2x a2y a2z
a3x a3y a3z
In this program, the reciprocal lattice vectors are deﬁned byai = (aix , aiy , aiz ),(37)

a3 × a1
a1 × a2
a2 × a3, b2 =, b3 =.
2πa1 · a2 × a3
2πa1 · a2 × a3
2πa1 · a2 × a3

(38)

and obtained asb1 =• kp_file : This line speciﬁes the ﬁle for the k-point list used in the sampling.
In this package, the list ﬁles are prepared and written in the input manually.
’DP10GM’ is for the cubic system using only the symmetry of the point groupO.
• E_cut is the cutoﬀ energy Ecut in the atomic unit, Ry.
• band_calc is a switch to change the calculation scheme. If it is ’metal’,
the occupation number around the Fermi level is smeared with the Fermi
distribution function. If ’bands’ is selected, the band structure calculation
using the mode ’bas’.
• f_crt : The convergence criteria for the interatomic force.
• f_tol : The convergence criteria for the CG iteration in the electronic statecalculation.

16

• calc_mode : The switch for the calculation mode. If it is ’bas’, optimizationof the atomic position is possible. If it is ’ext’, optimization of the unit cell
is possible with a constant pressure mode. ideb : A parameter controling theoutput (1 for the verbose mode, 0 for the ordinal mode.)+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++<<<<<<<<<<<<<< Mode dependant parameters >>>>>>>>>>>>>>>>>>>+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++==============<< Case 1;

band = yes or metal. >>=============______________________________________________________________| pwidth: Width of graduall cutoff for occupation number[Ry].|| incw:
Eigen states is sorted in every incw CG steps.
|| ncband: Number of un-occupied states.
|~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~&METALDpwidth=0.005,
incw=10,ncband=4&END
The part “Case 1” gives parameters for the mode ’metal’. pwidth determines thewidth for the Fermi distribution. incw is the number of an interval between theneighboring steps, when the program checks the order and the occupation of theenergy eigen states around the Fermi level.
====<<

Case 2;

cmode = ext. (Variable cell simulation.) >>====___[ Constant-cutoff is mimiced with following parameters. ]___||| vke:
Kinetic energy is scaled in vke where|||k+G|^2 is beyond E_cut.
|| sigma: Above scaling is faded in over sigma.
|| defmax: Maximum deformation of unit-cell allowed in this||simulation.
|~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Original
vke = 2.0,
sigma = 0.05,
defmax = 3*0.05
17

&mimicd
vke = 2.0,
sigma = 0.10,
defmax = 3*0.05&end
The part “Case 2” controls a function to reduce a possible sudden change in the listof the G vectors, when the constant-pressure molecular dynamics is done. This part
is not in use for the present package.

8.2

SOPT.CNTL“Sopt” performs the ﬁxed cell mode. The cell is not necessarily the same as theunit cell of a crystal, but may be a super cell. INPUT_DATA/SOPT.CNTL is written in
the namelist format of Fortran. The contents are the following.
&soptcntl
itmaxs=1,
f_crt=0.005,
ftol = 1.0d-9,
itmaxd=200,
delta = 0.4,
ideb = 0&end• itmaxs : The maximum allowed steps for the structural optimization.
• f_crt : The convergence criterion for the structural optimization. When the
maximum of the absolute value of an interatmic force becomes less than thisvalue, the convergence is assumed to be achieved. The unit is A.U.
• ftol : The convergence criterion for the electronic state calculation. Thisvalue is compared with the relative change in the total energy per a CG step.
• itmaxd : The maximum value of CG steps in the electronic state calculation.
• delta : The width of the optimization step for the atomic structure. This is
a proportional constant for the interatomic force.

18

8.3

ALLOPT.CNTL“Allopt” performs the variable cell mode. Using this mode, the full structural optimization and the constant-pressure molecular dynamics are realized. The contentsof INPUT_DATA/ALLOPT.CNTL are the following.
0. Chose a mode.
mode = ’md|op’: Molecular Dynamics or Optimization.
In case that MD is chosen, you decide whether "Nose’s constanttemparature formula" is used or not.
&modesw
mode=’md’,
nose = ’yes’,&end• Using functions of “Allopt”, refering the interatomic force and the internal
stress, we can make the simulator to perform the molecular dynamics by
mode=’md’, or the structural optimization by mode=’op’. When nose=’yes’
is speciﬁed, the Nosé thermostat is used to control the temperature. If nose=’no’,
opt performs a molecular dynamics without temperature control.
0-a. Mode dependant parameters.
nose:
unit:
temp:
q_inv:

In case that nose = ’yes’, constant temparature simulation
with following parameters is performed.
The unit of temparature.[au|K|eV]
Temperature.
This value corresponds to Q^(-1).
In case that q_inv = 0, then zeta = 0.
See, Nose, Solid State Physics, Vol. 24, pp. 98-106(1989).
Nose, Mol. Phys., Vol. 52, p 255(1984).
q_inv = 0.05??

&noseparam
unit = ’K’,
temp = 100.0,
q_inv = 0.01,&end
p_crt: Convergence parameter to judge convergence of pressure in [GPa].
19

gamma: Friction ratio. gamma is in [0-1]. If gamma =1, MD is done.
&optparam
f_crt=1d-10,
ftol=1d-8,
p_crt=0.0001,
gamma = 0.5,&end
• unit determines the unit of temperature, which may be in Kelvine or in electron volt (eV). The value is given by temp. Another parameter q_inv is theNosé parameter.
• f_crt The convergence criterion of the interatomic force for the structural
optimization. The unit is in the atomic unit.
• ftol The convergence criterion of the total energy of the electronic state.
• p_crt The convergence criterion of the pressure for the structural optimization. The unit is in Ry.
• gamma A damping factor of the velocity for both the atoms and the cells.
===========================================================================
1. Define external environment.
-------------------------------------------------------------------------1-1. Pressure
iso: A switch to select isotropic or anisotropic pressure.
Case; iso = ’no’, pext_3d_[i|f](1:3) is read.
Case; iso = ’yes’, pext_3d_[i|f] is not read.
pext_ini: initial pressure
pext_fin: final pressure&extpress
iso=’no’,pext_init=0.0,
pext_final=0.0,
pext_3d_i=0.0,0.0,0.0,
pext_3d_f=0.0,0.0,0.0,
&end

20

• iso determines whether the isotropic pressure is applied or not.
• A set of pext_init and pext_final or another set of pext_3d_i and pext_3d_f
we can specify the initial and ﬁnal values of the external pressure in the molecular dynamics.
--------------------------------------------------------------------------1-2. Parameters concernig how the pressure is introduced and how MD runs.
nsec_p:nstp_i:
nstp_g:External pressure is gradually enforced.
Sectioned by nsec_p steps between pext_fin and pext_ini.
nstp_i MD steps are performed at the pext_ini.
nstp_g MD steps are performed at each pressure
except for pext_ini.

&prscntl
nsec_p=10,
nstp_i=10,
nstp_g=1,&end• Using nsec_p, nstp_i, and nstp_g we can control the change in the externalpressure.
==========================================================================
3. Parameters for controlling evolution.
delt_t:ww:itmaxa:itmaxd:Time interval in femto second.
Fictitious mass ratio [w/mass]
Maximam iterations for MD.
Maximam iterations in diagw.

&mdparams
delt_t=5.0,
ww = 100000.0,
itmaxa=500,
itmaxd=200,&end• delt_t：The time step for the molecular dynamics.
21

• ww：The imaginative (artiﬁcial) mass for the unit cell.
• itmaxa：The total steps of the molecular dynamics.
• itmaxd：The maximum steps for the CG steps in the electronic structure
calculation.
=========================================================================
4. Initial condition.
ekin_init:
Initial kinetic energy is scaled in this value.
init_conf:
Switch how is the initial velocities are given.
Case; init_conf = ’given’, read from iveloc namelist in this file.
Case; init_conf = ’rando’, random vectors, generated with ’iseed’,is used.
ivpoint:
Initial kinetic energy is given at ivpoint step.
&initveloc
ekin_init = 10.0,
ckin_init = 0.0,
init_conf = ’given’,
iseed = 180,
ivpoint = 1,
&end-----------------------------------------------------------------------4-1. Define initial direction to move.
This parameters are read only in case init_conf = ’given’.
how:Case; how = ’rel’, final conf - initial conf is used.
Case; how = ’abs’, final conf. is used as a vector.

latv_i(3,3): initial primitive vector in real space [Ang]latv_f(3,3): final primitive vector in real space [Ang]coordi(3,natm): initial coordinate
coordf(3,natm): final coordinate
Example for structural transformation from cubic-diamond to beta-tin
of silicon.
&iveloc22

how = ’rel’,latv_i =2.715, 0.000, 2.715,
2.715, 2.715, 0.000,
0.000, 2.715, 2.715,
latv_f =2.715, 0.000, 2.715,
2.715, 2.715, 0.000,
0.000, 2.715, 2.715,
coordi =0.000000000, 0.000000000,
0.500000000, 0.500000000,
coordf =0.000000000, 0.000000000,
0.500000000, 0.500000000,&end

0.000000000,
0.000000000,
0.000000000,
0.000000000,• This part describes how to ﬁx the initial velocity. In ESopt, this ﬁnal part ofcoordi and coordf are skipped.

8.4

FILE IO.CNTL

The contents of INPUT_DATA/FILE_IO.CNTL are the followings.
======================================================================
1. Parameter(s) to controll output.

&file_io_cntl
read_file=’no’,
initial_file=1&end

intchg: interval of charge output&outcntl
intchg = 10,
intefl = 20,
intwf = 10,&end23

2. Debug Flags&debug_fl
ideb = 0&end• read_file : This parameter speciﬁes a way to read the result of the atomic
conﬁguration and the wavefunction determined in a former calculation. Oneof ’cont’, ’wyes’, ’cyes’, ’syes’ or ’no’ are selected. The action for each
option is the following.
– ’cont’: The structure given by STRUCT/strct and the wavefunction in
fort.98 are read.
– ’wyes’: The wavefunction in fort.98 is read.
– ’cyes’: Structural data in STRUCT/cell-XXX and STRUCT/latcrd-XXX
is read. Here, XXX means a number given in the format I3.3, which is
speciﬁed by initial_file. The wavefunction in fort.98 is read.
– ’syes’: Structural data in STRUCT/cell-XXX and STRUCT/latcrd-XXX
is read. Here, XXX means a number given in the format I3.3, which is
speciﬁed by initial_file.
– ’no’: The structure is given by CORD in INPUT_DATA and the wavefunction
prepared in the program is used.
• initial_file : When read_file=’cyes’ or read_file=’syes’, this value
gives XXX. If read_file has other value, this option is neglected.
• intchg : The interval for the steps, when the charge density is output.
• intefl : The interval for the steps, when MONIT/energy is output. If the user
wants to reduce the time cost for IO, this value should be increased.
• intwf : In an older version, this value speciﬁed the interval for output of
MONIT/wave000 storing the wave function. In the present version, the value
is given in fort.99. If the user wants to use this information to continue thesimulation, this ﬁle may be copied into fort.98.

8.5

kpoint

The contents of kpoint is the following.
24

number of line (i4) =7^^^^

G X K G L K W X
number of mesh (i4) = 10
1.000000000 0.000000000
1.000000000 0.500000000
number of mesh (i4) = 10
1.000000000 0.500000000
0.750000000 0.375000000
number of mesh (i4) = 100.750000000 0.375000000
0.000000000 0.000000000
number of mesh (i4) = 10
0.000000000 0.000000000
0.500000000 0.500000000
number of mesh (i4) = 100.500000000 0.500000000
0.750000000 0.375000000
number of mesh (i4) = 100.750000000 0.375000000
0.750000000 0.250000000
number of mesh (i4) = 100.750000000 0.250000000
0.500000000 0.000000000

0.000000000
0.500000000
0.500000000
0.375000000
0.375000000
0.000000000
0.000000000
0.500000000
0.500000000
0.375000000
0.375000000
0.500000000
0.500000000
0.500000000

The value of number of line determines the number of k lines in the Brillouinzone, on which the energy dispersion is obtained. The number is limited by 20.
Assume here that the number is 7 as above. These k lines have to form a continuous
loop in the zone. The next ﬁlled line below an empty line shows the name of thespecial points determining two ends of each k lines. If you want, you may use
conventional names owing to their symmetry. Next, we use three lines per each kline. The number of divisions of the k line and two end points in the reciprocallattice of b1 , b2 , b3 are given.

8.6

Outputs

In the calculation, a sub-directory named by MONIT is created in the working
directry, and in this directory there apperas temporal information on the energy
etc.. The total energy is in energy. The wavefunction may be given in intwf,
if properly speciﬁed in the code. (This option is switched of in ESopt.) When25

the structural optimization is ﬁnished, wave999 for the wavefunction is given. The
format of the wavefunction data is unformatted.
The charge density is stored in another subdirectory CHARGE. The ﬁle name ischgdnsXXX.
The major output is given in the following ﬁles.
• For Sopt, see the directory RESULTS.
– cord***: Atomic coordinates in the real space.

– eig***: The eigen values of the Kohn-Sham equation.

– gradient*** : The interatomic force.
– H***: The diagonal elements of a 3× 3 matrix given by the unit
vectors, the total energy and the volume of the cell.
– Htot***: Each contribution for the total energy: 1. the kinetic
energy of the electron system, 2. the Hartree term, 3. the exchangecorrelation term, 4. the potential energy of the local potential part, 5.
the potential energy given by the non-local pseudo-potential part, 6. theα term, 7. the Ewald summation, 8. and the summation from 1 to 7.
• For Allopt, see the directory STRUCT.
– cell***: The lattice vectors.

– latcrd***: The internal coordinates.

– press***: The external pressures.

– realcrd***.xyz: The atomic coordinates in the real space.
– result***: The total energy same as Htot*** for Sopt, the enthalpy, the stress tensor, and the atomic forces.
– strct***: The parameters of the unit cell, as cell***, latcrd,
and zeta in this order.
– strn***: The strain tensor.

26

9

Source codes of the program

In this section, we show some example of the implementation of each calculation
procedure found in the subroutines, in order to provide information for understanding the action of the program.

9.1

Kinetic energy

The kinetic energy is given by a function ekin(inw) in etotal.f90.
!**********************************************************************
double precision function ekin(inw)use parameters
use globalarray
implicit none
include ’commons/common.h’integer :: ik,ib,ipw,inw
real(8) frm,z,fd,x!

-----<< definition of fermi distribution function >>----frm(z) = 1d0 + vke_mim/( 1d0 + exp( -1d0*beta_mim*z ) )! ----------<< calc. of KE (mimicing) >>---------!* epw[Ry],gk[a.u.],nadd
ekin = 0d0
do ib = 1,mxband_g
do ik = 1,nk_g
fd = fdist_g(ib,ik)
do ipw = 1,mxpwk(ik)x = 0.5d0*( gk_g(ipw,ik) - ec_mim )ekin = ekin &+ gk_g(ipw,ik)*frm(x) &* ( dreal(wvfn_g(ipw,ib,ik,inw))**2 &+ dimag(wvfn_g(ipw,ib,ik,inw))**2 ) * fd
end doend doend do27

ekin = 2d0*ekin/dfloat(natm_g)/effnk
returnend• inw: To realize the CG steps, the wavefunction is stored in a few arrays. This
index speciﬁes the array for the wavefunction.
• ib, ik, ipw: These indices represent the band index, k, and G. The list of G
appearing in the planewave exp(i(k + G) · r)) is given for each k so that theenergy of the plane wave is below Ecut .
• gk_g(ipw,ik): The array to store |k + G|2 .
• wvfn_g(ipw,ib,ik,inw): The array to store the wavefunction ϕk,n (G).
In the variable ekin, we have
Ekin = 2∑|k + G|2 |ϕk,n (G)|2 .

k,n,G

The value of ekin is thus the energy per a unit cell in Ry. Here the electron mass
becomes m = 1/2. The spin degrees of freedom is just taken into account by multiplying the factor 2. Another factor frm(x) represents a smearing method to achievea constant-cutoﬀ simulation for the constant-pressure variable-cell simulation.

9.2

Hartree energy

The Hartree energy is given by a function eee(vcell) in etotal.f90.

!**********************************************************************
double precision function eee(vcell)use parameters
use globalarray
implicit none
include ’commons/common.h’integer :: igreal(8) :: vcell
eee = 0d028

do ig = 2,ng_g
eee = eee + (dreal(chgg_g(ig))**2+dimag(chgg_g(ig))**2)*gg_g(ig)end doeee = eee*4d0*pi/vcell/dfloat(natm_g)
returnend• chgg_g(ig): An array to store the Fourier transform of the charge n(G).
• gg_g(ig): An array storing 1/|G|2 .
We have the value of the Hartree energy
Eel−el =′4π ∑|n(G)|2,Ωcell G |G|2

in eee. The value of eee is also given per a unit cell in Rydberg (Ry). The electron
charge is thus e2 = 2. The reason for the factor of Ω−1cell is that Eel−el is originallygiven in the double integral in the real space.

9.3

Total energy

The total energy is given by a subroutine etotal(vcell,inw,etot) in etotal.f90.
!**********************************************************************
subroutine etotal (vcell,inw, etot)
use parameters
use globalarray
implicit none
integer :: lmx,inw,i,jreal(8) :: etot(8),vcell,vatm, &ekin,eee,exc,evl,evnl,t1,t2,dsecnd
parameter( lmx = 2 )include ’commons/common.h’
t1 = dsecnd(’cpu’)vatm = vcell/dfloat(natm_g)29etot(1) = ekin(inw)etot(2) = eee(vcell)etot(3) = exc(vcell)etot(4) = evl()
etot(5) = evnl(inw)etot(8) = etot(1)+etot(2)+etot(3)+etot(4)+etot(5)+etot(6)+etot(7)t2 = dsecnd(’cpu’)
cpu_time(1) = cpu_time(1) + t2-t1
returnend

The energy is given for each contribution and the meaning of these parts are thefollowings.
• ’ ekin= ’,etot(1) : The kinetic energy.

• ’ eee= ’,etot(2) : The Hartree energy.

• ’ exc= ’,etot(3) : The exchange-correlation energy.

• ’ evl= ’,etot(4) : The local-potential contribution.

• ’ evnl= ’,etot(5) : The contribution from the non-local potential.

• ’ ealp= ’,etot(6) : The αI term.

• ’ eewd= ’,etot(7) : The ion-ion interaction (The Ewald summation).

• ’ etot= ’,etot(8) : The total energy.

9.4

Sopt control routine

The main routine is given to switch the calculation mode by calling a subroutine.
• Allopt : The optimization of the lattice structure and the atomic structure.
The constant-pressure molecular-dynamics is possible.
• Sopt : The optimization of the atomic structure with a ﬁxed unit cell.
• DMD : The damped dynamics with a damping factor.
ESopt is provided by stopping some functions except for Sopt and limited Allopt.
Here, we consider the function of sopt. The code below is given by omitting someextra comments.
30

!**********************************************************************
subroutine sopt()
use parameters
use globalarray
implicit real*8(a-h,o-z)
include ’commons/lattice.h’include ’commons/common.h’logical(4) :: exist1, exist2
real(8) :: etot(8),delta,cpu0,cpu1
namelist / soptcntl / itmaxs, ftol, f_crt, itmaxd, delta, ideb!==================================================================
inquire(file=’SOPT.CNTL’,exist=exist2)
inquire(file=’../INPUT_DATA/SOPT.CNTL’,exist=exist1)
if(exist2) thenopen(10,file=’SOPT.CNTL’,form=’formatted’, &iostat=ios,status=’old’)read(10,soptcntl,iostat=ios)close(10)else
if(exist1) thenopen(10,file=’../INPUT_DATA/SOPT.CNTL’,form=’formatted’,&iostat=ios,status=’old’)read(10,soptcntl,iostat=ios)close(10)elsestop ’Can not open SOPT.CNTL’
endif
endif
!

----------< generate g-points >--------------call gptgen()!

-----< calc. some components of total energy >---------call vinit()!

-----< calc. ewald sum (independent part with coordinates of atoms ) >----calliniewd()31

! Main loop starts.
do its = 1, itmaxs! ----------< optimization for electronic state >---------call getll(etot,nloop)
call charge(0)! ----------< calc. Hellmann-Feynman force >---------call force(vcell)
call eforce(natm_g,stoa,grd,frc,fn)
call wfile(etot,its)! ----------< output data >---------call wfile0(000,its)! ----------<

judgement of convergence >-------------

if(ideb.eq.1) write(*,*)’fn, f_crt = ’,fn,f_crt
if(fn.lt.f_crt) exit
! ----------< move all atoms >-------------------------------call mvatm(delta)
end do! Main loop ends.
if(fn.lt.f_crt) thencall wfile0(999,its)
end ifwrite(*,*)’elapse time in sopt = ’, cpu1-cpu0
returnend subroutine sopt

We can see that the algorithm consists from the following steps.
1. Initialization.
(a) Input SOPT.CNTL.
(b) gptgen() : Creation of lists for G, etc.
32

(c) vinit() : Input the potential data.
(d) iniewd() : Initialization of parameters for the Ewald calculation.
2. The optimization loop.
(a) getll(etot,nloop) : The electronic structure determination process.
(b) charge(0) : The charge density determination step.
(c) force(vcell) : The force determination step.
(d) eforce(natm_g,stoa,grd,frc,fn) : An optimization of force for a
special usage.
(e) wfile(etot,its), wfile0(999,its), wfile0(000,its) : The output
routines.
(f) mvatm(delta) : The routine to obtain new atomic coordinates.

9.5

Subroutine for diagonalization

In the subroutine diagw, the CG method is applied to obtain the electronic structure.

subroutine diagw (nloop,vcell,alp,etot)
use parameters
use globalarray
implicit none
include ’commons/common.h’integer :: nloop, inw, itmax, ik, in, ig, ist, loop, its, nwf, nnn, &istc, ireal(8) :: etot(8), vcell, alp, eps, ang, d1, d2, dnf, step, ep, e_p,&fret, fp, gd, avm, dmm, t1, t2, dsecnd,e2,e1
complex(8) :: ggc,gam,s, chg1
parameter (eps = 1d-10)
parameter (inw = 2)
parameter (e_p = 5d0)
! -----< calc. charge density for G-space >----call schmit(inw)33call toreal(inw)chg1 = chgg_g(1)etot(6)=alp*dreal(chg1)! -----< calc. total energy(fp) >----call etotal (vcell,inw,etot)
fp = etot(8)! -----< calc. force >----call wforce(vcell,inw)
call precnd(e_p)
call schmi2(inw,1)
itmax = max(itmaxd*2,npw_g*mxband_g*nk_g*4)! nwf(= npw*nband*nk) is num. of wave func.
! -----< store old data of w.f. >----do in = 1, mxband_g
do ik = 1, nk_g
do ig = 1, mxpwk(ik)wvfn_g(ig,in,ik,4)=wvfn_g(ig,in,ik,1)
wvfn_g(ig,in,ik,5)=wvfn_g(ig,in,ik,1)
end doend doend do! -----< iteration start >----ist = 1
istc = 0loop = 0! Main loop starts
do its = 1, itmaxs = (0d0,0d0)
d1= 0d0
d2= 0d0! -----< summation of some value >----nwf = mxband_g*nk_g34

do in = 1,mxband_g
do ik = 1,nk_g
do ig = 1,mxpwk(ik)s = s + dconjg(wvfn_g(ig,in,ik,1))*wvfn_g(ig,in,ik,2)
d1= d1+dreal(wvfn_g(ig,in,ik,1))**2+dimag(wvfn_g(ig,in,ik,1))**2
d2= d2+dreal(wvfn_g(ig,in,ik,2))**2+dimag(wvfn_g(ig,in,ik,2))**2
end doend doend doang = acos(cdabs(s)/sqrt(d1*d2))/pi
dnf = sqrt( d1/nwf )if(ang.lt.0.46d0.and.ang.ge.0.385d0) thencall schmi2(inw,2)elseif(ang.lt.0.385d0)thencall schmi2(inw,1)
endif
call linmin(vcell,etot,fret,ist)if(ist.eq.0) istc=istc+1
step=2d0*dabs(fret-fp)/(dabs(fret)+dabs(fp)+eps)!

----------< judgement for convergency >---------if((step.le.ftol.and.ist.eq.1).or.its.gt.itmaxd.or.istc.gt.3) thenloop = loop + 1
if(loop.gt.2) thenetot(6)=alp*dreal(chg1)
call efile(etot,ang,dnf,its,ftol,step,loop)
nloop = its
write(*,*) ’Count (its,loop)’, its,loop
returnelse
ist = 0
endif
endif

!

---------------------------------------------------

35ep=fret-fp
fp = fretcall wforce(vcell,inw)
call precnd(e_p)
call eigenv(avm,dmm,its,nnn)
if(ist.eq.0.or.nnn.eq.1) thencall toreal(inw)chg1 = chgg_g(1)call wforce(vcell,inw)
call precnd(e_p)
endif
if(its.eq.0) ist = 0
if(ist.eq.1) thengd = 0d0
ggc = (0d0,0d0)
do in = 1, mxband_g
do ik = 1, nk_g
do ig = 1, mxpwk(ik)gd = gd + dreal(wvfn_g(ig,in,ik,4))**2 &+ dimag(wvfn_g(ig,in,ik,4))**2
ggc = ggc + dconjg(wvfn_g(ig,in,ik,1) &- wvfn_g(ig,in,ik,4))*wvfn_g(ig,in,ik,1)
end doend doend doif(gd.eq.0d0) thenwrite(*,*)’warning. really gd = 0 ?’returnendif
gam = ggc/gd
do in = 1, mxband_g
do ik = 1, nk_g
do ig = 1, mxpwk(ik)wvfn_g(ig,in,ik,4)=wvfn_g(ig,in,ik,1)
end doend doend dodo in = 1, mxband_g
do ik = 1, nk_g
do ig = 1, mxpwk(ik)36wvfn_g(ig,in,ik,5)=wvfn_g(ig,in,ik,4) &+ gam*wvfn_g(ig,in,ik,5)
end doend doend dodo in = 1, mxband_g
do ik = 1, nk_g
do ig = 1, mxpwk(ik)wvfn_g(ig,in,ik,1)=wvfn_g(ig,in,ik,5)
end doend doend doelseif(ist.eq.0) thendo in = 1,mxband_g
do ik = 1,nk_g
do ig = 1,mxpwk(ik)wvfn_g(ig,in,ik,4)=wvfn_g(ig,in,ik,1)
wvfn_g(ig,in,ik,5)=wvfn_g(ig,in,ik,1)
end doend doend doist = 1endif
nnn = mod(its,intefl)
if(nnn.eq.0) thenif(ideb.eq.1) write(*,*) ’efile in’call efile(etot,ang,dnf,its,ftol,step,loop)
if(ideb.eq.1) write(*,*) ’efile out’
endif
end do! Main loop ends
stop ’Maximum iterations exceeded (diagw).’
end subroutine diagw

37