量子化学計算の大規模化1
石村 和也(分子科学研究所 理論・計算分子科学研究領域)
ishimura@ims.ac.jp
2017年度計算科学技術特論A 第14回
2017年7月20日

講義の内容• 本日(第14回)• 分子科学分野におけるスパコン利用の歴史• 量子化学計算とは• 量子化学計算方法とコスト• 大規模計算を行うためには• 高速化例• 来週(第15回)• 高速化・並列化例• 新たなオープンソースソフトウェアの開発• メニーコア時代に向けた開発

1

分子科学分野スパコンの変遷１•

自然科学研究機構 岡崎共通研究施設 計算科学研究センター (旧分子科学研究所
電子計算機センター)におけるCPU能力の変遷 https://ccportal.ims.ac.jp/

演算性能は10年で数百倍のペースで向上
年

理論総演算性能 (TFLOPS)

1979

3.6x10-5

1989

0.002

1999

0.092

2009

13

2015

492

2017/104,0762

分子科学分野スパコンの変遷2•

自然科学研究機構 岡崎共通研究施設 計算科学研究センター (旧分子科学研究所
電子計算機センター)におけるCPU能力の変遷 https://ccportal.ims.ac.jp/• 並列計算は必要不可欠で、大規模並列計算も当たり前の時代になりつつある• 2017年10月には、演算性能で京コンピュータ並のマシンが稼働する予定
年

機種1979 HITACHI M-180 (2台)
2000 IBM SP2 (Wide24 台)
IBM SP2 (Thin24 台)
NEC SX-5 (8CPU)
Fujitsu VPP5000 (30PE)
SGI SGI2800 (256CPU)
合計
2017/10 NEC LX (Xeon 20コア x 2, 2.4GHz) (32560 コア)
NEC LX (Xeon 18コア x 2, 3.0GHz) (5724 コア)
NEC LX (Xeon 12コア x 2, 3.0GHz) (1536 コア)(+ NVIDIA Tesla P100 x 2)
合計
地球シミュレータ(2002)
京コンピュータ(2011)

理論総演算性能(TFLOPS)
3.6 x 10-5
0.007
0.003
0.064
0.288
0.153
0.514 ←Core i7 7700K(4.2GHz,4コア)
2,500
2台分
549
221
806
4,076 ← 京コンピュータの約半分
35
10,5103分子科学分野スパコン利用の変遷• 設立当初から多くの研究者がスパコンを利用• システムの一部を長時間占有する大規模計算が近年増加– 2017年7月時点では、クラスタ専有利用枠は最大4096コア、7日間x12回/3か月

4元素と分子
1

2

H

He

3

4

5

6

7

8

9

10

Li

Be

B

C

N

O

F

Ne

11

12

13

14

15

16

17

18

Na Mg

Al

Si

P

S

Cl

Ar

31

32

33

周期表

19

20

21

22

23

24

27

28

29

34

35

36

K

Ca Sc

Ti

V

Cr Mn Fe Co

Ni

Cu Zn Ga Ge As Se

Br

Kr

37

38

39

40

41

42

46

47

53

54

Rb

Sr

Y

Zr

55

56

*1

72

Hf
*2

104 105 106 107 108 109 110 111 112 113 114 115 116 117 118

Cs Ba

25
43

26
44

45

30
48

49

50

51

52

Nb Mo Tc Ru Rh Pd Ag CdIn

Sn Sb Te

I

Xe

73

74

75

76

77

78

79

80

81

82

83

84

85

86

Ta

W

Re Os

Ir

Pt

Au Hg

Tl

Pb

Bi

Po

At

Rn

87

88

Fr

Ra

Rf

Db Sg Bh Hs Mt Ds Rg Cn Uut Uuq Uup Uuh Uus Uuo*1Lanthanoid*2
Actinoid

57

58

•
•
•

59

60

61

62

63

64

65

66

67

68

69

70

71

La Ce Pr

Nd Pm Sm Eu Gd Tb Dy Ho

Er Tm Yb Lu

89

91

92

93

97

98

99

100 101 102 103

Ac Th PaU

Np Pu Am Cm Bk

Cf

