Data structure of OpenMX and ADPACK• Structure of ADPACK• Solving the 1D Dirac equation• Pseudopotentials• Localized basis functions• Structure of OpenMX• Input_std()• Total energy• Data structure for parallelization
Taisuke Ozaki (Univ. of Tokyo)
OpenMX developer’s meeting, Dec. 21-22, 2015

What is ADPACK?
ADPACK (Atomic Density functional program PACKage) is a software to
perform density functional calculations for a single atom
The features are listed below:• All electron calculation by the Schrödinger or Dirac equation• LDA and GGA treatment to exchang-correlation energy• Finite element method (FEM) for the Schrödinger equation• Pseudopotential generation by the TM, BHS, MBK schemes• Pseudopotential generation for unbound states by Hamann's scheme• Kleinman and Bylander (KB) separable pseudopotential• Separable pseudopotential with Blöchl multiple projectors• Partial core correction to exchange-correlation energy• Logarithmic derivatives of wave functions• Detection of ghost states in separable pseudopotentials• Scalar relativistic treatment• Fully relativistic treatment with spin-orbit coupling• Generation of pseudo-atomic orbitals under a confinement potential
The pseudopotentials and pseudo-atomic orbitals can be the input data for OpenMX.

ADPACK is freely available

ADPACK can be available under GNU-GPL.

http://www.openmx-square.org/

Programs of ADPACK
Programs: 65 C rounties and 5 header files (50,000 lines)
Link:
LAPACK and BLAS
Main routine: adpack.c
Input: readfile.c, Inputtool.c
Output: Output.c
All electron calculations: All_Electron.c, Initial_Density.c, Core.c
Numerical solutions for Schroedinger and Dirac eqs.: Hamming_I.c, Hamming_O.c
Density: Density.c, Density_PCC.c, Density_V.c
Exchange-Correlation: XC_CA.c, XC_EX.c, XC_PW91.c, XC_VWN.c, XC_PBE.c

Mixing: Simple_Mixing.c
Pseudopotentials: MBK.c, BHS.c, TM.c
Pseudo-atomic orbitals: Multiple_PAO.c
The global variables are declared in adpack.h.

1D-Dirac equation with a spherical potential
1-dimensional radial Dirac equation for the majority component G is given by

The mass term is given by

Minority component

By expressing the function G by the following form,
One obtain a set of equations:

Lの満たすべき条件

The charge density is obtained fromSolving the 1D-Dirac equation
By changing the varial r to x with, and applying a predictor andcorrector method, we can derive the following equations:

For a given E, the L and M
are solved from the origin anddistant region, and they arematched at a matching point.

LM

In All_Electron.c, the calculation is performed.

Scattering by a spherical potential
Incident wave

eiq r

1 2
 q
2
Scattered wave
eiqriqr
 sca ( , r )  e (2l  1)eil sinl Pl (cos  )qr lPhase shift  ld
jl (kr ) |r0  Dl ( ) jl (kr0 )drtan l ( , r0 ) 
dr0 nl (kr ) |r0  Dl ( )nl (kr0 )drLogarithmic derivative of ψdDl ( , r )  r ln l (r )drr0r0
222Dl ( , r ) |r 0  
drr|
(r)|ld
r0 l 2 (r0 ) 0

If the norm of pseudized wave is conserved within r0 and thelogarithmic derivative coincides with that for the all electron case,the phase shift coincides with the all electron case to first order.

Norm-conserving pseudopotential by Troullier and Matins
N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993 (1991).

For, the following form is used.

Putting ul into radial Schroedinger eq. and solving it with respect to V, we havec0 ～c12 are determined by the following conditions:• Norm-conserving condition within the cutoff radius• The second derivatives of V (scr) is zero at r=0• Equivalence of the derivatives up to 4th orders of ul at the cutoff radius

Ultrasoft pseudopotential by Vanderbilt
D. Vanderbilt, PRB 41, 7892 (1990).

