Numerical Methods for Large-scale Quantum
Material Simulations
Hiroto Imachi
Tottori University
January 2017

Contents
1

Introduction

7

2

Background
2.1 Basics of parallel computing . . . . . . . . . . . . . . . . . . . .
2.1.1 Hardware . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.2 Software . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.3 Performance evaluation of parallel program . . . . . . . .
2.1.4 Trends of large-scale parallel computing . . . . . . . . . .
2.2 Basics of numerical linear algebraic routines . . . . . . . . . . . .
2.3 Basics of electronic state calculations . . . . . . . . . . . . . . .
2.3.1 Theoretical foundation of electronic state calculation . . .
2.3.2 Derivation of generalized eigenvalue problem . . . . . . .
2.3.3 Numerical aspect of large-scale calculations . . . . . . . .
2.4 Notes on physical and methodological concepts . . . . . . . . . .
2.4.1 Participation ratio . . . . . . . . . . . . . . . . . . . . .
2.4.2 Mobility . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4.3 Diﬀusion constant . . . . . . . . . . . . . . . . . . . . .
2.4.4 Atomic unit . . . . . . . . . . . . . . . . . . . . . . . .

11
11
11
12
12
13
14
15
15
17
18
18
18
21
22
23

3

Hybrid parallel eigenvalue solver
3.1 Concept of hybrid solver . . . . . . . . . . . . . . . . . . . . . .
3.2 Parallel eigenvalue solvers of ScaLAPACK, EigenExa and ELPA
3.3 Mathematical formulation . . . . . . . . . . . . . . . . . . . . .
3.4 Benchmark result . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.1 Result with the matrix size of M = 430, 080 . . . . . . . .
3.4.2 Benchmark with the matrix sizes of M=90,000, 22,500 .
3.4.3 Benchmark for a million dimensional matrix . . . . . . .
3.5 Discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.5.1 Preparation of initial distributed data . . . . . . . . . . .
3.5.2 Data conversion overhead . . . . . . . . . . . . . . . . .
3.5.3 Reducer analysis and future outlook . . . . . . . . . . . .

25
25
27
28
30
31
35
39
39
39
40
41

3

CONTENTS

4
3.6

Details of the code . . . . . . . . . . . . . . . . . . . . . . . . .
3.6.1 Execution flow of EigenKernel . . . . . . . . . . . . . . .
3.6.2 Quick start of EigenKernel . . . . . . . . . . . . . . . . .
3.6.3 Use with real applications . . . . . . . . . . . . . . . . .
3.6.4 ELSES matrix library . . . . . . . . . . . . . . . . . . .

43
43
44
45
46

4 Extreme scalable organic material simulations
47
4.1 Problem and previous works . . . . . . . . . . . . . . . . . . . . 47
4.1.1 Ultra-flexible device and organic materials . . . . . . . . 47
4.1.2 Electronic wavefunction in organic materials . . . . . . . 48
4.1.3 Concept of order-N method . . . . . . . . . . . . . . . . 49
4.1.4 Highly parallelizable mathematical structure . . . . . . . 51
4.2 Ground algorithm design . . . . . . . . . . . . . . . . . . . . . . 52
4.2.1 Generalized shifted linear equations . . . . . . . . . . . . 52
4.2.2 Calculation of projected physical quantity . . . . . . . . 53
4.3 Methodological details for large scale calculations . . . . . . . . 54
4.3.1 Quantum molecular dynamics simulation for organic polymers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.3.2 Details in parallelization . . . . . . . . . . . . . . . . . . 55
4.3.3 Sparsity of matrices . . . . . . . . . . . . . . . . . . . . . 55
4.3.4 MPI Communication . . . . . . . . . . . . . . . . . . . . 56
4.3.5 Use of Extensible Markup Language (XML) in File I/O . 56
4.3.6 Performance tuning on the K computer . . . . . . . . . . 57
4.4 Benchmarks and their analysis on the K computer . . . . . . . . 58
4.4.1 Architecture . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.4.2 Conditions of benchmarks . . . . . . . . . . . . . . . . . 59
4.4.3 Preparation of disordered structures of condensed polymers 60
4.4.4 Analysis on strong scaling and time-to-solution . . . . . . 61
4.4.5 Detailed analysis . . . . . . . . . . . . . . . . . . . . . . 63
4.5 Wavepacket dynamics simulation . . . . . . . . . . . . . . . . . 66
4.5.1 Concept of multi-time-scale quantum wavepacket simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
4.5.2 Basic mathematical formulation of generalized eigenvalue
equation . . . . . . . . . . . . . . . . . . . . . . . . . . 68
4.5.3 Multi-time-scale simulation (I) . . . . . . . . . . . . . . 69
4.5.4 Multi-time-scale simulation (II) . . . . . . . . . . . . . . 71
4.5.5 Simulation with modeled thermal atomic motion . . . . . 72
4.5.6 Calculation of diﬀusion coeﬃcient . . . . . . . . . . . . . 73
4.6 Organic material simulations . . . . . . . . . . . . . . . . . . . . 74
4.6.1 Quantum molecular dynamics simulation . . . . . . . . . 74

CONTENTS
4.6.2

4.7
4.8

5

Transport calculations with wavepacket dynamics simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Data scientific research for prescreening . . . . . . . . . . . . . .
108 -atom simulation . . . . . . . . . . . . . . . . . . . . . . . .
4.8.1 Network analysis of electronic wavefunctions . . . . . . .
4.8.2 Quantum wavepacket dynamics simulation on polymer
networks . . . . . . . . . . . . . . . . . . . . . . . . . .

5

Summary and future outlook

76
77
80
80
83
85

Chapter 1
Introduction
Nowadays, quantum material simulations, electronic state calculations, form a rigorous theoretical foundation of material science and engineering [1]. These simulations treat each electron as a quantum ‘wave’. Several simulation softwares, like
Gaussian [2] and VASP [3], are the de facto standards and are used widely amongexperimental and industrial researchers, as well as theoretical ones. Small system,
such as benzene molecule (C6 H6 ), can be calculated even by a laptop computer.
These simulations, however, are not applicable to large-scale systems, for example, one-hundred-nanometer-scale systems, owing to heavy computational costs
and a strong need appears for numerical methods in such large-system calculations. On the other hand, the current and next-generation supercomputers work as
massively parallel computing and the numerical methods for these machines are
crucial issues in computational science.
The present thesis is devoted to numerical methods for large-scale quantum
material simulation, so as to satisfy the need, in particular, for industrial application [4, 5, 6]. We focused on fundamental numerical problems in quantum
material simulations and these problems appear the Schrödinger equations in thetime-independent or time-dependent forms. The thesis consists of three main issues; hybrid parallel generalized eigenvalue problem solver, quantum wavepacket
dynamics simulation solver, and extreme parallelism on the full system of the
K computer for an order-N electronic structure calculation code (N is the number of atoms in a system). The three parts are relevant as follows; an orderN electronic structure calculation code ELSES(=Extra Large Scale Electronic
Structure calculations) [7] enables one-hundred-nm-scale quantum molecular dynamics (MD) simulations. Electronic structure can be calculated via generalized
shifted linear equations or generalized eigenvalue equations. The order-N calculation is achieved when generalized shifted linear equations are solved. Generalized eigenvalue problem solver gives more reliable numerical result with larger(O(N 3 )) computational costs. The order-N solver and the hybrid parallel general7

8

CHAPTER 1. INTRODUCTION

ized eigenvalue problem solver are complement to each other. They are massively
parallel numerical solvers and show high scalability up to the full system of the K
computer, one of the fastest supercomputers in the world. Based on the quantum
MD results, quantum wavepacket dynamics simulations are executed to calculatetransport properties.
At first the background of the thesis such as overview of parallel computing
and organic materials is explained. Organic materials play a crucial role amongnext-generation IoT (Internet of Things) products, such as display, battery andsensor, since they form flexible atomic structures and enable ultra-thin, light, flexible (wearable) devices with a low fabrication cost. Because disorder (or randomness) in atomic structures is important for properties of organic materials, largescale (100nanometer or 108 atom scale) quantum material simulation methods are
required. Eﬃcient parallel computation is a key technique for realization of suchlarge-scale quantum material simulations. Algorithms and implementations must
be carefully chosen depending on a required physical quantity and parallel computer architecture.
Hybrid parallel generalized eigenvalue problem solver Optimal hybrid numerical solvers are constructed for massively parallel generalized eigenvalue problem (GEP). The strong scaling benchmark was carried out on the K computer andother supercomputers for electronic structure calculation problems in the matrix
sizes of M = 104 − 106 with up to 105 processor cores. The procedure of GEP
is decomposed into the two subprocedures of the reducer from the GEP to thestandard eigenvalue problem (SEP) and the solver of SEP. A hybrid solver is constructed, when a routine is chosen for each subprocedure from the three parallel
solver libraries of ScaLAPACK, ELPA and EigenExa. The hybrid solvers withthe two newer libraries, ELPA and EigenExa, give better benchmark results thanthe conventional ScaLAPACK library. The detailed analysis on the results implies
that the reducer can be a bottleneck in next-generation (exa-scale) supercomputers, which indicates the guidance for future research of parallel generalized eigenvalue problem solvers.
Quantum wavepacket dynamics simulation solver This part focuses ontransport calculations for condensed organic polymers. One-hundred-nm-scale
electronic structure calculations were carried out by ELSES on the K supercomputer. The transport calculations were carried out as a theoretical extension for thequantum (hole) wavepacket dynamics simulation. The calculation is based on atime-dependent Schrödinger type equation i dtd ψ = HWP ψ for a hole wavefunction
ψ(t) with a modelled Hamiltonian HWP . Analysis of time evolution of ψ(t) gives
mobility, which is an important parameter for electronic device performance. The
method was applied to a single polymer chain and condensed polymers. The resultof mobility calculation is consistent to the experimental trend.
Order-N electronic structure calculation code A novel parallel linear-algebraic

9algorithm was introduced to electronic state calculations. The benchmark shows
an extreme strong scaling (75% in parallel eﬃciency) and a qualified time-tosolution (less than 102 sec in elapsed time) on the full system of the K computer.
Their mathematical foundation is generalized shifted linear equations ((zS −H)x =b), instead of conventional generalized eigenvalue equations. The foundation has
a highly parallelizable mathematical structure and applicable to many scientificareas. The simulation of organic polymer devices was carried out in academicindustrial collaboration. Using the electronic state calculations, a network analysis on connected polymer networks was carried out. Small networks of several
polymers that electronically connected are extracted. The quantum wavepacket
dynamics simulation was employed for transport calculations for the extracted
networks. The simulation and data analysis reveal that electronic waves propagate
on connected polymer networks and contribute the electrical current. The present
simulation method will give the insights of next-generation electronic devices andtheir fabrication process.

Chapter 2
Background2.1

Basics of parallel computing

Parallel computing is a form of computational methods which exploits concurrency in a program to execute eﬃciently. Concurrency is a character of a programthat its partial tasks can be executed simultaneously at overlapping time intervals.
Here eﬃciency means mainly shortening of the elapsed time for a program. In another context, eﬃciency also means reduction of required memory on a machine,which leads to enable larger computation on whole machines. Transforming a
program to assign processors for the simultaneous partial tasks so that the elapsed
time is shorten is called parallelization. Parallelism, chance of parallelization, appears in various levels of a program. Parallelism is called coarse-grain when eachof the partial tasks is relatively large and consumes long elapsed time. Conversely,
it is called fine-grain when they are relatively small. Generally, coarse-grain parallelism requires smaller additional computational costs on parallelization. Parallelization needs support from both sides of hardware and software.

2.1.1 Hardware
A parallel computer is a computer which has multiple processors and parallel
programs can run on it. Most of recent computers are kind of parallel computerbecause they ordinarily have multiple cores on a CPU chip. A supercomputeris a parallel computer which has very high computational capacity and mainly
used for scientific computing. In general, computational capacity of supercomputers is compared by FLOPS (value of how many floating-point operations canbe executed per second). For example, the K computer, one of the fastest supercomputers in the world, has theoretical peak performance of 11.28 Peta FLOPS.
Parallel computers are classified into two classes, shared memory machines or11

12

CHAPTER 2. BACKGROUND

distributed memory machines, by their memory architecture. On a shared memorymachine all the CPU cores are able to access to the single memory unit. On a
distributed memory machine, there are multiple memory units and a CPU coreis able to access to a part of them. One should use communication between thememory units to share data on a distributed memory machine. Although it iseasier to implement a parallel program on a shared memory machine than on a
distributed memory machine, there is a limitation on the number of computationalcores on a shared memory machine due to costs for preserving memory coherence.
Therefore almost all of modern supercomputers are distributed memory machinesthat connect many shared memory computer nodes by communication network.
As stated above, a parallel computer consists of processors, memory units, andnetwork. Among these, costs of memory access and network access are relatively
large in modern supercomputers. Therefore it is important to reduce memory andnetwork access as much as possible for implementation of an eﬃcient parallel
program.

2.1.2 Software
Programming models for parallel computation are classified by their memory
model. On the shared memory model which a single memory space are sharedby all processors, APIs (application programming interface) like pthreads (POSIX
threads) or OpenMP [8] are commonly used. On the other hand, on the distributed
memory model which each processor refers local memory space, required data areshared by an explicit send and receive. This style of sharing data is called message passing. De facto standard API for the message passing programming modelis MPI [9]. Almost all of modern supercomputers provide an implementation of
OpenMP and MPI. An example of famous MPI implementation is MPICH [10].
It should be noted that the classification of programming models is independent from the classification of memory architecture of hardware (shared memorymachine or distributed memory machine). One can run MPI parallel program ona shared memory machine. Also, one can run a shared memory parallel program
written in such as a PGAS (Partitioned Global Address Space) language on a distributed memory machine [11].

2.1.3 Performance evaluation of parallel program
Performance of a parallel program is often evaluated by how the elapsed time T
decreases with increase of the number of processors n. It is called scaling property. There are two measures of scaling property, weak scaling and strong scaling.
In weak scaling, the problem size (computational complexity) is increased in proportion to the number of processors. In strong scaling, the problem size is fixed2.1. BASICS OF PARALLEL COMPUTING

13

and only the number of processors is changed. Ideally, T is independent of n in
weak scaling, and T is inversely proportional to n in strong scaling. However,in reality, T is larger than them because there are costs for parallelization itself.
Generally, it is more diﬃcult to improve strong scaling property than weak one
because the fraction of costs for parallelization itself increases as n increases.
Strong scaling property can be predicted by Amdahl’s argument (or Amdahl’slaw) [12]. Suppose that there are a program which is divided into two parts andan amount of a computational resource s (for example, the number of processors).
The elapsed time of one part decreases in proportion to s. That of the other partdoes not change while the amount of the resource changes. p is the percentage ofthe elapsed time for the former part when s = 1. T s is the total elapsed time as schanges. Under the above condition, speedup ratio r(s) ≡ T 1 /T s is determined by
1r(s) =.
(2.1)(1 − p) + ps
In context of parallelization, it indicates that benefit from parallelization is small
when the fraction of the non-parallelized part 1 − p is large. Also the upper limitof the speed up ratio is determined by1
lim r(s) =
.
(2.2)
s→∞
1− p
For example, when p = 0.9, the maximum speedup is ten times as fast as thebaseline.

2.1.4 Trends of large-scale parallel computing
Past research of parallel computing had a tendency to increase solvable problem
size, in other words, a tendency to improve weak scaling property, not strong
scaling property. Because recent supercomputers have very high computational
capacity, the maximum solvable problem size for them is sometimes out of therange of interest from a viewpoint of real application. Therefore strong scaling
property is becoming more and more important.
In recent supercomputers, although performance increase of computational
units is rapid, that of data transfer functions such as memory access and communication is relatively slow. The separation between them will expand and existing
parallel programs may suﬀer from low scalability on next generation supercomputers.
Those problems on modern large-scale parallel computing cannot be solved
by eﬀort in a single field. Related fields of hardware architecture, algorithm, implementation and application must be simultaneously considered and adjusted foreach other. Such a style of design with collaboration of related fields is often
called co-design.

CHAPTER 2. BACKGROUND

14

2.2 Basics of numerical linear algebraic routines
Numerical linear algebraic solvers are fundamental components of scientific computing. Examples of typical problems in numerical linear algebra are linear equation and eigenvalue problem. Elementary methods for the problems such as calculating inverse matrix for linear equation or solving characteristic polynomialfor eigenvalue problem are often not suitable to calculation on computers due tolarge computational cost and numerical error. Instead of them, algorithms andimplementations suitable to computers has been extensively studied [13, 14, 15].
Especially, numerical linear algebraic solvers for large matrices have strong needs
among various applications with the current and next-generation supercomputers.
Numerical linear algebraic solvers are classified into two categories. One isdense-matrix solver and the other is sparse-matrix solver. In the dense-matrix
solvers, all the matrix elements are treated as non-zero values and a routine requires O(M 2 ) memory cost and O(M 3 ) operation cost, typically, where M is thematrix dimension. In the sparse-matrix solvers, on the other hand, only non-zeroelements are treated explicitly and a routine typically requires O(N) memory costand cheaper operation cost, where N is the number of non-zero elements. Sparsematrix solvers are again classified into two categories. One is sparse direct solver
and the other is sparse iterative solver. Sparse direct solvers sequentially transforma matrix into some standard form to compute a solution. Sparse iterative solvers
do not transform a matrix and use it only for matrix-vector multiplications to compute a solution. The thesis does not treat sparse direct solvers. Sparse matrices
emerge from various areas of scientific computing. For example, discretizationof diﬀerential equation often generates a sparse matrix, reflecting spatial localityof diﬀerential operators. In the thesis, the hybrid generalized eigenvalue solver
in Chap. 3 is classified into dense eigenvalue problem solvers, and the electronic
state calculation solver in Chap. 4 is classified into sparse linear equation solvers.
Design of dense-matrix solvers does not relatively depend on property of input matrices. Therefore, there are de facto standard libraries for dense numerical linear algebra calculation, such as BLAS 1 , LAPACK 2 and ScaLAPACK 3 .
BLAS [16] is an API set of basic linear algebra operations like summation andmultiplication of matrices and vectors. LAPACK [17] is a solver collection ofnumerical linear algebra problems like linear equation, eigenvalue decompositionand singular value decomposition. LAPACK is constructed based on BLAS. This
separation of implementation level realizes both of high performance and portability. ScaLAPACK [18, 19] is the distributed parallel version of LAPACK. These
libraries are available on almost all of modern supercomputers and highly optiBLAS = Basic Linear Algebra Subprograms
LAPACK = Linear Algebra PACKage
3
ScaLAPACK = Scalable Linear Algebra PACKage

12

2.3. BASICS OF ELECTRONIC STATE CALCULATIONS

15

mized for each architecture. But several routines give severe bottlenecks in thecomputational speed with current massively parallel architectures.
On the other hand, design of sparse-matrix solvers highly depends on property
of input matrices. Moreover, parallelization and optimization are generally morediﬃcult than dense-matrix solvers due to lack of regularity in calculation. Therefore, it is important to consider specific structures in a problem when designingsparse-matrix solvers. In sparse iterative solvers, the most time-consuming procedure is often matrix-vector multiplication. CRS (Compressed Row Storage) [15]is a common sparse matrix format. In CRS format, a matrix data is stored by
two integer arrays colind, rowptr and a floating-point value array values. values
stores the values of the non-zero elements of the matrix in row-wise order. colind
has the same length with values and contains the column indices of the corresponding non-zero elements. rowptr has the length of the number of rows in thematrix and contains indices in values and colind where the data of each row starts.
Pseudo-code of matrix-vector multiplication y = Ax in CRS format is as follows:
1: procedure matrix-vector multiplication in CRS format(A, x)2:
for i = 1 . . . m do3:
for k = A.rowptr[i] . . . A.rowptr[i + 1] − 1 do
4:
j = A.colind[k]5:y[i]+ = A.values[k] ∗ x[ j]6:
end for7:
end for
8: end procedure

2.3

Basics of electronic state calculations

This section is devoted to a brief review of the theoretical foundation of electronic
state calculation. The explanation focuses on the mathematical structure of thetheory and ignores several physical aspects, such as spin freedom of electron, forsimplicity.

2.3.1 Theoretical foundation of electronic state calculation
The fundamental theory of quantum material simulation or electronic state calculation stems from Schrödinger’s equation for many electron wavefuntionsΨ ≡ Ψ(r1 , r2 , ...., rNelec ),(2.3)

a complex scalar function, where Nelec is the number of electrons and r j is the
position of the j-th electron. It is very diﬃcult, however, to solve the equation

CHAPTER 2. BACKGROUND

16

