計算科学技術特論B 第9回講義資料 2016年6月16日

オーダーN法2
第8回
2016年
6月9日

第9回
2016年
6月16日

• 低次スケーリング法の重要性• 計算のオーダーとは？
• DFT計算におけるオーダーN法• オーダーN Krylov部分空間法• 第一原理計算への拡張• 局在基底法• オーダーN法の超並列化の方法• オーダーN法の応用• 数値厳密な低次スケーリング法• O(N)準厳密交換汎関数

尾崎 泰助
東京大学 物性研究所

様々な局在基底法
第8回の講義で議論した様に、オーダーN法は基底関数の局在性もしくはKohn-Shamハミルトニアン行列の疎性を前提に定式化されている。

局在性を持つ基底関数や離散化方法として、以下の方法が知られている。
• ガウス基底法• 原子様数値基底法• 有限要素法• 有限差分法• Wavelet法• 特殊関数法(Sinc関数等)
本講義では原子様数値基底法を議論する。

原理的に完全系を張る系統的基底関数は有限要素法、有限差分法、Wavelet法、
特殊関数法である。準局所近似では超並列計算を行うことで、大規模計算を実行可能であるが、非局所汎関数を用いた場
合、行列の疎性が失われ、計算コストは大幅に増大する。一方、完全系ではない
ものの少数の最適化基底関数でほぼ収束解が得られる手法として、ガウス基底法、
原子様数値基底法がある。目的に応じて計算手法の使い分けが必要であろう。

LCPAO法(Linear-Combination of Pseudo Atomic Orbital Method)
一電子KS軌道を擬原子軌道の線形結合で表現する方法

特徴:
・原子様基底を用いているため、計算結果の物理的・化学的解釈が容易。
・基底関数系は厳密には完全系ではないため、完全な収束解を得ることは難しい。
しかし物理的起源に基づく基底のため、注意深く構成された場合には比較的少ない基底関数で実質的に収束解を得ることが可能。

・基底が局在しているため、硬い擬ポテンシャルの使用が可。高精度な擬ポテ
ンシャルを作成することが可能。
・基底の局在性のため、ハミルトニアン行列の計算量と保存メモリ量がO(N)となる。
・様々なO(N)法のアイデアと整合性が良い。

プリミティブ擬原子基底関数1.

閉じ込めポテンシャルの下で原子のKS方程式
を自己無撞着に解く。

s-orbital of oxygen

2. 次にこの閉じ込めポテンシャル下でノルム保存擬ポテンシャルを作成。
3. 各角運動量成分lに対して基底状態と励起状態をstep2で得られた擬ポテン
シャルに対して数値的に求める。

多くの系において、精度と計算効率は二つのパラメーターで制御される。

カットオフ半径基底関数の数PRB 67, 155108 (2003)
PRB 69, 195113 (2004)

基底関数に対する収束性
基底のカットオフ半径と基底を増やすとエネルギーと構造パラメーターは系統的に収束している。数学的には厳密ではないが、二つのパラメーターは変分パラメーターと見なし得る。

分子固体

二原子分子の基底状態の計算
擬原子基底による二原子分子の基底状態計算結果

局所密度近似(LDA)における全ての“成功”と“失敗”がDNP程度の基底サイズで再現できる。DNP=二つの価電子基底+一つの分極関数基底数値擬原子基底関数の変分最適化

LCPAO法の計算精度・計算速度は基本的に基底のカットオフ半径

基底の数で制御される。しかし基底関数の形自体を最適化出来るならば、計算量を増やすことなしに、計算精度を向上させることが出来る。

Ozaki, PRB 67, 155108 (2003)

基底関数の変分最適化
1電子KS波動関数

縮約軌道

縮約係数aを固定し、cに関して変分をとれば

a: 縮約係数→
変数cがaに対する従属変数であると見なし、さらに変数cの観点でKS方程式
が自己無撞着に解けていると仮定すると、次式が得られる。

Ozaki, PRB 67, 155108 (2003)

