GPU-accelerated large-scale excited-state
simulation based on divide-and-conquer
time-dependent density-functional tightbinding
Takeshi Yoshikawa
RISE Waseda University
8th ADACTokyo, Japan, October 30, 20191

Introduction• Waseda University
PhD, Department of Chemistry andBiochemistry, Nakai group (March, 2015)• April 2015- March 2016
Research Associate, Waseda University• April 2016- March 2018
Assistant Professor, Waseda University• April 2018- Present
Lecturer, Waseda University
Linear-scaling open-shell and excited-state theory
2

Contents• DC-DFTB-MD– Divide-and-conquer density-functional tightbinding molecular dynamics– Ground state, massively parallel• DC-TDDFTB-MD– Divide-and-conquer time-dependent densityfunctional tight-binding molecular dynamics– Excited state• Acceleration of DC-TDDFTB-MD withGPU
3

Development of DC-DFTB-MD program: DC-DFTB-MD•MD and geometry optimization for large systems•Written in Fortran 90/95•K computer, Fujitsu FX10, PC cluster (Linux), ...

Academic free

Main features
Energy and Gradient(Cluster and PBC models)

DFTB1/2/3,
Spin-polarized DFTB

Linear-scaling method

Divide-and-conquer (DC)Dispersion correction

Slater-Kirkwood, Lennard-Jones, DFT-D2, DFT-D3

SCC convergence

Broyden, Simple mixing, Anderson, DIISMD ensemble

NVE, NVT, NPH, NPT

Thermostat and [Barostat]

Velocity scaling, Nosé-Hoover chain, Berendsen,
Andersen, Langevin, [Berendsen]

MD option

RATTLE constraint, Soft potential, Simulated annealing,
Lagrange interpolation of initial guess charges

Geometry optimization

BFGS, Steepest descent, Conjugate gradient

http://www.chem.waseda.ac.jp/dcdftbmd/

4

Target: Chemical reaction dynamics of large Systems
Our recent illustrative applications
Chemical absorption of CO2
Grotthuss-type mechanism

Ion-pair mechanism

CO2 absorption

CO2 regeneration

Lithium ion devices

300 K
320 K
360 K
400 K

20

Potential [V]# of H+ hopping

Proton diffusion in bulk water

10
0

0

4

8

time [ps]

12

−1.3
−1.2
−1.1
−1.0
−0.9
−0.8

5.0 Å

Z

Quantum mechanical treatment of whole system with inexpensive computational cost5Our three pillars for chemical reaction dynamics

DC

DFTB

DivideandConquerO(N)DensityFunctional
TightBinding
SEQM-MD

Massively parallel
DC-DFTB-MD
computations

K computer•O(N): Linear-scaling•SEQM-MD: Semi-Empirical
Quantum MechanicsMolecular Dynamics

6Divide-and-Conquer (DC) method
One of fragmentation techniques to achieve linear-scaling[1,2]…

Subsystem
Buffer region

Divide…

…

Localization
regionSubsystem fragmentation
without overlap

Solve subsystems’ equations•Reduce diagonalization cost: O(N)
•Obtain subsystems’ orbitals…

Set common Fermi level (Conserve total number of electrons)
Conquer…
SCF•Total density matrix•Total energy•Other properties[1] M. Kobayashi, H. Nakai, in Linear-Scaling Techniques in Computational Chemistry and Physics, (2011), pp. 97–127.
[2] W. Yang, T.-S. Lee, J. Chem. Phys. 103, 5674 (1995).

7

Advantages of DC method in QM-MD simulations•Automatized fragmentation of subsystems at each time stept  t2t  t1

Fragmentation due to localityof density matrix
Accuracy and cost control bychanging buffer sizeNo need of chemical intuition

Applicable to dynamical bond formation and cleavage•Unnecessary to predefine electron & spin numbers of subsystems(Determine at each SCF step)
Conventional fragmentation method
S1  0 S 2  1 2 S3  0 S1  ? S 2  ? S3  ?

DC methodS 1 2S 1 2?
Applicable to electron and spin delocalized systems

8

Density-Functional Tight-Binding (DFTB) method
SEQM derived from density functional theory (DFT)[1,2]
DFTB1

DFTB2

DFTB3