with materials directly, because the computational cost increases exponentiallyas the function of Nelec . Walter Kohn, a Nobel Prize winner in Chemistry on
1998, called the situation ‘exponential wall’ [20]; If the continuum space of r j isdiscretized into p points, the number of elements for the discretized wavefunctionN
is M = pelec. In the case of Nelec = 36 electrons with p = 2, for example, thenumber M is given as
M = pNelec = 236 = 68, 719, 476, 736 ≈ 70G(2.4)

and the required memory size to store the wavefunction as complex double precision complex array is estimated to be
16B × M ≈ 1TB.

(2.5)

The above estimation indicates that the above calculation can not be carried out
with one hundred electrons.
Most of the present electronic state calculations are carried out with oneelectron wavefunction ϕ(r), instead of the many-electron wavefunction Ψ. See
textbooks for detail [1]. Eﬀective equations appear as eigenvalue equations incontinuum spaceĤeﬀ ϕ j = λ j ϕ j(2.6)

with the Hamiltonian operator Ĥ ofĤeﬀ ≡ −ℏ2∆ + Ueﬀ (r).
2me

(2.7)

Here me is the mass of electron and ℏ is Planck’s constant (me = 9.109534 ×
1031 kg, ℏ = 1.0545887 × 10−34 Js). The scalar function Ueﬀ (r) is the potentialfunction for electrons and diﬀers among materials. Since the potential function
includes the quantum mechanical interactions between electrons, Eq. (2.6) is anonlinear equation. The solution of Eq. (2.6) gives the eigenpairs of {λ j , ϕ j (r)}.
The eigenvalues of {λ j } are real and the index j is determined as λ1 ≤ λ2 ≤λ3 , .... The lowest Nelec eigenpairs are those for the electrons in material. The j-th
eigenvalue and eigenfunctions means the energy and state (wavefunction) of thej-th electron, respectively. The wavefunctions satisfy the orthogonality relation of∫ϕk (r)ϕl (r)dr = δkl ,(2.8)
unless the eigenvalues are degenerated (λk , λl ). The Nelec -th eigenvalue or eigenfunction is usually called, highest occupied (HO) energy or state, respectively.

2.3. BASICS OF ELECTRONIC STATE CALCULATIONS

17

The numbers of electrons Nelec in material and is usually proportional to the number of atoms N (Nelec ∝ N).
It is noteworthy that the one-electron wavefunction theory corresponds to an
approximation on the many-electron wavefuntion by the determinant form called
Slater determinant. The case of Nelec = 2, for example, is written as
1 ϕ1 (r1 ) ϕ2 (r1 )1
Ψ(r1 , r2 ) := √= √ (ϕ1 (r1 )ϕ2 (r2 ) − ϕ2 (r1 )ϕ1 (r2 )).
2 ϕ1 (r2 ) ϕ2 (r2 )2
The determinant form satisfies the mathematical relation ofΨ(r1 , r2 ) = −Ψ(r2 , r1 ),(2.9)(2.10)

known as a Fermion property of electrons.

2.3.2 Derivation of generalized eigenvalue problem
Here a generalized eigenvalue problem (GEP) is derived by discretization of real
space. The problem will be the central linear algebraic problem in this thesis.
When an electronic wavefunction ϕk (r) is expressed by the linear combination
with given basis functions of {χ j (r)} j=1,...M∑ϕk (r) =y jk χ j (r),(2.11)
ja generalized eigenvalue problem of
Hyk = λk S yk .
appears with M × M Hermitian matrices of H and S that are defined as∫
Hi j ≡
χ∗i (r)Ĥeﬀ χ j (r)dr∫
S ij ≡
χ∗i (r)χ j (r)dr.

(2.12)(2.13)(2.14)The matrix S is positive definite from the definition. The matrices H and S are
called Hamiltonian and the overlap matrices, respectively.
In the thesis, the real-space atomic-orbital representation is used and the matrix dimension is proportional to the number of atoms (M ∝ N). A basis function
χ j (r) is a localized function, called atomic orbital, of which center is located at
the position of one atom (nucleus). In addition, the basis functions are chosen as
real functions and the matrices of H and S are real symmetric. Since the basis
functions are normalized, all the diagonal elements of S are the unity (S ii = 1)
and all the oﬀ-diagonal elements satisfy the inequality of |S i j | < 1(i , j). The
orthogonality relation in continuum real space, Eq. (2.8), is transformed into
ykT S yl = δkl .

(2.15)

CHAPTER 2. BACKGROUND

18

2.3.3 Numerical aspect of large-scale calculations
As explained in the previous subsections, the GEP of Eq. (2.12) gives the mathematical foundation of electronic structure calculations or quantum mechanicalcalculations of materials, in which an electron is treated as a quantum mechanical‘wave’. Fig. 2.1(a) shows an example of the wavefunction. The number of therequired eigenvalues is, at least, on the order of the number of the electrons or theatoms in calculated materials.
When one think about large scale calculations, however, a potential diﬃcultyin parallelism is the orthogonality relation between eigenvectors. Since the eigenvectors should satisfy an orthogonality relation of Eq. (2.15), an explicit orthogonalization procedure is required. An orthogonalization procedure requires an operation cost of (O(N 3 )) and can be a bottleneck in parallelism. Since the orthogonality is inherent in quantum mechanics, the above potential diﬃculty appears
commonly among the large-scale electronic state calculations.
Here an order-N electronic state calculation code ELSES [7, 21] is explained
briefly, since the present research is a theoretical extension of the methods developed previously for ELSES. The simulations use novel ‘order-N’ linear-algebraic
methods in which the computational cost is ‘order-N’ (O(N)) or is proportional to
the number of atoms N. Their mathematical foundation is sparse-matrix (Krylovsubspace) solvers. The detailed theories are explained later in this thesis. Eﬃcientmassively parallel computation is found in Fig. 2.1, a strong scaling benchmark
on the K computer [22, 23] with one hundred million atoms or one-hundrednanometer scale materials. The simulated materials are a nano-composite carbon
solid with N = 103, 219, 200 or M = 412, 876, 800 [22] and an amorphous-likeconjugated polymer with N = 102, 238, 848 or M = 230, 776, 128 [23]. All the
calculations in this thesis were realized with first-principles-based modeled (tightbinding-form) theory. The details can be found in Ref. [21] and will be explainedlater in this thesis.

2.4 Notes on physical and methodological concepts
2.4.1 Participation ratio
Participation ratio (PR) was introduced in theoretical condensed matter physics,
mainly for the theory of localized electron wavefunctions [24, 25, 26, 27].
Originally the PR is defined for a normalized wavefunction ϕ(r) as
(∫
PR

(org)

(ϕ) ≡

)−1
|ϕ(r)| dr4.

(2.16)

2.4. NOTES ON PHYSICAL AND METHODOLOGICAL CONCEPTS

(a)

19

Elapse time (sec)

(b)

R

R

The number of processor cores

Figure 2.1: (a) The upper panel is a π-type electronic wavefunction in an
amorphous-like conjugated polymer (poly-((9,9) dioctyl fluorine)). The lower
panel shows the atomic structure (R≡C8 H17 ) [22]. (b) Strong scaling plot by
ELSES for one-hundred-million-atoms calculations on the K computer. [22, 23]
The calculated materials are a nano-composite carbon solid (the upper line) andthe amorphous-like conjugated polymer (the lower line). The number of used
processor nodes is from P = 4,096 to 82,944 (full nodes of the K computer).
Suppose D is a closed area whose volume is Ω and ϕ is uniformly distributed in
D as

1


 √Ω (r ∈ D)
ϕ(r) = (2.17)
0 (otherwise).
Then the PR of ϕ is

(

PR

(org)

1
(ϕ) =
Ω2

)−1

∫
dr

= Ω.

(2.18)

D

Namely, the PR indicates how the wavefunction spreads in space.
An analogue of the PR for a finite dimensional vector is defined. In context ofthe thesis, it is defined for an eigenvector y of the generalized eigenvalue equation
Eq. (2.12) as()2
yT S y.
PR(y) ≡ ∑(2.19)
4
j |y j |Note that the definition is invariant under scaling. The PR indicates a measure of
the number of non-zero elements, namely, how broadly the elements exist on theCHAPTER 2. BACKGROUND

20

indices. For example, the case of S = I and
(
)T
1 1 1
y ≡ √ , √ , √ , 0, 0, ..., 0 ,3 3 3(2.20)

gives the PR as−1  (
−1
)4 




∑




1



4
PR(y) = |y j | 
=
3 √ 



 = 3.

3 (2.21)

j