基底関数の変分最適化の手続き
前ページの変分計算において、“変数cがaに対する従属変数であると見なし”、という仮定が置かれていた。この数学的な構
造は構造最適化の際の座標Rと変数cの関係と同じものである。
この場合には変数cは座標Rの従属変数である。したがって、
構造最適化と同じ手続きで縮約係数aが最適化できる。
ステップ 1ステップ 2固定したaでKS方程式を自己無撞着に解く。
Eのaに関する勾配を求め、エネルギーが低下する方向にaを更新する。

勾配が閾値に達するまで、ステップ1と2をGDIIS法を用いて反復する。

Primitive vs. Optimized
エネルギーの収束性

炭素原子の動径関数

フッ素原子は炭素原子のp軌道から電子を引き抜く。
そのため、炭素原子のp軌道が大きく局在化する。

擬ポテンシャル及び基底関数のデータベース

高精度かつ高効率計算を実現するために、データベースを構築し、一般公開。

データベースを構築する上でのWishリスト• 擬ポテンシャルと全電子計算の対数微分は目視の範囲で完全に一致する。
• 固体の平衡格子定数がWien2kと比較して、0.005Å以下である。

• 固体の体積弾性率がWien2kと比較して、3GPa以下である。
• バンド構造がWien2kと比較し、ほぼ同等である。
• 分子の最適化構造がこれまでの報告例とほぼ同等である。
• 分子の原子化エネルギーがこれまでの報告例とほぼ同等である。
• 基底関数重なり誤差が0.5 kcal/mol以下である。

LCAO法における基底関数最適化の一般的考察 (1)NLCAO: LCAO法における基底関数の数
NPW:
平面波法における基底関数の数

NLCAO << NPW• LCAO法の基底数はPW法のそれと比較して、著しく小さい。ゆえにPW基底の張る空間に比べて、LCAO基底の張る空間は小さい。
したがって一般にはLCAO法はPW法より低精度であると考えられる。
• 一方で各元素は固体および分子中で、いくつかの類似した化学的環境に置かれている。したがって化学的環境のトレーニングセッ
トを適切に選択の上、それぞれの環境に対して基底関数を予め最適化し、一つの基底関数セットに集約しておけば、非常に高精度
なLCAO基底を作成できる可能性がある。
• 実際の応用計算ではトレーニングセットとは異なる化学的環境にあるが、近似的にはいくつかの選択した化学的環境の線形結合で表現可能であると考えられる。

LCAO法における基底関数最適化の一般的考察 (2)全ヒルベルト空間

ある元素の張る部分空間

各元素が張る空間は全ヒルベルト空間の部分空間である。この部分空間を張る
最適な基底セットを見つけることが出来るならば、可搬性の高い高精度基底が作成できる可能性がある。そこで以下の手続きで基底関数の最適化を試みる。

1.
2.
3.
4.

化学的環境のトレーニングセットの選択
各環境に対する基底関数の変分最適化
部分空間の回転による冗長成分の消去
異なる最適化基底関数の一つの基底セットへの集約

LCAO基底の最適化: 類似する試み
類似するアイデアでLCAO基底の最適化に関する試みが行われている。
• “Gaussian basis sets for accurate calculations on molecularsystems in gas and condensed phases”, J. VandeVondele
and J. Hutter, JCP 127, 114105 (2007). → CP2K code
• “Ab initio molecular simulations with numeric atom-centeredorbitals”, V. Blum et al., Comp. Phys. Comm. 180, 2175
(2009). → FHI-AIMS code

両グループは適切な化学的環境下で基底関数を最適化した場合、高精度なLCAO基底関数を構築できることを示した。

基底関数の最適化手順
1. 典型的な化学的環境を選択

2. 動径波動関数を変分最適化

3. 部分空間内で回転し冗長成分を取り除いた後に統合

××

変分最適化基底: P原子の場合

孤立原子とダイマーの最適化基
底はほぼ同一。BCCの最適化基
底は局在化している。

配位数が増加するに伴い、最適化基底は局在化している。

