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

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

第9回
2016年
6月16日

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

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

工学分野での第一原理計算へ向けて
103 – 106 atom

鉄鋼材料

時間スケールDFT calculations of thousands atoms
is still a grand challenge.

O(N3)102 atomLow-orderDNA

多くの適用例がある.
マテリアルデザインも試みられ、多くの成功事例がある。

系のサイズ

リチウム電池大規模計算の必要性-物性と材料特性 材料特性は完全結晶や分子が持つ固有の物性と不均一な組織構造 (粒界、異種物質界面、析出等)から生じる二次的な材料特性の複合的要因で決定される。
 工業的に使用される多くの材料は結晶構造と高次構造を適切に選択することで材料特性の向上がはかられる。
 例えば高性能磁石であるNd-Fe-B系磁石の保磁力はその結晶構造、粒径、界面構造で決定される。

http://ev.nissan.co.jp/LEAF/P
ERFORMANCE/ムーアの法則によれば…
10 Zeta Flops (Zeta = 1021)

Top 500
http://www.top500.org/

2035

20年後に計算できるサイズの見積もり #1
今後もムーアの法則が成り立ち、並列効率がシステムサイズや計算機の並列度に依存しないと仮定すると (理想的な条件下)、20年後に計算できるサイズの上限値を見積もることが可能である。

計算量cN

計算能力

b exp(aT )

p

c:
N:
p:
b:
a:

定数
システムサイズ
計算のオーダー定数定数系のサイズをα倍: N → α N の場合を考える。
計算能力の向上比が計算量の増大比と等しいとおけるので、c( N ) p b exp(a(T  T ))pcNb exp(aT )⇒   exp(aT / p)

14か月で計算能力が二倍とすれば a=ln(2)/14、20年後⇒ ΔT = 240 か月。

20年後に計算できるサイズの見積もり #2
方法

計算コスト
(N: 系のサイズ)

10倍の系を扱う
場合の計算コストの増大比

20年後に計算できるサイズ

CCSD(T)

N710000000

5.5CCSD, MRMP

N61000000

7.2

MP2

N410000

20HF, DFT

N3100052O(N) DFT

N110145000

今後の２０年間、Mooreの法則が続いたしても、低次スケーリング法を用いない限り、大規模計算はかなり困難である。

計算のオーダーとは？
計算の複雑さを特徴づける量は計算のオーダーとして知られている。取り扱う系のサイズをNとした場合に、計算量がこのNに対してどの様な依存関係があるのか、O表記を使って記述する。例えば、Nに対して計算量が3乗で増加する場合には3O(N )と書き、オーダーN三乗と読む。この場合、Nのサイズが2Nとなったら、(2N)3=8N3なので、計算量は8倍になることが分かる。
同じ目的の計算を行う場合に、計算のオーダーが小さい手法
が良い計算手法である。 (計算精度が保持されていることが前提である)

内積の計算ベクトルを一次元の配列によって表現し、内積計算を行う。
ベクトルA,Bの初期化
for ループによる内積計算
計算オーダーは?
forループはi=0からN-1までの計算を繰り返す。

N回の乗算
N回の加算

O(N)

つまりベクトルの内積計算はO(N)の計算である。

行列ベクトル積の計算
ベクトルを一次元の配列によって、行列を二次元の配列によって表現し、行列ベクトル積の計算を行う。
初期化

行列ベクトル積