In cases of S , I, the value of PR can be smaller than one (PR(y) < 1) and
the PR does not indicate the measure of the number of non-zero elements. In
calculation of PRs for the eigenvectors of a generalized eigenvalue problem, suchcases appear in high eigenvalue regions. For example, the case of
)
)
(
(
1 s
0 −t
, S =
H=
(2.22)
s 1−t 0
with t > 0 and 0 < s < 1. The eigenpairs are obtained as
( )
−t
1
1
λ1 =
, y1 = √1+s2(1 + s) 1
and
t
λ2 =
,
1−s

1
y2 = √
2(1 − s)

(

)
1
.
−1

The PR of the eigenvectors are given by
 (
−1
)4 
{
}−1




1
1


PR(y1 ) = == 2(1 + s)2
2 √



2


2(1
+
s)
2(1 + s)(2.23)(2.24)(2.25)

and
 (
−1
)4 
{
}−1




1
1


PR(y2 ) = == 2(1 − s)2 .
2 √



2


2(1
−
s)
2(1 − s)(2.26)√The latter value can be smaller than one (PR(y2 ) < 1), when s > 1 − 1/ 2.
However, only low eigenvalue regions are interesting in most cases and thevalue of PR in these regions usually larger than one and indicate the measure of
the number of non-zero elements.
It is noteworthy that the unit vector of y ≡ (1, 0)T satisfies the normalization
condition of y T S y = 1 and gives the PR of PR(y) = 1.

2.4. NOTES ON PHYSICAL AND METHODOLOGICAL CONCEPTS

21

2.4.2 Mobility
In this section, mobility of a carrier in semiconductors is derived by Einstein’s
relation as in many textbooks, for example, Ref. [28]. Mobility is an important
physical quantity of semiconductor and means the ability of a carrier to movethrough a material. Since a material with higher mobility value can be a better
device material, the material exploration for higher mobility is crucial for material
design.
Here the derivation is explained in the one-dimensional case for a carrier, holeor exited electron. The theoretical extension to two- and three-dimensional cases
is straightforward. The mass, charge, position of the carrier is denoted as m, q andx, respectively. The sign of q is denoted as s ≡ sign(q) (s = +1 for q > 0 ands = −1 for q < 0). The velocity is denoted as v ≡ dx/dt. The potential for the
carrier is denoted as V(x). The electric field is given as
E(x) ≡ −dV.
dx

The density of carriers, denoted as ρ(x), is given as()
sqV(x)
ρ(x) = ρ0 exp −kB T(2.27)(2.28)

by the Boltzmann distribution with the temperature T and Boltzmann’s constantkB .
The drift current is derived from Newton’s equation of motionmdvm= qE − v.
dtτ

(2.29)The second term of the right hand side is a modeled friction force with the time
constant of τ(> 0) called mean free time. A stationary state results in dv/dt = 0
and
v=

qτ
E.
m

(2.30)The drift current is defined as
jdrift ≡ ρqv =

ρq2 τ
E.
m

(2.31)

If the mobility µ is defined as
µ≡s

qτ
,
m

(2.32)

CHAPTER 2. BACKGROUND

22
one obtains
jdrift = sρqµE.

(2.33)

Hereafter, Eq. (2.33) is used as the definition of the mobility µ.
The diﬀusion current, on the other hand, is defined as
jdiﬀ ≡ −sqDdρdx

(2.34)

with the diﬀusion coeﬃcient D(> 0). The diﬀusion current appears from theparticle flow from a higher density region into a lower density region.
The total current is given as
jtot ≡ jdrift + jdiﬀ
= sqρµE − sqD

dρ
dx

dρ dV
dV
− sqD
dx
dV}dx
{
dV
dρ
= −sq
ρµ + Ddx
dV= −sqρµ(2.35)

If no voltage is imposed, the total current is zero ( jtot = 0) and one obtains
1 dρ
ρ dV
)
d (
−D
log ρ
dV(
)
d sqV
D
dV kB T
Dsq
kB T
D|q|
.
kB T

µ = −D====(2.36)

Eq. (2.36) is Einstein’s relation.

2.4.3 Diﬀusion constant
In diﬀusion phenomena, a characteristic diﬀusion speed called diﬀusion constantD, which has dimension of L2 T −1 , can be obtained. In the one dimensional diﬀusion equation∂f∂2 f= D 2,∂t∂x(2.37)

2.4. NOTES ON PHYSICAL AND METHODOLOGICAL CONCEPTS
the solution of the initial value problem from f (x, 0) = δ(x) is()
1x2
f (x, t) = √
exp −
,
4Dt
2 πDt

23

(2.38)

which is the normal distribution with variance of σ2 = 2Dt. Inversely the diﬀusion
constant D is determined from the solution f (x, t) as⟨x2 (t)⟩
D=
2t
where

(2.39)∫⟨x (t)⟩ ≡2

x2 f (x, t)dx.

(2.40)

2.4.4 Atomic unit
The atomic unit (a.u.) is widely employed in electronic state theory papers and isused throughout this thesis, except where indicated. The atomic unit is defined sothatme = 1,ℏ = 1,

e2
= 1,
4πϵ0

(2.41)

where me , e, ℏ and ϵ0 are the mass of electron, the charge of electron (e > 0),
Planck’s constant, and the permittivity of vacuum, respectively. The length atomic
unit
(
)−1
ℏ2 e2
aB ≡
≈ 0.52918 × 10−10 m(2.42)
me 4πϵ0is called Bohr radius. The energy atomic unit
(
)2
me e2
EH ≡ 2≈ 4.3598 × 10−18 Jℏ 4πϵ0(2.43)

is called Hartree unit. The time atomic unit is given by
T au ≡

1ℏ≈ 2.4189 × 10−17 s ≈ fs.
EH
40(2.44)

Several energy conversion relations among the atomic unit, electron volt (eV) and
Kelvin (K) are noted below;
1 a.u. = 27.211 eV,

1 eV = 11, 604 K.

(2.45)

Chapter 3
Hybrid solver for massively parallel
eigenvalue computation
3.1

Concept of hybrid solver

The present chapter presents a hybrid solver collection, EigenKernel [29], formassively parallel eigenvalue computation, as a foundation of large-scale electronic state calculations. Nowadays ScaLAPACK[18, 19] 1 is the de facto standard solver library for parallel computations but several routines give severe bottlenecks in the computational speed with current massively parallel architectures.
Novel solver libraries were proposed so as to overcome the bottlenecks. Since
the performance of numerical routines varies significantly with problems and architectures, the best performance is achieved, when one constructs an optimal‘hybrid’ among the libraries.
The concept of hybrid solver is illustrated in Fig. 3.1. It is a numerical middleware and has a unique data interface to real applications. One can choose theoptimal workflow for each problem without any programming eﬀort.
The present thesis focuses on dense-matrix solvers for generalized eigenvalue
problems (GEPs) in the form of
Ayk = λk Byk(3.1)

with the given M×M real-symmetric matrices of A and B. The matrix B is positive
definite. The eigenvalues {λk } and the eigenvectors {yk } will be calculated. The
computational cost is O(M 3 ) or is proportional to M 3 . The present hybrid solvers
are constructed among ScaLAPACK and the two newer libraries of ELPA [30, 31,
32] 2 , and EigenExa [33, 34, 35, 36]. The ELPA and EigenExa libraries are written
1
2

ScaLAPACK = Scalable Linear Algebra PACKageELPA = Eigenvalue soLvers for Petascale Applications

25

26

CHAPTER 3. HYBRID PARALLEL EIGENVALUE SOLVER

Figure 3.1: Concept of hybrid solver; Structure of the program code (a) withoutand (b) with hybrid solver or numerical middleware.
in Fortran and appeared in 2000’s for eﬃcient massively parallel computations.
The present dense-matrix solvers are complementary methods to the order-N
calculations, because the order-N calculation gives approximate solutions, while
the dense-matrix solvers give numerically exact ones with a heavier (O(M 3 )) computational cost. The use of the two methods will lead us to fruitful researches.
The exact solutions are important, for example, when the system has many nearly
degenerated eigen pairs and one would like to distinguish them. The exact solutions are important also as reference data for the development of fine approximate
solvers.
The matrices of A and B in the present benchmark appear on ‘ELSES Matrix Library’. [37] See Sec. 3.6.4 for details of ELSES Matrix Library. The
benchmark was carried out with the data files of ‘NCCS430080’, ‘VCNT22500’‘VCNT90000’ and ‘VCNT1008000’ for the matrix sizes of M=22,500, M=90,000,
M=430,080, M=1,008,000, respectively. A large matrix data (> 0.5GB) is uploaded as a set of split files for user’s convenience.
The physical origin of the matrices is explained briefly. The files in the present
benchmark are carbon materials within modeled tight-binding-form theories basedon ab initio calculations. The matrix of ‘NCCS430080’ appears in the material
research on a nano-composite carbon solid (NCCS) [38]. An sp-orbital form [39]is used and the system contains N = M/4 = 107, 520 atoms. The other files are
generated for thermally vibrated single-wall carbon nanotubes (VCNTs) within asupercell. An spd-orbital form [40] is used and each system contains N = M/9
atoms. The VCNT systems were prepared, so as to generate matrices systematically in diﬀerent size with similar eigenvalue distributions. These matrices were3.2. PARALLEL EIGENVALUE SOLVERS OF SCALAPACK, EIGENEXA AND ELPA 27

Figure 3.2: Workflow of the hybrid GEP solver.
used for the investigation on π-electron materials with the present dense-matrix
solver and the order-N solver. 3

3.2

Parallel eigenvalue solvers of ScaLAPACK, EigenExa
and ELPA

A hybrid solver is constructed, when a routine is chosen for each subprocedurefrom ScaLAPACK, EigenExa and ELPA. The code was developed as a general
middleware that can be connected not only to ELSES but also to any real application software, as in Fig. 3.1. A mini-application was also developed and usedin the present benchmark. In the benchmark, ScaLAPACK was used as a built-in
library on each machine. EigenExa in the version 2.2a 4 and ELPA in the version
2014.06.001 were used. ELPA and EigenExa call some ScaLAPACK routines.
3

The present matrices are sparse, which does not lose the generality of the benchmark, since
the cost of the dense matrix solver is not dependent on the number of non-zero elements of thematrix.
4
The present EigenExa package does not include the GEP solver. The GEP solver routine for
EigenExa in the present thesis is that of the version 2.2b of KMATH EIGEN GEV [41] that shares
the SEP solver routine with the EigenExa package.

28

CHAPTER 3. HYBRID PARALLEL EIGENVALUE SOLVER

3.3 Mathematical formulation
The GEP of Eq. (3.1) can be written in a matrix form of
AY = BYΛ,(3.2)

where the matrix Λ ≡ diag(λ1 , λ2 , . . . ) is diagonal and the matrix Y ≡ (y1 y2 · · · )satisfies Y T BY = I. In the solvers, the GEP of Eq. (3.1) is reduced to a standard
eigenvalue problem (SEP) of
A′ Z = ZΛ,(3.3)

where the reduced matrix A′ is real symmetric [14] and the matrix of Z ≡ (z1 z2 · · · )contain eigenvectors of A′ . The reduction procedure can be achieved, when the
Cholesky factorization of B gives the Cholesky factor U as an upper triangle matrix:
B = U T U.

(3.4)

The reduced matrix A′ is defined by
A′ = U −T AU −1 .

(3.5)

The eigenvectors of the GEP, written as Y ≡ (y1 y2 · · · ), are calculated from thoseof the SEP by
Y = U −1 Z.

(3.6)

This procedure is usually called backward transformation.
The GEP solver is decomposed into the two subprocedures of (a) the solver of
the SEP in Eq. (3.3) and (b) the reduction from the GEP to the SEP ((A, B) ⇒ A′ )and the backward transformation (Z ⇒ Y). The subprocedures (a) and (b) arecalled ‘SEP solver’ and ‘reducer’, respectively, and require O(M 3 ) operations.
Figure 3.2 summarizes the workflows of the possible hybrid solvers. A hybrid
solver is constructed, when one choose the routines for (a) the SEP solver and (b)the reducer, respectively.
For (a) the SEP solver, five routines are found in the base libraries; One routine is a ScaLAPACK routine (routine name in the code : ‘pdsyevd’) that uses theconventional tridiagonalization algorithm. [42] The ELPA or EigenExa library
contains a SEP solver routine based on the tridiagonalization algorithm. The routine in ELPA is called ‘ELPA1’ (routine name in the code : ‘solve evp real’) inthis thesis, as in the original paper [31], and the one in EigenExa called ‘Eigen s’or ‘EIGS’ (routine name in the code : ‘eigen s’). ELPA and EigenExa also contain the novel SEP solvers based on the narrow-band reduction algorithms without

3.3. MATHEMATICAL FORMULATION

29

the conventional tridiagonalization procedure. The solvers are called ‘ELPA2’(routine name in the code : ‘solve evp real 2stage’) for the ELPA routine and‘Eigen sx’ or ‘EIGX’ (routine name in the code : ‘eigen sx’) for the EigenExa
routine in this thesis. See the papers [34, 32] for details.
For (b) the reducer, three routines are found in the base libraries and are called
ScaLAPACK style, ELPA style, and EigenExa style reducers in this thesis. In the
ScaLAPACK style, the Cholesky factorization, Eq. (3.4) is carried out and thenthe reduced matrix A′ , defined in Eq. (3.5), is generated by a recursive algorithm(routine name ‘pdsygst’) without explicit calculation of U −1 nor U −T . Details of
the recursive algorithm are explained, for example in Ref. [43]. In the ELPA style,
the Cholesky factorization (routine name: ‘cholesky real’) is carried out, as in the
ScaLAPACK style, and the reduced matrix A′ is generated by the explicit calculation of the inverse (triangular) matrix R ≡ U −1 (routine names : ‘invert trm real’)and the explicit successive matrix multiplication of A′ = (RT A)R (routine names:‘mult at b real’) [31] 5 . In the EigenExa style, the Cholesky factorization is not
used. Instead, the SEP for the matrix B
BW = WD,(3.7)

is solved by the SEP solver (Eigen sx), with the diagonal matrix of D ≡ diag(d1 , d2 , ...)and the unitary matrix of W ≡ (w1 w2 ....). A reduced SEP in the form of Eq. (3.3)
is obtained by
A′ = (D−1/2 W T )A(WD−1/2 )Y = WD−1/2 Z,

(3.8)
(3.9)

because of Z = D1/2 W T Y and W −T = W. Equation (3.9) is solved by the SEP
solver (Eigen sx).
Though the SEP solver of Eq. (3.3) requires a larger operation cost than the
Cholesky factorization (See Fig.1 of Ref. [44], for example), the elapsed time can
not be estimated only from the operation costs among the modern supercomputers.
The benchmark of the hybrid GEP solvers was carried out for the eight workflows listed in Table 3.1. In general, a potential issue is the possible overheadof the data conversion process between libraries. This issue will be discussed inSec. 3.5.2.
5The benchmark was carried out in an ELPA style reduction algorithm. The ScaLAPACK
routine of ‘pdtrmm’ is used for the multiplication of the triangular matrix R from right, while asample code in the ELPA package uses the ELPA routine (‘mult at b real’). The diﬀerence isignored, since the elapsed time of the above procedure is not dominant.

CHAPTER 3. HYBRID PARALLEL EIGENVALUE SOLVER

30

Table 3.1: List of the workflows in the benchmark. The routine names for the SEP
solver and the reducer are shown for each workflow. Abbreviations are shownwithin parentheses.
Workflow
SEP solver
ReducerA
ScaLAPACK (SCLA) ScaLAPACK (SCLA)
B
Eigen sx (EIGX)
ScaLAPACK (SCLA)C
ScaLAPACK (SCLA)
ELPA
ELPA2
ELPA
D
E
ELPA1
ELPA
F
Eigen s (EIGS)ELPAGEigen sx (EIGX)
ELPA
H
Eigen sx (EIGX)
Eigen sx (EIGX)

3.4 Benchmark result
Strong scaling benchmarks are investigated for the hybrid solvers. The elapsed
times were measured for (i) the full eigenpair calculation (T full ) and (ii) the ‘eigenvalueonly’ calculation (T evo ). In the latter case, the elapsed time is ignored for the calculation of the eigenvectors. The two types of calculations are important among
electronic structure calculations. [31] The present benchmark ignores small elapsed
times of the initial procedure for distributed data and the comments on them willappear in Sec. 3.5.1.
The benchmark was carried out on three supercomputers; the K computer at
Riken, Fujitsu FX10 and SGI Altix ICE 8400EX. The K computer has a single
SPARC 64 VIIIfx processor (2.0GHz, 8-core) on node. The FX10 is Oakleaf-FX
of the University of Tokyo. Fujitsu FX10 is the successor of the K computer andhas a single SPARC64 IXfx processor (1.848 GHz, 16-core) on each node. 6 SGI
Altix ICE 8400EX of Institute for Solid State Physics of the University of Tokyo
was also used. It is a cluster of Intel Xeon X5570 (2.93GHz, 8-core). The byteper-flop value (B/F) is B/F=0.5, 0.36 or 0.68, for the K computer, FX10 or SGI
Altix, respectively. The numbers of used processor nodes P are set to be squarenumbers (P = q2 ) except in Sec. 3.4.3, since the ELPA paper [31] reported thatthe choice of a (near-)square number for P can give better performance.
When the non-traditional SEP solver algorithm of ELPA is used on Altix, one
can choose an optimized low-level routine using SSE instructions(‘REAL ELPA KERNEL SSE’) and a generic routine
6

Additional options of the K computer and FX10 are explained; An MPI process shape on
the Tofu interconnect was not specified. The rank directory feature to alleviate I/O contention wasused.

3.4. BENCHMARK RESULT

31(‘REAL ELPA KERNEL GENERIC’). [31] The optimized code can run only on
the Intel-based architectures compatible to SSE instructions and was prepared so
as to accelerate the backtransformation subroutine. Among the results on Altix,
the ‘ELPA2’ solver and the workflow D on Altix are those with the optimized
routine, while the ‘ELPA2′ ’ solver and the workflow D′ are those with the genericroutine.

3.4.1 Result with the matrix size of M = 430, 080
The benchmark with the matrix size of M = 430, 080 was carried out for up to
P = 10,000 nodes on the K computer. The elapsed times for P=10,000 nodes
is shown in Table 3.2. The elapsed times for all the cases are shown in Fig. 3.3
for the (a) full (T full ) or (b) eigenvalue-only (T evo ) calculations. The decomposed
times are also shown in Fig. 3.3 (c) for the SEP solver (T SEP ) and the reducer(T RED ) (T full = T SEP + T RED ). Table 3.3 shows the decomposed time of the SEP
solvers for P=10,000. A SEP solver routine is decomposed into three subroutinesof (i) the tridiagonalization or narrow-band reduction (‘TRD/BAND’), (ii) the divide and conquer algorithms for the tridiagonal or narrow-band matrices (‘D&C’)
so as to compute the eigenvalues, and (iii) the backtransformation of eigenvectors(‘BACK’) so as to compute the eigenvectors of the GEP. One can observe several
features on the results; (I) In the full calculation benchmark (Fig. 3.3(a)), the best
data, the smallest elapsed time, appears in the workflow G for P=10,000. The
workflow G is the hybrid solver that uses the ‘Eigen sx’ SEP solver in EigenExa
and the ELPA style reducer, since these routines are the best among the SEP
solvers and the reducers, respectively, as shown in Fig. 3.3(c) and Table 3.3. In−1
Table 3.2, the speed (T full
) of the workflow G is approximately four times fasterthan that of the conventional workflow A (11,634 sec) / (2,734 sec) ≈ 4.3).
(II) Fig. 3.3 (c) shows that the ELPA style reducer gives significantly smaller
elapsed times than those of ScaLAPACK and those of EigenExa. The elapsed
time for P=10,000 is T RED = 1,261 sec with the ELPA style reducer and is T RED= 2,157 sec with the EigenExa reducer. The elapsed time with the EigenExa
reducer is governed by that of the SEP solver for Eq. (3.7) (T SEP = 1,473 sec inTable 3.3). (III) In the eigenvalue-only calculation (Fig. 3.3(b)), the best data, thesmallest elapsed time, appears in the workflow D for P=10,000. The workflow
D is the solver that uses the ‘ELPA2’ SEP solver and the ELPA style reducer
and the eigenvector calculation consumes a large elapsed time of T vec ; T vec ≡T full − T evo = (4,242 sec) - (2,227 sec) = (2,015 sec) in Table 3.2. The time T vec
is contributed mainly by the backward transformation subroutine (T BACK =1,892
sec) in Table 3.3, because the backward transformation subroutine in ELPA2 uses
a characteristic two-step algorithm (See Sec. 4.3 of Ref. [31]).

32

CHAPTER 3. HYBRID PARALLEL EIGENVALUE SOLVER

Table 3.2: Selected results of the benchmark. The elapsed time for the full (eigenpair) calculation (T full ) and that for the eigenvalue-only calculation (T evo ) with theworkflows. The recorded time is the best data among ones with diﬀerent numbersof the used nodes. The number of used nodes (P) for the best data is shown withinparentheses. The best data among the workflows are labelled by ‘[B]’. The saturated data are labelled by ‘[S]’. The workflow D′ on Altix is that without the SSE
optimized routine of the ‘ELPA2’ SEP solver. See the text for details.
Size M/Machine WF
1,000,080/FX10 G
430,080/K
A
B
C
D
E
F
G
H
90,000/K
A
B
C
D
E
F
G
H
90,000/FX10
A
B
C
D
E
F
G
H
90,000/Altix
A
B
C
D
D′
E
F
G
H

T full (sec)
T evo (sec)39,919 (P = 4,800)
35,103 (P = 4,800)
11,634 (P = 10,000)
10,755 (P = 10,000)
8,953 (P = 10,000)
8,465 (P = 10,000)
5,415 (P = 10,000)
4,657 (P = 10,000)
4,242 (P = 10,000)
2,227 (P = 10,000)[B]
2,990 (P = 10,000)
2,457 (P = 10,000)
2,809 (P = 10,000)
2,416 (P = 10,000)
2,734 (P = 10,000)[B]
2,355 (P = 10,000)
3,595 (P = 10,000)
3,147 (P = 10,000)
590 (P = 4,096)
551 (P = 4,096)493 (P = 1,024)[S]449 (P = 1,024)[S]318 (P = 4,096)
298 (P = 4,096)
259 (P = 4,096)
190 (P = 4,096)[B]
229 (P = 4,096)[B]194 (P = 4,096)
233 (P = 4,096)
210 (P = 4,096)
258 (P = 4,096)
240 (P = 4,096)
253 (P=4,096)
236 (P=4,096)
1,248 (P = 1,369)
1,183 (P = 1,369)
691 (P = 1,024)[S]648 (P = 1,024)[S]835 (P = 1,369)
779 (P = 1,369)
339 (P = 1,369)
166 (P = 1,024)[B][S]
262 (P = 1,369)
233 (P = 1,024)[S]250 (P = 1,369)[B]
222 (P = 1,369)
314 (P = 1,024)[S]283 (P = 1,024)[S]484 (P=1,369)
456 (P=1,369)
1,985 (P = 256)
1,675 (P = 256)
1,883 (P = 256)
1,586 (P = 256)
1,538 (P = 256)
1,240 (P = 256)
1,621 (P = 256)594 (P = 256)
2,621 (P = 256)585 (P = 256)[B]
1,558 (P = 256)
1,287 (P = 256)
1,670 (P = 256)
1,392 (P = 256)
1,453 (P = 256)[B]
1,170 (P = 256)
2,612 (P=256)
2,261 (P=256)

3.4. BENCHMARK RESULT

33(Continuation of Table 3.2)
Size M/Machine WF
22,500/K
A
B
C
D
E
F
G
H
22,500/FX10
A
B
C
D
E
F
G
H
22,500/Altix
A
B
C
D
D′
E
F
G
H

T full (sec)
65.2 (P = 1,024)45.8 (P = 1,024)[S]41.7 (P = 2,025)
28.4 (P = 2,025)
28.3 (P = 2,025)[B]
28.8 (P = 1,024)[S]29.7 (P = 1,024)[S]39.3(P=1024)[S]126.2 (P = 256)71.3 (P = 256)[S]103.5 (P = 256)[S]30.5 (P = 529)[B]
34.3 (P = 256)[S]32.1 (P = 529)
45.3 (P = 529)
74.9(P=529)
51.4 (P = 256)70.0 (P = 256)45.6 (P = 256)41.8 (P = 256)59.6 (P = 256)32.3 (P = 256)[B]
48.5 (P = 256)57.2 (P = 256)71.2 (P=256)

T evo (sec)
59.6 (P = 256)43.2 (P = 1,024)[S]37.8 (P = 2,025)
22.6 (P = 1,024)
22.6 (P = 1,024)[B]
26.9 (P = 1,024)[S]27.8 (P = 1,024)[S]37.5(P=1024)[S]118.1 (P = 256)67.1 (P = 256)[S]96.3 (P = 256)[S]24.4 (P = 529)[B]
31.2 (P = 256)[S]29.4 (P = 529)
42.5 (P = 529)
72.2 (P=529)
42.1 (P = 256)50.7 (P = 256)35.5 (P = 256)22.3 (P = 256)[B]
21.8 (P = 256)[B]
26.7 (P = 256)37.3 (P = 256)39.6 (P = 256)64.1 (P=256)

Table 3.3: Decomposition of the elapsed time (sec) of the SEP solvers with
M=430,080 and P = 10, 000. See the text for the subroutine names of‘TRD/BAND’ , ‘D&C’ and ‘BACK’.
SEP solver TRD/BAND D&C
SCLA
3,055
465
ELPA2
966
141
1,129
138
ELPA1
EIGS
1,058
196
828
390
EIGX

BACK Total (T SEP )
633
4,152
1,892
2,999
400
1,667
265
1,521
255
1,473

34

CHAPTER 3. HYBRID PARALLEL EIGENVALUE SOLVER

Figure 3.3: Results with M=430,800 on the K computer. The elapsed times areplotted with the workflows for the (a) full (T full ) and (b) eigenvalue-only (T evo )calculations. (c) The decomposed times for the SEP solver (T SEP ) and for the
reducer (T RED ) are plotted. The routines for the reducers is labeled by ‘(RED)’.
Detailed decomposed times for subprocedures of the ELPA style reducer and theCholesky decomposition in the ScaLAPACK style reducer are also plotted in (c).
The ideal speedup in parallelism is drawn as a dashed gray line.

3.4. BENCHMARK RESULT

35

3.4.2 Benchmark with the matrix sizes of M=90,000, 22,500
The benchmark with the smaller matrix sizes of M = 90, 000 and 22,500 are also
investigated. The maximum number of used processor nodes is Pmax = 4,096,
1,039 and 256, on the K computer, FX10, and Altix, respectively. 7 Figures 3.4
and 3.5 show the data with M=90,000 and with M=22,500, respectively. The
decomposed times are shown in Fig. 3.6. Table 3.2 shows the best data for each
workflow among the diﬀerent numbers of used nodes.
The results will help general simulation researchers to choose the solver andthe number of used nodes, since the elapsed times in Table 3.2 are less than ahalf hour and such calculations are popular ‘regular class’ jobs among systematicinvestigations. 8Here, the results are discussed; (I) Table 3.2 shows that the smallest elapsed
time in the full calculation appears among the workflows with the ELPA style reducer (the workflows D, E, F, and G) and that in the eigenvalue-only calculation
appears with the workflow D. The above features are consistent to the resultsin the previous subsection. (II) Unlike the result in the previous subsection, thespeed up is sometimes saturated. An example is observed in Fig. 3.5 (a), in thefull calculation with M=22,500 on the K computer, because the elapsed time in
the workflow F gives a minimum as the function of P at P=1,024. The decomposition analysis of Fig. 3.6(b) indicates that the saturation occur both for the SEP
solver and the reducer, which implies that the improvement both on the SEP solver
and the reducer is desirable. The saturated cases are marked in Table 3.2 with thelabel of ‘[S]’. 9 (III) Finally, the SSE-optimized routine in the workflow D is compared with the generic routine in the workflow D′ in the case of M = 90, 000 on
Altix with P =256. The SSE-optimized routine is prepared only in the backward
transformation process. Since the process with the SSE-optimized routine or thegeneric one gives the elapsed time of T BACK = 929 sec or T BACK = 1,872 sec, respectively, the process is accelerated with the SSE-optimized routine by 1, 872 sec/ 929 sec ≈ 2.02. As shown in Table 3.2, the full calculation is accelerated with
the SSE-optimized routine by 2, 621 sec / 1, 621 sec ≈ 1.62.
It is observed on Altix that the ‘ELPA2’ and ‘ELPA2′ ’ SEP solver required nonblocking communication requests beyond the default limit number of NMPI MAX = 16, 384
and the job stopped with an MPI error message. Then the limit number was increased to
NMPI MAX = 1, 048, 576, the possible maximum of the machine by the environment variable‘MPI REQUEST MAX’ and the calculations were completed.
8
One should remember that supercomputers are usually shared by many researchers who run
many calculations in similar problem sizes successively and/or simultaneously.
9
No saturation is found on Altix, unlike on the K computer and FX10, partially because themaximum number of used nodes (Pmax = 256) is smaller.
7

36

CHAPTER 3. HYBRID PARALLEL EIGENVALUE SOLVER

Figure 3.4: Benchmark with the matrix size of M=90,000, (I) on the K computer
for the (a) full (eigenpair) and (b) eigenvalue-only calculation, (II) on FX10 for
the (c) full and (d) eigenvalue-only calculation, (III) on Altix for the (e) full and(f) eigenvalue-only calculation. The ideal speedup in parallelism is drawn as adashed gray line.

3.4. BENCHMARK RESULT

37

Figure 3.5: Benchmark with the matrix size of M=22,500, (I) on the K computer
for the (a) full (eigenpair) and (b) eigenvalue-only calculation, (II) on FX10 for
the (c) full and (d) eigenvalue-only calculation, (III) on Altix for the (e) full and(f) eigenvalue-only calculation. The ideal speedup in parallelism is drawn as adashed gray line.

38

CHAPTER 3. HYBRID PARALLEL EIGENVALUE SOLVER

Figure 3.6: Decomposition analysis of the elapsed time into those of the SEP
solver and the reducer (I) on the K computer with (a) M=90,000 and (b)M=22,500, (II) on FX10 with (c) M=90,000 and (d) M=22,500, (III) on Altix
with (c) M=90,000 and (d) M=22,500. The routines for the reducers is labeledby ‘(RED)’. The ‘ELPA2′ ’ SEP solver is that without the SSE optimized routine.
The ideal speedup in parallelism is drawn as a dashed gray line.

3.5. DISCUSSIONS

39

3.4.3 Benchmark for a million dimensional matrix
Finally, the benchmark for a million dimensional matrix is discussed. A press
release at 2013 [45] reported, as a world record, a benchmark of a million dimensional SEP carried out by EigenExa, in approximately one hour, on the full
(82,944) nodes of the K computer. An eigenvalue problem with a million dimensional matrix (M=106 ) seems to be the practical limitation of the present supercomputer, owing to the O(M 3 ) operation cost.
We calculated a million dimensional GEP at Dec. 2014 on the full (4,800)
nodes of Oakleaf-FX. Since computational resource was limited, only one calculation was carried out with the workflow G, because it gives the best data among
those with M = 430, 080 in Table 3.2. The calculation finished in approximately
a half day, as shown in Table 3.2 (T full = 39,919 sec and T evo = 35,103 sec). The
elapsed time of the reducer (T RED = T full − T SEP = 15,179 sec) is smaller than but
comparable to that of the SEP solver (T SEP = 24, 740). The benchmark provedthat the present code qualifies as software applicable to massively parallel computation with up to a million dimensional matrix.
The million dimensional GEP was also solved with the workflow G by the K
computer. The elapsed time was T elaps = 9, 939 sec with nnode = 41, 472 nodes
and T elaps = 5, 516 sec on the full system (with nnode = 82, 944 nodes).

3.5

Discussions

3.5.1 Preparation of initial distributed data
In the benchmark, the initial procedures including file reading are carried out forthe preparation of distributed data. Its elapsed time is always small and is ignoredin the previous section. 10 These procedures, however, may consume significant
elapsed times, when the present solver is used as a middleware with real applications. The discussions on such cases are beyond the present scope, since they
depend on the program structure of the real applications. Here, several comments
are added for real application developers; In general, the matrix data cost is, at
most, O(M 2 ) and the operation cost is O(M 3 ) in the dense-matrix solvers and one
should consider a balance between them. In the case of M = 430, 080, for example, the required memory size for all the matrix elements is 8 B × M 2 ≈ 1.5 TB,
which can not be stored on a node of the K computer. Therefore, the data shouldbe always distributed. In the real application (ELSES), the initial distributed data
10

In the case of the workflow G on the K computer with M=430,080 and P=10,000, for example, the elapsed time of the initial procedures is T ini =123sec and is much smaller than that of thetotal computation (T tot = 2, 734sec. See Table. 3.2). It is noteworthy that the present matrices aresparse, as explained in Sec. 3.1.

CHAPTER 3. HYBRID PARALLEL EIGENVALUE SOLVER

40

Table 3.4: The elapsed times for data conversion; ‘(b → 1)’, ‘(1 → b)’ and ‘T RED ’
are the times in seconds for, the conversion process from block cyclic into cyclic
distributions, the inverse process and the whole reducer procedure, respectively.
The saturated data of T RED are labelled by ‘[S]’. The ‘ratio’ is ((b → 1) + (1 →b)) / T RED .
Size M
1,008,000
430,080
90,00022,500Machine(P)FX10(4,800)
K(10,000)K(4,096)
FX10(1,369)
Altix(256)K(2,025)
FX10(529)
Altix(256)(b → 1)
51.4
13.4
6.89
1.89
2.01
0.571
0.328
0.120

(1 → b) T RED
51.7
8,208
6.48
1,261
0.797
124[S]
0.973 84.0[S]2.02394
0.610 11.3[S]
0.176
9.20
0.279
11.9

ratio[%]
1.26
1.58
6.21
3.41
1.02
10.4
5.48
3.35

is prepared, when only the required elements are generated and stored on eachnode.

3.5.2 Data conversion overhead
As explained in Sec. 3.3, several workflows require data conversion processes between distributed data formats, since ScaLAPACK and ELPA use block cyclic
distribution with a given block size nblock (> 1) and EigenExa uses cyclic distribution (nblock ≡ 1). In the present benchmark, the block size nblock in ScaLAPACK
and ELPA was set to be nblock = 128, a typical value. Consequently, the workflows
B, F, G require data conversion processes. In the present thesis, the elapsed timeof the conversion procedures is included in the reducer part (T red ).
Table 3.4 shows the elapsed time for the data conversion. The elapsed times
are shown in the cases with the maximum numbers of used nodes (P = Pmax )among the present benchmark. Two data conversion procedures are required. One
is the conversion from the block cyclic distribution into the cyclic distribution,
shown as ‘(b → 1)’ in Table 3.4 and the other is the inverse process shown as ‘(1→ b)’. The two procedures are carried out, commonly, by the ‘pdgemr2d’ routine
in ScaLAPACK.
Table 3.4 indicates that the overhead of the data conversion procedures is always small and is not the origin of the saturation. In general, the conversion
requires an O(M 2 ) operation cost, while the calculation in a dense-matrix solver
requires an O(M 3 ) operation cost. The fact implies the general eﬃciency of hybrid
solvers, at least, among dense-matrix solvers.

3.5. DISCUSSIONS

41

3.5.3 Reducer analysis and future outlook
The decomposition analysis of the ELPA-style reducer is focused on, since the
ELPA-style reducer is fastest among the three libraries. Figure 3.3 (c) shows thecase on the K computer with M=430,080. The elapsed times of the subproceduresof the ELPA-style reducer are plotted; ‘ELPA(R1 )’ is the Cholesky factorizationof Eq. (3.4), ‘ELPA(R2 )’ is the explicit calculation of the inversion R = U −1 of theCholesky factor U, ‘ELPA(R3 )’ and ‘ELPA(R4 )’ are the successive matrix multiplication of Eq. (3.5) and ‘ELPA(R5 )’ is the backward transformation of eigenvectors by matrix multiplication of Eq. (3.6). The elapsed times of the Cholesky
factorization in the ScaLAPACK style reducer is also plotted as ‘SCLA(R1 )’ asa reference data. The same decomposition analysis is carried out also for othercases, as shown in Fig. 3.7. One can observe that the Cholesky factorization of
the ELPA-style reducer does not scale and sometimes is slower than that of the
ScaLAPACK reducer. In particular, the saturation of the ELPA-style reducer iscaused by that of the Cholesky factorization in Fig. 3.7 (a)(b)(c).
The above observation implies that the reducer can be a serious bottleneck in
the next-generation (exa-scale) supercomputers, though not in the present benchmark. One possible strategy is the improvement on the Cholesky factorization for
better scalability and another is the development of a reducer without the Cholesky
factorization, as in the EigenExa-style reducer.
As a future outlook, the present code for the hybrid solvers is planned to beextended by introducing the solvers with diﬀerent mathematical foundations. A
candidate is the parallel block Jacobi solver [46, 47]. Since the solver is applicable
only to standard eigenvalue problems, the hybrid solver enables us to use thesolver in generalized eigenvalue problems.

42

CHAPTER 3. HYBRID PARALLEL EIGENVALUE SOLVER

Figure 3.7: Decomposition analysis of the elapsed time of subprocedures of the
ELPA style reducer and the Cholesky factorization in the ScaLAPACK style reducer (I) on the K computer with (a) M=90,000 and (b) M=22,500, (II) on FX10
with (c) M=90,000 and (d) M=22,500, (III) on Altix with (c) M=90,000 and (d)M=22,500. The ideal speedup in parallelism is drawn as a dashed gray line.

3.6. DETAILS OF THE CODE

3.6

43

Details of the code

3.6.1 Execution flow of EigenKernel
The code of the hybrid eigenvalue problem solver collection, named EigenKernel,
is released at GitHub [29]. EigenKernel can be built to a static library and astandalone executable. The executable form of EigenKernel is called eigbench.
execution flow of eigbench is as follows. The list mostly corresponds to the source
file src/main.f90.
1. Read commandline arguments
Read commandline arguments in the standard Unix style. eigbench cansolve both of SEP and GEP in almost the same manner. In the minimum,
eigbench takes one option (-s) to choose the solver name and one (in SEP)or two (in GEP) arguments of matrix file paths. Names of GEP solvers
have the prefix of ’general ’ like ’general scalapack’. In commandline arguments, output file paths of eigenvalues (-o), inversed participation ratios(see Sec. refSEC-PR-INTRO) of eigenvectors (-i), eigenvectors (-d) and
execution log file (-l) can be specified. The range of the indices of eigenvectors to be output can also be specified (-p <start>,<end>).
2. Read matrix files in MatrixMarket format and broadcasting
Input file format in eigbench is the MatrixMarket format [48]. An input fileis read only on the master node (0-th node in MPI) and broadcast to all thenodes. After broadcasting, matrix data are converted into a format for themain calculation. As explained in Sec. 3.5.2, EigenKernel uses block cyclicor (non-block) cyclic format in the main calculation.
3. Main solver routine call
Call the specified SEP or GEP solver. The data conversion process runs
when needed as explained in Sec. 3.5.2. When EigenKernel is used in thestatic library form, users can call only the main solver routines.
4. (Partially optional) output eigenvalues, inversed participation ratios and eigenvectors
Output calculated eigenvalues, inversed participation ratios of eigenvectors
and eigenvectors itself. Output of eigenvectors is optional process executed
only when the -p option is specified because in most cases the eigenvector
matrix is dense even if the input matrices are sparse and therefore requires
large disk space. Also, output of eigenvectors is parallelized by MPI. Output
eigenvector indices are partitioned into the compute nodes evenly in cyclic

44

CHAPTER 3. HYBRID PARALLEL EIGENVALUE SOLVER
manner. This parallelization is especially eﬀective in systems which have a
local file system for each compute node.
5. (Optional) compute residual norms and orthogonality of calculated eigenpairs
EigenKernel has two accuracy checking routines. One is to compute residual norms
rk ≡ ∥Ayk − λk Byk ∥ (k = 1, . . . , M),(3.10)and the other is to compute orthogonality of eigenvectorst ≡ ∥Y T BY − I∥F(3.11)

where ∥·∥F is the Frobenius norm.

3.6.2 Quick start of EigenKernel
In minimum configuration, EigenKernel can be built with an MPI Fortran compiler and ScaLAPACK (without ELPA and EigenExa). Settings depend on a machine environment should be written in Makefile.in. The sample Makefile.inin the repository supposes there are mpif90 over gfortran as a fortran compiler andlibgomp as an OpenMP library. Following commands will test solving a generalized eigenvalue problem with the matrix size of M = 30.
tar zxvf eigenkernel-*.tar.gz
cd eigenkernel-*
cp Makefile.in.gfortran.noext Makefile.inmakempirun -np 4 bin/eigbench -s general_scalapack \matrix/ELSES_MATRIX_BNZ30_A.mtx matrix/ELSES_MATRIX_BNZ30_B.mtx
After executing eigbench, there are output files named eigenvalues.dat, ipratios.dat,
log.json. eigenvalues.dat and ipratios.dat contain computed eigenvalues and inversed participation ratios of computed eigenvectors respectively. log.json contains execution information such as given commandline options in the JSON format and elapsed time for some routines. In the execution command -s <solver> isa mandatory option to specify the solver routine. The general scalapack solver is,of course, a pure ScaLAPACK solver. The last two arguments are paths to input
matrix files in the MatrixMarket format. Note that the input matrices must be real
symmetric and moreover the latter one must be positive definite. You can also
solve standard eigenvalue problems.

3.6. DETAILS OF THE CODE

45

mpirun -np 4 bin/eigbench -s scalapack \matrix/ELSES_MATRIX_VCNT400std_A.mtx
To utilize full functions of EigenKernel, it should be built linking ELPA andEigenExa. After installing them, edit Makefile.in and set LIBS variable properly to indicate the paths where .a and .mod are installed. Then EigenKernel with
ELPA and EigenExa should be rebuilt likeemacs Makefile.in # Edit $LIBS properly
make clean
make WITH_EIGENEXA=1 WITH_ELPA=1
If the rebuild is succeeded, truly hybrid solvers can be selected by the -s option.
mpirun -np 4 bin/eigbench -s general_elpa_eigensx \matrix/ELSES_MATRIX_VCNT900_A.mtx \matrix/ELSES_MATRIX_VCNT900_B.mtx

3.6.3 Use with real applications
EigenKernel can be built into a Fortran library as well as the standalone miniapplication eigbench. For example, one of optional solvers in ELSES uses EigenKernel. The library is automatically generated in default build setting. A minimal
example of calling EigenKernel library to solve a positive symmetric generalized
eigenvalue problem is as follows.
integer :: ierrtype(process) :: proc
type(eigenpairs_types_union) :: eigenpairs
call mpi_init(ierr)call setup_distribution(proc)
call reduce_generalized(dim, A, desc_A, B, desc_B)
call eigen_solver_scalapack_all(proc, desc_A, A, eigenpairs)
call recovery_generalized(dim, dim, B, desc_B, &eigenpairs%blacs%Vectors, eigenpairs%blacs%desc)
call mpi_finalize(ierr)
Call of the subroutines ‘mpi init’ and ‘mpi finalize’ is mandatory in MPI programs. ‘setup distribution’ determines the processor grid layout in ScaLAPACK.
‘reduce generalized’ and ‘recovery generalized’ execute the reduction and the inverse reduction between GEP and SEP in the ScaLAPACK style, respectively.
‘eigen solver scalapack all’ solves the reduced SEP. After calculation, eigenpairs
are stored in the variable ‘eigenpairs’. Variables ‘A’, ‘B’, ‘desc A’ and ‘desc B’

46

CHAPTER 3. HYBRID PARALLEL EIGENVALUE SOLVER

are the standard distributed matrices and their descriptors in ScaLAPACK. There
are also subroutines of solvers using external libraries like ELPA or/and EigenExa.
They have a unified interface and users can easily exchange solver routines.

3.6.4 ELSES matrix library
ELSES Matrix Library [37] is a collection of matrix data generated by ELSES
for quantum material simulations. They appears in standard or generalized eigenvalue problems. It contains from small (dimension M = 30) to very large (M =1, 008, 000) matrices. The matrices are sparse, real-symmetric or Hermitian. The
matrix data are recorded in the MatrixMarket format [48]. The atomic unit isused for energy unit. Some matrix packages contain a list of all eigenvalues andinversed participation ratios of eigenvectors.

Chapter 4
Extreme scalability of 108-atomsimulation for organic materials
4.1

Problem and previous works

4.1.1 Ultra-flexible device and organic materials
Ultra-flexible devices are based on organic materials and play a crucial role amongnext-generation IoT products, such as display, battery and sensor. Organic materials form flexible atomic structures and enable ultra-thin, light, flexible (wearable) devices with a low fabrication cost [49, 50, 51]. Alan J. Heeger, Alan G.
MacDiarmid and Hideki Shirakawa won the Nobel Prize in Chemistry in 2000
for their pioneering research on organic polymer devices in 1970’s. The atomic
structure of organic materials is disordered [52, 53, 54, 55] and the simulationof huge, 100-nanometer-scale, disordered systems are required. Nowadays manyelectronic state calculation codes are available but such large-scale calculations
are far beyond their computational limit.
The present thesis reports that a new electronic state calculation code ELSES[7] is based on novel linear algebraic algorithms suitable to massive parallelism[56, 57, 58, 21, 59, 60, 23, 4, 5], and was applied to the above challenging industrial problem. Since the algorithm has a highly parallelizable mathematicalstructure, the code shows excellent results both in the strong scaling and the timeto-solution. The code realized calculations with one-hundred-million (108 ) atoms
or 100-nm-scale systems, the largest system among electronic state calculations.
Since the algorithms are mathematical, the method is applicable to various materials, such as semiconductor and metal. The algorithm is applicable also to othercomputational science area. The code has been developed for large-scale and/or
fast electronic state calculations by academic-industrial collaboration. A recent
collaboration research is one for organic device materials (Refs. [21, 60, 5] and47

48CHAPTER 4. EXTREME SCALABLE ORGANIC MATERIAL SIMULATIONS
device(schematic)

simulation ofcondensed
polymers< 1,000 nm( 108 atoms or
105 polymers )close-up( ≒200nm )close-up:connected
polymers

close
-up

70 nm

C

C

n

2 nm

Figure 4.1: Overview of the present material research on organic polymer device;
The comparison of the length scales in real device and the present simulation ofpoly-(phenylene-ethynylene) (PPE).

the present work) with Sumitomo Chemical Co., Ltd, and another is one for battery materials [61, 62] with Toyota Motor Corporation. The other application
researches can be found in the reference lists of the above papers.
Figure 4.1 indicates the overview of the present material research. The lengthscales of real device and the present simulation are compared.

4.1.2 Electronic wavefunction in organic materials
The electric conductivity in organic devices and related materials like graphenestems from the characteristic electronic wavefunction called π wavefunction. A π
wavefunction lies, typically, on a benzene ring. Figure 4.2(a) shows a π wavefunc-

4.1. PROBLEM AND PREVIOUS WORKS

(a)

49

(b)

Figure 4.2: Organic materials and π electronic wavefunction. (a) Atomic structure of benzene and an example of π electronic wavefunction. (b) A schematic
figure of condensed organic polymers in a disordered structure. The vectors of
n1 , n2 , n3 , n4 are defined as ones perpendicular to the benzene ring planes.

tion in benzene (C6 H6 ). The charge of the π wavefunction of |ϕ(r)|2 is depicted
as isosurfaces. The wavefunction extends to the direction perpendicular to theplane of benzene, which is the origin of the strong anisotropy in electric conductivity. Figure 4.2(b) shows schematically the structure of condensed polymers inreal devices. Disorder appears both in intra- and inter-chain structures. Here thevectors of n1 , n2 , n3 , n4 are depicted as ones perpendicular to the benzene ringplanes. When neighboring two vectors of ni and n j are almost parallel (ni ∥ n j ),the π wavefunctions on their benzene rings connect easily with each other and theelectrical current can propagate between them.
In conclusion, quantum simulations with π-type wavefunctions are crucial forthe organic electronic devices, which is a challenging problem in computational
science.

4.1.3 Concept of order-N method‘Order-N’ method is a general name of the large-scale calculation methods inwhich the computational cost is O(N) or proportional to the number of atoms
N. The present code is one of them. Other order-N calculation codes were also
developed [63, 64, 65, 66].
A key concept for the order-N method was stated by Walter Kohn. His paperin 1996 shows that the above potential diﬃculty in electronic state calculation can

50CHAPTER 4. EXTREME SCALABLE ORGANIC MATERIAL SIMULATIONS
be avoided, when the theory is not based on an eigenvalue equation [67].
The theory [67] focuses on a physical quantity defined as∑⟨X⟩ ≡f (εk ) ykT Xyk ,

(4.1)

k

with a given sparse real-symmetric matrix X. Equation (4.1) is found in elementary textbooks of electronic state calculations. Here the function of f (ε) isa weight function, called Fermi function, and is defined as{ε − µ }−1f (ε) ≡ 1 + exp(
) .
τ

(4.2)

The weight function is a ‘smoothed’ step function with a smoothing parameterτ(> 0), because the Heaviside step function will appear in the limiting case ofτ → +0 ( f (ε) = 1(ε < µ) and f (ε) = 0(ε > µ)). The smoothing parameter τ hasa physical meaning of temperature of electrons. The parameter µ is the chemical
potential and the value should be determined so that the number of electrons is thecorrect one. The case in X = H, for example, gives the electronic energy∑⟨H⟩ ≡f (εk ) εk .
(4.3)
kOther quantities, such as the force on atom and the density of states (eigenvalue
histogram), are also described in the above form.
A quantity in Eq.(4.1) is transformed into the trace form of⟨X⟩ = Tr[ρX](4.4)

with the density matrixρ≡∑f (εk ) yk ykT .

(4.5)

k

The order-N property can be found, as follows; A density matrix element ρ ji
is not required when Xi j = 0, because the element ρ ji does not contribute to thephysical quantity∑
Tr[ρX] =ρ ji Xi j ,
(4.6)
i, j

even if its value is nonzero (ρ ji , 0). Consequently, the number of the required
density matrix elements ρ ji is O(N). The above principle is called ‘quantum locality’ or ‘nearsightedness principle’ [67].

4.1. PROBLEM AND PREVIOUS WORKS(a)51

iteration step
matrixgeneration

matrix solver

physicalquantities

conventional

G. E
Eig
Eigenvalue
igenvalu
lue eqn.

( wavefunction )physicalquantities

present

G. Shifted Linear eqns.

( Propagation (Greenʼs) fn. )(b) conventional method

present method・potential diﬃculty in parallelism(ex. orthogonalization procedures)・O(N3) operation costs

・highly parallelizable
mathematical structure (See (c))

・max. target size: 105 atoms

・O(N) operation costs
・max. target size: 108 atoms,
within the elapsed time < 102 s(c) highly parallelizable mathematical structure

Figure 4.3: Overview of the algorithm for large-scale electronic state calculationsin the extreme strong scaling; (a) The ground algorithm design. (b) Comparison between the conventional method with eigenvalue equation and the present
method with generalized shifted linear equations. (c) Schematic figure of thehighly parallelizable mathematical structure.

4.1.4 Highly parallelizable mathematical structure
The above formulation has a highly parallelizable mathematical structure and theoriginal problem is decomposed mathematically into parallel subproblems. The

52CHAPTER 4. EXTREME SCALABLE ORGANIC MATERIAL SIMULATIONS
trace in Eq. (4.4) can be decomposed as⟨X⟩ = Tr[ρX] =M∑

eTj ρXe j ,

(4.7)

j

with the j-th unit vector of e j ≡ (0, 0, 0, ....1 j , 0, 0, ..., 0)T . Here the quantity ofeTj ρXe j is called ‘projected physical quantity’, because the quantity is defined by
the projection onto the vector of e j . The essence of the parallelism is the factthat the projected physical quantity of eTj ρXe j is calculated almost independently
among diﬀerent indices of j.

4.2 Ground algorithm design
4.2.1 Generalized shifted linear equations
The ground algorithm design of the present method is summarized in Fig. 4.3(a).
The method is based not on the eigenvalue equation of Eq. (2.12) but on linear
equations in the form of(zS − H)x = b.

(4.8)

Here z is a (complex) energy value and the matrix of (zS − H) can be nonHermitian. The vector b is an input and the vector x is the solution vector. A set of
linear equations in the form of Eq. (4.8) with diﬀerent energy values (z = z1 , z2 , ...)
is called generalized shifted linear equations. The equations in the case of S = I
are called shifted linear equations.
The generalized shifted linear equations are solved on an iterative Krylovsubspace solver. Such solvers are sometimes called shifted Krylov-subspace solvers.
Krylov subspace is defined as the linear space ofKν (A; b) ≡ span[b, Ab, A2 b, ..., Aν−1 b],(4.9)

with a given vector b and a given square matrix A. Krylov subspace solver is thesolver in which the solution is calculated within a Krylov subspace. An example
is Conjugate Gradient method and the dimension of the subspace ν is the numberof iterations.
The use of Eq. (4.8) results in the Green’s (propagation) function formalism,
since the solution x of Eq. (4.8) is written formally as
x = Gb

(4.10)

4.2. GROUND ALGORITHM DESIGN

53

with the Green’s function G ≡ (zS − H)−1 . The Green’s function and the eigenvectors holds the relationship ofG(z) =

M∑
yk y Tkkz − εk.

The density matrix is also given by the Green’s function as
∫
−1 ∞
ρ=
f (ε) Im[G(ε + i0)] dε.
π −∞(4.11)(4.12)The use of the (generalized) shifted linear equations, instead of eigenvalue
equations, is a powerful strategy among large scale calculations. Such methodologies have been investigated, in particular, from 2000’s, partially because thestrategy is suitable to parallelism. Since the solver algorithms are mathematical,
they are applicable to many scientific areas, such as, QCD [68], large-scale electronic state calculation [56, 57, 58, 21], quantum many-body electron problem[69], nuclear shell model problem [70], first-principle electronic excitation problem [71], and first-principle transport calculation [72].

4.2.2 Calculation of projected physical quantity
In the present method, the formalism has a highly parallelizable mathematicalstructure and the projected physical quantity of eTj ρXe j , should be calculated as
explained in Sec. 4.1.4. Figure 4.3(b) shows the comparison between the conventional method with eigenvalue equation and the present method with generalized
shifted linear equations. The parallel procedure is illustrated in Fig. 4.3(c).
In the code, the projected physical quantity of eTj ρXe j is obtained from thegeneralized shifted linear equations of(zS − H)x( j) = e j(4.13)

within an iterative Krylov-subspace solver, called multiple Arnoldi solver [21]. In
short, Eq. (4.13) is solved within the subspace ofLν (e j ) ≡ Kν/2 (H; e j ) ⊕ Kν/2 (H; S −1 e j )(4.14)

with an even number of ν. The number ν is typically, ν = 30 − 300 and thecalculations in the present thesis was carried out with ν = 30 as in the previous
one [21]. The second term in the right hand side of Eq. (4.14) appears so as tosatisfy several conservation laws [21]. A reduced (small) ν×ν eigenvalue equation
is solved and the solution vector of Eq. (4.13) is given byx( j) := G( j) (z)e j(4.15)

54CHAPTER 4. EXTREME SCALABLE ORGANIC MATERIAL SIMULATIONS
with
G (z) ≡( j)

ν∑v ( j) v ( j)T
m

m

m
.
z − ε(mj)(4.16)

Here ε(mj) and vm( j) is an eigenvalue and eigenvector of the reduced equation (m =1, 2, ....ν). When the Green’s function of G in Eq. (4.12) is replaced by G( j) (z) in
Eq. (4.16), the projected physical quantity with the index of j is given by∫−1 ∞T
e j ρXe j :=f (ε) Im[eTj G( j) (ε+i0)Xe j ] dε
π −∞
ν
∑
=
f (ε(mj) )eTj vm( j) vm( j)T Xe j .
(4.17)m

Equation (4.17) will be exact, if the subspace dimension of ν increases to theoriginal matrix dimension (ν = M). As an additional technique in large-scalecalculations, the real-space projection method [21] was also used. In short, theKrylov subspace is generated by a Hamiltonian matrix projected in real spaceH ( j) ≡ P( j) HP( j) , where the projection matrix P( j) projects a function onto thespherical region whose center is located at the atomic position of the j-th atomic
basis function (χ j (r)). The radius of the spherical region is determined with an
input integer parameter κ, so that the region contains κ atoms or more. The same
technique is used also for the overlap matrix S . The value of κ is set to κ = 100 in
the present thesis as in the previous one [21]. As results, numerical problems in
the form of Eq. (4.13) are solved with the matrix size of, typically, M ′ = 200−400
among the simulations of the present thesis.

4.3 Methodological details for large scale calculations
4.3.1 Quantum molecular dynamics simulation for organic polymers
The present code is written in Fortran 90 with the MPI/OpenMP hybrid parallelization and the main purpose is quantum molecular dynamics simulation. Herean electron is treated as a quantum mechanical wave in the Green’s (propagation)function formalism, while an atom (a nucleus) is treated as a classical particle in
Newtonian equation of motionMI

d 2 RI= FI ,
dt2(4.18)

4.3. METHODOLOGICAL DETAILS FOR LARGE SCALE CALCULATIONS 55
where MI and RI are the mass and the position of the I-th atom and FI is theforce on the I-th atom. Other variables, such as the electronic charge on each atom{qI }I , can be calculated also in the dynamical simulation. As technical details, thequantum molecular dynamics simulations are realized with first-principles-basedmodeled (tight-binding-form) theory. See Refs. [21, 5] and the references thereinfor details. In the code, these variables that consume O(N) memory cost are stored
redundantly among nodes. All the matrix elements of H and S are generated fromthese variables on each node without any inter-node communications. Several
matrix elements are generated redundantly among nodes. The force and charge
on each atom is calculated in the trace form of Eq. (4.4)(See Ref. [21] for details).
The time evolution was realized numerically by the velocity-Verlet algorithm. The
finite-temperature simulations were realized by a thermostat method.

4.3.2 Details in parallelization
The details in parallelization are explained here, so as to clarify the condition of aneﬃcient parallelism. In the present parallel computation, the parallel subproblems
for calculating the projected physical quantity of eTj ρXe j should be carried out,
as in Fig. 4.3(c). In the code, the basis index j is treated as a composite index.
When a wavefunction is expressed by atomic orbitals in Eq. (2.11), a couple ofatomic orbitals χ j (r) are prepared at each atom. These atomic orbitals are diﬀerent
in shape. Among the present calculations, one (s-type) orbital is prepared at eachhydrogen (H) atom, while four (s-, p x -, py -, pz -type) atomic orbitals at each carbon(C) atom. Consequently, the basis index j is the composite indices of the atomindex J and the orbital index β ( j ≡ j(J, β)). The projected physical quantityof eTj ρXe j is calculated under the double loop of J and β. In the present code,
the outer loop, the loop with the atom index (J) is parallelized both by the MPI
and OpenMP methods, and a projected physical quantity of eTj ρXe j is calculated
within a thread (a CPU core). An eﬃcient parallel computation is realized on thecondition that the number of atoms is larger than that of cores (N > ncore ).

4.3.3 Sparsity of matrices
Sparsity of the matrices of Hi j and S i j are explained briefly. See Ref. [21] andreference therein for details. As explained in the previous subsection, the indices iand j are the composite indices of the atom indices I and J and the orbital indices αand β, respectively (i ≡ i(I, α), j ≡ j(J, β))). Therefore, an element of the matrices
H and S is expressed by the four indices as HIα;Jβ and S Iα;Jβ , respectively. Since amatrix element value decreases quickly and monotonically as the function of theinter-atomic distance between the I-th and J-th atoms (rI J ), a cutoﬀ distance rcut
was introduced. A matrix element, HIα;Jβ or S Iα;Jβ , is ignored, if rI J > rcut , which