Es Fm Md No

90

94

95

96

Lr

原子の種類(元素)は110程度
現在1億3000万以上の分子が確認されている (例：H2O,
CH4)
原子間の結合は電子が担っている
5

量子化学計算とは• 量子化学計算：分子の電子分布を計算し、分子の構造、反応性、物性などを解析・予測する• 入力：原子の電子分布 (原子軌道)、原子座標• 出力：分子の電子分布 (分子軌道)• 計算量は計算方法により異なり、原子数の3乗から7乗(もしくはそれ以上)に比例して増加する
原子軌道 (C,H)
原子座標

分子軌道 (ベンゼン(C6H6))

H
H

H

H

H

HartreeFock計算

H6

量子化学計算で得られるもの
分子のエネルギー
安定構造、遷移状態構造
化学反応エネルギー
光吸収、発光スペクトル
振動スペクトル
NMR(核磁気共鳴)スペクトル
各原子の電荷溶媒効果
結合軌道解析など

化学反応 A+B → C+D
遷移状態
触媒

エネルギー

•
•
•
•
•
•
•
•
•

活性化エ
ネルギー

反応前
A+B

生成熱

反応後
C+D
反応座標
7

kernel and pentameric Au5(SeR) 6 staple motif areunprecedented and provide a structural basis for understanding the gold−selenolate nanoclusters.

T

量子化学計算適用例

hiolate-stabilized gold nanoclusters have attracted wide
research interest in recent years. To date, a number ofze-discrete gold nanoclusters have been identified,1−14 and aw of them have been structurally characterized by singleystal X-ray crystallography.15−23 In parallel with the thiolateotected gold nanoclusters, recent works have revealed that,changing theligand of thegold nanoclustersfromthiolatetoenolate (HSeR), more-stable gold nanoclusters can be
oduced, and the related properties have been studied.24−27
hese studies also found that selenolate-protected Aun(SeR)mnoclusters possess characteristics different from those of theun(SR)m counterpartsand thushaveconsiderable potential asew functional nanomaterials. However, there have been no
ports thus far on the successful crystallization of Aun(SeR)mnoclusters. In order to clarify theprecisecorrelation between
e ligand and cluster stability, the structure of nanoclusters
otected by selenolate should be pursued. Herein we report
e first structure of selenolate-stabilized Au24(SeC6H5) 20
noclusters.
Details of the synthesis are provided in the Supporting
formation. Briefly, HAuCl4·3H2O wasdissolved in water anden phase-transferred to CH2Cl2 with the aid of tetraoctymmonium bromide (TOAB). Then, both C6H5SeH andaBH4 were added simultaneously to convert Au(III) intou(I) or Au(0) by co-reduction. After reaction overnight, theueousphasewasremoved. Themixturein theorganicphaseas rotavaporated, and then washed several times with
H3OH/hexane. Dark brown crystals were crystallized from
H2Cl2/ethanol over 2−3 days. The crystals were thenllected. The structure of Au24(SeC6H5)20 was determined
X-ray crystallography. The optical absorption spectrum of
u24(SeC6H5)20 nanoclusters(dissolved in tolueneor CH2Cl2)
owsthreestepwisepeaksat380, 530, and 620nm(Figure1).
note, the optical spectrum of the thiolate counterpart, i.e.,u24(SC2H4Ph) 20 nanoclusters (dissolved in toluene or

Figure 1. Optical absorption spectrum of Au24(SeC6H5)20 nanoclusters.

CH2Cl2), shows a distinct band at 765 nm and a shoulder
band at 400 nm.28
The total structure of the Au24(SeC6H5)20 nanocluster isshown in Figure2. Asimilar structurewasdiscussed in previous

Se-Protected Au24(SePh)20 Nanocluster

O

N

貴金属

Figure 2. Crystal structure of a selenophenol-protected
Au24(SeC6H5)20 nanocluster. (Color labels: yellow = Au, violet = Se,gray = C; all H atoms are not shown).