•Dμν: Density matrix
1 atom
1 atom2
EDFTB   D H    E    AB q AqB   AB q AqB •ΔqA: Induced Mulliken
2 AB
3 AB

A B
charge on atom A
①
②
③
④
AO

atom

0

rep
AB

Evaluate using parameters determined by DFT calculations①: Charge independent termPrecomputed Hamiltonian and overlap matricesTwo-center approximation②: Short-ranged two-body repulsive termCore-core repulsion
DFT double-counting contribution③, ④: Charge dependent term
ΔqA is determined self-consistently (SCC)Monopole approximation[1] M. Elstner, G. Seifert, Phil. Trans. R. Soc. A 372, 20120483 (2014).
[2] M. Gaus, Q. Cui, M. Elstner, WIREs Comput. Mol. Sci. 4, 49 (2014).

No integral evaluation
at runtime
Key ingredient forreasonable accuracy
and transferability
Proper description ofchemical bonds in MDwith small cost9

Parallel procedure of DC-DFTB energy calculation
Subsystem α

Buffer•Hybrid of MPI/OpenMP…
…

MPI: Assign subsystems to processes
OpenMP: Operations for localization regions
DFTB1

DFTB3

DFTB2
atom

Evaluate H0,α, Sα, γα, Γα, Erep1
 AC   BC qC
2
C

0
H   H 
 S  

atom
q
1

 S    q AAC  qB BC   C CA  CA qC
6

C 3

Update Hamiltonian HαDiagonalize to get cα & εαMPI comm.

Set common Fermi level εF Solve equations:  ci H    i S   0Constraint of total number of electrons ne:
subsy stem MO( )ne  2

MPI comm.

Construct density DαSCCf       p c c S




i

i

i

i




D

MO( )

 2  f   F   i ci ci *i

Evaluate Mulliken charge qα



F



qA    DS A L ( )fβ(x): Fermi functionf  x  11  exp(  x)pα: partition matrix

MPI comm.

Evaluate Mulliken charge q

SCC Convergence check: max q  qold  

10Parallel performance of DC-DFTB energy calculation
Size dependence of computational time
Computational time [sec]

300.0

200.0

(H2O)nn = 4,000–500,000# of cores
8,000
16,00032,00064,000128,000

H2O × 500,000

100.0

0.0
0.0

O(N1.4)108000256000 364500

24.66 nm

500000

Number of H2O

Computational time of (H2O)500000# of cores

16,00032,00064,000128,000

Time [sec]

193.4

107.0

61.1

46.5

Speed-up

16000.0

28908.7

50673.972935.4•Density: 1.0 g/cm3•DC-DFTB2•Subsystem: 1 H2O•Buffer radius：6 Å•K computer

11

Contents• DC-DFTB-MD– Divide-and-conquer density-functional tightbinding molecular dynamics– Ground state, massively parallel• DC-TDDFTB-MD– Divide-and-conquer time-dependent densityfunctional tight-binding molecular dynamics– Excited state• Acceleration of DC-TDDFTB-MD withGPU
12

DC-based excited-state theory
DC-based time-dependent density-functional tight-binding theory
Ground-state(DC-DFTB)
Select EC region
Excited statefor EC region(DC-TDDFTB)

Excitation center (EC)in effect outside EC

Casida’s formalism
A
 s*
B
s

B  X 
1 0  X    s 
s*  s 
A  Y  0 1  Y sss[1] T. Yoshikawa, M. Kobayashi, A. Fujii, and H. Nakai,
J. Phys. Chem. B, 117, 5565 (2013).
[2] N. Komoto, T. Yoshikawa, J. Ono, Y. Nishimura,and H. Nakai, J. Chem. Theor. Comput., 15, 1719(2019).

13

Time-dependent DFT(B) energy/energy gradient• Excitation energy: linear-response theory– Hermit eigenvalue problem: Davidson algorithm A  B

 A  B  A  B  cI  I2cI
Aia, jb  ijab  a  i   Kia, jb Bia , jb  K ia ,bj
1/2

1/2

 : orbital energyi,j: occupied orbital

a,b: unoccupied orbital• Excitation energy gradient0H K 
S

1

D     WR x  R x
2  R x
 R x– Z-vector method• Evaluate one-body, energy-weighted, and two-body density matrices• Avoid the explicit evaluation of MO coefficient derivative for 3×Natom14

