GPGPUによるフラグメント分子軌道法の高速化/* フラグメント分子軌道法(FMO)およびエネルギー表示法(ER)を活用した自由エネルギー計算のGPGPUによる高速化 */Acceleration of fragment molecular orbital method by GPGPU

Ryota Koga 古賀 良太
President of X-Ability Co.,Ltd. ㈱クロスアビリティ# Collaborative Researchers
Yuki Furukawa(X-Ability Co.,Ltd.), Koji Yasuda(Nagoya University)
Nobuyuki Matsubayashi, Shun Sakuraba(Kyoto University)
30 May 2012
第１７回計算工学講演会＠京都教育文化センター

X-Ability Co.,Ltd. ㈱クロスアビリティ• 5 Office (foundation : 15th Jan 2008)
Japan : Hongo, Tokyo + The K-comp, Kobe
Thailand, Indonesia, and U.S. branch• 2 Business(1) Scientific Computing (mainly Chemistry)
GPGPU computation(2) Sensor Networking• By only 3 Regular Board Members

Drugs using Computational Chemistry©wikipedia
Erlotinib hydrochloride©wikipediaZanamivir(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) 2012 X-Ability Co.,Ltd.

3

Energy Representation (ER)• Free energy computation method which is ofcomparable accuracy to exact method but is fasterby two orders of magnitude• Ermod http://sourceforge.net/p/ermod/wiki/Home/– Product at a Japanese national project– Applicable to supercritical fluid, ionic liquid, micelle, andlipid membrane– All-atom computation of free energy of protein
solvation is now feasible.
– Current theoretical challenge is the ligand-binding intoprotein.
Copyright (C) 2012 X-Ability Co.,Ltd.

Trial of free energy calc on GPGPU• (a) FMO based MD # time consuming– MD using forcefield with QM(FMO) charge• (b) MD/ER– Free energy evaluation by ER using MD trajectory• (a) + (b) FMO based MD/ER -> Free energy– a combination FMO based MD with MD/ER– Protein solvation as a whole is possible, but ligand
binding scheme is now in development.

Focus on (a)FMO using GPU in the present talk
Copyright (C) 2012 X-Ability Co.,Ltd.

5Plan for implementation
GAMESS-US(FMO free application)

Replace HF-SCF, FMO-ESP
with CUDA scratch code

GAMESS-US+ GPGPU add-on

GPGPU add-on is developed by us.
1. Use best tuned algorithm as far as GPU’s resource permits.
2. Shared memory is used as L1 cache.
3. Use CPU threads (pthread) to harness GPU calculations
for complicated thread control (synchronization, queuein..)
4. Use OpenMP for CPU thread parallelization.
5. Single/double mixed-precision calculation.

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)

spz

dyz

 i (r )   Ck(i )  k (r )kKs ：  (r )   d k ek 1F (C )C  SC k ( r  A ) 2

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

Matrix form of eigenvalue problem(Self consistent)

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

7

Hartree-Fock SCF (HF-SCF) 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/ N

corea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.

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 toavoid expensive
memory access.

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

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

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

Test ModelInsulin (PDBID:2HIU)
Small Protein
44 amino acids

System Configuration(1) Intel Core i7-3930K @3.2GHz (6 core) ,GTX580 x 2, 32GB, CUDA 4.0(2) Intel Core i7-3930K @3.2GHz (6 core),GTX580 x 4, 32GB, CUDA 4.0(3) Intel Core i7-3930K @3.2GHz (6 core) ,GTX680 x 2, 32GB, CUDA 4.0+ Intel Composer XE 12.0 + MKL 10.3

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

Total Energy of Insulin
Time [sec]Total Energy [a.u.]

Original GAMESS

3279.944-21635.4488652520

Our GPGPU work

790.527-21635.4488649211(1) Intel Core i7-3930K @3.2GHz (6 core) , GTX580 x 2, 32GB, CUDA 4.0(2) Intel Core i7-3930K @3.2GHz (6 core), GTX580 x 4, 32GB, CUDA 4.0(1) + (2)-Error of total energy of 44 fragments is smallenough.
Energy of each amino acid is highly consistent.

FMO2-ESP and HF-SCF calculation
time of Insulin in total computation
GAMESSThis work

speedup

ESP part(GPU)2571.490

170.897

x 15.0

HF-SCF part(host)708.454619.630x 1.14

Total

3279.944

790.527x 4.15(1) Intel Core i7-3930K @3.2GHz (6 core) , GTX580 x 2, 32GB, CUDA 4.0(2) Intel Core i7-3930K @3.2GHz (6 core), GTX580 x 4, 32GB, CUDA 4.0(1) + (2)

ESP part is much better than total acceleration ratio.
HF-SCF part is on-memory calculation (faster than GPGPU).

Comparison of architectures(FMO2 HF-SCF)
Time [sec]Total Energy [a.u.]

Original GAMESS

7026.264-21635.4488652592

Our work [GTX580 x 2](1)1670.931-21635.4488649623

Our work [GTX680 x 2](3)1708.522-21635.4488649862(1) Intel Core i7-3930K @3.2GHz (6 core) , GTX580 x 2, 32GB, CUDA 4.0(3) Intel Core i7-3930K @3.2GHz (6 core) , GTX680 x 2, 32GB, CUDA 4.0

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

FMO2 HF-SCF 6-31G* (d-orbital)
ESP [sec]

Total time [sec] Total energy [a.u.]

Original

19915.745

25360.740 -21643.9995562825

Our GPGPU
work

1142.092

5839.737 -21643.9995952946

speedupx 17.4x 4.3(3) Intel Core i7-3930K @3.2GHz (6 core) , GTX680 x 2, 32GB, CUDA 4.0

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

FMO2 HF-Gradient 6-31GESP [sec]Total time [sec] Total energy [a.u.]

Original

13234.570

20222.516 -21635.4478224023

Our GPGPU
work

1497.2567917.169 -21635.4478217928

speedupx 8.8x 2.5(3) Intel Core i7-3930K @3.2GHz (6 core) , GTX680 x 2, 32GB, CUDA 4.0

Gradient is important for structure optimization.

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

18

FMO3 HF-SCF 6-31G
ESP [sec]

Total time [sec] Total energy [a.u.]

Original

53898.693

72271.047 -21635.6075667898

Our GPGPU
work

2589.916

19302.927 -21635.6075665173

speedupx 20.8x 3.7(3) Intel Core i7-3930K @3.2GHz (6 core) , GTX680 x 2, 32GB, CUDA 4.0

FMO3 is fmo calculation with many-body effect.

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

19

Thank you for your attention.