DFT calculations by Pei et al. on thiolate-capped Au24(SR)20
nanoclusters.29 To find out details of the atom-packing
structure, we focus on the Au24Se20 framework without the
carbon tails (Figure 3A). The Au24Se20 can be divided into a
prolate Au8 kernel (Figure 3A, highlighted in green), twotrimeric Au3Se4 staple-like motifs (Figure 3B, labeled i andhighlighted with blue curves), and two pentameric Au5Se6
staple motifs (Figure 3B, labeled ii and highlighted with blue
curves). Following this anatomy, the Au24Se20 framework can

反応分子と触媒表面のくっつき ナノサイズの分子・クラス
Received: December 25, 2013
やすさ（相互作用）を計算して、Published:ターの計算(反応性、安定
February 18, 2014より反応が進む金属を調べる
2963
性)が可能になり、取り扱
える対象が大幅に拡大© 2014 American Chemical Society

dx.doi.org/10.1021/ja4131142 | J. Am. Chem. Soc. 2014, 136, 2963−2965

新たな理論や並列アルゴリ
ズムの開発によりタンパク質の計算も可能になった8

量子化学計算方法とコスト•

演算内容から分類
計算方法

Hartree-Fock（SCF), DFT法

N: 基底(or電子)数
演算量

O(N4)

データ量

通信量

O(N2)

2電子積分計算(キャッシュ内演算)
密対称行列の対角化(カットオフでO(N3)程度)O(N2)

摂動(MP2,MP3,..)法,結合クラスター(CCSD, CCSD(T),..)法

O(N5～)

O(N4～) O(N4～)

O(N5～)

O(N4～) O(N4～)

密行列-行列積

配置間相互作用(CIS, CISD,...)法
疎行列の対角化

基底の数が 2倍になると、計算量はN3：8倍、N5：32倍基底の数が10倍になると、計算量はN3：1,000倍、N5：100,000倍••

複素数演算: 重原子を含むときは相対論効果が重要になる。複素数の演算が必要になる場合がある。
高精度演算 : 将来の巨大系において、倍精度では有効桁数が不足する
可能性あり。ポスト京の次では4倍精度が必要になるかも。
9Hartree-Fock法
HartreeFock計算

原子軌道の線形結合係数(分子軌道係数)を求める

原子軌道Gauss関数

FC  εSC
Fock行列F: Fock行列, C: 分子軌道係数
S: 基底重なり行列, e: 分子軌道エネルギーF  H    Ci Ci 2 |     |  i , ,

原子軌道(AO)2電子積分

1𝜇𝜈|𝜆𝜎 = න 𝑑𝒓1 න 𝑑𝒓2 𝜙𝜇 𝒓1 𝜙𝜈 𝒓1
𝜙 𝒓 𝜙 𝒓𝑟12 𝜆 2 𝜎 2

𝜙𝜇 𝒓1 : 原子軌道Gauss関数

初期軌道係数C計算

AO２電子反発積分計算+
Fock行列への足し込み (O(N4))
Fock行列対角化 (O(N3))
分子軌道C収束
計算終了

分子軌道
収束せず

10

2次の摂動(MP2)法•••

Hartree-Fock計算で分子のエネルギーの約99%を求めることができるが、定量的な議論を行うためには
残り1%の電子相関エネルギーが重要
MP2法は最も簡便な電子相関計算方法
積分変換(密行列-行列積)計算が中心

Hartree-Fock計算(|)計算 (O(N4))

𝑜𝑐𝑐 𝑣𝑖𝑟(i|)計算 (O(N5))

𝑖𝑗(i|j)計算 (O(N5))

𝑎𝑖|𝑏𝑗 2 𝑎𝑖|𝑏𝑗 − 𝑎𝑗|𝑏𝑖
𝐸𝑀𝑃2 = ෍ ෍𝜀𝑖 + 𝜀𝑗 − 𝜀𝑎 − 𝜀𝑏𝑎𝑏 𝐴𝑂

𝑎𝑖|𝑏𝑗 = ෍ 𝐶𝜇𝑎 𝐶𝜈𝑖 𝐶𝜆𝑏 𝐶𝜎𝑗 (𝜇𝜈|𝜆𝜎)
分子軌道(MO) 𝜇𝜈𝜆𝜎
2電子積分