56CHAPTER 4. EXTREME SCALABLE ORGANIC MATERIAL SIMULATIONS
makes the matrices to be sparse. In the present thesis, the cutoﬀ distance rcut is set
to be rcut = 5au(≈ 0.2646nm) for diamond crystal and rcut = 10au(≈ 0.5292nm)
for condensed polymers. A longer cutoﬀ distance is used for condensed polymers,
so as to include the interaction between polymers.

4.3.4 MPI Communication
The communication among nodes is required, only when a summation is performed in the form of Eq. (4.7), as shown in Fig. 4.3(c). See Ref. [21] for detailed
formulations. The code is written in the MPI/OpenMP hybrid parallelization andthe summation is carried out hierarchically; First, the summation is carried out
on each node by OpenMP directives and then the summation is carried out between nodes by MPI_Allreduce(). The pure MPI parallelism is also possible
but consumes larger memory costs.

4.3.5 Use of Extensible Markup Language (XML) in File I/O
In our simulations, the main input and output files are written in the format ofExtensible Markup Language (XML), since the XML format is simple, flexible
and widely used on the Internet [23]. For example, the minimum information for
an atom is written as follows;<atom element ="C"><position unit =" angstrom "> 1.0 d0 0.0 d0 0.0 d0 </position ></atom >The above description means that a carbon atom is located at the position of(x, y, z) = (1, 0, 0), where Angstrom unit (1 Angstrom = 10−10 m) is used. The
method for reading XML files should be chosen properly, according to the purpose. In our simulations, the XML files are read by two methods, Document
Object Model (DOM) method and Simple API for XML (SAX) method. In general, the DOM method is easier in programming and results in huge memory andtime consumption for large-size data, while the SAX method is more diﬃcult inprogramming and results in tiny memory and time consumption. Two input files,
configuration and structure XML files should be prepared for each simulation andthey are quite diﬀerent in their file size. (i) The configuration XML file describes
calculation conditions, such as temperature of the system. A typical file consistsof several tens of lines and the file size does not depend on the system size N. The
configuration file is read by the DOM method, since its file size is always tiny. (ii)
The structure XML file describes the atomic structure data, as shown in the above
example. The structure XML file is read by the DOM method, since the file size is
proportional to the system size N and can be huge. Since the atomic structure data

