Towards Large Scale Quantum ChemistryCalculations on the K and Fugaku Supercomputer
William Dawson
Takahito Nakajima
RIKEN Center for Computational Science

ADAC8 Workshop, University of Tokyo, Kashiwa, Japan.
1Computational Molecular Science Research Team

• Also in collaboration with Luigi Genovese (French AlternativeEnergies and Atomic Energy Commission).

2

Theme: Chemical Properties from First Principles• The goal of our research is to develop theoretical methods,

numerical algorithms, and software to study the properties ofmolecules and materials from the basic laws of the universe (inparticular: quantum mechanics).

Atomic Positions

Electronic Structure
Software

System Properties• Electronic structure methods allow for a quantum mechanicaltreatment of systems, giving scientists tools to go beyond thetechnical limits of experiment in understanding the natural
world and designing new materials.

• A broad field with many diﬀerent methods and algorithms.
3

Example Application: Bioremediation• Bioremediation: the use of naturally

occurring organisms to break down
environmental pollutants.

• Example: using enzymes to degrade
toxins, which can contaminate foodsupplies.

• Requirement 1: computationalmethods and software that canperform accurate calculations on suchlarge systems (10,000 atoms).

• Requirement 2: methodologies thatcan extract scientific insight fromthese calculations.

Zaccaria, Marco, WD, Viviana Cristiglio, Massimo Reverberi, Laura E. Ratcliff, Takahito Nakajima, Luigi Genovese, and Babak
Momeni. "Designing a bioremediator: mechanistic models guide cellular and molecular specialization." Current opinion inbiotechnology 62 (2020): 98-105.

4

Requirement 1: Methods for Large Systems• One of the major computational bottlenecks for quantum chemistry
software is computing the density matrix D from the hamiltonian H.

• This is typically done by solving an eigenvalue problem.
• HV = 𝝀V (where V is a nbasis x nelectrons matrix).
• Then compute the density matrix as D = VVT.
• Computational cost scales with the third power of the systemsize.