原子軌道(AO)
2電子積分

𝜀𝑖 : 軌道エネルギー, 𝐶𝜇𝑎 : 分子軌道係数例)・ C60(6-31G(d))の場合：O:120, V:720, N:840
・データ量1/2O2V2×8Byte = 25GB
・Xeon E5-2698(32コア, 2.3GHz)1ノード,
SMASH program,対称性無しで約771秒(ai|j)計算 (O(N5))(ai|bj)計算 (O(N5))
MP2エネルギー計算 (O(N4))

11配置間相互作用(CI)法••

主な計算内容は、疎行列の行列要素とその対角化
CIS法：1電子励起配置の線形結合
𝑜𝑐𝑐 𝑣𝑖𝑟

Ψ𝐶𝐼𝑆 = ෍ ෍ 𝑐𝑖𝑎 Φ𝑖𝑎𝑖•𝑎Hartree-Fock
電子配置

+

+…

+

+…

CISD法：HF,1,2電子励起配置の線形結合
𝑜𝑐𝑐 𝑣𝑖𝑟

𝑜𝑐𝑐 𝑣𝑖𝑟𝑎𝑏
Ψ𝐶𝐼𝑆𝐷 = 𝑐𝐻𝐹 Φ𝐻𝐹 + ෍ ෍ 𝑐𝑖𝑎 Φ𝑖𝑎 + ෍ ෍ 𝑐𝑖𝑗𝑎𝑏 Φ𝑖𝑗
𝑖

𝑎

𝑖𝑗

𝑎𝑏

例) C60(6-31G(d))の場合：O:120, V:720, N:840
CISの配置の数：O×V= 1.0×105
CISDの配置の数： O2×V2 = 1.0×1010
配置の数が疎行列の次元数12大規模計算を行うためには• 近似の導入– FMO, DC, ONIOM, QM/MMなど分割法– ECP, Frozen core, 局在化軌道など化学的知見の利用• 高速化(青字：プログラミングが特に重要になる内容)– 演算量の削減– 収束回数の削減– 実行性能の向上• 並列化– 計算機間の並列化– 計算機内の並列化– データの分散13近似の導入1• Fragment MO(FMO)法K. Kitaura, E. Ikeo, T. Asada, T. Nakano, and M.
Uebayasi, Chem. Phys. Lett., 313 (1999) 701-706.

– 巨大な分子(例えばタンパク質)をフラグメント(アミノ酸1残基もしくは複数残基)ごとに分割して、1量体と2量体のエネルギーから全体のエネ
ルギーを計算𝐸𝐹𝑀𝑂2 = ෍ 𝐸𝐼 + ෍ 𝐸𝐼𝐽 − 𝐸𝐼 − 𝐸𝐽
𝐼

𝐼>𝐽𝐸𝐼 : 𝐼番目の1量体エネルギー
𝐸𝐼𝐽 : 𝐼番目と𝐽番目の2量体エネルギー– より正確なエネルギーを求める場合は、3、4量体の計算を行う– エネルギー計算がアミノ酸残基ごとの相互作用エネルギー解析にもなっている– 現在、京数万ノードを使った計算が行われている• 他に、分割統治(DC)法、ONIOM法、QM/MM法などの分割法がある

14

近似の導入2• Effective core potential(ECP)法– 原子間の結合は価電子が重要であるため、内殻電子をポテンシャルに置き換え– Hay-Wadt(LANL2), SBKJC, Stuttgart, ....

• Frozen core近似– 電子相関計算において、内殻軌道からの励起配置を考慮しない

非占有軌道
価電子軌道

+

+…

内殻軌道
Hartree-Fock計算

電子相関計算15

近似の導入3• Localized MO法– 電子相関計算において、通常分子全体に広がっている分子軌道を局在化させて近くの軌道間の相関のみ考慮

H2O
通常のCanonical MOs

Localized MOs• Resolution of the identity(RI)法, 密度フィッティング法– 補助基底を導入して、4中心積分を3中心積分などの積で表現することで計算量やデータ量を削減• Fast Multipole法(FMM)など他にも様々な方法がある

16高速化1(演算量の削減)••