The phase shift is reproduced around multiple referenceenergies by the following non-local operator.

If the following generalized norm conserving condition is fulfilled,
the matrix B is Hermitian. Thus, in the case the operator VNL is also
Hermitian.
If Q=0, then B-B*=0

How the non-local operator works?
Operation of the non-local operator to pseudized wave function
Note thatIt turns out that the following Schroedinger equation is satisfied.

The matrix B and the generalized norm conserving condition
The matrix B is given byThus, we have

By integrating by parts
・・・(1)
By performing the similar calculations, we obtain for the all electron wave functions
・・・(2)
By subtracting (2) from (1), we have the following relation.

Norm-conserving pseudopotential by MBK
I. Morrion, D.M. Bylander, and L. Kleinman, PRB 47, 6728 (1993).

If Qij = 0, the non-local operator can be transformed to a diagonal form.

The form is exactly thesame as that for the Blochl
expansion, resulting in noneed for modification ofOpenMX.
To satisfy Qij=0, the pseudized wave function is written by

The coefficients can be determined by matching up to the third derivatives tothose for the all electron, and Qij=0. Once c’s are determined, χ is given by

The form of MBK pseudopotentials
The pseudopotential is given by the sum of a local term Vloc and non-local term VNL.

V(ps) Vloc (r)  VNL

The local term Vloc is independent of the angular channel l. On the other hand,the non-local term VNL is given by projectors

VNL   i | i i |i
The projector consists of radial and spherical parts, and depends on species,
radial and l-channels.

Optimization of pseudopotentials(i) Choice of parameters
1.
2.

3.
4.

Choice of valence electrons (semi-core included?)
Adjustment of cutoff radii by monitoring shape ofpseudopotentials
Adustment of the local potential
Generation of PCC charge

Optimization of PP
typically takes a halfweek per a week.

(ii) Comparison of logarithm derivatives
If the logarithmic derivatives
for PP agree well with thoseof the all electron potential,
go to the step (iii), or returnto the step (i).

good

No good(iii) Validation of quality of PP by performinga series of benchmark calculations.

goodNo good

Good PPDatabase (2013)

Optimized VPS and PAO
are available, which
includes relativistic effectwith spin-orbit coupling.

http://www.openmx-square.org/

Close look at “vps” files #1
Input fileIn the header part, the input filefor the ADPACK calculations
are shown, which maybe
helpful for the next generationof pseudopotentials.

Close look at “vps” files #2
Eigenvalues for all electron calculation

The eigenvalues with
j=l±1/2 for the allelectron calculations by
the Dirac equation areincluded, which can beused to estimate thesplitting by spin-orbit
coupling

Close look at “vps” files #3
Information for pseudopotentials
The specification for thepseudopotentials ismade by vps.type,number.vps, andpseudo.NandL.
The project energies λ is shown as follows:Close look at “vps” files #4
The generated pseudopotentials are output by Pseudo.Potentials

1st column: x
2nd column: r=exp(x) in a.u.
3rd column: radial part of local pseudopotential
4th and later columns: radial part of non-local pseudopotentials.

Close look at “vps” files #5
Charge density for partial core correction

1st column: x
2nd column: r=exp(x) in a.u.
3rd column: charge density for PCC

Primitive basis functions
1. Solve an atomic Kohn-Sham eq.
under a confinement potential:
s-orbital of oxygen

2. Construct the norm-conserving
pseudopotentials.
3. Solve ground and excited states for thethe peudopotential for each L-channel.

In most cases, the accuracy and efficiency can be controlled byCutoff radius
Number of orbitals

PRB 67, 155108 (2003)
PRB 69, 195113 (2004)

Variational optimization of basis functions
One-particle wave functions

Contracted orbitals

The variation of E with respect to c with fixed a gives→Regarding c as dependent variables on a and assuming KS
eq. is solved self-consistently with respect to c, we have

