Priority Issue 9
to be Tackled by Using Post K Computer“Elucidation of the Fundamental Laws
and Evolution of the Universe”
KAKENHI grant 17K05433, 25870168

第2回HPC-Phys 勉強会＠理研和光
2018/12/01

大規模並列計算用原子核殻模型コードの開発
清水 則孝
東京大学大学院理学系研究科
附属原子核科学研究センター

Outline• 原子核物理における殻模型計算とは• KSHELLコードの開発– M-scheme 表現– コード開発・基本アルゴリズム・ランチョス法– 並列性能– 行列・ブロック化による高速化• まとめポスト京重点課題9サブ課題B 原子核研究
複雑な核力から出発しつつ、大規模量子多体計算により、原子核の多体構造を明らかにし、その性質を計算する。素粒子、宇宙、エネルギーなどの問題へ応用。
核力（+3体力）量子多体計算 ⇒ 殻模型計算
軽い核の第一原理計算
有効相互作用の構築
摂動論(Extended KK)
3体力効果
モンテカルロ殻模型

中性子過剰核のエキゾチック核構造の解明
基礎的な核データの計算によりさまざまな応用を目指す。
中性子陽子

・ 素粒子（ニュートリノレス二重ベータ崩壊）
・ 宇宙（重元素合成過程の解明）
・ 原子力(長寿命核分裂生成物の構造計算)

C. Petrache et al., in prep.

原子核実験では、例えば、、、
1. 加速器実験
136Xe

820MeV Beam

3. 準位を決定・議論
136Ba Exp.

170Er etc.

Target136BaHigh-spin state

2. 脱励起のガンマ線を測定

理論（殻模型）との比較・議論
実験研究者に強い需要

原子核殻模型計算における配位混合• 殻模型計算における波動関数+++ ・・・

配位の数（Mスキーム次元）は中重核領域で急激に増加する。

シュレディンガー方程式を解く

   vm mH  E 

大次元疎行列の固有値問題に帰着
ランチョス法が効率的

mm H m ' v m  Ev m
m''

現在のスパコンでは、１０００億次元程度の計算が可能。

殻模型計算のスキーム
核力、(NN, NNN)

有効相互作用

模型空間

不活性コア
多体摂動計算(EKK)、
ユニタリ模型演算子法、
etc.

厳密対角化殻模型計算、
モンテカルロ殻模型

コード開発

136

56Ba80 : 殻模型計算 vs. 実験値

励起エネルギー (MeV) … 行列の固有値

arrow width : B(E2)arrow width : intensity

4316

11+(nh11/2)-2

with C. Petrache et al.

密度汎関数法

殻模型計算
第一原理計算Ref. UNEDF SciDAC collaboration

原子核殻模型のハミルトニアン
特定の質量領域用にチューンされたハミルトニアンが様々に公開されている。
⇒ 実験解析用にも多用される。
使いやすい殻模型計算コードの開発が重要。

核子は平均場中における有効相互作用を通して運動する– 一粒子エネルギー (SPE)– 二体行列要素 (TBME)

H    a na a†V(abcd;J)A

JM ( a , b) AJM ( c, d )a b ,c  d , JMa  (na , la , ja ) 一粒子軌道を指定するインデックス

na … number operators of orbit ana   cma

A†J ( a, b) †
a ,ma a ,mac11††[c†a  c†b ]( J ) jmjmJMcc
 a ab ba ,ma b ,mb
1   ab1   ab ma ,mbどのくらいの配位をとりあつかえるか？

（配位の数）

東大グループ
＠京計算機

世界最大級の大次元対角化計算（1011次元)に成功！
N. Tsunoda, T. Otsuka, N. Shimizu, M. H.-Jensen, K. Takayanagi,and T. Suzuki, Phys. Rev. C 95, 021304R (2017)

モンテカルロ殻模型
重要な「状態」だけを選び出す東大グループ独自の計算手法
旧来型計算
モンテカルロ殻模型
モンテカルロ殻模型＋京計算機