演算量の少ない1,2電子積分計算アルゴリズム開発
Schwarzの不等式を用いたカットオフ– AO2電子積分(|)の計算実行を基底シェルごと(M,N,L,S)に判断– O(N2)のデータ(|)でO(N4)の計算を判断– 内殻電子の基底exp(-ar2)のaが大きい場合、異なる原子間では重なりが小
さくなり、大半がスキップされる– 密度行列 (𝐷𝜇𝜈 = 2 σ𝑜𝑐𝑐
𝑖 𝐶𝜇𝑖 𝐶𝜈𝑖 )の差も掛けた上で判定すると、収束に近づくにつれてより多くの計算がスキップされる
exp(-ar2)
Fock行列作成4重ループ
do M1, Nbasis
do N=1, M
do L=1, M
do S=1, L
if({DD((MN|MN)(LS|LS))1/2}≧threshold) |ブロック計算
enddo
enddo
enddo
enddo

a：大

a：小
17

高速化2(演算量の削減)•

対称性の利用– 基底が実数の場合のAO2電子積分： (|) = (|) = (|) = (|) = ...
– ⇔、⇔、⇔の入れ替えが可能で、計算量は1/8に1𝜇𝜈|𝜆𝜎 = න 𝑑𝒓1 න 𝑑𝒓2 𝜙𝜇 ∗ 𝒓1 𝜙𝜈 𝒓1
𝜙𝜆 ∗ 𝒓2 𝜙𝜎 𝒓2
𝑟12
𝜙𝜇 𝒓1 : 基底関数(原子軌道Gauss関数)– 具体的には、Fock行列計算において、(𝜇𝜈|𝜆𝜎)1つで複数のAO2電子積分をカ
バー
1F  H    D  |     |  
2


 ,
𝐹𝜇𝜈 = 𝐹𝜇𝜈 + 𝐷𝜆𝜎 𝜇𝜈|𝜆𝜎 + 𝐷𝜎𝜆 𝜇𝜈|𝜎𝜆𝐹𝜆𝜎 = 𝐹𝜆𝜎 + 𝐷𝜇𝜈 𝜆𝜎|𝜇𝜈 + 𝐷𝜈𝜇 𝜆𝜎|𝜈𝜇1
𝐹𝜇𝜆 = 𝐹𝜇𝜆 − 𝐷𝜈𝜎 𝜇𝜈|𝜆𝜎2
1𝐹𝜈𝜎 = 𝐹𝜈𝜎 − 𝐷𝜇𝜆 𝜆𝜎|𝜇𝜈2…

18

高速化3(収束回数の削減)• SCFの収束回数削減– Direct inversion of the iterative subspace(DIIS)法– 近似Second-Order SCF(SOSCF)法– Quadratically convergent SCF法
– Level shift法 (HOMO-LUMOギャップが小さい時に有効)– 小さな基底で分子軌道計算 -> それを初期軌道にして大きな基底での計算 (DFT計算で重要)• 構造最適化回数削減– Newton-Raphson法– Hessianのアップデート: BFGS法– 効率的な座標系 : Redundant coordinate, Delocalized coordinate (分子の結合、角度、二面角の情報から初期Hessianを改良)

19

高速化4(計算機の効率的な利用)••

1990年代ごろまでは、SCF計算においてAO2電子積分を一度だけ計算してディスクに保存し、何度もディスクから読み込んでいたCPU性能が大幅に向上したため、AO2電子積分を何度も計算するDirect
SCF法が主流になった
Conventional SCF法

Direct SCF法

初期軌道係数C計算

初期軌道係数C計算

AO２電子反発積分計算+
ディスクへの書き込み

AO２電子反発積分計算+
Fock行列への足し込み

AO２電子反発積分読込+
Fock行列への足し込み

Fock行列対角化

Fock行列対角化

分子軌道
収束せず

分子軌道
収束せず

分子軌道C収束

計算終了

分子軌道C収束

計算終了

20