最適化基底における冗長性
Mn原子の擬ポテンシャルにはFull space
3s, 3p, 3d, 4s 電子が含まれている。
3s電子と3p電子は化学的環境依存性が小さいため、全ての最適化された3s基底と3p基底を含めると基底が冗長になる。
これらの最適化基底関数の冗長性は取り除かれるべきである。
しかし最適化の際に3s (3p)および4s(4p)軌道の混合が生じるため、冗長性を取り除くことが困難となる。
この冗長性は最適化基底セットを部分空間内で回転することで、取り除くことが出来る。

最適化基底の部分空間内での回転(1) 例えば、s5>3 を考える。この記号は5つのプリミティブs基底から3つの最適化s基底を作成することを意味する。最適化によって得られた縮約係数は右図の行列で与えられる。
(2) 一般化ガウス消去法を適用することで、行列Aは行列Bに変形できる。
(3) さらに行列Bに対してグラム-シュミット直交化を課すことで、行列Cが得られる。

この変換により、最初の最適化基底は元のプリミティブ基底とほぼ同形となり、基底の冗長性を取り除くことが出来る。基底最適化の効果は二番目以降に押し込まれる。

DFT計算の再現性に向けた試み
Science 351, aad3000 (2016)

15コード
69名が参加
71 elemental bulks
GGA-PBE
スカラー相対論

Siの平衡格子定数 by PBEの推移
誤差の原因 • 基底関数

• 擬ポテンシャル• 実空間積分• 逆格子空間積分

Δ-ゲージ
異なるコードの精度を比較する方法

GGA-PBEと実験の比較: Δ-factor

GGA-PBEと実験との比較において、58元素の平均Δ-factor
は23.5meV/atomである。擬ポテンシャル、基底関数、数値積分から生じる誤差は1/10程度、2meV/atom程度である
ことが望ましい。

Δ-ゲージを用いた各コード間の比較

OpenMXの平均Δ-ゲージは2.0meV/atomである。

オーダーN法の超並列化手法• 原子の3次元空間分割
修正再帰二分法, 慣性モーメントテンソル法• グリッドの領域分割
4種のグリッド構造による空間分割• 高次元FFTの並列化
References:
A three-dimensional domain decomposition method for large-scale DFT
electronic structure calculations
T.V.T. Duy and T. Ozaki, Comput. Phys. Commun. 185, 777-789 (2014).
A decomposition method with minimum communication amount forparallelization of multi-dimensional FFTs
T.V.T. Duy and T. Ozaki, Comput. Phys. Commun. 185, 153-164 (2014).

O(N)法の並列化
ほぼ独立に各原子の計算が実行可能。各MPIプロセスに原子を割り当て、計算する。

原子の3次元空間分割法
通信量・メモリ量を軽減させるため、原子集団を空間的に局所性を保ちつつ、分割する必要がある。

望ましい条件:• 空間局所性• ほぼ等価な計算量• 任意の系に適用可能• 単純で計算コストの小さなアルゴリズム

原子の空間分割(1) 修正再帰二分法により、MPIプロセスを分割(2) 慣性テンソル行列から決定された主軸上で原子を分割する。

原子の空間分割
修正再帰二分法により、MPIプロセスを分割。例えば19プロセスの場合には以下の様に分割する。

2のべき乗以外の任意のプロセス数に適用可能。

慣性モーメントテンソルを用いた原子の並び替え

再帰的な原子の空間分割

主軸への射影による空間
分割を再帰的に行う。分割の比は予め作成してあるMPIプロセスのバイナ
リーツリーの比で決まる。
分割比の重みは各原子の1ステップ前の計算時間を用いる。
空間局所性とほぼ等価な計算量が実現できる。

原子の空間分割: 実例
ダイヤモンド 16384 atoms, 19 processes

枝分かれCNT, 16 processes

グリッドデータの領域分割
行列要素の計算、Poisson方程式の解法等で、ユニフォームメッ
シュが使用される。数万原子系ではその自由度は1億を超える。
グリッド分割構造A,B,C,Dを構築
し、それぞれの計算で最適な構造を選択。
分割構造間の通信が最小化される様に設計。

3D-FFTの2D並列化の分割方法
Row-wise型