4.3. METHODOLOGICAL DETAILS FOR LARGE SCALE CALCULATIONS 57
contains three (x, y, z) components in the double precision (8 B) value for each
atom, the required data size with 107 (= 10 M) atoms is estimated to be 3 × 8 B
×10 M = 240 MB. A typical size of the structure XML file with 107 atoms is one
GB. In addition, the parallel file reading is used for large-scale calculations withsplit XML files for the structure file and gives a significant acceleration [38]. The
K computer and FX10 support the parallel file IO, called ’rank directory’ function,at the hardware level and are suitable to the parallel file reading with split XML
files.

4.3.6 Performance tuning on the K computer
Performance tuning on two bottlenecks, matrix-vector multiplication (matvec) in
the iterative Krylov-subspace solver and frequent access to shared disk space, wascarried out on the K computer as follows.
For old implementation, performance profiling using the basic profiler on the
K computer showed that matvec routine in the iterative Krylov-subspace solver
took 39% execution time in the whole program. Because parallelism is alreadyutilized by higher levels of the algorithm (see Sec. 4.3.2), only single node tuning
was enough. Because the matvec routine in old implementation had a loop overfour indices of matrix elements (see Sec. 4.3.3), loop length was short and localityof reference was low, which led to low eﬃciency. Accordingly, the matvec routine
was split into two subroutines. The one is to generate matrix data in compressedrow storage (CRS) format, and the other is to compute matvec with the CRS
format matrix. Now the latter subroutine only has a loop over two indices.
Moreover, matrix data formats other than CRS are implemented for matvec
routine, because the most eﬃcient matrix data format depends on sparsity of matrices. In the present code, the conventional quadruple loop (called ‘quad’), CRS
format (‘crs’) and dense format (‘dense’) are implemented. For CRS and dense
format, matvec routines that utilize symmetry of matrices are also implemented(‘crs-sym’ and ‘dense-sym’, respectively). The dense matvec routines use LAPACK library for eﬃciency. When sparsity of matrices is high and low, the CRS
matvec routines and the dense ones will be eﬃcient, respectively. The user canselect an optimal matvec routine with a configuration file.
The elapsed times per step T elaps by the five matvec routines (quad, crs, dense,
crs-sym and dense-sym) were measured, for two systems that have the same number of atoms N=1,228,800 and diﬀerent sparsity. The systems are a disorderedsystems of condensed polymers called ‘P1’ and the ideal diamond solid called‘D1’. Details of the systems and definition of T elaps will be described in Sec. 4.4.2.
The matrix generated from ‘D1’ is more sparse than ‘P1’ case. The numberof nodes nnode was fixed to 2, 592. Table 4.1 shows comparison of T elaps by thematvec routines. When the best matvec routines for each system are chosen, T elaps

58CHAPTER 4. EXTREME SCALABLE ORGANIC MATERIAL SIMULATIONS
is 1.3 − 1.6 times faster than cases of the conventional quadruple loop code. It wasalso found that a dense matvec is suitable for the dense system ‘P1’, and a sparse
matvec is suitable for the sparse system ‘D1’, as expected. However, utilizationof symmetry of matrices did not contribute much to speed up because the matrix
size for each node was small.
Table 4.1: The measured elapsed times T elaps (sec) by the five matvec routines forthe ideal diamond solid (‘D1’) and the condensed polymer system with (‘P1’).
The best elapsed time for each system is written in bold. Note that the matrices
generated from the two systems have diﬀerent size though the number of atoms N
is the same, due to diﬀerence of elements.
SystemP1
D1

quad crs
16.5 11.2
15.5 11.6

crs-sym dense
11.3
10.6
11.5
22.8

dense-sym
10.4
20.8

On the K computer, access to a disk space that is shared by all the computenodes happens when, for example, nameless file in Fortran is used. Frequent
access to the shared disk space can cause unexpected performance decrement.
To avoid this, an environment variable TMPDIR which specifies the directory
to create temporary files can be set to ’.’, which means the local filesystem ofcompute nodes.

4.4 Benchmarks and their analysis on the K computer
The present benchmarks were carried out so as to show (i) the extreme strong scaling, on the full system of the K computer, and (ii) the time-to-solution qualified
for a practical research. Our target value of the time-to-solution is T elaps = 102 s
for the elapsed time per step in a quantum molecular dynamics simulation, sothat a dynamical simulation of nstep = 103 steps can be executed within one day(T elaps nstep = 105 s ≈ one day).

4.4.1 Architecture
The architecture of the K computer is explained briefly. The K computer achieved
10 PetaFLOPS Linpack performance for the first time in the world.
The K computer consists of 82,944 compute nodes which have one CPU, one
Inter Connect Controller (ICC) , which is responsible for interconnection between

4.4. BENCHMARKS AND THEIR ANALYSIS ON THE K COMPUTER

59

compute nodes, and 16GB oﬀ-chip memory. Each CPU has 8 cores and 6MB L2
cache shared by cores and each core has 4 double precision FMA units, which canbe driven by 2 × 2way SIMD operations, and frequency is 2GHz. Therefore, total
double precision peak performances are 16 GigaFLOPS per core, 128 GigaFLOPS
per CPU and 10.62 PetaFLOPS per whole system, respectively [73].
The interconnection between compute nodes of the K computer is named‘Tofu’ which constructs physical 6 dimensional mesh/torus network topology.
User can use it as logical 3 dimensional torus network with redundant links which
contribute robustness of the connection and flexibility of the operation. Each compute node has 6 links for 3 dimensional directions to connect neighboring nodes
for each direction. Each link has 5GB/s full-duplex bandwidth for each direction.
The oﬃcial operation of the K computer started at September 2012, the K
computer has been kept one of the largest and fastest super computer in the world.
Actually, the K computer keeps top level position in major benchmark rankings.
For example, according to the TOP500 list NOVEMBER 2015 [74], the K computer is ranked in No.4 as well as No.1 in Graph500 list November 2015 [75]which is a benchmark for capability of graph analysis and No.2 in November
2015 HPCG Results [76] which is a benchmark for applications with iterative
matrix solver in Conjugate Gradient method.
The results show that the K computer is applicable to general purpose computing.
In general massive parallel computations, collective communication uses synchronizations among compute nodes, which requires extremely high overheadcaused by increased number of compute nodes. To address the problem, theICC of the K computer provides hardware barrier function which is applied to
MPI_Barrier(),
MPI_Allreduce(), and so on. In particular, the feature is strongly eﬀective in
reduction of MPI_Allreduce().

4.4.2 Conditions of benchmarks
The benchmark jobs were executed in the MPI/OpenMP hybrid parallelization.
The number of the MPI processes is set to that of the compute nodes and thenumber of the OpenMP threads is set to be eight, the number of cores per computenode. All the numerical calculations were carried out in double precision. The
jobs were executed by specifying the three-dimensional node geometry on the(y)(z)
K computer (nnode ≡ n(x)node × nnode × nnode ) for optimal performance or minimum
hop count. The number of used nodes (node geometry) is listed below; nnode= 2, 592(= 12 × 12 × 18), 5, 184(= 12 × 18 × 24), 10, 368(= 18 × 24 × 24),
20, 736(= 24 × 27 × 32), 41, 472(= 27 × 32 × 48), and 82, 944(= 32 × 48 × 54,
the full system). The MPI_Allreduce() optimized on the K computer were used

60CHAPTER 4. EXTREME SCALABLE ORGANIC MATERIAL SIMULATIONS[77].
The simulations were carried out for three disordered systems of condensed
polymers, poly-(phenylene-ethynylene) (PPE). The three systems are called ‘P100’,‘P10’ and ‘P1’ and contain N=101,606,400 (≈ 108 or 100M), N=10,137,600 (≈10M) and N=1,228,800 (≈1M) atoms, respectively. The periodic boundary condition is imposed. The size of the periodic simulation box is 134 nm × 134 nm ×
209 nm for the ‘P100’ system. Among them, the systems with N = 108 , 107 atoms
are the main targets for the full-system calculation. As another reference data foran ‘ideal’ system without structural disorder, the simulations were carried out also
for the ideal diamond solid called ‘D100’ that contains N = 106, 168, 320 (≈ 108
or 100M) atoms in the ideal periodicity, because the structural disorder can be asource of load imbalance.
A calculation job contains the initial preparation procedure including the reading of the input (atomic structure) file and the simulation procedure for five steps,
as shown in Fig. 4.3(a). In the benchmark, a finite-temperature dynamics simulation was carried out with N = 106 , 107 atoms, while a ‘snapshot’ simulation, a
successive calculation of the given atomic structures with slightly diﬀerent valuesof volume, was carried out with N = 108 atoms. The snapshot was carried out
with N = 108 atoms, because of the limitation of built-in memory size (16GB
per node). For example, the position data RI ≡ (xI , yI , zI ) with N = 108 atoms
occupies the memory size of 8B ×3N = 2.4GB. The present snapshot calculation
consumes 9 GB per node. A dynamical simulation should store additional variables such as velocity and force, and requires a larger memory size. The memory
size is estimated to be 60 GB per node, which exceeds the limitation of builtin memory. The above diﬀerence between the simulations does not change theconclusions of the present discussions on elapsed time, since the elapsed time isdominated by that for electronic state calculation.
The elapsed time per step (T elaps ) is recorded by averaging the elapsed times
among the four steps except the first step, because the first step includes several
initial procedures, unlike the other steps. Since the elapsed time for the initial
preparation procedure and the first step is in the same order of that for one step, it
is negligible in practical dynamical simulations with ntime = 103 steps or more.

4.4.3 Preparation of disordered structures of condensed polymers
The initial atomic structures for the disordered systems were generated from classical molecular dynamics simulations by GROMACS [78]. GROMACS is free
software for various kind of classical molecular dynamics simulations and highly
parallelized with MPI. Classical simulations work faster but do not treat elec-

4.4. BENCHMARKS AND THEIR ANALYSIS ON THE K COMPUTER

61

tronic (quantum) waves that are responsible for the device property explained inSec. 4.1.2.
The detailed procedure for generating the three initial disordered structures ofcondensed polymers ‘P100’, ‘P10’ and ‘P1’ was as follows. At first a straightpolymer chain of PPE with 100 monomers (N = 1200 atoms) is generated. Thencopies of the chain are bundled in parallel on a lattice. The number of copies is84672, 8448 and 1024 for ‘P100’, ‘P10’ and ‘P1’ respectively. For further details, when the direction of the chain axis is set to z-axis, the number of copies foreach axis (n x , ny , nz ) is (168, 168, 3), (44, 48, 4) and (32, 32, 1) for ‘P100’, ‘P10’
and ‘P1’ respectively. The numbers are determined so that the entire bundle approximates a cubic region. In the bundles each chain is randomly rotated aroundthe z-axis and random displacement along the z-axis is imposed. By the method,
structures which have moderate density are generated without MD simulations
with pressure coupling.
After that molecular dynamics simulations by GROMACS were carried out forthe bundle structures to generate totally disordered systems. Generalized AMBER
Force Field (GAFF) [79, 80] was used for the force field. In the simulations
the periodic boundary condition is imposed. Total simulation time is 10 ps andtemperature is 1, 200 K. Simulation time step is 10−4 ps and therefore the numberof time steps is 105 . This simulation condition is enough for relaxation of theinitial lattice bundle configuration. All the simulations were carried out on SGI
ICE XA whose each compute node consists of two Intel Xeon 2.5 GHz (12 cores)
processor. The elapsed time was 12,109 sec on 36 nodes for ‘P100’, 4,024 sec on
9 nodes for ‘P10’ and 737 sec on 4 nodes for ‘P1’.

4.4.4 Analysis on strong scaling and time-to-solution
The elapsed time per step was analyzed so as to observe the strong scaling property and the qualified time-to-solution. The measured elapsed time is summarizedin Table 4.2. Here the parallel eﬃciency ratio α is defined byα≡

T elaps (n0 )/T elaps (nnode )nnode /n0

(4.19)

with n0 ≡ 2592. For example, the parallel eﬃciency ratio α with 108 atoms andthe maximum number of nodes (nnode = 82, 944) is α = 0.92 for ‘D100’ andα = 0.75 for ‘P100’.