高速化5(実行性能の向上)• SIMD演算の利用– 1つの演算命令で複数のデータを処理– do、forループを多用• 時間のかかる演算回数を削減– 割り算、組み込み関数の利用回数削減 (まとめて演算する、ループ
内で同じ演算はループの外に出す)• データアクセスを効率的にしたアルゴリズム・プログラム開発– 連続したデータアクセスでキャッシュミスの削減– キャッシュ内のデータを再利用し、メモリアクセスを削減– 多次元配列の取り扱い 例：A(x, y, z) or A(z, y, x)
演算器
アクセススピード
容量L1キャッシュ速小メモリ

L2キャッシュ遅大21

高速化6(実行性能の向上)• コンパイラの最適化オプションの設定– 基本的にはコンパイラが最適化しやすいように、わかりやすいシンプ
ルなコードを書く必要がある– データの依存関係や多重ループ内にIF文があるなど複雑になると、最適化されない場合が多い• BLAS, LAPACKなど数学ライブラリの利用– BLASライブラリはCPUの性能を引き出してくれるが、小さい配列(100次元程度)の場合、サブルーチンコールのオーバヘッドの方が大きくなる可能性がある– BLAS2を数多く使うより、BLAS3でまとめて実行
演算量

データ量

BLAS2 (行列-ベクトル)

O(N2)

O(N2)

BLAS3 (行列-行列)

O(N3)

O(N2)
22

並列化
• ノード間(MPI)、ノード内(OpenMP)それぞれで並列化– 式の変形– 多重ループの順番の工夫• 均等な計算負荷分散– ノード間で分散、さらにノード内でも分散• 大容量データの分散– 京のメモリは1ノード16GB, 8万ノードでは1PB以上– ノード間では分散、ノード内では共有– 中間データ保存用としてディスクはあまり期待できない• 通信量、通信回数の削減– 並列計算プログラムのチューニングで最後に残る問題は通信関係(特にレイテンシ)が多い

23

並列化率と並列加速率(Speed-Up)
計算負荷が均等に分散されている場合
数十コア(PCクラスター)

数千コア(スーパーコンピューター)並列化率99%でもある程度速くなる

並列化率99.9%でも不十分
8192.064.0
99.999%99.99%
99.9%
99.0%

6144.0Speed-up

Speed-up

48.0

32.0

4096.016.0

2048.0

0.0

0.0

0

16

32
CPUコア数99.999%99.99%
99.9%
99.0%

48

64

0

2048

4096
6144
CPUコア数8192

24

並列化手法OpenMP(ノード内)
使い始めるのは簡単だが、性能を引き出すのは大変
今後ノードあたりのコア数はさらに増えると予想され、ますます重要になる

MPI(ノード間、ノード内)

難易度

使い始めは大変だが、性能を出すのはOpenMPよりも容易

BLAS,LAPACKライブラリ利用(ノード内)
行列・ベクトル演算を置き換え、行列×行列は効果大Intel MKL
AMD ACML

コンパイラの自動並列化(ノード内)
コンパイルするときにオプションを付けるだけ
ループの計算を分散、あまり効果が得られないことが多い
pgf90 -Mconcur
ifort -parallel25

MPI+OpenMPハイブリッド並列化
ノード間をMPI、ノード内をOpenMPで並列化
MPI

MemoryMPIOpenMP

core core

core core

core core

core core

core core

core core

core core

core core

core core

core core

MemoryMemoryOpenMPメリット
並列化効率の向上

メモリの有効利用

ノード内ではOpenMPで動的に負荷分散
MPIプロセス数削減による計算負荷バランスの向上
MPIプロセス数削減による通信の効率化
OpenMPによるノード内メモリの共有

デメリット
アルゴリズム、プログラムが複雑になり、開発コスト上昇
26

MPI通信の最適化••

通信時間はデータサイズによって主要因が異なり、対応策が異なる– 小さいデータ：レイテンシ(遅延)時間 → 送受信回数削減– 大きなデータ(約1MB以上)：バンド幅 → 送受信データ量削減
小さいデータの場合：何度も送受信するより、配列にデータを集めて一度に送受信– 現在のMPI通信のレイテンシは0.1秒以上 (Xeon E5-2698(32コア,
2.3GHz)では約10万回の演算に相当)
MPI_Send,Bcast3回
A(10)
B(20)
C(30)
node 0