新しい領域分割法を提案

Pencil型

これまで用いられてきた領域分割法

Duy and Ozaki, CPC 185, 153 (2014).

3D-FFTの2D-並列化(row-wise型)

1D-並列化と比較し、
分割数Nまでは通信量は増加せず。
分割数N2でも2倍の増加。

通信量の比較：三つの並列化手法
Compared to 1D-parallelization, no increase of MPI communication up to N. Even at
N2, just double communication.

OpenFFTと他のライブラリとの比較
新しい並列化手法により高並列FFT計算が実現。
Cray-XC30上で実施

1283 data points
OpenFFTはGNU-GPLコードとして以下のwebよりdownloadできる。http://www.openmxsquare.org/openfft/

「京」上でのO(N)法の並列化効率

系: 131072原子からなるダイヤモンド格子

13万1072コアを使用しておよそ68％の並列化効率

オーダーN法の応用例• bcc-鉄と炭化物の界面構造の計算
H. Sawada, S. Taniguchi, K. Kawakami, and T. Ozaki,
Modelling Simul. Mater. Sci. Eng. 21, 045012 (2013).

bcc-鉄中の炭化物の析出
By 新日鐵住金・澤田英明 氏

純鉄は構造材料としては強度が十分でない。炭化物の析出を利用し、 強度・靱性の制御が行われる。

整合析出

HRTEM image
準整合析出

不整合析出

析出物:

TiC, VC, NbC

析出物直径と抵抗力
Y. Kobayashi, J. Takahashi and K. Kawakami, Scripta Mater. 67 (2012) 854

TiC

TiC

TiC

粒子1個あたりの抵抗力 F (N)

粒子1個あたりの抵抗力 (N)

3.E-092.E-09TiC析出モデル鋼
析出強化量
0MPa
Ti

1.E-093D アトムプローブ像

TiC 球状183MPa10nm

TiC 板状317MPaCu 球状

0.E+00012

粒子直径 R (nm)

3

4

5

界面の面積当たりのエネルギー

整合・準整合析出の関係

歪場

準整合界面エネルギー
整合界面エネルギー

析出物直径

部分整合格子の単位胞
部分整合格子の格子不整合の大きさ
Fe/NbCの場合計算に必要な原子数 a0 ( Fe)  n  a0 ( NbC) （Å）
mabs(a(NbC)-a(Fe))

1.6
4

1.4
1.2

3

1.0

0.8

2

6

number of
Fe unit, m

14
12

20

0.2
0.0

1

10

8

ＶＣ
Fe: 36
VC: 35

ＮｂＣ
Fe: 9
NbC: 8

3

2118

11238

627

5

3530

18730

1045

7

4942

26222

1463

16

0.6
0.4

ＴｉＣ
Fe: 16
TiC: 15

15

11

7

# of
layer

9

19

17
18

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

number of NbC unit, nStructures of coherent andsemi-coherent precipitates
Baker-Nutting orientation:(001)precipitate//(001)iron and [100]precipitate//[110]ironCoherent caseSemi-coherent case

NbC layer
Nb

Nb

C

C

Fe

Fe

c

Fe layer

21 atomsab1463 atoms

整合・部分整合析出の最適化計算

c0

21

14

13

20

19

12

11

18

17

10

9

16

15

8

interface energy (J/m2)

5.0

Fe7/(NbC)7

5
4
3

13

20

19

12

11

18

17

10

Fe-Nb

9

16

Fe-C

15

8

semi-coherent

7

3.0

2.0
1.0

6

0.0

5

24.0

25.0

26.0

27.0

28.0

4

29.0
3

2
1

14

4.0

7
6

21

c0 （ Å ）2
1

 界面エネルギー： 整合（Fe-C） ＜ 部分整合 ＜ 整合（Fe-Nb）
 部分整合でのFe相とNbC相の層間距離はFeとNbが近接する整合界面に近い

19

12

11

18

17

10

9

16

15

8

-8

-4

4
3
2

DOS (staes/eV/atom)

5

-2

0

2

4

Nb
C

c15

C