次元の制限を超えられるが、とり扱いが難しい。
今回はとりあげない。

殻模型計算コードKSHELLの開発• 大次元・実対称・疎行列の固有値問題を解く。最小固有値近傍のみが必要– ランチョス法、（クリロフ部分空間法）• 大次元の問題にチャレンジ– ハミルトニアン行列の行列要素はオンザフライで生成、ベクトルのみメモリーに保持• PCから京・ポスト京まで単一のコード、使い勝手を実現。

https://sites.google.com/a/cns.s.u-tokyo.ac.jp/kshell/
Fortran 90/95, 20,000行

KSHELLコードの公開・普及中• 公開済みの殻模型計算コードは複数あるが、並列計算可能なものはなかった。
– c.f. OXBASH 1985, ANTOINE, NuSHELL, BigSTICK, ...

• 並列性能と使い勝手を両立したコードKSHELLを開発、公開。
• 講習会を随時開催– 東大本郷（清水）、理研仁科センター（清水）、 埼玉大（清水）、ケープタウン（角田佑）、U. Oslo (J. Midtbo), TRIUMF(J. Midtbo), …

– 実験研究者の需要が強い。

ハミルトニアン行列のMスキーム表現
M=0  u1

M=0

M=0

 u3

 u2

 u1 m1  u2 m2  u3 m3  3
1 1 3m   ,  , ,
22 2 2

ハミルトニアン行列の固有値問題を解く
H ij  mi H m j

ハミルトニアン行列は疎行列。なぜ？
ハミルトニアンは2体相互作用：† †H   Tij c†cVc ijkl i c j cl cki jij

例： 2体相互作用では

m1 H m3  0ijkl

323j
2j



殻模型におけるハミルトニアン行列⇒ 実対称・疎行列

44Ti in pf-shell (4000次元)

行列要素全てをメモリーに保持することは困難
オンザフライ生成による行列・ベクトル積の実現、最低2本のランチョスベクトルをメモリーに保持する。

ランチョス法：アルゴリズム• 初期ベクトル u0• 反復する k  0, 1, 2, ...
行列ベクトル積による部分空間の拡張

ak  u Hu k
Tku 'k 1  Huk  ak uk  bk 1uk 1
bk  u 'Tk 1 u 'k 1uk 1  u 'k 1 / bk

• 3重対角行列

Gram-Schmidt直交化

規格化

T  uuiT iTHuAuj jを対角化する。

Tu
注意： 数学的には、 i u j   ij が保たれているはずである。

しかしながら、反復を繰り返していくと、丸め誤差によって直交性が失われていく。

再直交化が必須。 Thick restartが有効。

ランチョス法の収束性
56Ni pf殻の殻模型計算： 109次元の対角化計算

10個の最小固有値

8GB ランチョスベクトル... 241回の反復

Excitation energies of 56Ni
Ref. M. Horoi et al. Phys. Rev. C73 061305R (2006)

code “KSHELL”
9 sec/iteration @FX10 240 nodes (SPARC 64 IXfx, 3840 cores), total 35min.

KSHELLのアルゴリズム：ビットによる配位の表現
M=0  u1

M=0

 u2

M=0

 u3

 u1 m1  u2 m2  u3 m3  4 5 6 70 1 2 33
1 1 3m   ,  , ,
22 2 2† †m1  c†0 c†
c
1 3 c6 m

0 1 2 3 4 5 6 7 Decimal num.

|m1>

1 1 0 1 0 0 1 075|m2>

1 0 0 1 0 1 1 0105|m3>

0 0 0 0 1 1 1 1

240

323j
2j



「オンザフライ」行列要素生成による行列ベクトル積
ハミルトニアン行列要素 H ij  mi H m j の数は多すぎてメモリーに保存できない。
ランチョス法では、行列ベクトル積を実現すれば行列自体は不要。
4 5 6 70 1 2 3m3  c c c c m2
† †4 7 3 0m