a(10)
b(20)
c(30)
node ***MPI_Send,Bcast1回
A(10)
B(20)
X(60)
C(30) コピーnode 0

a(10)
x(60)
b(20)
コピー c(30)node ***

27

OpenMPの変数(Fortran)•

OpenMP並列領域では、すべての変数をノード内で共有する変数(shared)と各スレッドが別々の値を持つ変数(private)に分類– 既存のコードにOpenMPを導入する場合、privateにすべき変数の指定
忘れによるバグに注意– デフォルトはshareduse module_A: cmsi
real*8 : tcci
– common、module変数はshared!$OMP parallel do– do変数はprivate
do j = 1,n– OpenMP領域で呼ばれた関数や
call test(j, tcci)
サブルーチン内で新規に使われる
enddo
変数はprivate...
Subroutine test(j, tcci)
use module_A: cmsi
integer : ims
28

OpenMPのオーバーヘッド•••!$OMP parallelや!$OMP do (特にschedule(dynamic))のコストは、OpenMP領域の計算量が少ないと無視できない→ できるだけ多くの計算(領域)をまとめてOpenMP並列化→ 多重ループの分散では、外側のループを並列化
排他的処理criticalやatomicを多用すると、待ち時間が増加し効率が低下する→ できるだけ上書きをしないコードにする→ スレッドごとに変数を用意(private, reduction)common、module変数をprivate変数にする場合threadprivateは便利だが、OpenMP領域で頻繁に呼ばれる関数やサブルーチンで利用するとそのコストが大きくなる場合がある→ common、module変数からサブルーチン、関数の引数に変更

29高速化と並列化の重要性と難しさ
•
•
•

•
•
•
•
•

ノードあたりのCPUコア数はますます増加すると予想される
高速化と並列化は、スパコンだけではなく研究室レベルのPCクラスタでも必須になりつつある
スカラCPU搭載スパコンもPCクラスターも、開発方針、ソースコードはほぼ同じ
どの程度力を入れるかは、目的によって異なる分野によっては、式の導出からやり直す必要がある
京コンピュータはスカラCPUであるため、京だけのチューニングは多くない
既存のソースコードの改良では性能を出すのが難しい場合がある
開発コストは増える一方なので、オープンソースでの公開と共有で分野全体でのコスト削減を考える段階に来ている

30

原子軌道(AO)2電子積分計算アルゴリズム
𝜇𝜈|𝜆𝜎 = න 𝑑𝒓1 න 𝑑𝒓2 𝜙𝜇∗1
𝒓1 𝜙𝜈 𝒓1
𝜙𝜆 ∗ 𝒓2 𝜙𝜎 𝒓2
𝑟12
𝜙𝜇 𝒓 : 基底関数(原子軌道Gauss関数)

•
•
•
•
•
•
•
•

Rys quadrature (Rys多項式を利用)
Pople-Hehre (座標軸を回転させ演算量削減)
Obara-Saika (垂直漸化関係式(VRR))
McMurchie-Davidson (Hermite Gaussianを利用した漸化式)
Head-Gordon-Pople (水平漸化関係式(HRR))
ACE (随伴座標展開)
PRISM(適切なタイミングで短縮(contraction)を行う)
他にも数多くのアルゴリズムが提案されている

31

AO2電子積分計算のコスト𝜇𝜈|𝜆𝜎 = න 𝑑𝐫1 න 𝑑𝐫2 𝜙𝜇∗𝐾1
𝐫1 𝜙𝜈 𝐫1
𝜙𝜆 ∗ 𝐫2 𝜙𝜎 𝐫2𝑟12基底関数 𝜙𝜇 𝐫 = ෍ 𝑐𝑖 𝑥 − 𝐴𝑥 𝑛𝑥 𝑦 − 𝐴𝑦𝑛𝑦𝑧 − 𝐴𝑧 𝑛𝑧 exp −𝛼𝑖 𝐫 − 𝐀 2𝑖

水素のSTO-3G
基底関数(K = 3)HS 3
3.42525091 0.15432897
0.62391373 0.53532814
0.16885540 0.44463454