• Instead, we can compute the density matrix directly using the fermi
function:β(H-μ) + I).
D=1/(e•• This function here is a matrix function. Is it possible to use thisframework to reduce the computational cost?

5Introduction to Matrix Functions• Simple definition: We are all familiar with functions of a single
variable f(x). In the matrix function case, just replace thevariable x with a matrix A.
Standard Function

Matrix Function

Interpretationf(x) = x2

f(A) = A2

Matrix Productf(x) = 1/x

f(A) = A-1

Matrix Inverse

f(x) = exf(A) = eA

Matrix Diﬀerential Equationf(x) = sign(x)f(A) = sign(A)Projection on to Subspace• Matrix functions have many diﬀerent applications including
solving diﬀerential equations and the study of complex
networks.

Higham, Nicholas J. Functions of matrices: theory and computation. Vol. 104. Siam, 2008.

6Motivating Matrices• In many domains, the problem of interest can be represented using asparse, hermitian matrix.

• Under certain conditions, not only is the matrix A sparse, but also somematrix functions f(A) are sparse.

• Estrada matrix exponential eβA contains a scaling factor which might beinterpreted as a unit of edge weight.

• For certain values of β, the matrix exponential of small world matrices isalso sparse.

102

Sparsity(%)

101

10=1= 120= 14= 181= 161= 32

29210

211

212
Nodes

213

214

7

Motivating Matrices - Chemistry• For insulating systems (and metals at high temperature), it is

known that the density matrix is sparse when represented in a
localized basis.

• Example: the Hamiltonian and Density Matrix of 1920 watermolecules computed using BigDFT.

Hamiltonian

Density Matrix

8

Sparsity Aware Matrix Function Calculation• From the list of methods for computing matrix functions, we will
select calculations based on matrix polynomials.

• e.g. Chebyshev polynomials: f(A) ≅ ∑ciTi(x).
T2(A) = 2A2 - I
T4(A) = 8A4 - 8A2 + I• T0(A) = I
3 - 3A
T(A)=AT(A)=4A
13•• Computing a matrix polynomial requires only two core routines:
matrix addition, matrix multiplication.

• Easy to parallelize.
• Many functions can be tuned through just two routines.
• In the case of sparse matrices, we replace these kernels withsparse matrix addition, and sparse matrix multiplication.

9Parallel Implementation: Topological Aware• Approach: a three dimensional data distribution maps well tothe high level 6D Tofu Network topology of K/Fugaku.

• Map a 2D matrix onto a 3D topology, and exploit the extra 3
dimensions for eﬃcient collective communication.

• For each dimension, we are free to choose whether to partition
or replicate data structures.

Proc 2,1,1

X-Dimension

Proc 1,2,1

Proc 2,2,1

Proc 1,1,2

Proc 2,1,2

Z-Dim

Proc 1,2,2

Proc 2,2,2

ensio

Y-Dimension

Proc 1,1,1

X-Dimension

Y-Dimension

Y-Dimension

X-Dimension

Proc 1,1,3

Proc 1,2,3

Proc 2,1,3

Proc 2,2,3n10

Parallel Implementation: Matrix Multiplication• Core routine of computing matrix powers.
• Overdecomposition of the matrix enables strong scaling performance.
• Gather+Reduce implemented as collective communication.
Y-Dimension

Partial Gather (X/Y)

Local MultiplyReduce (Z)(1,1,1) (1,2,1)(2,1,1) (2,2,1)

X-Dimension

Y-Dimension

Z-Dimension

X-Dimension= A2(1,1,2) (1,2,2)(2,1,2) (2,2,2)[1] Solomonik, Edgar, and James Demmel. "Communication-optimal parallel 2.5 D matrix multiplication and LU factorization algorithms." In European Conference on Parallel Processing,pp. 90-109. Springer, Berlin, Heidelberg, 2011.
[2] Dawson, William, and Takahito Nakajima. "Massively parallel sparse matrix function calculations with NTPoly." Computer Physics Communications 225 (2018): 154-165.

11

Parallel Implementation: Thread Parallel• OpenMP loop parallelism doesn’t work well with overlapping
communication. Instead we use OpenMP task framework.

• Creating a task manager, and dependency graph.
• Block levels tasks help mitigate load imbalance from sparsematrix and communication in diﬀerent directions.
Local Gather Row Ai

Transpose Column Bj

Transpose Column Bk

Global Gather Row Ai

Local Gather Column Bj

Local Gather Column Bk

Transpose Row Ai

Global Gather Column Bj

Global Gather Column Bk

Compute Block Cij

Compute Block Cik

Sum Z Direction Cij

Sum Z Direction Cik
12

NTPoly Features• General Polynomials• Standard Polynomials• Chebyshev Polynomials• Hermite Polynomials• Transcendental Functions• Trigonometric Functions• Exponential and Logarithm• Matrix Roots• Square Root and Inverse

Square Root
• Matrix pth Root and Inverse
pth root• Matrix Inverse (and MoorePenrose Inverse)• Sign Function/PolarDecomposition• Quantum Chemistry• Density Matrix Purification• Chemical PotentialCalculation• Density Matrix Extrapolation• Other• Parallel File I/O• MIT License• Fortran/C/C++/Python

Interface• Automatic data redistribution• Real and Complex Matrices

WD, and Takahito Nakajima. "Massively parallel sparse matrix function calculations with NTPoly." Computer Physics Communications 225 (2018): 154-165.

13

Performance: Quantum Chemistry• Standard eigensolvers can make limited use of the sparsity of amatrix, but will be outperform by matrix function basedapproaches.

• Calculation of water clusters of various sizes, 6-31G basis set,using the TRS2 density matrix method to approximate the fermifunction.

14

Going Beyond Quantum Chemistry• By building a general framework, we can apply NTPoly to many diﬀerent
fields.

• Example: computing the matrix exponential eβA to study social networks.

CaroleArthur
Patrick
Bob
KimPaul

Alice
Mary

AmandPeterJames

Estrada, Ernesto, Naomichi Hatano, and Michele Benzi. "The physics of communicability in complex networks."
Physics reports 514, no. 3 (2012): 89-119.

15

Going Beyond Quantum Chemistry• By building a general framework, we can apply NTPoly to many diﬀerentfields.

• Example: computing the matrix exponential eβA to study social networks.

CaroleArthur
Patrick
BobKimPaul

AliceMary

AmandPeterJames

Estrada, Ernesto, Naomichi Hatano, and Michele Benzi. "The physics of communicability in complex networks."
Physics reports 514, no. 3 (2012): 89-119.

16

Going Beyond Quantum Chemistry• By building a general framework, we can apply NTPoly to many diﬀerent
fields.

• Example: computing the matrix exponential eβA to study social networks.

CaroleArthur
Patrick
Bob
KimPaulAliceMaryAmandPeterJames

Estrada, Ernesto, Naomichi Hatano, and Michele Benzi. "The physics of communicability in complex networks."
Physics reports 514, no. 3 (2012): 89-119.

17

Going Beyond Quantum Chemistry• By building a general framework, we can apply NTPoly to many diﬀerent
fields.

• Example: computing the matrix exponential eβA to study social networks.

CaroleArthur
Patrick
Bob
KimPaulAliceMary

AmandPeterJames

Estrada, Ernesto, Naomichi Hatano, and Michele Benzi. "The physics of communicability in complex networks."
Physics reports 514, no. 3 (2012): 89-119.

18

Going Beyond Quantum Chemistry• By building a general framework, we can apply NTPoly to many diﬀerentfields.

• Example: computing the matrix exponential eβA to study social networks.

CaroleArthur
Patrick
Bob
KimPaul

Alice
Mary

AmandPeterJames

Estrada, Ernesto, Naomichi Hatano, and Michele Benzi. "The physics of communicability in complex networks."
Physics reports 514, no. 3 (2012): 89-119.

19

Going Beyond Quantum Chemistry• By building a general framework, we can apply NTPoly to many diﬀerentfields.

• Example: computing the matrix exponential eβA to study social networks.

CaroleArthur
Patrick
Bob
Kim

Paul

AliceMary

AmandPeterJames

Estrada, Ernesto, Naomichi Hatano, and Michele Benzi. "The physics of communicability in complex networks."
Physics reports 514, no. 3 (2012): 89-119.

20

Example: The Resilience of Social Networks• Network Resiliency calculations:• Data Set: Israeli Social Network “TheMarker Cafe”• Nodes: 69413. Sparsity: 0.04%.
• Procedure:• Remove a node from the graph• Random Node• Node with largest degree• Compute the matrix exponential• Compute the sparsity• Repeat50
45

Average Number of Edges

40
35
30
25
20
15Random
Largest Degree

100500

1000

1500
2000
Nodes Removed

2500

WD, and Takahito Nakajima. "Massively parallel sparse matrix function calculations with NTPoly." Computer Physics Communications 225 (2018): 154-165.

3000

21Conclusion• Electronic structure calculations enable us to predict theproperties of molecules using the laws of quantum mechanics.
• Using the world’s largest supercomputers to push quantum
mechanical methods into the 10,000 atom regime, enabling
calculations of previously inaccessible systems and properties.
• Use of matrix functions to enable large scale quantum chemistrycalculations.
• A highly scalable library that can be applied to many diﬀerent
disciplines.
• Complexity reduction approach to fragment large systems.

• NTPoly: https://github.com/william-dawson/NTPoly• BigDFT: https://gitlab.com/l_sim/bigdft-suite• PyBigDFT: https://bigdft-suite.readthedocs.io/projects/PyBigDFT/22