Ozaki, PRB 67, 155108 (2003)

Optimization of basis functions
1. Choose typical chemical environments

2. Optimize variationally the radial functions

3. Rotate a set of optimized orbitals within the subspace, anddiscard the redundant funtions

××

Close look at “pao” files #1
Contraction coefficients

The optimized PAO are
generally obtained by two or
three calculations for orbital
optimization. The contraction
coefficients are attached in theheader of the pao file.

Close look at “pao” files #2
Input file

The input file for the ADPACK
calculations is attached which
can be helpful for checking thecomputational condition.

Close look at “pao” files #3
e.g., Si7.0.pao

14825376The eigenvalues of pseudized and confined states are shown.
The information can be used to choose basis functions step by step.

One may find that thestructural properties canbe obtained by Si7.0s2p2d1.
Also, it is found that thecutoff of 7.0 (a.u)
reaches the convergent
result from thecomparison between
Si7.0-s3p2d2f1 andSi8.0-s3p2d2f1.

Close look at “pao” files #4
Valence density

1st column: x
2nd column: r=exp(x) in a.u.
3rd column: valence density for the chosen states

Close look at “pao” files #4
The generated radial functions are output bypseudo.atomic.orbitals.L=0,1,…

1st column: x
2nd column: r=exp(x) in a.u.
3rd and later columns: radial part of basis functions

Δ-factor
The delta factor is defined as difference of total energy between Wien2k(FLAPW+LO) and a code under testing, which is shown as shaded
region in figure below, where the volume is changed by plus and minus
6 % taken from the equilibrium V0 .

Lejaeghere et al., Critical Reviews in Solid State and Materials Sciences 39, 1-24 (2014).

Comparison of codes in terms of Δ-factor

http://molmod.ugent.be/deltacodesdft

Structures of OpenMX••••••••

Language: C, fortran90
265 sub files, about 1000 sub routines
21 header files
About 300,000 lines
Compilation by makefile
Eigenvalue solver: ELPA included
Linking of LAPACK, BLAS, FFTW3
Hybrid parallelization by MPI,OpenMP

Basic structure

openmx_common.hCommon functions andcommon variables aredeclared in the header file.

openmx()Main computational flow

Input_std()

Reading input filetruncation()

Analysis of geometrical structure,
memory allocation, analysis ofcommunication patternDFT()

SCF, total energy, forces

MD_pac()Molecular dynamics andgeometry optimization

OutData()

OutputInput_std.c
The input file is analyzed in Input_std() which employs Inputtool.c. After searching
keywords, the value after the keyword is set for the keyword. The variables for thekeywords are declared in openmx_common.h.

Inputtool.c allows us to write the input file in arbitrary order for the keywords.

Implementation: Total energy (1)The total energy is given by the sum of six terms, and a proper integration
scheme for each term is applied to accurately evaluate the total energy.

Kinetic energy
Coulomb energy with external potential

Hartree energy

Exchange-correlationenergyCore-core Coulomb energy
TO and H.Kino, PRB 72, 045121 (2005)

Implementation: Total energy (2)The reorganization of Coulomb energies gives three new energy terms.
sThe neutral atom energy

Short range and separable to twocenter integrals

Difference charge Hartree energy

Long range but minor contribution

Screened core-core repulsion energy

Short range and two-center
integrals
Difference charge

Neutral atom potentialImplementation: Total energy (3)So, the total energy is given by

Each term is evaluated by using a different numerical gridwith consideration on accuracy and efficiency.
The relevant subroutines
Set_OLP_Kin.c, Total_Energy.c}}

Spherical coordinate inmomentum space

Set_ProExpn_VNA.c, Total_Energy.c
Set_Nonlocal.c, Total_Energy.cReal space regular mesh

Poisson.c, Total_Energy.c
Set_XC_Grid.c, Total_Energy.cReal space fine mesh

Total_Energy.c

