実空間差分法による大規模第一原理電子状態計算
岩田潤一, 白石賢ニ, 押山淳
筑波大学計算科学研究センター,筑波大学数理物質科学研究科
iwata@comas.frsc.tsukuba.ac.jp
概要：物性物理とは原子が多数寄り集まった系の多体問題の研究である。密度汎関数理論により現在数百原子系の電子の振る舞いを量子力学の第一原理に基づいて記述することが可能となっているが、世の中には多種多様な物質が存在し、より大規模な系を扱いたいという要求が常にある。本発表では実空間差分法による超並列計算機上での大規模第一原理計算手法の開発とその応用を紹介する。

1

はじめにためには従来の計算手法を、テラフロップス、ペ物性物理学の研究対象は原子核と電子からなるタフロップス、・・・と今後益々進歩して行くであ多体系である。電子は通常量子力学的に扱う必要ろう計算機の能力をフル活用する形でコーディンがあり、その運動は Schrödinger 方程式によってグを行っていくことが重要と考えられる。本発表規定される。しかしながらＮ体系の Schrödinger

では、今後の計算機の発展において必須のアーキ

方程式を直接数値的に解くことはＮが非常に小さテクチャと考えられる超並列計算機上での実空間い場合を除いてほぼ不可能であり、実際には何ら差分法による第一原理計算手法の開発とその応用かの近似を導入する必要がある。

について紹介する。

1964 年、Hohenberg-Kohn-Sham らによって「Ｎ体系の Schrödinger 方程式の基底状態を求め2

実空間差分法る問題を、Ｎ本の１体系の Schrödinger 方程式をDFT に基づく第一原理計算は、数学的には解く問題に焼き直す」という密度汎関数理論

Kohn-Sham 方程式と呼ばれる３次元偏微分方程(DFT)が提唱され、これにより O(N3)という比較的式の固有値問題の解を固有値の小さいものから順少ない計算コストで大規模系を扱う道が開かれたに N 個求めるという問題である。

[1,2]。もちろん DFT においても実際の計算を行う

 1 2

 − 2 ∇ + v[ ρ ](r ) + vˆNL  φn (r ) = ε nφn (r )にはある種の近似が必要であるが、現在広く用いられている局所密度近似および一般化密度勾配近似はその単純さに比して非常に定量性のある結果を与えることが知られており、これによって多体系の電子状態を、経験的なパラメータ等を一切使用せずに求める、第一原理計算と呼ばれる物性研究の手法が確立した。現在では物性物理学のみならず量子化学、材料科学、分子生物学などにおいても第一原理計算による物性予測は高い信頼性を持って受け入れられるようになっている。
近年、ナノサイエンスという言葉が大きく取り沙汰されるようになり、物性物理においてもナノ

ここでNρ (r ) = ∑ φn (r )(2)n =1は電子密度と呼ばれるものである。(1)式の左辺の演算子（ハミルトニアン）が電子密度に依存しているので通常の固有値問題と異なり、(1)式はハミルトニアンと固有ベクトル（波動関数）が自己無撞着となるように解かなければならない。ハミルトニアンの第３項目は非局所擬ポテンシャルと呼ばれるもので、波動関数への作用はNk上という系‐が注目されるようになってきている。

vˆNLφn (r ) = ∑ vk (r ) ∫またバイオサイエンスの分野においてもタンパクとなる[3,4]。

構造‐原子数にして数千〜数万個あるいはそれ以(1)k =1r ′ ≤ Rcdr′vk* (r′)φn (r′)(3)

質その他の巨大分子が興味の対象であり、やはり従来、物性の第一原理計算では固体のような周数千〜数万原子オーダーの系を扱う必要が出てく期構造を扱う関係で平面波基底を用いるのが一般る。このような超大規模系を第一原理計算で扱う的であった。しかし近年では周期的でない系を扱う場合も多く、また平面波基底を用いる場合ハミ

表１要となり、並列化した場合に非常に多くの全対全部分対角化

Time (sec.)1046

通信を行うことになるので、並列化の効率という

Gram-Schmidt
CGTOTAL