c18

-6

-4

-2

0

2

-8

-6

nb18

-4

Fe4
Fe7

2.0

1.5
1.0
0.5

-2

0

2

4

Energy (eV)

2.0
1.8
1.6
1.4
1.2
1.0
0.8
0.6
0.4
0.2
0.0

c11

-6

-4

-2

0

2

4

Energy (eV)

-8

-6

-4

-2

0

Energy (eV)

2

4

13

20

19

12

11

18

17

10

9

16

15

8

Fe4

6

Fe7

2.0

5

1.5

4

1.0

3

0.5

2

0.0

0.0

14

7

Fe

2.5

21

c8

C

-8

4

nb15

Nb

3.0

Fe

2.5

2.0
1.8
1.6
1.4
1.2
1.0
0.8
0.6
0.4
0.2
0.0

Fe

Energy (eV)

3.0

6

1

Nb11

Energy (eV)

2.0
1.8
1.6
1.4
1.2
1.0
0.8
0.6
0.4
0.2
0.0
-8

7

-6

DOS (staes/eV/atom)

20

Nb8

Nb

DOS (staes/eV/atom)

13

2.0
1.8
1.6
1.4
1.2
1.0
0.8
0.6
0.4
0.2
0.0

DOS (staes/eV/atom)

14

DOS (staes/eV/atom)

21

DOS (staes/eV/atom)

界面における局所状態密度

-8

-6

-4

-2

0

Energy (eV)

2

4

1

部分整合析出における界面構造

H. Sawada et al., Modelling Simul. Mater.
Sci. Eng. 21, 045012 (2013).

48

歪みエネルギーの推定
Model potential method： Finnis-Sinclair
Fe atoms：432,000

Precipitation

Mother phase

8.94Å

Precipitation

Mother phase

Coherent：
10% expansion
Semi-coherent：0% expansion

4% expansion（Due to dislocation）

整合/部分整合構造の遷移

計算値
実験

2.3 nm
2～3 nm

TEM像と硬さ試験より見積もられた値
Kobayashi et al., Scripta Materialia 67, 854 (2012).

低次スケーリング数値厳密対角化法
T. Ozaki, PRB 75, 035123 (2007).
T. Ozaki, PRB 82, 075131 (2010).

Density functionals as a functional of ρ
Density functionals can be rewritten by the first order reduced density matrix: ρ

where the electron density is given by ρ

All we need to calculate

1. Non-zero matrix elements of H and S
2. Density matrix corresponding to the non-zero matrix elements

Main difficulty: ‘diagonalization’O(N3) method - Numerically exact diagonalization
Householder+QR method
Conjugate gradient (CG) method
Davidson method
Even if basis functions are localized in real space, Gram-Shmidt (GS) type method is needed tosatisfy orthonormality among eigenstates, which results in O(N3) for the computational time.

O(N) method - can be achieved in exchange for accuracy.
O(N) Krylov subspace method,DC, DM, OM methods, etc..

O(N2~) method
Is it possible to develop O(N2~) methods without
introducing approximations? → No more GS process.

対角化: 直交化を避けることは可能か？
1. Green関数法

2. 密度行列法
LNV, PRB 47, 10891 (1993)

3. 反復法

低次スケーリング数値厳密対角化法
1. 必要な密度行列要素のみを留数定理で計算

2. 必要なGreen関数成分のみを“入れ子分解”法とLDLT分解から導出した漸化式で計算

特徴
計算コスト: ３乗以下
計算精度:
数値厳密
適用性:
絶縁体&金属を区別しない次元性:
1D, 2D, 3D全てに適用可能
並列性:
高い並列性を有している基底関数:
局在基底(PAO, FEM, FD等)PBR 75, 035123 (2007).
PRB 82, 075131 (2010).

Fermi関数の連分数表示
PBR 75, 035123 (2007).

Hu et al., JCP 133, 101106 (2010)
Karrasch et al., PRB 82, 125114 (2010).
Lin Lin et al., Chinese Annals of
Mathematics (CAM), Ser. B.

