Acceleration of computational quantum chemistryby heterogeneous computer architectures

Ryota Koga
President of X-Ability Co.,Ltd.
# Collaborative researchers
Yuki Furukawa(X-Ability Co.,Ltd.), Koji Yasuda(Nagoya University)
14 Dec 2011
GTC Asia 2011, BeijingAbout the X-Ability Co.,Ltd.
• 3 Office (foundation : 15th Jan 2008)
Hongo, Tokyo / The K-comp, Kobe / U.S. branch• 3 Business(1) Scientific Computing (mainly Chemistry)developing XA-CHEM-SUITE(2) Sensor Networking(3) ChemBioinformatics (machine-learning etc.)• By only 4 Regular Board Members

Drugs using Computational Chemistry©wikipedia
Erlotinib hydrochloride©wikipedia
Zanamivir(trade name Tarceva)

trade name Relenza
Influenza preventive medicine

Lung cancer, pancreatic cancer
Steve Jobs might have used
this to extend his life…

Copyright (C) 2011 X-Ability Co.,Ltd.

3

Copyright (C) 2011 X-Ability Co.,Ltd.

4

Abstract

Benchmark : multicore CPUs

• Vectorized PRISM by AVX on Sandybridge CPU• HF (Hartree-Fock) : x 3 faster• DFT by CUDA on NVIDIA GPU• J-matrix formation : x 245• Exchange correlation : x 4
NVIDIA GPU• ERI J-matrix (Coulomb potential) by CUDA on NVIDIA GPU• FMO Environmental Electrostatic Potential : x 17
GAMESS-US(free application)

Replace HF, DFT, FMO ESP
with CUDA/AVX scratch code
Acceleration so muchGAMESS-US+ XA-AVX/CUDA-QM

CPU[AVX]/GPU[CUDA] hybrid programming• No GPU Program run without CPUs.
– CPU+GPU architecture is general.

• SSE/AVX should be appropriate if the process issuitable for SIMD parallelization and not for GPU.
• AVX computes 4 double precision arithmetics at once.
• AVX is theoretically twice faster than SSE.

– AVX is available on latest CPU (SandyBridge).
• A normal Desktop PC is used in this research.

AVX : Pros and Cons
Basically, it is like a vector processor.
Pros

Cons

- Twice faster than SSE- SIMD parallelization- utilizes main memorysame as FPU computing- Poor for algorithms thatuse conditional branching
and/or random memory
access- Automatic vectorization isefficient only for simple
casesPlan for implementation
AVXGPU- Use best tuned algorithm as far as- Define overloaded
GPU’s resource permits.
operators of C++ to avoid- Shared memory of Fermi core isusing AVX intrinsic functions
used as L1 cache.
explicitly.
- Use CPU threads (pthread) to- Use function pointers of
C/C++ to eliminate branches
in most inner loop.
-harness GPU calculations forcomplicated thread control(synchronization, queueing…) .
Use OpenMP for CPU threadparallelization.
Single/double mixed-precisioncalculation.

System Configuration• Core i5 2500 @ 3.30GHz• NVIDIA Geforce GTX580– 512 cuda cores• DDR3 8GB• CUDA 4.0• Intel Composer XE 2011

Copyright (C) 2011 X-Ability Co.,Ltd.

9

How to solve Schrӧdinger equation?
H (1 N )(1 N )  E(1 N )Expand wave function with“contracted Gaussian basis set”
(Linear Combination of Atomic Orbitalapproximation)

spzdyz i (r )   Ck(i )  k (r )kKs ：  (r )   d k ek 1

F (C )C  SC k ( r  A ) 2

Kpx：  (r )   d k ( x  Ax )e k ( r  A ) 2

k 1

Matrix form of eigenvalue problem(Self consistent)

Copyright (C) 2011 X-Ability Co.,Ltd.

10

Summary of Hartree-Fock Procedure
F (C )C  SC

：SCF(ab | cd )  

 a (r )  b (r )  c (r ' )  d (r ' )dr ' dr ：ERI| r  r '|F  H    Ca C*a 2      
2/ Ncorea



1


core
 H   D       2



core
 H 
 G
Density Matrix (D) : is updated at every SCF cycle to solve nonlinear equation
after initial guess.
ERIs (ab|cd) : may be calculated only once and store them on memory inprinciple, but it needs large O(N4) memory. They are recalculated at every cycle(direct SCF) to reduce expensive disk I/Os. This step is bottleneck.
Density Matrix×ERI = J-matrix : Coulomb potential matrix
Density Matrix×ERI = K-matrix : HF exchange matrix. GPU implementation is
not easy because it needs a lot of registers.