62CHAPTER 4. EXTREME SCALABLE ORGANIC MATERIAL SIMULATIONS
Table 4.2: The measured elapsed times T elaps (sec) for ideal diamond solid with
108 atoms (‘D100’) and condensed polymer systems with 108 atoms (‘P100’),
with 107 atoms (‘P10’) and with 106 atoms (‘P1’). The ideal or measured speedup ratio is shown inside the parenthesis.
nnode
D100
P100
P10
P1
2,592 (1)1001.4 (1)741.1 (1)81.4 (1)
10.3 (1)5,184 (2)502.2 (1.99) 378.5 (1.96) 43.7 (1.86) 5.95 (1.73)
10,368 (4) 252.6 (3.96) 195.2 (3.80) 24.3 (3.35) 3.28 (3.14)
20,736 (8) 127.9 (7.83) 103.0 (7.19) 11.4 (7.14) 1.96 (5.26)
41,472 (16) 65.6 (15.3)
57.1 (13.0) 6.32 (12.9) 1.25 (8.23)
82,944 (32) 34.1 (29.4)
30.9 (24.0) 3.60 (22.6) 0.84 (12.2)
Table 4.3: Communication time T comm and barrier time T barr of the elapsed time
T elaps . See the caption of Table 4.2 for notations. The values are listed as T comm /T barr (sec).
nnode
D100
P100
P10
P1
2,592 1.04 / 7.16 1.60 / 28.14 0.382 / 6.69 0.0617 / 2.68
5,184 1.04 / 3.34 1.60 / 20.75 0.378 / 4.86 0.0689 / 1.85
10,368 1.05 / 2.22 1.61 / 14.97 0.384 / 4.34 0.0734 / 1.13
20,736 1.05 / 1.34 1.61 / 9.05 0.218 / 3.31 0.0712 / 0.674
41,472 1.06 / 1.03 1.64 / 6.87 0.215 / 2.18 0.0727 / 0.409
82,944 1.06 / 0.485 1.65 / 5.96 0.218 / 1.28 0.0613 / 0.227
Table 4.4: Communication time T comm and barrier time T barr of the elapsed time
T elaps (sec) for the molecular dynamics simulation in the ‘P10’ case.
nnode T elaps
2,592 81.8
5,184 43.9
10,368 24.5
20,736 14.4
41,472 9.38
82,944 6.62

T comm
0.635
0.637
0.674
0.676
0.696
0.677

T barr
6.84
4.95
4.42
3.44
2.37
1.45

The data of Table 4.2 is plotted in Fig. 4.4. One can find that the target timeto-solution of T elaps ≈ 102 s is fulfilled by nnode = 20,736 and 2,592 for the polymer
systems with 108 and 107 atoms, respectively. The above results indicate that thepresent simulation method is qualified for practical device material research with
107 -108 atoms.

63

Elapsed time (s)4.4. BENCHMARKS AND THEIR ANALYSIS ON THE K COMPUTER

Number of nodes

Figure 4.4: Strong scaling benchmarks for ideal diamond solid with 108 atoms(‘D100’) and condensed polymer systems with 108 atoms (‘P100’), with 107
atoms (‘P10’) and with 106 atoms (‘P1’). Dashed lines are drawn for ideal scaling.

Since the computational cost is proportional to the number of atoms in theorder-N method, the number of atoms per cores NA/C ≡ N/nnode is an importantparameter for qualified simulations. The above two cases with the qualified timeto-solution (T elaps ≈ 102 s) give, commonly, the value of NA/C ≡ N/nnode ≈500.
The condition of NA/C = 500 will be a typical condition for qualified simulations.
Figure 4.4 indicates also that the strong scaling property is found in all cases,
because the elapsed time T elaps decreases monotonically as the function of thenumber of used nodes. As explained in Sec. 4.3.2, an eﬃcient parallelism appearson the condition that the number of atoms per core is larger than one (NA/C ≡N/ncore > 1). In the case of ‘P1’, for example, the number of atoms per core is
NA/C ≈ 1.85 with the maximum number of nodes (nnode ≡ 82, 944). The value of
NA/C ≈ 1.85 is close to the lower limit and the simulation shows a lower value ofthe parallel eﬃciency ratio (α = 0.38).

4.4.5 Detailed analysis
Detailed analysis was carried out for the eﬀect of communication time and load
imbalance. In the simulations, not only the total elapsed time T elaps , but also theaccumulated MPI communication time T comm , and the accumulated barrier time
T barr , are recorded at every step on all nodes. The barrier time includes the time
to wait for other processors. Table 4.3 shows the average value among the nodes

64CHAPTER 4. EXTREME SCALABLE ORGANIC MATERIAL SIMULATIONS
for T elaps and T comm and the maximum value for T barr . The communication time
T comm is consumed by inter-node data communications, while the barrier time T barr
appears from a load imbalance.

Figure 4.5 plots the data in Table 4.3. Three points are discussed; (i) The
communication time T comm is almost constant among all the four cases, because
the communication time is consumed by MPI_Allreduce() within a tree algorithm. (ii) Among all the cases, the barrier time T barr seems to include a term
(1)
(1)
T barrthat is proportional to the total elapsed time T elaps (T barr
∝ T elaps ), though the(1)
ratio of T barr /T elaps is not serious. (iii) Among the condensed polymer systems(Figs. 4.5(b)-(d)), the barrier time T barr is always much larger than the communication time T comm and can aﬀect the eﬃciency in parallelism with the maximum
number of nodes (nnode = 82, 944).

To end up this section, two comments are addressed; (I) The further tuning
should be focused mainly on single-core calculations, since the most routines are
executed as single-core calculations as in Fig.4.3(c). The profiler reported that theperformance is 2.3 % of the peak for the ‘P100’ case with nnode = 82, 944 in Table
4.2. The severest limitation in the present calculations is the memory size of the
K computer (16GB per node) and the present code was written in the memorysaving style, in which the memory cost should be minimized and the time costis sometimes sacrificed. Since the situation can diﬀer among materials and/or architectures, a possible way is to add another workflow in the time-saving style.
The routines can be classified into those for the generation of matrix elements
and for the Krylov subspace solver as in Fig.4.3(a). The matrix-vector multiplication gives a large fraction of the total elapsed time, as usual in a Krylov-subspacesolver, and a typical fraction is 21 % among the present condensed polymer systems. The result suggests that the matrix generation part gives a larger fraction.
(II) Table 4.4 shows the detailed elapsed time for the molecular dynamics simulation in the ‘P10’ case, the possible maximum size (See Sec. 4.4.2). For example,(MD)
the elapsed time per molecular dynamics time step is T elaps= 81.8 sec or 6.62
sec in nnode =2,592 or 82,944, respectively. Fig. 4.6 plots the data in Table 4.4
in the same manner of Fig. 4.5(c). For comparison, Fig. 4.6 also shows the data
in Fig. 4.5(c), the data with the electronic structure calculation part. The elapsed
time is much smaller than the target time-to-solution (102 s) and the method isqualified well for a real research. However, non-negligible time costs appear in(MD)
the total elapsed time (T elaps) among the cases with nnode > 2 × 104 , because of
the additional routine for MD simulation. Therefore, further tuning of the codefor faster MD simulations is needed.

4.4. BENCHMARKS AND THEIR ANALYSIS ON THE K COMPUTER

65

(b)

Elapsed time (s)Elapsed time (s)(a)

Number of nodes(c))

Number of nodes

Elapsed time (s)Elapsed time (s)(d)

Number of nodes

Number of nodes

Elapsed time (s)Figure 4.5: Decomposition analysis of the elapsed times. The total elapsed time
T elaps , the communication time T comm and the barrier time T barr are plotted for (a)ideal diamond solid with 108 atoms (‘D100’) and condensed polymer systems (b)with 108 atoms (‘P100’), (c) with 107 atoms (‘P10’) and (d) with 106 atoms (‘P1’).
Dashed lines are drawn for ideal scaling.

Number of nodes

Figure 4.6: Decomposition analysis of the elapsed time for the MD simulation in
(MD)
(MD)
the ‘P10’ case. The total elapsed time T elaps(Tot(MD)), the barrier time T barr
(MD)
(Barr(MD), and the communication time T comm(Comm(MD)) are plotted per MD
step in the same manner of Fig. 4.5(c). The data for the electronic state calculation(Tot, Barr, Comm) are also plotted for comparison.

66CHAPTER 4. EXTREME SCALABLE ORGANIC MATERIAL SIMULATIONS

4.5 Wavepacket dynamics simulation
To investigate transport property of organic device materials, quantum wavepacket
dynamics simulations are carried out. The simulation is based on a Schrödingertype equation∂t Ψ = −iHeﬀ Ψ(4.20)

for a hole wavefunction Ψ(r, t) with a modelled Hamiltonian Heﬀ . In general,
the time interval of simulation step is much smaller than that in the molecular
dynamics simulation (hWP ≪ hMD ) and is hWP = 0.1 fs or smaller. A typicalpost-simulation analysis is the estimation of mobility explained in Sec. 2.4.2. The
estimation procedure is explained; After the wavepacket simulation, the meansquare displacement (MSD) of the wavepacket is calculated at each time step(⟨r2 (t)⟩). Then the MSD behavior is plotted as the function of time. See Fig. 4.7
for typical graphs. One can find a linear behavior in Fig. 4.7 and the diﬀusion
constant D can be fitted as⟨r2 (t)⟩.
t→∞2t(4.21)De
kB T(4.22)

D = lim
Finally mobility µ is given asµ=

under Einstein’s relation in Sec. 2.4.2.
In a previous paper, [55] a wavepacket dynamics simulation was carried out
for an organic polymer, in which a π-orbital-only model is constructed for disordered atomic structures in dynamical simulations. The analysis of the result
gives mobility. The present method is a theoretical generalization of the previous one. The methodological details will be explained in rest of this section anda brief outline is explained here. The initial wavepacket Ψ(r, t = 0) is set as an
eigenstate, such as the HO state. Two kinds of dynamical simulation methods areconstructed; The first simulation method is called ‘multi-time-scale simulation’ or
the combined simulation between molecular dynamics and wavepacket dynamics
simulations. In the simulation, a double loop structure is used in the code for thetime evolution; The outer loop is used for the atomic thermal motion, a slower dynamics, while the inner loop is the quantum wavepacket dynamics, as explainedin Sec. 4.5.1. The second simulation method is called ‘simulation with modeled
atomic motion’, in which the atomic motion is introduced as a modeled stochastic
term in the wavepacket dynamics. The latter method gives a faster simulation.

4.5. WAVEPACKET DYNAMICS SIMULATION67[ A 2]

mean squaredisplacementabcTime [ps]Figure 4.7: Examples of the mean square displacement during the wavepacket
dynamics simulations. The three simulated samples of (a)(b)(c) are disordered organic polymer (poly-(phenylene-ethynylene);PPE) with disordered structures by
the thermal motion.

4.5.1 Concept of multi-time-scale quantum wavepacket simulations
Rest of this section is devoted to the investigation on the practical methods ofquantum wavepacket dynamics, as a real-time non-stationary simulation, in which
the atomic position and the wavepacket are updated iteratively. The method should
be beyond the adiabatic theory and is not yet well established. The purpose of
the present method is to construct a modeled theory for a practical large-scale
device material simulation within a moderate computational cost. The simulated
system is a material with a single carrier, hole or excited electron, and the carrier
is described as a wavepacket. Such simulations were carried, for decades, on sitemodels, in which a (small) molecule is treated as a site. See Ref. [81] for example.
In the present thesis, a multi-time-scale simulation is adopted, since the atomic
motion is much slower than the wavepacket dynamics. In other words, the situation is assumed to be nearly adiabatic. Therefore, the time evolution is realized
by a double loop structure in the code; The atomic position is updated in the outerloop, while the wavepacket is updated in the inner loop. The time interval of theouter loop is set to be h = 40au(≈ 1fs), a typical value for molecular dynamics

68CHAPTER 4. EXTREME SCALABLE ORGANIC MATERIAL SIMULATIONS
simulation, while the time interval of the inner loop should be set to be hWPD ≪ h.

4.5.2 Basic mathematical formulation of generalized eigenvalue
equation
The basic mathematical formulation is summarized for the generalized eigenvalue
equation
H0 yk = λk S yk .

(4.23)The eigenvectors satisfy the orthogonality relationyk† S yl = δkl .

(4.24)

Eq. (4.23) is rewritten as∑

Hi j Y jk = λk∑j

S i j Y jk

(4.25)jor
HY = S YΛ,(4.26)

where yk ≡ (Y1k , Y2k , ...., Y Mk ), Λ ≡ diag(λ1 , λ2 , ...λ M ). Eq. (4.24)) is rewritten as
Y †S Y = I(4.27)

Y −1 = Y † S .

(4.28)

orThe wavepacket vector on the LCAO basis u ≡ (u1 , u2 , ..., u M )T is expressedby the linear combination of the eigenvectors∑ui =αk Yik
(4.29)koru = Yα.

(4.30)

with α ≡ (α1 , α2 , ..., α M )T . Eqs. (4.30) and (4.28) give
α = Y −1 u = Y † S u.

(4.31)

4.5. WAVEPACKET DYNAMICS SIMULATION

69

4.5.3 Multi-time-scale simulation (I)Our multi-time-scale simulation is explained. The wavepacket dynamics, the inner loop, is considered for the time period ofnh ≤ t ≤ (n + 1)h,(4.32)

where n indicates the MD time step or the iteration number of the outer loop. The
wavepacket is written as the linear combination of atomic orbitals {χ(n)j (r)} as∑(n)Ψ(r, t) =(4.33)u(n)j (t)χ j (r).
j

The suﬃx of ‘(n)’ indicates that the used atomic orbital basis set is that for theatomic position at the n-th MD time step (t = nh). The atomic orbitals {χ(n)j (r)}
are real functions. Here we consider an eﬀective Schrödinger-type equation of the
wavepacket∂t Ψ = −iĤeﬀ Ψ

(4.34)

with an eﬀective Hamiltonian operator Ĥeﬀ . Eqs. (4.33) and (4.34) lead us to theequation of
S (n) u̇(n) = −iH (n) u(n) ,with the definitions of(4.35)

∫
Hi(n)
j

≡

S i(n)
j

≡

∫

(n)
χ(n)
i (r) Ĥeﬀ χ j (r)dr(4.36)(n)χ(n)
i (r)χ j (r)dr.

(4.37)As a practical methodology, the matrix of H (n) is set to be that in Eq. (4.23) (H (n) =H0(n) ). When the generalized eigenvalue equation(n) (n)H (n) yk(n) = λ(n)k S yk

(4.38)

is solved, the wavepacket vector u(n) (t) is expressed by the eigenvectors of Eq. (4.38)
as∑(n)u(n) (t) =α(n)(4.39)
k (t)ykk(n)(n) Twith the coeﬃcient vector of α(n) ≡ (α(n)
1 , α2 , ...α M ) . Since Eqs. (4.35) and(4.39) are reduced to(n) (n)α̇(n)k (t) = −iλk αk (t),(4.40)

70CHAPTER 4. EXTREME SCALABLE ORGANIC MATERIAL SIMULATIONS
the time evolution of α(n)k (t) is solved analytically asα(n) ((n + 1)h) = D(n) α(n) (nh)(4.41)

with the diagonal matrix of
(n)

(n)

(n)

D(n) ≡ diag(e−iλ1 h , e−iλ2 h , ..., e−iλM h ).

(4.42)At the time of t = (n + 1)h, the wavepacket Ψ should be expressed as the linear
combination for the updated atom position (R(t = (n + 1)h))∑Ψ(r, t) =u(n+1)(t)χ(n+1)(r)(4.43)jjj

for the time period of(n + 1)h ≤ t ≤ (n + 2)h(4.44)

and a relation between u(n+1) ((n + 1)h) and u(n) ((n + 1)h) is needed for the simulation. A reasonable relation isu(n+1) ((n + 1)h) := u(n) ((n + 1)h),(4.45)

because the relation will be exact, when all the atoms are not mobile and the basis
set is unchanged (χ(n+1)(r) = χ(n)jj (r)).
Since Eqs. (4.30) and (4.31) giveu(n) ((n + 1)h) = Y (n) α(n) ((n + 1)h)α(n+1) ((n + 1)h) = Y (n+1)† S (n+1) u(n+1) ((n + 1)h),(4.46)(4.47)

the relation between α(n+1) ((n + 1)h) and α(n) ((n + 1)h) is obtained asα(n+1) ((n + 1)h) = A(n) α(n) ((n + 1)h)(4.48)

A(n) ≡ Y (n+1)† S (n+1) Y (n) .

(4.49)

with

Eqs. (4.41) and (4.48) result inα(n+1) ((n + 1)h) = A(n) D(n) α(n) (nh).

(4.50)

Finally, the normalization procedure is imposed on the resultant vector of α(n+1) ((n+1)h), since Eq. (4.45) introduces a deviation in the normalization condition, when4.5. WAVEPACKET DYNAMICS SIMULATION

71

S (n+1) , S (n) . All the multi-time-scale simulations in Sec. 4.6.2 were carried out
by the above formulation.
In Eq. (4.30), the eigenvectors for the expansion of the wavepacket vector canbe partially selected. This approximation is called eigenvector filtering. Namely,a matrix of filtered eigenvectors
Yf = (yl yl+1 ... yl+m−1 ) ∈ C M×m(4.51)

is used instead of Y. Small m is advantageous in computational cost, and unphysical artifacts can be removed when the filtering region l, . . . , l + m − 1 is adequately
selected. Although Yf−1 is no longer exists in this case, the relations Yf† S Yf = Im
and Yf+ = Yf† S hold instead of Eqs. (4.27) and (4.28), respectively, where Yf+ isthe pseudo-inverse of Yf . In general, for any (not necessarily square) matrix A, its
pseudo-inverse A+ solves the least square problem
minx ||Ax − b||2

(4.52)by x = A+ b. In this case, A = Yf . By these relations, almost the same formulation
with the original case is valid. Hereafter the subscript f is dropped.

4.5.4 Multi-time-scale simulation (II)
Although the methodology in the previous subsubsection can realize a multi-timescale simulation code, the present subsubsection will introduce several corrections for better description of the wavepacket dynamics. The corrections wereintroduced from the energetic aspect based on the fact that, since the wavepacket
dynamics is motivated from the atomic thermal motion, the energy of wavepacketuT Hu should fluctuate only within the energy scale of the atomic thermal motion.
The energy discontinuity can appear in Eq. (4.48) and the correction is introducedbyα(n+1) ((n + 1)t) = Ãα(n) ((n + 1)t)Ãi j ≡ δ̃(λi(n+1) − λ(n)j )Ai j(4.53)(4.54)

with the ‘smoothed’ delta function ofδ̃(x) ≡ √1
2πbe−(1/2)(x/b) .
2(4.55)The energy parameter b(> 0) should be given in the order of the thermal fluctuation. The correction is introduced so that the transition between diﬀerent eigenstates ( j → i) can occur only when the two eigenenergies is so close that|λ(n+1)− λ(n)
ij | ≤ b.

(4.56)

72CHAPTER 4. EXTREME SCALABLE ORGANIC MATERIAL SIMULATIONS
Another correction is introduced by the fact that the thermal relaxation eﬀect
should induce an energy flow from the wavepacket into the atomic motion. The
standard equilibrium state theory gives the carrier energy distribution in the Boltzmann distribution (∝ exp(−|λ − λF |/kB T )). If the k-th eigenenergy is far from theFermi level (|λk − λF | ≫ kB T ), its component αk tends to be relaxed into otherlevels nearer the Fermi level. Such relaxation eﬀect is realized, when the diagonal
matrix D(n) is replaced by
(n)

(n)

(n)

(n)

D̃(n) ≡ diag(e−iλ1 h−γ1 h , e−iλ2 h−γ2 h , ...e−iλk h−γk h , .., e−iλM h−γM h ).

(4.57)

with the definition ofγk ≡|λk − λF |.
kB T τ

(4.58)The relaxation eﬀect is governed by the time parameter τ(> 0), since the k-th
eigenstate component αk decays as()
h |λk − λF |∝ exp −
(4.59)
τ kB Twith the time interval of h. In short, the parameter τ means a relaxation time,since the k-th eigenstate component with |λk − λF | ≈ kB T decays in the time orderof τ. It might be noteworthy that Eq. (4.57) can be understood, formally, as theintroduction of the imaginary part of the eigenenergy
( (n)
)
(
)
(n)
D̃(n)
(4.60)kk ≡ exp −iλk h − γk h = exp −i(λk − iγk )h .
The above two corrections modify Eq. (4.50) intoα(n+1) ((n + 1)h) = Ã(n) D̃(n) α(n) (nh).

(4.61)The method with the corrections is now under test, in which the test simulations
are running with several diﬀerent values of the parameter set. The computational
cost of the corrections is negligible. Quantitative discussion will be a future topic.

4.5.5 Simulation with modeled thermal atomic motion
The present subsection explains the simulation method with modeled atomic motion, in which the atomic thermal fluctuation is treated as a modeled stochastic
term in the wavepacket dynamics. In the simulation, only the input matrices are
required only for the initial structure (H (0) , S (0) ). When Eq.(4.35) is solved, theHamilton matrix at the n-th MD step, H (n) , is determined by(n)
H (n) := H (0) + Htherm(4.62)