TDDFTB energy: Davidson algorithm• Diagonalization for sparse matrix– Hermit eigenvalue problem0Ωc I   I2c I0Ω   A  B

1/2

c †I

 A  B  A  B 

Large matrix

1/2

cI

Small matrix

Diagonalizationc I , I

15

TDDFTB energy: Davidson algorithm• Coupling matrix K (for example DFTB2)
Atomμ,ν: AO
α,β: Atom

Kia , jb   qia  qjb
 

• Mulliken atomic transition charges
1q     Ci S  C  C i S C 
2  ia•

Construction of Ω   A  B 

1/2

C: MO coefficent
S: overlap matrix A  B  A  B 

1/2– Standard: AO-MO transformation for each Davidson iterationia, I  ia2 Fia  4  ia Ci C K    jb C j C bc jb, Iia

jb 

– In-core: store q in memory   F  4 q    
ia , I

2
ia

ia
A

ia

A

ia

 B



jb

q
F
AB  
jb A
jb  

 jb


16

Algorithm of (DC-)TDDFTB
Step 1: Set initial trial vector
Sort in ascending order for difference of orbital energies
ia   a   i  : orbital energy i,j: occupied orbital a,b: unoccupied orbital
Step 2: Construct hermit matrix ΩStandard
Step 2-1: MO-AO transformation
for trail vector
Step 2-2: Construct AO based Ω
Step 2-3: AO-MO transformation Ω

In-coreStep 2-1: Construct Mulliken
atomic transition charges
Step 2-1: Construct Ω directlyfrom trial vector

Step 3: Diagonalize hermit matrix Ω
SCF

Step 4: Update trial vector

SCF
17

TDDFTB energy gradient: Z-vector method• Lagrangian of excitation energy in TDDFTBL  X  Y  A  B X  Y  X Y  A  B X Y

   X  Y X  Y  1   Z ia H ia   W pqTD ( S p q   p q )pq , p  qia• Z-Vector equation

L
0
C p

 ( A  B)

Linear equationsRia   ( X  Y )ib H ab[( X  Y)]  H ia [T]

Z jb   Ria

jb

Stationary conditionbia , jb  ( X  Y ) ja H [( X  Y)] .

ji


H pq
[V ]   K pq ,ruVru
ruj• One-body, two-body, and energy weighted DM by Z-vector0H K 
S

1

D     WR x  R x
2  R x
 R x
18

TDDFTB energy gradient: Z-vector method• Preconditioned conjugate gradient (PCG) method– Solve iteratively symmetric and positive linear equations

Step 1: Construct Z-vector
Step 1-1: Set initial Z vector (Jaccobi preconditioner) Z ia   a   i
Step 1-2: Construct residual vector0

riak 1  riak   k   A  B ia , jb diakjb

SCF

Step 1-3: Update Z vector

Z iak 1  Z iak   k diak

Step 2: Construct density matrix

Step 3: Evaluate gradient19Performance of (DC-)TDDFTB method• Wall-clock time [second]– Molecular: p-cumaric acid + 228H2O– Method: TDDFTB2/3ob
EnergyEnergy gradient

TDDFTB DC-TDDFTBTDDFTB DC-TDDFTB

Step 1

0.367

0.137

Step 1

19.565

4.222

Step 2

1.943

0.293

Step 2

0.200

0.066

Step 3

0.000

0.000

Step 3

0.240

0.091

Step 4

0.127

0.037

DC calculations reduce the computational cost,especially for largest computational step in each algorithm.
CPU: Intel Xeon Gold 6136 (3.00 GHz) 12 cores

20

System-size dependency: total TDDFT calculation• CPU time and required memory– Molecular: p-cumaric acid + nH2O
– Method: TDDFTB2/3ob

In core algorithm reduce the wall times by large memory region
DC calculations reduce the CPU time and required memory
CPU: Intel Xeon Gold 6136 (3.00 GHz) 12 cores

21

Contents• DC-DFTB-MD– Divide-and-conquer density-functional tightbinding molecular dynamics– Ground state, massively parallel• DC-TDDFTB-MD– Divide-and-conquer time-dependent densityfunctional tight-binding molecular dynamics– Excited state• Acceleration of DC-TDDFTB-MD withGPU
22