178.8
350
1661ルトニアンの演算に高速フーリエ変換(FFT)が必面で不利となる。そこで我々は（１）周期系、孤立系など様々な境界条件に柔軟に対応できる、

MFLOPS/node
40002600
701300（２）FFT の必要がなく並列性に優れている、という２つの長所を併せ持つ実空間差分第一原理計を用いて行った。DFT 全体の計算量は部分対角化

算法[5]に着目しコード開発を進めている。

および Gram-Schmidt 直交化の O(N3)でスケール実空間法では(1)式を３次元格子上で離散化し、されるが、これらのルーチンはともに大部分の演波動関数は各格子点上の値で表現する。ラプラシ算が BLAS3 レベルの行列×行列積の形で処理でアンは高次差分で近似し、通常 4 次〜６次（9 点きるため非常に高い演算性能を引き出すことがで〜13 点）公式を用いる。並列化は空間格子をノーき、実際これらの部分の計算時間は O(N2)の CGド数分の小領域に分割することで行う。必要なノと大差ない程度に納まっている。

ード間通信は、主に差分および非局所擬ポテンシ
ャルの演算のための隣接通信と総和演算のための4まとめ通信である。また部分対角化（後述）ルーチンに実空間法による超並列計算機上での物性第一原おいて Scalapack の実対称行列（またはエルミー理計算コードの開発を行い、現在様々なテスト計ト行列）対角化ルーチンを使用している。プログ算を経て、具体的な物性物理の問題への応用を始ラム言語には Fortran90 を使用し、並列化は MPIめている。DFT 計算の演算量は O(N3)でスケールを用いて行っている。

されるが、その最も計算の重い部分は、最も計算機の性能を引き出せる演算（BLAS3）で処理でき3実際の計算とパフォーマンスるため、並列化および BLAS3 を使う等、計算機実空間法による第一原理計算の流れは次のよう能力をフル活用するコーディングにより、当初想になる。まず適当な初期密度を与えてハミルトニ像されていたよりも遥かに巨大な系に対する第一アンを作り、(1)式に従って共役勾配法(CG)あるい原理計算の適用可能性が見えてきた。ペタフロッは最小残差法(RMM)により波動関数の最適化をプス級の計算機が完成すれば数万原子からなる超行う。ただしハミルトニアンと波動関数が自己無大規模系の物性予測が可能になると思われる。

撞着となるように解を求める必要があるので、CGの反復は４〜５回程度で強制的に打ち切る。その後 Gram-Schmidt 法によりＮ本の波動関数の規格直交化を行い、さらにＮ本のベクトルだけで張られる空間内でハミルトニアンの対角化を行う（部分対角化）。このようにして 少しだけ 更新された波動関数で電子密度を作り直し、ハミルトニアンを作り直し、再び CG に戻り、この過程を自己無撞着解に到達するまで繰り返す。この反復計算を「自己無撞着なポテンシャル場を求める」という意味で SCF(Self-Consistent Field)計算と呼ぶ。
表１は Si 原子 4096 個からなるバルク固体のDFT 計算の SCF 反復１ステップにかかる時間を示したものである。収束解を得るまでにこの反復を約 50 回程度繰り返すことになる。計算規模は格子点数 884736 点、求める波動関数の本数は 8196本である。計算は筑波大学計算科学研究センター
PACS-CS（LV Xeon EM64T 2.8GHz）256 ノード

参考文献[1] P. Hohenberg and W. Kohn,「InhomogeneousElectron Gas」, Phys. Rev. 136, B864, (1964).
[2] W. Kohn and L. J. Sham,「Self-ConsistentEquations
Including
ExchangeandCorrelation Effects」, Phys. Rev. 140, A1133,(1965).
[3] N. Troullier and J. L. Martins,「Efficientpseudopotentials
forplane-wave
calculations」, Phys. Rev. B 43, 1993, (1991).
[4] L. Kleinman and D. M. Bylander,「EfficaciousFormforModelPseudopotentials」, Phys. Rev. Lett. 48, 1425,(1982).
[5] J. R. Chelikowsky, et al., 「 Higher-orderfinite-difference pseudopotential method」,
Phys. Rev. B 50, 11355, (1994).