0 1 2 3 4 5 6 7 Decimal num.

|m1>

1 1 0 1 0 0 1 075|m2>

1 0 0 1 0 1 1 0105

|m3>

0 0 0 0 1 1 1 1240得られたビット表現の位置を2分法で探索するc†4 c†7 c3 c0
off 0,3 bits
andon 4, 7 bits

高速化手法：陽子空間・中性子空間への分解   vi mi   v p ,n m (p ) mn(n )ip ,nm (p )mn(n )
Ref. C. W. Johnson et al., Comp. Phys. Comm. 184, 2761 (2013)

高速化手法： 陽子・中性子相互作用部分の行列・ベクトル積（ボトルネック部分）
†( )
(n )
†( ) †(n )( )(n )c†cccmm
ccmccmmmi j l kpni kpj lnp'n'c†j cl mn(n )

( )
c†c
mi kp

†
(n )
†( )
ccmccm陽子側の1体作用 i k pと中性子側の1体作用 j lをnあらかじめ求めておくことで、陽子・中性子相互作用演算を高速化する。

KSHELLにおける並列計算• 全てのランチョスベクトルをメモリーにのせられれば、再直交化の計算時間は無視できる。行
列・ベクトル積がほぼすべての計算時間を占める。
• 行列・ベクトル積をロードバランスよく並列化したい。
– 対角ブロックは計算量が大きい– 分散保持されたベクトルの通信量も無視できない

KSHELLにおける行列・ベクトル積の並列化
9プロセスの場合: 各プロセスの担当• 対角ブロックが均等に分散されている。
• 通信量が(V:ベクトルサイズ、N:プロセス数）に抑えられる。
• 再直交化の並列化は自明。

並列性能（ストロングスケーリング）
56Ni, pf-shell

109 dimension

105 スレッド並列までの良好な並列性能

計算速度の向上に向けて• オンザフライ生成による行列・ベクトル積の実行
時間がボトルネック– より詳細には、2分法によるビット表現の探索• ブロック化による高速化を試みる。
（ブロックランチョス法、ブロック櫻井・杉浦法）

ブロックランチョス法• あるいは、リスタート法を組み合わせて Thickrestart block Lanczos method
行列・ブロックベクトル積による部分空間の拡張とはのブロック行列

ブロック化による性能向上
行列・ベクトル（ブロック）積の計算時間

行列・ブロック積の計算時間

オンザフライ生成

1ベクトルあたりの計算時間

48Cr in pf-shell, 1,963,461 dimension, Xeon 20 coresブロックランチョス法による性能向上

48Cr in pf-shell, 1,963,461 dimension, Xeon 20 cores

近似固有値と厳密値との差32番目の固有値の収束性

ブロックランチョス法による性能向上

トータルの計算時間収束までの反復回数

48Cr in pf-shell, 1,963,461 dimension

112Sn w/ 100Sn core

6,210,638 dimension

ブロック櫻井杉浦法の導入： 原子核準位密度の計算
中性子捕獲反応の理解に重要な物理量。

z-Pares libraryによる櫻井ー杉浦法の導入・
確率的固有値密度推定法
200億次元行列の固有値密度推定を可能に( ( E )   tr (z-H 1 dz

Ref. N. Shimizu Y. Utsuno, Y. Futamura, T. Sakurai, T. Mizusaki andT. Otsuka, Phys. Lett. B 753, 13 (2016)

Summary• 原子核殻模型計算– 原子核の微視的構造を記述、第一原理計算にも用いられる。
– 殻模型計算・コードは実験研究者含め多岐にわたる需
要• 殻模型計算コードKSHELLの開発・公開– 大規模並列計算＋オンザフライ行列要素生成による大次元計算– ブロック化による性能向上• Future work

– 2体カレント、3体力対応などの機能拡張– ポスト京にむけて