Hermite Gaussian Basis Set| p]  H t ( x  Px )e ( x  Px ) 2

 ( y factor )  ( z factor )t-th Hermite polynomial
A Product of two Gaussian functions |ab] is expanded exactly by HermiteGaussians |p].
(Cartecian) Gaussian

Hermite Gaussian

J ab  Dcd [ab | cd ]~Dcd | cd ]  Dq | q ]~~J p  Dq [ p | q ]~J ab | ab]  J p | p ]ERI calculation in Hermite Gaussians space is faster than that in thenormal Gaussians.
Because of the communication cost we can’t get each ERI [p|q] from GPU.
Copyright (C) 2011 X-Ability Co.,Ltd.

Procedure of J-matrix (Coulomb)
Constructshell-pair data|ab], |p], Dq
expand themin terms ofprimitive HGs

Sort |p] and |q]
for faster exitof loops.

J p   [ p | q]Dq|p]: [p|p]1/2|q]: [q|q]1/2 Dq

descending orderconverted to Jabp, q, Dqq
1. Each thread
between eachpair of [p|q].

Accumulated Jp

2. No Integral
symmetry to
avoid expensive
memory access.

Less than 10% of total cost
Copyright (C) 2011 X-Ability Co.,Ltd.

Density Functional Theory (DFT)
Wave function of electron (orbital)[ 2 / 2  veff (r )] i (r )   i i (r )Kineticenergy

Coulomb potential from nuclear (Hcore) andelectron (J) + exchange-correlation potential

Hartree-Fock :coreF  H  J   K DFT :core
F  H 
 J   v

Exchange correlation vkl    k (r )  l (r ) f (  (r ),  (r ))dr

Acceleration of J-matrix + Exchange Correlationis necessary.

Procedure of Exchange-Correlation
E xc   wi  (  (ri ),  (ri )),  (r )   (r )   (r )i  (ri )2 (ri ) l (ri )  (ri )   l (ri )

 


 k v xc  l   wi  k (ri )
i

Quadrature point ri andweights wi are generated.

ri, wi,DExc, potential on
quadrature point fi, gi ( ri )f i  wi(ri)∇(ri)

Electron density (ri),∇(ri)on quadrature points
parallelization

 (ri )   Dkl  k (ri )  l (ri )kl ( ri )g i  2 wi ( ri )fi, gi
vkl matrix

Vxc matrix：vkl    f i  k (ri )  g i   k (ri ) l (ri )Copyright (C) 2011 X-Ability Co.,Ltd.

i

DFT acceleration [GPU]Model, Basis set

Time [sec]Total energy [a.u]GAMESS

929.407-2912.2041896614

This work

305.709-2912.2041830108

GAMESS

1296.509-2927.4589680121

This work

370.807-2927.4589838167

GAMESS

2186.225-3772.6098820622

This work

651.099-3772.6098692643

GAMESS

3010.225-3792.2248655337

this work

800.743-3792.2248881864

Paclitaxel(C47H51NO14), 3-21G

Paclitaxel(C47H51NO14), 6-31G

Valinomycin(C54H90N6O18), 3-21G

Valinomycin(C54H90N6O18), 6-31G

Timing details of DFT calculation [GPU]Valinomycin,BLYP/3-21G unit:secGAMESSThis work

speedup

HF Fock Matrix
formation

248.92776.241x 3.27DFT J matrix
formation

619.465

2.529

x 244.94(※)

DFT exchange
correlatin matrix
formation

972.364218.387

x 4.45

total

2186.225

651.099x 3.36# original GAMESS J-matrix is very slow because ERIs are explicitly
calculated in Cartesian Gaussian basis.

PRISM implementation (1) [AVX]• PRISM Algorithm
One of the fastest algorithms to evaluate ERIs.
adopted by Gaussian (famous QM software)Constantly fast for various contraction length /angular momentum.
Convert Boys function [0](m) ∝∫u^m exp(-Tu)duto ERI using recurrence with a guide of “PATH”.
Chooses the best path for different contraction
length / angular momentum.

PRISM implementation (2) [AVX]
All recurrence formulae are(Integral A) = a * (Integral B) + b *(Integral C)

The same recurrence formula isapplied to the intermediate
integrals of the same shell type(angular momentum, bra/ketcontraction length).
It would be suitable for SIMD
parallelization.

Gill and Pople(1991)

PRISM implementation (3) [AVX]
class double4# operator overload{
public:double4& operator=(double d){m_d = _mm256_set1_pd(d);return *this;}double4 operator+(const double4& dd){
return double4(_mm256_add_pd( m_d, dd.m_d) );}// … define other operators__m256d m_d;}

Difference between PRISM of ours and Gaussian‘sGaussian

X-Ability