↑
𝛼𝑖

↑
𝑐𝑖

2電子積分の演算量=xK4+yK2+zK(基底関数の短縮数)に依存している(sp,sp|sp,sp)の計算コスト
Method
PH
MD
HGP

DRK

x

220

1500

1400

1056

y

2300

1700

30

30

z

4000

0

800

800

32

高速化例(新規アルゴリズム開発)
K. Ishimura, S. Nagase, Theoret Chem Acc, (2008) 120, 185-189.



Pople-Hehre (PH)法
座標軸を回転させることにより演算量を減らす
xAB=0
yAB=0
yCD=0

xPQ=一定yPQ=一定B
Pz(AB|CD)
D

y

x
Q



McMurchie-Davidson (MD)法C
漸化式を用いて(ss|ss)型から角運動量を効率的に上げる[0](m)(=(ss|ss)(m))-> (r) -> (p|q) -> (AB|CD)A

2つの方法の組み合わせ
座標軸回転 → 漸化式を用いて軌道角運動量を上げる → 座標軸再回転

33

新規アルゴリズムの特徴[0](m)(=(ss|ss)(m))[A+B+C+D][A+B|C+D][AB|CD]xPQ=一定yPQ=一定



xAB=0
yAB=0
yCD=0
演算量=xK4+yK2+z (K:基底関数の短縮数)(sp,sp|sp,sp)型の場合
MethodPH

PH+MD

Method

PH

PH+MD

x

220

180

K=1

6520

6583

y

2300

1100

K=2

16720

12490

z

4000

5330

K=3

42520

29535

6-31G(d)やcc-pVDZなど適度な短縮数の基底関数で性能発揮

34

プログラム開発
•
•
•
•
•

•
•

(ss|ss)から(dd|dd)までの21種類の積分計算ルーチンを作成
座標回転のメリット(xAB=0, yAB=0など)を考慮した漸化式を導出、その計算コードを作る自動生成プログラムをFortranとPerlで作成
文字列をFortranで出力し、Perlで整形
約2万行のコードを自動生成し、デバッグを含めた開発時間を短縮
自動生成以外のコードは、基本的にdoループで書けるよう配列データの並びを工夫し、さらに割り算やsqrt、expなど組み込み関数はできるだけまとめて演算詳細は、SMASHプログラムのint2sp.F90, int2spd[1-4].F90、もしくはGAMESSプログラムのint2b.src (sp関数), int2[r-w].src (spd関数)を参照
SMASHでは2電子積分ルーチンをライブラリ化し、引数でのみデータのやり取りを行っているため、容易に他のプログラムへ移植することができる

35

計算結果
GAMESSに実装して、Fock行列計算時間(sec)を測定
計算機：Pentimu4 3.0GHz
分子

•
•
•
•
•

Taxol(C47H51NO14)

Luciferin(C11H8N2O3S2)

基底関数

STO-3G(361 AOs)

6-31G(d)(1032 AOs)

aug-cc-pVDZ(550 AOs)

Original GAMESS (PH)85.7

2015.2

2014.9

PH+MD

69.9

1361.81154.5

2 – 4割計算時間を削減
d関数を含む場合、削減の割合は大きくなる
2005年からGAMESSにデフォルトルーチンとして正式導入
演算の約8割はdoループ内で行われるため、現在のCPUに適した方法
座標軸を元に戻す変換行列に6dから5dへの変換を組み込むことが可能

36

まとめ
•
•
•
•
•

•

分子科学分野では35年以上前からスパコンを使い、それぞれの時代で可能な大規模計算を行ってきた量子化学計算は、分子のサイズが大きくなると計算量、データ量、通信量が急激に増加する
高速化・並列化はスパコンだけではなく、研究室レベルのPCクラスタでも不可欠になりつつある
京コンピュータレベルの計算機が全国のスパコンセンターに徐々に導入され始めている
計算機はますます複雑になり(CPU内部もコア数も)、アルゴリズムとプログラム開発はさらに大変になるため、分野全体での開発コスト削減を考える必要がある
プログラミングだけではなく、演算量の削減、収束回数の削減、適切な近似の導入も重要である

37