Atomic 3D atomic partitioning
How one can partition atoms tominimize communication and memory
usage in the parallel calculations ?
Requirement:• Locality• Same computationalcost• Applicable to anysystems• Small computational
overhead
T.V.T. Duy and T. Ozaki, Comput. Phys.
Commun. 185, 777-789 (2014).

Modified recursive bisection
If the number of MPI processes is 19, then the following binary
tree structure is constructed.

In the conventional recursive bisection, the bisection is made so thata same number can be assigned to each region. However, themodified version bisects with weights as shown above.
This is performed in Set_Allocate_Atom2CPU.c

Reordering of atoms by an inertia tensor

This is performed in Set_Allocate_Atom2CPU.c

Recursive atomic partitioning
The method guarantees
 Locality of atomicpartitioning Balanced
computational cost
 Applicability to anysystems
 Small computational
cost

T.V.T. Duy and T. Ozaki, CPC 185, 777 (2014).

Allocation of atoms to processess
Diamond 16384 atoms, 19 processes

Multiply connected CNT, 16 processes

Definition of neighboring atoms
Each atom has strictly localizedbasis functions. Thus, the nonzero overlap between basis
functions occurs if r12 < (rc1 + rc2).

rc2

The analysis is performed in
Trn_System() of truncation.c, andrelevant variables for the
information areFNAN[]: # of neighboring atoms
natn[][]: global index of theneighboring atomncn[][]: cell index of theneighboring atom

r12
rc1

Three indices for atoms
Global index: if there are N atoms in the unit cell, then each atom has a global
index which is within 1 to N.
interMediate index: a set of atoms (Np atoms) are assigned to each
MPI process. Then index within the MPI process each atom has anintermediate which is within 1 to Np.

Local index: The second index of natn[][] and ncn[][] runs for a local
index which is within 1 to FNAN[].

Example: the three indices of atoms
Only the non-zero Hamiltonian matrix elements are stored in H. An example isgiven below to show how the conversion among the three indices is made.

where Mc_AN is the intermediate index, and Gc_AN is the global index, h_AN isthe local index, respectively.
M2G, F_G2M, and S_G2M can be used to convert the indices.

Cutoff energy for regular mesh
The two energy components Eδee + Exc are calculated on real space
regular mesh. The mesh fineness is determined by plane-wave cutoff
energies.

The cutoff energy can be related to the meshfineness by the following eqs.

Partitioning of grids

Uniform grid is used to calculate
matrix elements and solve Poisson’sequation. A hundred million grid
points for a few dozen thousandatoms.
A proper one of four data structures
for grid is used for each calculation.
They are designed to minimize theMPI communication.
These data structure are all constructed in Construct_MPI_Data_Structure_Grid() and
Set_Inf_SndRcv() of truncation.c.

Case study #1
Values of basis functions are calculated on grids in the structure A.

In Set_Orbitals_Grid.c

Nc runs over grids inthe structure A.

Case study #2:
2D-parallelization of 3D-FFT in Poisson.c→→ccabaCompared to 1D-parallelization, no increaseof MPI communication up to N. Even at N2,
just double communication.

bCase study #3
Each MPI process calculates a part of the energy terms evaluatedon the regular mesh using the structure B(ABC) as follows:

In Total_Energy.c
BN runs over grids inthe structure B(ABC).

Case study #4
Each MPI process calculates exchange-correlation potential andenergy density using the structure D as follows:

In Set_XC_Grid.c

MN runs over grids
in the structure D.

Summary
The data structure of OpenMX is carefully designed so that the sizeof memory and MPI communication can be minimized.
The first step to construct the data structure is how atoms areallocated to each MPI process by using the modified bisection andinertia tensor projection methods.
The second step to construct the data structure is how four sorts ofgrid structure are constructed.
Each of calculations are performed by choosing one of the four gridstructures, and changing the grid structure requires the MPI
communication, of which pattern is determined beforehand intruncation().