4.5. WAVEPACKET DYNAMICS SIMULATION73(n)
is introduced for the atomic thermal fluctuationand S (n) := S (0) . The matrix Htherm
term. The matrix was chosen to be the diagonal and the element is( (n) )kB T
Htherm =cos 2πxi(4.63)
ii
2

with the temperature of T and a random number of xi . The random numbers of{xi }i are updated with the time interval of MD (h = 1fs) or the time scale of theatomic thermal motion. The same random number is applied for the bases on
the same atom. The above method will be exact in the zero-temperature limit, inwhich the matrices of (H (n) , S (n) ) are unchanged throughout the simulation.
This method runs faster than the multi-time-scale simulation in the previous
subsections. We believe that the two methods are complementary. Since our
purpose is the systematic survey among diﬀerent samples, so as to obtain insights
for better device materials and their fabrication process, the problem is the balance
between the computational cost and the accuracy. If the simulation runs faster, one
can investigate a dynamics simulation among many sample in diﬀerent structures.
The number of the samples is crucial, in particular, among data scientific researchand the present method should be fruitful.

4.5.6 Calculation of diﬀusion coeﬃcient
As explained in Sec. 2.4.2, the mobility µ is proportional to the diﬀusion coeﬃcient D (Eq. (2.36)). In the present thesis, D is determined from time evolutionof the wavepacket Eq. (4.33) as follows. At first the Mulliken charge [82] for j-th
atomic orbital m j (t) is defined as[](n)(n) (n)m j (t) ≡ Re u j (t)S u (t) .
(4.64)
Note that a normalization condition∑m j (t) = 1(4.65)

jis satisfied due to the normalization condition of α(t), Eqs. (4.30) and 4.27). Denoting the coordinate of the atom which the j-th atomic orbital belongs to by r j (t),the mean square displacement (MSD) ⟨r(t)2 ⟩ is determined as∑⟨r(t)2 ⟩ ≡m j (t)∥r j (t) − ⟨r(t)⟩∥2(4.66)jwhere ⟨r(t)⟩ is the mean coordinate of the atoms weighted by the Mulliken charge∑⟨r(t)⟩ ≡m j (t)r j (t).
(4.67)
j

74CHAPTER 4. EXTREME SCALABLE ORGANIC MATERIAL SIMULATIONSIn the same manner as Sec. 2.4.3, diﬀusion constant for the wavepacket simulations is determined as the slope ofd(t) ≡ ⟨r(t)2 ⟩/2.

(4.68)

Because d(t) is generally not a linear function of time, the slope is calculated byleast-square fitting. Moreover, trajectories of d(t) by wavepacket simulations from
diﬀerent initial atomic structures should be averaged for reliable values of D.

4.6 Organic material simulations
Organic materials play a crucial role among next-generation IoT products, such as
display, battery and sensor, since they form flexible atomic structures and enableultra-thin, light, flexible (wearable) devices with a low fabrication cost. Since theatomic structure of organic device materials is disordered, the simulation of suchmaterials requires 100-nanometer-scale systems, which is beyond the computational limit of the present standard electronic structure calculations.
The transport of organic device is governed by π-type wavefunctions that appears typically on benzene-ring regions. Figure 4.8(a) shows an example of π-type
wavefunctions in an organic polymer, poly-((9,9)-dioctyl fluorene) with n = 2.
The transport property stems from both ballistic and diﬀusive conduction mechanisms and can be simulated by wavepacket dynamics. In an ideally orderedsample, for example, several wavefunctions are extended throughout the sample
and the sample shows ballistic conduction. In a fairly disordered sample, on theother hand, wavefunctions are localized and the sample shows diﬀusive conduction motivated by the atomic thermal motion.

4.6.1 Quantum molecular dynamics simulation
Quantum molecular dynamics simulations were carried out for condensed organic
polymers using the order-N method. Figure 4.9 shows an example of bundle-likepoly-(phenylene vinylene) (PPV), The sample consists of 169 PPV polymers andeach polymer has n = 50 monomers.
The length of a polymer is L ≈ 40nm. The total number of atoms is N =117, 962. The van der Waals force [83] is included in the simulation. The parallel
computation was carried out by 960 nodes (15,360 cores) of the supercomputer
Oakleaf-FX of the University of Tokyo. The time interval of simulation step is
hMD = 1 fs. A finite-temperature simulation in T = 600 K was performed fornstep = 5000 iteration steps or the period of hMD nstep = 5ps. The total elapsed
time is approximately 10 hours.

4.6. ORGANIC MATERIAL SIMULATIONS

75

Figure 4.8: (a) Example of organic polymer: poly-((9,9) dioctyl-fluorene) (PFO)with n = 2. The HO state is drawn as an example of π wavefunction. (b) Quantum
wavepacket dynamics on a polymer of poly-(phenyleneethynylene). The numberof monomer unit of n = 1000 and the length of the polymer chain is L ≈ 700 nm.
The figure is a close-up. (c) Quantum wavepacket dynamics on pentacene thinfilm with single layer. (d) A condensed polymer of poly-((9,9)-dioctyl fluorene).
The three polymers in a periodic cell form an amorphous-like structure. The threepolymers are painted in diﬀerent colors for better understanding. (e)-(f) Quantum
wavepacket dynamics on the condensed polymers. The initial state at t = 0 (e)and the final state at t = 1 ps (f) are drawn.

76CHAPTER 4. EXTREME SCALABLE ORGANIC MATERIAL SIMULATIONS

Figure 4.9: Quantum molecular dynamics of poly-(phenylene vinylene) (PPV) at
the initial time (t = 0) and the final time (t = 5ps).

4.6.2 Transport calculations with wavepacket dynamics simulations
Several results for transport calculation of organic materials are discussed below. The calculations were carried out with the multi-time-scale simulation, except where indicated. Figure 4.8(b) shows the wavepacket dynamics simulationin organic polymers of poly-(phenyleneethynylene), where the maximum polymer lengths is L = 700 nm the figure is shown as a close-up. A characteristic semi-localized wavepacket was observed, since the wavepacket spreads over
several monomer units but does not over the whole region of the polymer. Preliminary results on calculated mobility by the present method are consistent to
the experimental trend [55] in which the mobility of meta (zig-zag) type polymers is larger than that of para (linear) type polymers. The mobility values were
measured as the average of three samples in the simulation method with modeled atomic motion. The calculated mobility values of the meta- and para-type
polymer are (µmeta , µpara ) = (1.41, 0.98) [cm2 / V s], which are comparable tothe experimental values of (µmeta , µpara ) = (2.1, 0.75) [cm2 / V s]. Figure 4.8(c)shows the wavepacket dynamics simulation in a disordered pentacene thin filmwith a single layer. The periodic boundary condition is imposed on in the upper and lower area in Fig. 4.8(c). The wavefunction propagates from a molecule
into another. Figure 4.8(d)-(f) shows condensed polymers of poly-((9,9)-dioctyl
fluorene). The initial structure was generated by a classical molecular dynamics
simulation, so as to save the total computational time in research. Three polymers
in a periodic cell form an amorphous-like structure, as shown in Fig. 4.8(d). The
number of atoms in a polymer is Npolymer = 692 and the total number of atoms is
N = 692 × 3 = 2076. Figure 4.8(e)-(f) shows a result of wavepacket dynamics.

4.7. DATA SCIENTIFIC RESEARCH FOR PRESCREENING77

The initial wavepacket is localized in one polymer and extends over diﬀerent polymers. More quantitative discussions are under development and will be a future
topic.

4.7

Data scientific research for prescreening

In this section, a machine learning approach to the research of organic polymer device materials is presented 1 . The purpose is a pre-screening procedure of atomic
structures before large-scale quantum wavepacket dynamics calculations. Atomic
structures of condensed polymers are highly varied due to their flexibility and disorder. Therefore we should filter hopeful structures in the sense of high mobility
before we actually execute quantum wavepacket dynamics calculations. The previous simulation research [5, 6] implies that the relative rotation angles between
benzene rings strongly aﬀect locality of π-electron wavefunctions and then mobility.
As a method for the pre-screening procedure, clustering of disordered polymer structures by locality of π-electron wavefunctions was carried out. K-means
method, which is a simple clustering method [84], is used. K-means method canbe found in standard textbooks of machine learning, such as Ref. [85]. At first
PPE polymer chains with 240 atoms (20 monomers) whose relative rotation angles between adjacent benzene rings are given from the normal distribution withstandard deviation of 20 / 60 degree are generated 100 times for each class (total
N = 200 samples). After the rotations, small fluctuations taken from the normal
distribution with standard deviation of 0.01 Åis given to all the coordinates of all
atoms to remove degeneracy. For each structure (i = 1, . . . , N), a singleshot calculation by ELSES generates Hamiltonian H (i) and overlap matrix S (i) , and thegeneralized eigenvalue problem is solved numerically as(i) (i) (i)H (i) y (i)j = λ j S y j ( j = 1, . . . , M).

(4.69)The matrix size M, which is determined by the number of each element in theatomic structure, was M = 714. Then the K-means clustering method is applied
for the atomic structures using the list of the participation ratios of the eigenvectors
defined in Sec. 2.4.1(i) Tp(i) := (PR(y1(i) ), . . . , PR(y M))(4.70)

as a feature vector for the i-th structure.
1

H. Imachi, Y. M. Tsai, T. Hoshi, and W. Wang, A data scientific research on organic polymer
device materials, The 19th Asian Workshop on First-Principles Electronic Structure Calculations,
National Chiao Tung University, Hsinchu, Taiwan, 31. Oct.- 2. Nov., 2016.

78CHAPTER 4. EXTREME SCALABLE ORGANIC MATERIAL SIMULATIONS

Figure 4.10: Eigenvector participation ratio diagrams of diﬀerent atomic structureclasses. The x-axis is the eigenvalue and the y-axis is the participatin ratio of thecorrespoiding eigenvector. (a) and (b) are diagrams of a structure with adjacentrotation angles in standard deviation of 20 and 60 degree, respectively. (c) and (d)are the same diagrams in partial eigenvalue range (−0.6 ≤ λ ≤ −0.4 au).
Figure 4.10 shows typical examples of the eigenvector participation ratios forthe two structure classes. It indicates that participation ratios of π-electrons arestrongly aﬀected by the disorder in rotation angles. When the disorder in rotation
angles is small (Fig. 4.10(a)), π-electrons are delocalized and their participation
ratios are high, and vice versa, because the participation ratio of an eigenvector indicates its spatial extent. For each structure classes, a typical shape of π electronic
wavefunction is visualized in Fig. 4.11. One can confirm the relation between thedisorder in rotation angles and the participation ratio, because in Fig. 4.11(b), thewavefunction is more localized than in Fig. 4.11(a) due to the high rotation angles
at the edges of the wavefunction.
The number of clusters, a parameter for the K-means algorithm, was two.

4.7. DATA SCIENTIFIC RESEARCH FOR PRESCREENING

79

Figure 4.11: Typical atomic structures and π electronic wavefunctions of diﬀerent
atomic structure classes. Structures of (a) and (b) have adjacent rotation angles instandard deviation of 20 and 60 degree, respectively. The participation ratios of
the wavefunctions (a) and (b) are 122.7 and 48.9, respectively.
Figure 4.12 shows the result of the K-means clustering. Each atomic structure isplotted on the plane by two statistical quantities with markers corresponding to
the cluster labels given by the K-means algorithm. As a result, the two clusters
perfectly matched to the structure classes. Namely, all the structures with adjacentrotation angles in standard deviation of 20 degree are clustered into one group, andall the structures with adjacent rotation angles in standard deviation of 60 degree
are clustered into the other group.

Figure 4.12: K-means clustering of polymer structures by their eigenvector participation ratios. Each structure is plotted by two statistical quantities. The x-axis
is the mean of the rotation angles between adjacent benzene rings. The y-axis is
the resulting standard deviation of the small fluctuation to the coordinates. Twomarkers, red circle and blue triangle, correspond to the two clusters labeled by theK-means algorithm.
The result of K-means clustering indicates that the PRs of the eigenvectors in-80CHAPTER 4. EXTREME SCALABLE ORGANIC MATERIAL SIMULATIONS(b)= 41 nm= 7 nm(a)= 9 nm
(c)

(d)

time

= 53 nm

Figure 4.13: Visualization of partial regions in a 108 -atom system of condensed
polymers (PPE). The system has the periodic cell lengths of (265nm, 206nm,
239nm). (a) The region has the edge lengths of (L x , Ly , Lz )=(53 nm, 41 nm, 48nm) and contains approximately 8×105 atoms. (b) The region has the edge lengthsof (L x , Ly , Lz )=(9 nm, 7 nm, 8 nm) and contains approximately 4 × 103 atoms. (c)(d) Schematic figures of a local network of connected polymers.

herit the disorder property in atomic structures of polymers, and inversely, atomic
structures can be classified using the PRs before executing wavepacket simulations for them.

4.8 108-atom simulation
This section is devoted to the scientific investigation with the ‘P100’, 108 atom
(largest) system of condensed organic polymers. A smaller system with 105 atoms
was calculated as a finite temperature simulation by 104 cores and the wholeelapsed time is 10 hours for 5000 iteration steps, as reported in Sec. 4.5.1
Such a dynamical simulation is, however, impractical with 108 atoms at thepresent day and this section demonstrates a part of the possible future research.

4.8.1 Network analysis of electronic wavefunctions
A large-scale data analysis was carried out for a fundamental research of electronic devices, so as to characterize the propagation of electronic wave through adisordered structure.
Figures 4.13(a)(b) show partial regions of a 108 atoms and one can observe that

4.8. 108 -ATOM SIMULATION81

the structure is disordered. In general, electronic wavefunctions are localized ina disordered structure. The electrical current can propagate among polymers thatare ‘connected’ locally by π electronic wavefunctions, as explained in Sec. 4.1.2.
Therefore, the local network of connected polymers should be investigated.
Figures 4.13(c)(d) show schematically small local networks of connected polymers. Three polymers are drawn and atoms are depicted as filled circles. The
figures include a small local network that consists of two polymers connected bya dashed line. The networks are dynamical, as seen later in this section.
A network analysis on connected polymer networks was carried out fromelectronic states, since a connection between polymers is formed by the electronic waves that spread over the polymers. The analysis was carried out with theGreen’s function G obtained by the parallel order-N simulation, as follows; Stage
I: The order-N simulation gives a ‘connectivity’ matrix of BI J∑∑
ρIα;Jβ H Jβ;IαBI J ≡(4.71)
αβ

where I and J are the atom indices. The connectivity matrix is called integrated
crystal orbital Hamiltonian population (ICOHP) among physics papers [86, 59].
The quantity is a partial sum of the electronic energy ⟨H⟩ in Eq. (4.3) (⟨H⟩ =∑I J BIJ ). Since the matrix elements are calculated always during the parallelorder-N simulation, the elements can be obtained independently among nodes,
without any additional operation or communication cost. A matrix element BI J
has a physical meaning of a local bonding energy between the I-th and J-th atoms;
If the value of |BIJ | is significantly large, the two atoms are ‘connected’ by electronic wave. The matrix B is sparse, since only atoms within a short distancecan be ‘connected’ with each other. Stage II: Since every atom belongs to one ofpolymers, the connectivity matrix for polymers is defined byB(poly)
≡
PQ

∈Q
∈P ∑
∑
I

BIJ ,

(4.72)

J∑
where the summation of ∈P
I , for example, means the summation among the atoms
that belong to the P-th polymer. If an element B(poly)
PQ shows a meaningful non-zerovalue, the P-th and Q-th polymers are connected by electronic wave. The matrix
B(poly) is sparse. The dimension of B(poly) is equal to the number of polymers
N (poly) = 83, 349 and is much smaller than that of B (N = 108 ). Stage III: Asa coarse grained analysis, the eigenvalue equation of B(poly) z = λz in the matrix
dimension of N (poly) was solved by the parallel eigenvalue solver [4]. Then theparticipation ratio PR(z) is computed for each eigenvector z.

82CHAPTER 4. EXTREME SCALABLE ORGANIC MATERIAL SIMULATIONS

Figure 4.14: Participation ratio histogram of the eigenvectors {zk } of the connectivity matrix for the PPE polymers B(poly) . The highest participation ratio
PR(z61731 ) = 2.814 is achieved for the eigenvector z61731 .

Table 4.5: The number of eigenvectors zk of the connectivity matrix for the PPE
polymers B(poly) in some characteristic ranges of the participation ratio PR(zk ).
Note that PR(zk ) ≥ 1.0 holds for any k because the standard eigenvalue problem
is considered here (See Sec. 2.4.1).
Range of PR(zk )1.0 ≤ PR(zk ) < 1.05
1.05 ≤ PR(zk ) < 1.5
1.5 ≤ PR(zk ) < 2.0
2.0 ≤ PR(zk ) < 2.5
2.5 ≤ PR(zk )Count of eigenvectors zk
80,521
2,135
654
36
3

4.8. 108 -ATOM SIMULATION

83

Fig. 4.14 shows the participation ratio histogram for the eigenvectors of B(poly) ,and Table 4.5 summarizes the histogram by selecting some characteristic ranges.
Several eigenvectors with participation ratios PR(z) ∼ 3.0 are found, which means
the presence of small local polymer networks with several (about 3) polymers, as
illustrated schematically in Fig. 4.13(c). Such networks are responsible for theelectrical current, which will be confirmed in Sec. 4.8.2. Fig. 4.14 and Table 4.5
also show the almost all of eigenvectors have the participation ratios PR(z) ∼ 1.0,
which means that many other polymers are almost isolated or are not connected
to others.

4.8.2 Quantum wavepacket dynamics simulation on polymer
networks
The quantum wavepacket dynamics simulation was carried out, so as to confirm
that electronic wave can propagate on the small polymer networks detected in theprevious subsection. From the eigenvectors z with high PR value PR(z), the small
polymer networks were extracted as isolated systems by selecting several polymer
indices which have the largest absolute values. Then the wavepacket simulations
were carried out for them.
Several tens of the wavepacket simulations were carried out in diﬀerent atomic
structures and initial wavepackets (Ψ(r, t = 0)). Figure 4.15 shows a typical result,in which the dynamical simulation was carried out approximately for 1 ps. The
initial structure is extracted from the eigenvector z61731 with the highest participation ratio shown in Fig. 4.14. Because the participation ratio is PR(z61731 ) =2.814 ∼ 3, 3 polymers are selected by the 3 indices with the largest absolute values
in z61731 . The simulation consumes six hours with 128 nodes. The wavepacket aswell as the atom positions change dynamically, though the atomic motion is muchslower than that of the wavepacket. In Fig. 4.15, the initially localized wavepacket
propagates first through one polymer (Fig. 4.15(a)→(b)) and then propagates toothers (Fig. 4.15(b)→(c)). The above behavior is reasonable, because the connection between polymers is much weaker than that within a polymer. The above
observation confirms that the network analysis is fruitful for the investigation of
charge propagation in organic device materials.

84CHAPTER 4. EXTREME SCALABLE ORGANIC MATERIAL SIMULATIONS

70 nm

(a)

(b)

(c)

(d)

(e)

(f)

15 nm

Figure 4.15: Quantum wavepacket dynamics of the condensed polymers. The
charge density |Ψ(r, t)|2 is plotted at (a) t = 0, (b) t = 50 and (c) t = 948fs,
respectively. Figures (d), (e) or (f) is a close up of (a), (b) or (c), respectively.

Chapter 5
Summary and future outlook
Numerical linear algebraic methods were developed for large-scale electronic
state calculations on modern massively parallel supercomputers, for example, thefull system of the K computer.
The optimal hybrid generalized eigenvalue solver was constructed with thethree parallel dense-matrix solver libraries of ELPA, EigenExa and ScaLAPACK.
The hybrid solvers with ELPA and EigenExa give better benchmark results thanthe conventional ScaLAPACK library. The code was developed as a middlewareand a mini-application and is available online. The benchmark was carried out
with up to a million dimensional matrix on the K computer and other supercomputers. The decomposition analysis of the elapsed time reveals a potential bottleneck part on next-generation (exa-scale) supercomputers, which indicates theguidance for future development of the algorithms and the codes.
An extreme scalability was realized by the code optimization based on thenovel Krylov-subspace solver. The simulation was carried out with 108 -atom simulations for flexible (organic) electronic devices, key devices of next-generation
IoT products. The mathematical foundation is generalized shifted linear equations((zS − H)x = b), instead of conventional eigenvalue equations. Since the foundation has a highly parallelizable mathematical structure, the benchmark shows
an extreme strong scaling on the full system of the K computer and a qualified
time-to-solution for material researches in industrial purpose. The method wasapplied to the condensed polymer materials in academia-industrial collaboration.
A combined research of the quantum simulation and the data analysis reveals thatpolymers form small local networks and electronic waves propagate on the networks. The present research will give insights of flexible devices and their fabrication process. Since the present method is general, the combined research of huge
(100-nm-scale) quantum material simulations and the post-simulation data analysis will be fruitful not only for organic devices but also for many other industrial
materials.
85