Implementation with GPU: TDDFT energy
Overlap

Small matrix
Small matrix

Low efficiency
on GPUIn-core algorithm because of implementation with dgemm
ia , I   Fia  4 q
2
ia

A

ia
A




jb
ia    AB    jb qA Fjb  
 jb

 B

23

Implementation with GPU: TDDFT energy gradientAcceleration of most procedure with cublas library
CublasDGEMM, CublasDAXPY, Cublasddot
DC method: small require memory and communications
with CPU and GPU

24

Performance of (DC-)TDDFTB method• Wall-clock time [second]:
TDDFTB energy– Molecular: p-cumaric acid + 228H2O– Method: TDDFTB2/3ob
TDDFTBCPU
Step 1

0.367

TDDFTBGPU
0.013

DC-TDDFTBCPU

DC-TDDFTB
GPU

0.137

0.005

Step 2

1.943

0.075

0.293

0.010

Step 3

0.000

0.002

0.000

0.002

Step 4

0.127

0.003

0.037

0.002

Step 1, 2, and 4: computational cost ↓Step 3: computational cost ↑CPU: Intel Xeon Gold 6136 (3.00 GHz) 12 cores, GPU: V100

25

Performance of (DC-)TDDFTB method• Wall-clock time [second]:
TDDFTB energy gradient– Molecular: p-cumaric acid + 228H2O– Method: TDDFTB2/3ob
TDDFTBCPU
Step 1

19.565TDDFTBGPU
0.630

DC-TDDFTB
CPU
4.222

DC-TDDFTB
GPU
0.066

Step 2

0.200

0.023

0.066

0.009

Step 3

0.240

0.014

0.091

0.005

All steps: computational cost ↓Step 4: Required communication cost ↑ but small especially DC
CPU: Intel Xeon Gold 6136 (3.00 GHz) 12 cores, GPU: V100

26

System-size dependency: total TDDFT calculation• Wall-clock time and required memory– Molecular: p-cumaric acid + nH2O
– Method: TDDFTB2/3ob

DC method: required memory small GPU good!!
Multi-GPU: about 3.12 times with 4 GPUs on same bord
CPU: Intel Xeon Gold 6136 (3.00 GHz) 12 cores, GPU: V100

27

Application: photoactive yellow protein
Photoisomarization
p-cumaric acid(pCA)• Photoreaction cycle

hν
bR
pG

O

bR hν

pG*
K

N
pB

L
pR

M2

M1

O transN

K cis

Hydrogen bondfor around amino acidInLPYPIn water

284 nm

446 nm

M2 redshift
M1
pCA

Glu46

Tyr42

28

Application: photoactive yellow protein
MethodExcitation energy [eV]

CPU time

pCA+461H2O(diff.)

PYP(diff.)

DC-CIS

4.64(0.64)

4.22(1.44)~30 min.

DC-TDBLYP

3.76 (-0.24)

2.33(-0.45)~3 hours

DC-TDLCBLYP

4.10 (-0.27)

2.85(0.07)~3 hours

DC-TDDFTB

3.68 (-0.32)

2.17(-0.61)

1.27 sec.

DC-TDLCDFTB

4.07 (-0.30)

2.80(0.02)

2.52 sec.

Exptl.

4.00

2.78

―

―

―

TDBLYP: 6-31G**,TDDFTB/3ob, Excited center: pCA, Subsystem: 1H2O,
1 amino acid, pCA, Buffer size: 5.0 Å (pCA+461H2O), 6.0Å(for PYP),
16 cores of an Intel Xeon Gold 6134 (3.20 GHz)

JCTC, 15,
TDDFTB method could reproduce the excitation energy
of 2019(front cover)
TDDFT calculations within seconds in single point calculations.
29

Application: excited-state single proton transfer
Excited state proton transferof 2-acetylindan-1,3-dione (AID)

S1enolGround state~160 fs

S1ketoS0ketoS0enol

AID-keto
AID-enol 1 ps DC-DFTB-MD simulationin ground state Solvent: acetonitrile
 NVE ensembles
 3ob parameter set Buffer size: 5.0 Å

30

Application: excited-state single proton transfer
Excited state proton transferof 2-acetylindan-1,3-dione (AID)

S1enol

S0enolExcited state~160 fs