連分数展開と一般化固有値問題
連分数展開は一般化固有値問題に変換できる。→

極の計算は固有値問題に

極の分布と留数積分

連分数展開の極は虚数軸にある。

極の分布は非常に特徴的。最初は2πの間隔(松原極と同じ)。
実軸から離れるに従い、急激に間隔が大きくなる。→ グ
リーン関数の構造と一致し、収束が高速化されることが期待される。

Convergence of ρ w.r.t. poles
The calculation of ρ can be expressed by a contour integration:

Lehmann rep.

The analysis shows that the number of poles for each eigenstate for a sufficientconvergence does not depend on the size of system if the spectrum radius does notchange. → The scaling property is governed by the calculation of G.

Convergence property of the contour integration
Total energy of aluminum as a function of the number of poles
by a recursion method at 600 K.
Nicholson et al., PRB
50, 14686 1994

The energy completely
converges using only 80
poles within double
precision.

How can Green’s funtion be evaluated ?
・ The Green’s function is the inverse of a sparse matrix (ZS-H).

・ Selected elements of G(Z), which correspond to non-zero elementsof the overlap matrix S, are needed to calculate physical properties.

Our idea
1. Nested dissection of (ZS-H)
2. LDLT decomposition for the structured matrix→ a set of recurrence relations
TO, PRB 82, 075131 (2010)

入れ子分解(nested dissection)法
George, SIAM J. Numer. Anal. 10, 345 (1973).

相互作用は二元ツリー構
造にマップされる。

Nested dissection of a sparse matrix(i)-(v)のプロセスをドメイン毎に再帰的に行い、系のnested dissectionを行う。

(i) オーダリング:
各軸方向で基底座標に基づき基底関数を整列。

(ii) スクリーニング:
長いテールを持った基底関数をセパレーターに押し込む。

(iii) 出発核:
最も重なり数が少ない基底関数を出発核とする。

(iv) 核成長:
出発核から徐々に核を成長させ、|N0-N1|+ Nsを最小化する核の成長段階を探す。

(v) Dissection:

N0: ドメイン0内の基底数
N1: ドメイン1内の基底数
Ns: セパレーター内の基底数(i)-(iv)の過程を全ての軸方向に行い、 |N0-N1| + Nsを最初化する軸方向でdissectionを実施。

Square lattice for the nested dissection

Schur補元による逆行列計算 #1
A matrix X can be factorized using a Schur complement into a LDLT form.

Then, the inverse of X is given by

Schur補元による逆行列計算 #2
例えば次の構造化された行列に対してLを計算してみる。
Lは漸化式で計算出来ることが分かる。

Schur補元による逆行列計算 #3
Lを漸化式により計算

Lを利用してさらに漸化式により逆行列を計算

LDLT分解から導出された漸化式
PRB 82, 075131 (2010)→ O(N7/3)
次元性により計算オー
ダーが変わる。

1D
2D
3D

O(N(log2N)2)O(N2)O(N7/3)

計算アルゴリズム
loop for finding μ {loop of poles for (1) contour int. { → MPI
parallelization(2) Inverse calculation by the recurrence
relations.
→ OpenMP parallelization}(3) Calculate the total number of electrons(4) Correct μ by the Muller method (準ニュートン法)}

計算時間の実測
1D
2D
3D

For the inverse calculation ofmodel TB hamiltonians

従来法との計算精度の比較
Total energy (Hartree)
Conventional -4130.938589871644
New method -4130.938589871899

DNA, 650 atoms
700K, 80 poles
SZ basis set

並列効率

1 SCF on Cray-XT5

まとめ• オーダーN法の定式化は通常、KSハミルトニアンの疎行
列性に基づいている。疎行列構造を与える様々な方法が知られているが、本講義では原子様数値基底法を紹介した。
• オーダーN法の超並列化を実現する手法として、修正再帰二分法と慣性モーメントテンソル法に基づく領域分割法を議論した。
• オーダーN法は近似計算であるため、応用計算の際に精度に対する十分な検証が必要である。検証方法をユー
ザー間で共有していくことで、今後、応用計算が進展していくであろう。