CHAPTER 5. SUMMARY AND FUTURE OUTLOOK

86

Finally two points are discussed as the future outlook; The first point is todevelop eﬃcient numerical method for the real-time quantum dynamics simulations with better computational speed and higher reliability. Very recently, forexample, we have developed a mini-application for a novel structure-preserving
numerical solver [87] with the two-dimensional non-linear Schrödinger equation
1. The second point is the convergence between simulation and machine learning.
Preliminary results are reported on organic polymer device materials in this thesis. The convergence will realize the industrial application of a high-throughput
material design.

1

S. Kudo, H. Imachi, Y. Miyatake, Y. Yamamoto, and T. Hoshi, Application of themathematical-structure-preserving method to computational material science, The 2nd CDSMI
Symposium : Creation of new functional Devices and high-performance Materials to Support
next-generation Industries, The University of Tokyo, 6-7. Dec. 2016.

Acknowledgment
The author thanks to Takeo Hoshi, the supervisor of this work, and colleagues fortheir devoted guidance and discussions. The author also thanks to Toshiyuki Imamura in RIKEN Advanced Institute of Computational Science (AICS) and Takeshi
Fukaya in Hokkaido University for fruitful discussions on EigenExa. The author
has greatly benefited from Kiyoshi Kumahata, Masaaki Terai, Kengo Miyamoto,Kazuo Minami and Fumiyoshi Shoji in RIKEN AICS for helpful advice on performance tuning of ELSES. Tomofumi Tada in Tokyo institute of Technology gives
the author invaluable comments about quantum wavepacket simulation methods.
The author is deeply grateful to Yusaku Yamamoto and Shuhei Kudo in The University of Electro-Communications, and Shao-Liang Zhang, Tomohiro Sogabe
and Dongjin Lee in Nagoya University, and Takafumi Miyata in Fukuoka instituteof technology for insightful comments on numerical linear algebra algorithms.
The author appreciates the suggestions on structure-preserving methods by Yuto
Miyatake in Nagoya University and Takayasu Matsuo in The University of Tokyo.
The author also appreciates precise feedback by Wei-Chung Wang and Yuhsiang
Tsai in National Taiwan University. The author also thanks to Masaya Ishida(Sumitomo Chemical Co.) for providing several atomic structure data.
The K computer of RIKEN was used in the research projects of hp150144,
hp150281, hp160066 and hp160222.
The supercomputer Oakleaf-FX of The University of Tokyo was used in the
research project of 14-NA04 in ‘Joint Usage/Research Center for Interdisciplinary
Large-scale Information Infrastructures’ in Japan, in the ‘Large-scale HPC Challenge’ Project, Information Technology Center, The University of Tokyo and Initiative on Promotion of Supercomputing for Young or Women Researchers, Supercomputing Division, Information Technology Center, The University of Tokyo.
The supercomputer of Kyoto University was used in Collaborative Research
Project for Young, Women Scientists, Institute for Information Management andCommunication, Kyoto University.
The author also used the supercomputers SGI altix ICE 8400EX and SGI ICEXA/UV hybrid system at the Institute for Solid State Physics of The University ofTokyo and the supercomputers at the Research Center for Computational Science,87

88
Okazaki.

CHAPTER 5. SUMMARY AND FUTURE OUTLOOKBibliography[1] R. M. Martin, “Electronic Structure: Basic Theory and Practical Methods,”
Cambridge University Press, Cambridge, 2004.
[2] http://www.gaussian.com/[3] http://www.vasp.at/[4] H. Imachi and T. Hoshi, “Hybrid numerical solvers for massively parallel
eigenvalue computation and their benchmark with electronic structure calculations,” J. Inf. Process. 24, pp. 164–172, 2016.
[5] H. Imachi, S. Yokoyama, T. Kaji, Y. Abe, T. Tada, and T. Hoshi, “Onehundred-nm-scale electronic structure and transport calculations of organic
polymers on the K computer,” AIP Conf. Proc. 1790, 020010, 4pp., 2016.
[6] T. Hoshi, H. Imachi, K. Kumahata, M. Terai, K. Miyamoto, K. Minami, andF. Shoji, “Extremely scalable algorithm for 108 -atom quantum material simulation on the full system of the K computer,” In Proceedings of ScalA16:
Workshop on Latest Advances in Scalable Algorithms for Large-Scale Systems, held as a part of SC16, Nov. 13, 2016, Salt Lake City, Utah, USA.
[7] http://www.elses.jp/[8] http://openmp.org/[9] http://mpi-forum.org/[10] W. Gropp, E. Lusk, N. Doss, and A. Skjellum, “A high-performance,
portable implementation of the MPI message passing interface standard,”
Parallel Comput. 22, 6, pp. 789–828, 1996.
[11] A, J. Van der Steen, “Overview of recent supercomputers,” Publication of
the NCF, 2008.
89

90

BIBLIOGRAPHY[12] G. M. Amdahl, “Validity of the single processor approach to achieving large
scale computing capabilities,” In Proceedings of the April 18-20, 1967,
spring joint computer conference (AFIPS ’67 Spring), 1967.
[13] L. N. Trefethen and D. Bau, Numerical Linear Algebra, Society for Industrial and Applied Mathematics, Philadelphia, 1997.
[14] G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed., Johns Hopkins University Press, Baltimore, 2013.
[15] Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd ed., Society for
Industrial and Applied Mathematics, Philadelphia, 2003.
[16] L. S. Blackford, J. Demmel, J. Dongarra, I. Duﬀ, S. Hammarling, G. Henry,
M. Heroux, L. Kaufman, A. Lumsdaine, A. Petitet, R. Pozo, K. Remington,
and R. C. Whaley, “An Updated Set of Basic Linear Algebra Subprograms(BLAS),” ACM Trans. Math. Soft., Vol. 28, Issue 2, pp. 135–151, 2002.
[17] E. Anderson, Z. Bai, C. Bischof, L. S. Blackford, J. Demmel, J. Dongarra, J.
Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen,“LAPACK Users’ Guide,” Society for Industrial and Applied Mathematics,
Philadelphia, 1999.
[18] L. S. Blackford, J. Choi, A. Cleary, E. D’Azevedo, J. Demmel, I. Dhillon,
J. Dongarra, S. Hammarling, G. Henry, A. Petitet, K. Stanley, D. Walker,
and R. C. Whaley, “ScaLAPACK Users’ Guide,” Society for Industrial andApplied Mathematics, Philadelphia, 1997.
[19] http://www.netlib.org/scalapack/[20] W. Kohn, “Nobel Lecture: Electronic structure of matter - wave functions
and density functionals,” Rev. Mod. Phys. 71, 1253, 1999.
[21] T. Hoshi, S. Yamamoto, T. Fujiwara, T. Sogabe, and S.-L. Zhang, “Anorder-N electronic structure theory with generalized eigenvalue equations
and its application to a ten-million-atom system,” J. Phys. Condens. Matter
21, 165502, 2012.
[22] T. Hoshi, K. Yamazaki, and Y. Akiyama, “Novel Linear Algebraic Theory
and One-Hundred-Million-Atom Electronic Structure Calculation on The K
Computer,” JPS Conf. Proc. 1, 016004, 2014.

BIBLIOGRAPHY

91[23] T. Hoshi, T. Sogabe, T. Miyata, D. Lee, S. L. Zhang, H. Imachi, Y. Kawai,
Y. Akiyama, K. Yamazaki, and S. Yokoyama, “Novel linear algebraic theory and one-hundred-million-atom quantum material simulations on the K
computer,” PoS(IWCSE2013) 065, 13pp, 2014.
[24] R. J. Bell and P. Dean, “Atomic vibrations in vitreous silica,” Disc. Faraday
Soc. 50, pp. 55–61, 1970.
[25] R. J. Bell, “The dynamics of disordered lattices,” Rep. Prog. Phys. 35, 1315,
1972.
[26] D. J. Thouless, “Electron in disordered systems and the theory of localization,” Phys. Rep. 13, 93, 1974.
[27] F. Wegner, “Inverse Participation Ratio in 2 + ε Dimensions,” Z. Physik B
36, 209, 1980.
[28] B. L. Anderson and R. L. Anderson, “Fundamentals of Semiconductor Devices,” McGraw-Hill, New York, 2005.
[29] EigenKernel; https://github.com/eigenkernel/eigenkernel/[30] http://elpa.rzg.mpg.de/[31] A. Marek, V. Blum, R. Johanni, V. Havu, B. Lang, T. Auckenthaler, A. Heinecke, H. J. Bungartz, and H. Lederer, “The ELPA Library – Scalable Parallel Eigenvalue Solutions for Electronic Structure Theory and Computational
Science,” J. Phys. Condens. Matter 26, 213201, 2014.
[32] T. Auckenthaler, V. Blum, H. J. Bungartz, T. Huckle, R. Johanni, L. Krämer,
B. Lang, H. Lederer, P. R. Willems, “Parallel solution of partial symmetriceigenvalue problems from electronic structure calculations,” Parallel Computing, Vol. 37, Issue 12, pp. 783–794, 2011.
[33] http://www.aics.riken.jp/labs/lpnctrt/index_e.html[34] T. Imamura, S. Yamada, and M. Yoshida, “Development of a highperformance eigensolver on a peta-scale next-generation supercomputer system,” Prog. Nucl. Sci. Technol., Vol. 2, pp. 643–650, 2011.
[35] T, Imamura, Y. Hirota, T. Fukaya, S. Yamada, and M. Machida,“EigenExa: high performance dense eigensolver, present and future,” 8thInternational Workshop on Parallel Matrix Algorithms and Applications(PMAA14), Lugano, Switzerland, 2014; http://www.aics.riken.jp/
labs/lpnctrt/index_e.html

92

BIBLIOGRAPHY[36] T. Imamura, “The EigenExa Library – High Performance & Scalable Direct
Eigensolver for Large-Scale Computational Science,” ISC 2014, Leipzig,
Germany, 2014.
[37] http://www.elses.jp/matrix/[38] T. Hoshi, Y. Akiyama, T. Tanaka, and T. Ohno, “Ten-million-atom electronic
structure calculations on the K computer with a massively parallel order-N
theory,” J. Phys. Soc. Jpn. 82, 023710, 4pp, 2013.
[39] G. Calzaferri and R. Rytz, “The Band Structure of Diamond,” J. Phys. Chem.
100, pp. 11122–11124, 1996.
[40] J. Cerdá and F. Soria, “Accurate and transferable extended Hckel-type tightbinding parameters,” Phys. Rev. B 61, pp. 7965–7971, 2000.
[41] http://www.aics.riken.jp/labs/lpnctrt/KMATH_EIGEN_GEV_e.
html[42] F. Tisseur and J. Dongarra, “Parallelizing the divide and conquer algorithm
for the symmetric tridiagonal eigenvalue problem on distributed memory architectures,” SIAM J. Sci. Comput., Vol. 20, Issue 6, pp. 2223–2236, 1999.
[43] J. Poulson, R. van de Geijn, and J. Bennighof, “Parallel algorithms for reducing the generalized hermitian-definite eigenvalue problem,” FLAME Working Note #56, The University of Texas at Austin, Department of Computer
Science, Tech. Rep. TR-11-05, 2011.
[44] M. P. Sears, K. Stanley, and G. Henry, “Application of a High Performance
Parallel Eigensolver to Electronic Structure Calculations,” In Proceedings ofthe 1998 ACM/IEEE conference on Supercomputing (SC ’98), Orlando, FL,
1998.
[45] RIKEN, Press Release at 5. Dec. 2013 (in Japanese); available from http://www.riken.jp/pr/press/2013/20131205_1/[46] Y. Takahashi, Y. Hirota, and Y. Yamamoto, “Performance of the block Jacobi
method for the symmetric eigenvalue problem on a modern massively parallel computer,” In Proceedings of ALGORITMY 2012, pp. 151–160, 2012.
[47] Y. Yamamoto, L. Zhang, and S. Kudo, “Convergence analysis of the parallel classical block Jacobi method for the symmetric eigenvalue problem,”
JSIAM Letters, Vol. 6, pp. 57–60, 2014.
[48] http://math.nist.gov/MatrixMarket/

BIBLIOGRAPHY

93[49] M. Lefenfeld, G. Blanchet, and J. A. Rogers, “High-performance contacts inplastic transistors and logic gates that use printed electrodes of DNNSAPANI doped with single-walled carbon nanotubes,” Adv. Mater. 14, pp.
1188–1191, 2003.
[50] M. C. Choi, Y. Kim, and C. S. Ha, “Polymers for flexible displays: frommaterial selection to device applications,” Prog. Polym. Sci. 33, pp. 581–
630, 2008.
[51] Sony Corporation, “Sony develops a ‘Rollable’ OTFT-driven OLED displaythat can wrap around a pencil,” News Release, 26. May. 2010; http://www.
sony.net/SonyInfo/News/Press/201005/10-070E/index.html[52] H. Bässler, “Localized states and electronic transport in single component
organic solids with diagonal disorder,” Phys. Stat. Solidi B 107, pp. 9–54,
1981.
[53] S. H. Chen, H. L. Chou, A. C. Su, and S. A. Chen, “Molecular packing in crystalline poly(9,9-di-n-octyl-2,7-fluorene),” Macromolecules 37, pp.
6833–6838, 2004.
[54] R. Noriega, J. Rivnay, K. Vandewal, F. P. V. Koch, N. Stingelin, P. Smith,
M. F. Toney, and A. Salleo, “A general relationship between disorder, aggregation and charge transport in conjugated polymers,” Nature Mater. 12, pp.
1038–1044, 2013.
[55] J. Terao, A. Wadahama, A. Matono, T. Tada, S. Watanabe, S. Seki, T. Fujihara, and Y. Tsuji, “Design principle for increasing charge mobility of
π-conjugated polymers using regularly localized molecular orbitals,” Nat.
Commun. 4, 1691, 2013.
[56] R. Takayama, T. Hoshi, and T. Fujiwara, “Krylov subspace method formolecular dynamics simulation based on large-scale electronic structure theory,” J. Phys. Soc. Jpn. 73, pp. 1519–1524, 2004.
[57] R. Takayama, T. Hoshi, T. Sogabe, S. L. Zhang, and T. Fujiwara, “Linear
algebraic calculation of the Green’s function for large-scale electronic structure theory,” Phys. Rev. B 73, 165108, 2006.
[58] H. Teng, T. Fujiwara, T. Hoshi, T. Sogabe, S. L. Zhang, and S. Yamamoto,“Eﬃcient and accurate linear algebraic methods for large-scale electronic
structure calculations with nonorthogonal atomic orbitals,” Phys. Rev. B 83,
165103, 2011.

94

BIBLIOGRAPHY[59] T. Hoshi, Y. Akiyama, T. Tanaka, and T. Ohno, “Ten-million-atom electronic
structure calculations on the K computer with a massively parallel order-N
theory,” J. Phys. Soc. Jpn. 82, 023710, 2013.
[60] T. Hoshi, K. Yamazaki, and Y. Akiyama, “Novel linear algebraic theory andone-hundred-million-atom electronic structure calculation on the K computer,” JPS Conf. Proc. 1, 016004, 2014.
[61] S. Nishino, T. Fujiwara, H. Yamasaki, S. Yamamoto, and T. Hoshi, “Electronic structure calculations and quantum molecular dynamics simulationsof the ionic liquid PP13-TFSI,” Solid State Ionics 225, pp. 22–25, 2012.
[62] S. Nishino, T. Fujiwara, and H. Yamasaki, “Nanosecond quantum molecular
dynamics simulations of the lithium superionic conductor Li4−x Ge1−x P x S4 ,”Phys. Rev. B 90, 024303, 2014.
[63] J. M. Soler, E. Artacho, J. D. Gale, A. Garcı́a, J. Junquera, P. Ordejön, and D.
Sänchez-Portal, “The SIESTA method for ab initio order-N materials simulation,” J. Phys. Condens. Matter 14, 2745, 2002.
[64] C. K. Skylaris, P. D. Haynes, A. A. Mostofi, and M. C. Payne, “Introducing
ONETEP: Linear-scaling density functional simulations on parallel computers,” J. Chem. Phys. 122, 084119, 2005.
[65] T. Ozaki, “O(N) Krylov-subspace method for large-scale ab initio electronic
structure calculations,” Phys. Rev. B 74, 245101, 2006.
[66] M. J. Gillan, D. R. Bowler, A. S. Torralba, and T. Miyazaki, “Order-N firstprinciples calculations with the conquest code,” Comp. Phys. Commun. 177,
14, 2007.
[67] W. Kohn, “Density functional and density matrix method scaling linearlywith the number of atoms,” Phys. Rev. Lett. 76, pp. 3168–3171, 1996.
[68] A. Frommer, “BiCGStab(l) for Families of Shifted Linear Systems,” Computing 70, pp. 87–109, 2003.
[69] S. Yamamoto, T. Sogabe, T. Hoshi, S. L. Zhang, and T. Fujiwara, “Shifted
Conjugate-Orthogonal-Conjugate-Gradient Method and Its Application to
Double Orbital Extended Hubbard Model,” J. Phys. Soc. Jpn., 77, 114713,
2008.
[70] T. Mizusaki, K Kaneko, M. Honma, and T. Sakurai, “Filter diagonalizationof shell-model calculations,” Phys. Rev. C 82, 024310, 2010.

BIBLIOGRAPHY

95[71] F. Giustino, M. L. Cohen, and S. G. Louie, “GW method with the selfconsistent Sternheimer equation,” Phys. Rev. E 81, 115105, 2010.
[72] S. Iwase, T. Hoshi, and T. Ono, “Numerical solver for first-principles transport calculation based on real-space finite-diﬀerence method,” Phys. Rev. E
91, 063305, 2015.
[73] K. Yamamoto, A. Uno, H. Murai, T. Tsukamoto, F. Shoji, S. Matsui, R. Sekizawa, F. Sueyasu, H. Uchiyama, M. Okamoto, N. Ohgushi, K. Takashina, D.
Wakabayashi, Y. Taguchi, and M. Yokokawa, “The K computer Operations:
Experiences and Statistics,” In Proceedings of the International Conference
on Computational Science (ICCS2014), pp. 576–585, 2014.
[74] TOP500 list; http://www.top500.org/[75] The Graph 500 List; http://www.graph500.org/[76] HPCG Benchmark; http://www.hpcg-benchmark.org/index.html[77] T. Adachi, N. Shida, K. Miura, S. Sumimoto, A. Uno, M, Kurokawa, F.
Shoji, and M. Yokokawa, “The design of ultra scalable MPI collective communication on the K computer,” Comput. Sci. Res. Dev. 28, pp. 147–155,
2013.
[78] B. Hess, C. Kutzner, D. van der Spoel, and E. Lindahl, “GROMACS
4: Algorithms for Highly Eﬃcient, Load-Balanced, and Scalable Molecular Simulation,” J. Chem. Theory Comput., 4 (3), pp. 435-447, 2008.;
http://www.gromacs.org/[79] J. Wang, W. Wang, P. A. Kollman, and D. A. Case, “Automatic atom type andbond type perception in molecular mechanical calculations,” J. Mol. Graph.
Model. 25, 247260, 2006.
[80] J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman, and D. A. Case, “Development and testing of a general AMBER force field,” J. Comp. Chem. 25,
pp. 1157–1174, 2004.
[81] A. Troisi and G. Orlandi, “Charge-transport regime of crystalline organic
semiconductors: diﬀusion limited by thermal oﬀ-diagonal electronic disorder,” Phys. Rev. Lett. 96, 086601, 2006.
[82] R. S. Mulliken, “Electronic Population Analysis on LCAO-MO Molecular
Wave Functions. I,” J. Chem. Phys. 23, pp. 1833–1840, 1955.

96

BIBLIOGRAPHY[83] F. Ortmann, F. Bechstedt, and W. G. Schmidt, “Semiempirical van der Waals
correction to the density functional description of solids and molecular structures,” Phys. Rev. B 73, 205101, 2006.
[84] J. MacQueen, “Some methods for classification and analysis of multivariate
observations,” In Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability, Volume 1: Statistics, pp. 281–297, 1967.
[85] C. M. Bishop, Pattern Recognition and Machine Learning, Springer-Verlag
New York, Inc., Secaucus, 2006.
[86] R. Dronskowski and P. E. Blöchl, “Crystal orbital Hamilton populations(COHP): energy-resolved visualization of chemical bonding in solids basedon density-functional calculations,” J. Phys. Chem. 33, 97, 1993.; http://www.cohp.de/[87] Y. Miyatake and J. C. Butcher, “A characterization of energy-preserving
methods and the construction of parallel integrators for Hamiltonian systems,” SIAM J. Numer. Anal., 54, 1993, 2016.