S1ketoS0ketoAID-keto
AID-enol1 ps simulations
CPU: ~11 hours GPU:~20 mins

 1 ps DC-DFTB-MD simulationin ground state Solvent: acetonitrile
 NVE ensembles
 3ob parameter set Buffer size: 5.0 Å

31

Application: excited-state single proton transfer
Excited state proton transferof 2-acetylindan-1,3-dione (AID)
Bond length [Å]~160 fs

S1ketoS0ketoS0enolr(O2‐H8)r(O1‐H8)
time [fs]Excitation energy

AID-keto
AID-enol[1] P. K. Verma, A. Steinbacher, A. Schmiedel, P. Nuernberger,
T. Brixner, Struct. Dyn., 3, 023606 (2015).

Excitation energy [eV]S1enolBond length between O-H atoms

Absorbance: 4.01 eV [1]Fluorescence: 2.32 eVtime [fs]32Application: excited-state single proton transfer
Excited state proton transferof 2-acetylindan-1,3-dione (AID)

0 .2

E n e rgy d iffe ren ce [h a r tree ]S1enolGround state~160 fs

S1ketoS0ketoS0enolT o taMl D en e rgy
K in e tic en e rgy
E lec tron ic en e rgy

0 .1

0

-0 .1

-0 .2
0

200

400

600

800

1000

800

1000

tmi e [fs ]Excited state

AID-keto
AID-enolE n e rgy d iffe ren ce [h a r tree ]0 .2

T o ta M
l D en e rgy
K in e tic en e rgy
E lec tron en e rgy

0 .1

0

-0 .1

-0 .2
0

[1] P. K. Verma, A. Steinbacher, A. Schmiedel, P. Nuernberger,
T. Brixner, Struct. Dyn., 3, 023606 (2015).

200

400

600

tm
i e [fs ]33

Excited-state intramolecular proton transfer (ESIPT)• Application for organic electroluminescence and biosensor• Large red-shift fluorescence spectrum

ESIPT

Excitedstate
Absorption

Fluorescence

Groundstate• 2,2’-Bipyridine-3,3’-diol (BP(OH)2)– Excited-sate intramolecular double proton transfer in polar solvent
O1

N2

H 1 N1H 2 O2

dienolmonoketo

diketo34

Excited-state intramolecular proton transfer (ESIPT)• Computational details
Model
BP(OH)2+ 185 acetonitrile molecules (droplet) after
equilibration

DC-(TD)LCDFTB metaDynamics










NVT run for 9 ps at 300 K in ground state (DC-LCDFB)
NVE run for 5 ps in excited state (DC-TDLCDFTB)
Δt = 0.5 fs
Parameter: OB2[1] Corrective variable: O-H distance
Subsystem: BP(OH)2, acetonitrile
 Gaussian potential height: 1.88 kcal・mol-1
Excitation canter: BP(OH)2
 Gaussian potential width: 0.025 Å
Buffer size: 5.0 Å Frequency of gaussian potential addition
Soft potential: r = 17.0 Å: 20 fs
Dispersion correction: Grimme’s DFT-D3  Bias temperature: 9000 K[1] V. Q. Vuong, J. A. Kuriappan, M. Kubillus, J. J. Kranz, T. Mast, T. A. Niehaus, S. Irle, and M. Elstner, J. Chem. Theory Comput., 14, 115 (2018).

35

Excited-state intramolecular proton transfer (ESIPT)• Free energy surface for O-H distances0

O-H distance [Å]O2-H2 distance [Å]

S0

0
-1

S1

-4

-2

-8

Free energy difference[kcal・mol-1]

-3
-12

3.0

0.0

-16

-5

-20

-6

O1-H1 distance [Å]O1-H1 distance [Å]Time series of O-H distances
系列1
O1-H1
系列2
O2-H2

2.0
1.0

-4

H2

2.0

O2

S1

0.0

0 1 2 3 4 5 6 7 8 9Time [ps]

H1

3.0

1.0

S0

O1

0

1

2
3
Time [ps]45

Proton transfer occurred easily in S1 state.

36

Summary• Development of Divide-and-conquer-based
time-dependent density-functional tightbinding theory for large excited-statesystems: DC-TDDFTB-MD• Compatibility of DC and GPU• Further work: Ground state for GPU
more large and long-time simulation
GPU

37

38