What is itclassifies shell quartets (ab | cd )by shell type and number ofbra/ket contraction so that
“Driver” is determined on eachshell quartet on the fly.

generates source code of all PRISM
paths and sets appropriate function
pointer for each calculation target
dynamically.

Pros

A recurrence formula is applied tolong-enough batch of shellquartets so that it is suitable fortraditional vector processors.

Only 4 shell quartets are processedat once, resulting to avoid cache
misses.

ConsNot suitable for scalar processors
because it might results in cash
misses. A compiler could not
optimize the program efficiently.

The code has tens of thousands linesand 300-400MB volume.

Copyright (C) 2011 X-Ability Co.,Ltd.

Calculation Result of Hartree-Fock [AVX]Time [sec]Total energy [a.u]GAMESS

184.048-2895.7814570171

This work

78.271-2895.7814570169

GAMESS

324.386-2910.6633340322

This work

153.632-2910.6633340179

GAMESS

476.829-3750.9205018138

This work

155.481-3750.9205017267

GAMESS752.839-3770.0595984968

this work

323.098-3770.0595984236

Paclitaxel, 3-21G

Paclitaxel, 6-31G

Valinomycin, 3-21G

Valinomycin, 6-31G

10-7hartree is enough accuracy. We replaced GAMESS ERI with XA PRISM.

What’s FMO(Fragment MO)

Ab initio forlarge Insulators(protein etc.)

Divide a molecule into N fragments
Generate initial density matrices for all fragments (monomers)
Bottleneck
Prepare environmental electrostatic potentials (ESP)
using previous density matrices
NOSolve Fock equations for all monomers

FICI = SICIeI, for I=1 to N

Are all monomer
densities converged?
YES
Prepare environmental electrostatic potentials for all
fragment pairs (dimers) using monomer densities
Solve Fock equations for all monomers FICI = SICIeI, for I=2 to N, J=1 to I-1
Calculate total energy and properties of molecule

Test ModelInsulin (PDBID:2HIU)
44 amino acids

Acceleration of Environmental Electrostatic
Fragment A Fragment B
Potential (ESP)• Utilize ERI J-matrix

J ab   ab cd Dcd

cd: basis on neighboring fragments

cd• Decompose protein into many small fragments
such as amino acid molecules.
Conventional SCF (ERIs saved on disk) is effective forusual amino acids (# of basis < 180).
ERI is not bottleneck in FMO.
Environmental electrostatic potential dominates.

• Utilize J-matrix acceleration program byTotal Energy of Insulin
Time [sec]Total Energy [a.u.]

Original GAMESS

14158.0-21635.4488653044

Our work

3807.5-21635.4488649218-Error of total energy of 44 fragments is smallenough.
Energy of each amino acid is highly consistent.

FMO ESP[GPU] and SCF[AVX]
calculation time of Insulin
GAMESSThis work

speedup

ESP part(GPU)10696.2

640.0

x 16.71

SCF part(AVX)

2475.1

2136.8x 1.16

Total14158.03807.5

x 3.72

ESP part is much better than total acceleration ratio.

Summary of performance results & discussion
AVXGPU

Summary of Hartree-Fock x 3 byperformanc Vectorized PRISM
e results
Discussion

Unpredictable
implementation cost

DFT
J-matrix formation: x 245
Exchange correlations: x 4
FMO
ESP : x 17
Difficulty of simple parallelization
by register shortage.
Few cost of the transfer speed
between host and device using
contracted Gaussian basis set.

GPGPU is better regarding price performance ratio.
Copyright (C) 2011 X-Ability Co.,Ltd.

Further acceleration and improvement• HF [AVX]
Improvement of evaluation part of Boys function [0](m)• DFT [GPU]Grid generation of DFT• FMO [AVX+GPU]
HF acceleration by AVX brings FMO acceleration, andinvestigate extra overhead time.

• Others, technically speaking
OpenCL, cluster tuning by hardware, etc.

Quntum Pharmaceutical System(funded by JST A-STEP)# JST is the biggest Japanese Government Funding Agency.

Batch process using combination ofsoftwares in XA-CHEM-SUITE

Visualization of PIO
by Winmostar/Tencube

Collaboration with Prof.Kitaura
Molecularmodeling software
for batch
Fastest algorithm on
fastest architecturefrom electronic statecalculation toaffinity evaluation

FMOorFMO based MD/ER
batch solver
Charge assignmentby GAMESS-US

PIO solver

XA-PIOPIE(D)ABinding free energy (ER)Molecular Dynamics(Amber / Gromacs)+ XA-SSE/AVX/CUDA-QM+ XA-ER

Start from 1st Oct 2011

Thank you for your attention.