行列ベクトル積
計算オーダーは？
jのループ内の計算は内積計算と同じで、この計算をN回繰り返しているので、N×N回の乗算
2)
O(NN×N回の加算

行列ベクトル積計算はO(N2)の計算である。

行列積の計算
行列を二次元の配列によって表現し、行列積の計算を行う。
この計算は以下の様な三重forループで書くことが出来る。
計算オーダーは？
kのループ内の計算は内積計算と同じで、この計算をN×N回繰り返しているので、N2×N回の乗算
N2×N回の加算

O(N3)
上記は主要部分のみ記載

行列積計算はO(N3)の計算である。

極端な事例: 連立一次方程式の解法
科学技術計算において連立一次方程式の解法は至る所に現れる。計算アルゴリズムによって極めて大きな計算量の差異が生じる。

Cramerの方法

O(N!)

Gaussの方法

O(N3)上記の二つの方法は良く知られた手法であるが、その計算量の違いは、極めて大きなものである。
極端な例であるが、アルゴリズムの重要性を示す一つの例である。

オーダーN第一原理電子状態計算• 準局所近似である局所密度近似(LDA)や一般化勾配近似(GGA）に基づく密度汎関数理論(DFT)計算の計算オーダーはO(N3)である。
• 構造材料、電池材料、磁石材料などは数十nmの複合的な組織構造が材料特性を決めており、可能であれば数千から数十万原子系の第一原理計算を実行することが望ましい。
• これらの大規模計算を実現するためには、従来のO(N3)法を超並列化するか、もしくは計算量が系のサイズに比例するオーダーN法を使用する場合がある。

本講義では後者のオーダーN法に関して議論を行う。
オーダーN法に関する総説論文
 ” Linear scaling electronic structure methods”, S. Goedecker, Rev. Mod. Phys. 71, 1085
(1999).
 “O(N) Methods in electronic structure calculations”, D.R. Bowler and T. Miyazaki, Reports
on Progress in Physics 75, 036503 (2012).

密度汎関数理論
基底状態のエネルギーは密度の汎関数で記述できる.
Hohenberg and Kohn (1964), Levy (1979) (N-表示可能条件).

基底状態の多体問題は有効ポテンシャルに対する一体問題に帰着される. Kohn-Sham, (1965).

KS方程式の数学的構造
3次元の連立非線形偏微分方程式を自己無撞着に解く必要がある。
O(N3)
O(N2)

O(N log(N))

O(N)
局所密度近似(LDA)や一般化勾配近似(GGA)の場合における計算オーダーを赤字で示す。
最も計算オーダーが大
きい部分は固有値問題の解法であり、系のサイズが大きくなるにつれ、計算オーダーはO(N3)に漸近する。

一次縮約密度行列nの関数
電子密度 ρ(r) は一次縮約密度行列nから次式で計算される。

 (r )   nij  j (r ) i (r )i, j密度汎関数は一次縮約密度行列nの関数として書き直すことが出来る。

Etot [  , n]  Tr(nH kin )   dr  (r )vext (r )1
 (r )  (r ')  drdr ' Exc [  ]2| r r'|

もし基底関数 χ が実空間において局在しているならば、全エネルギーの計算に必要なnの要素数はO(N)となる。
密度行列nの必要な要素のみが計算可能であれば、計算オーダーを低減させることが出来る。

オーダーN DFT計算への二つの道筋
密度汎関数の通常の表示形式は密度とKS軌道を用いたものである。
近似を導入することなしに変数変換を導入することで、密度行列もしくはWannier関数を用いた表示形に変換できる。
通常の表示

密度行列表示ψ: KS電子軌道
ρ: 電荷密度 : Wannier関数
n: 密度行列
Wannier関数表示

密度行列及びWannier関数の実空間の局所性を利用することで、計算オーダーを低減させることが可能となる。

Wannier関数と密度行列
Wannier関数 はBloch関数ψのunitary変換から得られる。
occV|   dk U  |  k  exp(ik R)
3 (2 ) BZ m

バンドギャップを有する場合

密度行列はBloch関数ψの密度演算子への射影から得られる。

n(r , r ')   nij , Rn i (r   i )  j (r '  j  Rn )n

離散化された次式を利用して上記の連続関数が得られる。
occ
1nij , Rn   dk  exp(ik Rn )ci ,k c j ,kVB BWannier関数の局所性
O-2px in PbTiO3

Orbital in Aluminum

指数関数減衰J.Battacharjee andU.W.Waghmare, PRB73, 121102 (2006)

べき減衰

半導体及び絶縁体においてはWannier関数は指数関数減衰し、金属ではべき減衰する。1D系に対しては数学的に詳細な解析(He and Vanderbilt,PRL 86, 5341)が行われている。一般の場合に対する条件付きの数学的証明はBrouder et al., PRL 98, 046402で議論されている。

密度行列の局所性

有限ギャップ系
指数関数減衰金属
T=0 べき減衰0<T 指数関数減衰
D.R.Bowler et al.,Modell.Siml.Mater.Sci.Eng.
5, 199 (1997)

T=0Kにおいて、半導体及び絶縁体においては密度行列は指数関数減衰し、金属ではべき減衰する。有限温度では金属においても指数減衰する。
減衰特性に対する数学的な議論は Ismail-Beigi et al, PRL 82, 2127を参照のこと。

様々なオーダーN法
Wannier 関数 (WF)

変分法 (V)

密度行列 (DM)

摂動法 (P)

二つの最適変数とその最適化手法の組み合わせから少なくとも4種類の方法論が考えられる。

WF+V

WF+P

Orbitalminimization

Hoshi
Mostofiby Galli, Parrinello,Ordejon, Tsuchida

DM+V

DM+P

Density matrix Krylov subspaceby Li and Daw
Divide-conquer
Recursion
Fermi operatorオーダーN・DFTコードの開発状況
Conquest (DM) Bowler(ロンドン), Gillan(ロンドン),
宮崎(物材機構), 大野(物材機構)
Siesta(OM) Ordejon et al.(スペイン)

ONETEP (DM) Hayne et al.(イギリス)

OpenMX (Krylov) 尾崎 et al. (東大)FEMTECK (OM) 土田(産総研)
AkaiKKR (screened KKR) 赤井 et al.(東大)

電子の近視性に基づくオーダーN法

仮定
各原子の“電子状態”は化学的環境をもたらす周辺領域の原子配置によって大部分は決定されている。

分割統治法による収束特性切り出したクラスターを解く。 →分割統治法

W.Yang, PRL 66, 1438 (1991)

絶縁体、半導体金属金属系は収束に対して大きなクラスターサイズを必要とする。
→ 直接的な分割統治法の適用は困難。

Krylov部分空間に基づくオーダーN法• Lanczos変換に基づく方法
R. Haydock, V. Heine, and M. J. Kelly, J. Phys. C 5, 2845(1972); R. Haydock, Solid State Phys. 35, 216 (1980).
T. Ozaki, Phys. Rev. B 59, 16061 (1999); T. Ozaki, M.
Aoki, and D. G. Pettifor, ibid. 61, 7972 (2000).

• Two-sided block Lanczos変換に基づく方法
T. Ozaki and K. Terakura, Phys. Rev. B 64, 195126 (2001).
T. Ozaki, Phys. Rev. B 64, 195110 (2001).

• Arnoldi変換に基づく方法
T. Ozaki, Phys. Rev. B 74, 245101 (2006).

量子力学における固有値問題(1) シュレディンガー方程式の解ψを(3) 代入し、展開すれば一電子エネルギーは次式で与えられる。

(2) 基底関数で展開する。

(4) 係数cに対して正規化の条件の下、最小化すると、(5) 行列要素は積分によって求められる。

,

エルミート行列の固有値問題
二つの固有値εi,εjに属する固有ベクトルをci, cjとすると、・・・ (1)・・・ (2)(1)の左辺から<cj|を作用させ、・・・ (3)(2)の複素共役を取り、転置すると・・・ (4)(4)の右辺から|ci>を作用させて、(5)-(3)で次式が得られる。
i=jであるとき、つまりεは実数である。

・・・ (5)・・・ (6)i≠jであるとき、 異なる固有値に属する固有ベクトルは直交する。

固有値問題の“図解”行列をベクトルに掛けるとは?
cからHcに変わる時、cはHによって回転と伸張を受けている。

cΘ

ε|c|
Hc
一般の場合にはn次元ベクトル空間に拡張される。

この様な見方をすると固有値問題とはHを掛け算しても回転しないベクトルを探すことに対応している。また固有値とは伸張の度合いになる。

行列のブラケット表示
固有値問題の解を用いて、・・・ (1)行列は次式で書くことが出来る・・・ (2)

また単位行列は次式で与えられる。
・・・ (3)(2)を(1)の左辺に代入すると確かに右辺になることが確認できる。

固有値問題とグリーン関数グリーン関数は次式で定義される。

・・・ (1)
Hの固有値と固有ベクトルを用いてグリーン関数は次式で書くことが出来る。

・・・ (2)(2)を(1)に代入すると確かにグリーン関数の定義が満たされている。

べき乗法
ランダムベクトルv0に行列Hを繰り返し掛け算していくと収束解が得られるだろうか?
初期ベクトルv0は固有ベクトルの線形結合で書ける。

v1 = Hv0
v2 = Hv1
・・・
vn = Hvn-1
v∞ → ???
つまり絶対値最大の固有ベクトルに収束していく。ただし絶対値最大の固有値の近くに別の固有値があると収束が遅くなることも右の考察から分かる。

v0にHをn回、掛け算する。

ここでε0は絶対値最大の固有値である。

Lanczos法
べき乗法は絶対値最大の固有値付近に他の固有値が存在する場合、収束が遅くなる。このような場合に対して、探索空間としてこれまでに得られた部分空間に対する直交ベクトルを選ぶことで、 Lanczos法は速やかに収束解を与える。

アイデア
行列を三重対角化し、その三重対角行列の問題を解く。

Cornelius Lanczos,
1893-1974

いかにしてUは見つけられるのだろうか?

Quoted from http://guettel.com/lanczos/

Lanczos法の導出 #1
HTD=U†HUを明示的に書き下すとこれをさらに各列毎に書く。
1

この方程式系は一般にunとun-1が分かっているとun+1が分かる形に書き換えることが出来る。
+1

Lanczos法の導出 #2
したがってu0が与えられると、漸化式を解くことで、unが計算できる。
これをアルゴリズムとして書くと次のようになる。

=+1+1

Lanczos法の導出 #3
この手続きで計算したun+1が部分空間{u0, u1,・・, un}に直交していることは以下のように数学的帰納法を用いて証明できる。
ランチョス法の手続きから+1・・・ (1)まずn=0の場合に<u0|u0>=1。またu1はu0に直交していることが確認できる。
この時、nまでは正規直交系であると仮定する。

(1)の左から<uk|を掛けて(k=0～n)、+1

k=n+1の場合にはであることが分かる。

ゆえにun+1は部分空間{u0, u1,・・, un}に直交している。したがってLanczos法によって正規直交系 (Krylov部分空間) を作り出すことが可能である。

Lanczos法とGreen関数の関係 #1
ランチョス変換によって得られた三重対角行列を用いて、グリーン関数の有用な表式が得られる。

Lanczos法とGreen関数の関係 #2
行列式は三重対角行列の場合、漸化式で記述される。

一般化するとこの漸化式を用いてGreen関数の対角項を評価する。

ラプラス展開最終的に次式を得る。

Green関数と物理量Green関数の虚数部を計算してみる。

虚数部を積分し、虚数部をプロットすると以下の様になる。
したがって、Green関数の対角項の虚部は状態密度を与える。

Lanczos法に基づくリカージョン法
アルゴリズム
1. 原子サイトiを初期ベクトルとする。
2. Lanczos変換を行い、H行列を三重対角化する。
3. 繰り返しステップNLでLanczosステップを打ち切る。
4. Green関数の連分数表示を用いて、原子iの局所状態密度を計算する。
5. ステップ1に戻り、次の原子サイトに対して同一のプロセスを繰り返す。

計算オーダー
1. Hが密行列で、Lanczosステップが原子数に比例2. Hが疎行列で、Lanczosステップが原子数に比例
3. Hが疎行列で、Lanczosステップが定数
4. 各原子毎にLanczosプロセスの飛び移り領域を制限→ O(N4)→ O(N3)→ O(N2)→ O(N)

解析的に解ける例
3-foldのベーテ格子

K-foldのベーテ格子: オンサイトのエネルギー ε, 最近接ホッピング積分 h.
あるサイトからスタートすると、n番目めのLanczosベク
トルは次式で与えられる。

Lanczos基底表示での行列要素

したがって連分数表示より着目するサイトの局所グリーン関数が得られる。

入れ子構造に着目すると次式が得られる。

Lanczosベクトルの発展
Lanczosベクトルの発展は電子のホッピング過程として理解できる。下記の例は正方格子におけるベクトルの発展過程である。

|L0>

|L1>

|L5>

周辺環境に関する情報は近接サイトから順番に含まれていく。その結果として、部分空間の次元が低減される。

Krylov部分空間とモーメントの関係
Lanczosベクトル(Krylov部分空間ベクトル)は次のKrylov部分空間を張っている。

U K  {| u0 , H | u0 , H 2 | u0 , H 3 | u0 ,....., H q | u0 }
Krylov部分空間ベクトルでハミルトニアンを表示する。

H mn  um | H | un K u0 | H m n 1 | u0  u0 |  ( m n1) | u0 モーメントの定義 ( p )  c p c† cc† Hcc† Hc  c† Hcc†Hp

q次のKrylov部分空間は(2q+1)次までのモーメントの情報を含んでいる。

一方、|ε/Z|<1の条件下で、グリーン関数はモーメントで表示可能である。
つまりクリロフ部分空間の次元はグリーン関数の精度に1対1対応する。したがって、局所電子構造は周辺の原子配置でほぼ決まるという化学的直観は数学的にも正当化される。

クリロフ部分空間法の第一原理計算への拡張
一般の局在基底を用いてKS軌道を展開し、KS方程式を解く場合には次の一般化固有値問題に帰着する。強結合モデルと異なり、重なり積分を明示的に考慮する必要がある。

Hc    Sc
本講義では二つの方法を紹介する。

• Two-sided block Lanczos変換に基づく方法
T. Ozaki and K. Terakura, Phys. Rev. B 64, 195126 (2001).
T. Ozaki, Phys. Rev. B 64, 195110 (2001).

• Arnoldi変換に基づく方法
T. Ozaki, Phys. Rev. B 74, 245101 (2006).

Two-sided block Lanczos変換に基づく方法 #1

G(Z )(ZS  H )  I
1を次式の様に変形し、G( Z ) S ( ZI  S H )  G (Z )(ZI  H ')  IL

非エルミート行列であるH’に対してblock化したtwo-sidedLanczos変換を行う。

~

ここで双対基底関数|iα>は次式で定義され、本来の基底と双直交性を満たす。

Two-sided block Lanczos変換に基づく方法 #2
Two-side block Lanczos変換により行列Hはブロック三重対角化され、グリーン関数のブロック対角要素とブロック非対角要素は次式で与えられる。
ブロック対角要素

ブロック非対角要素

下線はブロック要素であることを、LはLanczos基底表示であることを示す。

密度行列は次式から計算される。

nL   1 Im  dEGL (E  i 0 ) f ( E )本来の基底表示での密度行列nは逆変換から得られる。

nij ,   n0Ln , U nk , S jk1n,kArnoldi変換に基づく方法 #1
ハミルトニアンのブロック三重対角性を保持するのではなく、クリロ
フ部分空間ベクトルのS-直交性だけを保持する。
Green関数のモーメント表示

G( Z )  ( ZS  H )1  c( ZI   )1 c† ( p)p 1
Zp 0

したがってS-1Hを掛け算することでクリロフ部分空間が生成される。

アルゴリズムQ  S 1

Arnoldi変換に基づく方法 #2• 数値安定性を向上させるために、Arnoldi変換によるKrylov部分空間ベクトルの生成は第一SCFステップでのみ行う。
• 直交性を高精度で保持するために対角化による直交変換を行う。

• 各原子毎にSとHの実空間領域での打ち切りを導入することで、計算オーダーはO(N)となる。

Arnoldi変換に基づく方法 #3
埋め込まれたクラスター問題
Krylov部分空間表示において、標準固有値問題となる。

ここで HK は短距離と長距離の寄与からなる。

Green: コア部分
Yellow: バッファー部分• 埋め込まれたクラスターは他からのクーロン相互作用を受けている。
• 共通化学ポテンシャルを課すことで、電荷移動が許される。

Arnoldi変換に基づく方法 #4
アルゴリズム1. 有限クラスターを構築
2. クリロフ部分空間ベクトルを生成
3. 埋め込みクラスター問題を解く
4. 共通の化学ポテンシャルを決定
5. 全電荷密度を生成6. 全系に対するポアソン方程式の解法
7. SCF収束までステップ3に戻る

SCF収束性の比較

fcc Alの場合

TO, PRB 74, 245101 (2006).

Proposed (Arnoldi法)はSCF計算においてもロバストであるが、Recursion法(two-sided Lanczos法)は数値誤差の影響を受けやすく、SCF収束が困難。

全エネルギーの収束特性
精度と効率はクラスターサイズとKrylov部分空間の次元で制御。

全エネルギーの収束性(半導体、絶縁体、金属、分子系)

収束特性は単純ではないが、クラスターサイズで決まる次元の30%程度で収束。

分割統治法とKrylov部分空間法の比較

カッコ内の数字
1番目: 緑の部分内の原子数2番目: 黄色の部分の原子数• 金属に対してはKrylov部分空間法はDC法より高速である。
• 共有結合性の強い系では両者は同等である。

計算時間の比較
各原子に割り当てられたクラスターのサイズは系のサイズには依存しない。
したがって計算オーダーは厳密にO(N)となる。

Carbon diamond(a) 8 cpus(b) 32 cpus

Al表面上に担持されたカーボンナノチューブ

電荷分布に対し、通常の対角法とオーダーN法は良く一致している。

Al表面からCNT側面に電子移動し、電気二重層を形成。
Al表面の深さ方向にフリーデル振動的な振る舞いが見られる。

マリケン(スピン)電荷の比較

電荷及びスピン共に、通常の対角法とオーダーN法は良く一致している。

構造最適化の精度検証 #1
O(N)Krylov部分空間法で計算される全エネルギーの一次微分は近似的なものであり、厳密には全エネルギーの微分となっていない。近似一次微分で構造最適化が可能か検証する。
# of atoms in the truncated cluster

159275381525

ベンチマーク計算は62炭素原子、2ボロン原子からなるダイヤモンド格子。2つのボロン原子は最隣接原子として置換型で導入。

O(N)Krylov部分空間法の計算精度を向上させるに従い、力の収束性も向上していることが分かる。

構造最適化の精度検証 #2
ベンチマーク計算は62炭素原子、2ボロン原子からなるダイヤモンド格子。2つのボロン原子は最隣接原子として置換型で導入。

Time (s)/MD

139926

1890

2841

1421

計算コストと計算精度を考慮し、実際の応用計算を実施する際にはKrylov3の計算条件が推奨される。

固有エネルギーと波動関数の計算
1.
2.

O(N)法を用いてSCF計算を実施。
O(N)法で得られた電荷密度を利用し、one-shotの対角化計算を行う。

HO (guanine)

LU (cytosine)

状態密度の比較から、実用上、本手法は十分な精度を持っていると判断される。

まとめ• 密度汎関数理論に基づく第一原理電子状態計算の計算
オーダーは通常、O(N3)である。
• より現実を反映したシミュレーションを実現するためには数千原子以上からなる複雑で大規模な系の取り扱いが必要である。

• その様なシミュレーションを実現する手法として、計算コストが系のサイズに比例したオーダーN法がある。
• 講義ではKrylov部分空間に基づくオーダーN法に関して紹介し、さらに第一原理計算への拡張を議論した。