スーパーコンピュータによる量子力学第一原理シミュレーションの大規模化・高精度化の試み
広瀬喜久治

後藤英和

稲垣耕司

小野倫也

大阪大学 大学院工学研究科

１．はじめに量子力学に基づくコンピュータシミュレーションは、化学反応素過程の解析法として実用化の域に達しつつある。特に「第一原理分子動力学シミュレーション」は、スーパーコンピュータの能力向上によりモデルの大規模化が可能になり、安定構造探索など静的現象の解析にとどまらず、化学反応素過程など動的現象の取り扱いも可能となってきた。これまでに我々が開発してきた、スーパーコンピュータ上で動作する高速シミュレーションプログラムにより、原子数十個程度のモデルならば実用的な時間で分子動力学シミュレーションを実行できるまでになってきた1)2)。これによって、我々の研究グループの目的である「原子レベルの表面創成」の研究に活用できるものとなりつつある。しかしながら、計算能力の向上と共に、より大きなモデルに対しての解析が期待されるため、未だに計算能力は十分とはいえない状況である。
我々は、これまでに開発したプログラムを、平成 15 年 1 月に東北大学情報シナジーセンターに導入された最先端スーパーコンピュータである SX-7 に対応させ、より大規模な計算を高速に実行するためのプログラム開発を行うとともに、まったく新しい原理に基づく計算方法の開発を開始した。本稿では、SX-7 においてその計算能力を発揮させるために必要であったプログラム改良の要点について述べるとともに、従来法による大規模高速シミュレーションの結果について述べ、上記の新しい計算法であるオーダーN プログラム開発の現状もあわせて述べる。
２．SX-7 における高速計算の要点
平成 15 年 1 月に導入された新スーパーコンピュータ SX-7 は導入当初より安定した高いパフォーマンスを発揮している。本項では、この新コンピュータの性能を十分発揮させる上で重要となる項目について述べる。
スーパーコンピュータの特徴は、もちろん数値計算速度が高いことであるが、最近のパソコンレベルのCPUでもGFLOPSオーダーの計算が可能で、これらとの明確な違いは計算能力そのものよりもむしろメモリアクセス性能の桁違いな高速性にある。
一般に大量の計算が必要になりやすい、密に数値が詰まった行列×ベクトル積や行列×行列積の計算には、ブロック化の手法が容易に適用できるため、ベクトルレジスタ（あるいはキャッシュメモリ）が十分あれば、あまりメモリアクセスは重要にはならない3)。一方、われわれが扱うようなスパース行列やスパースベクトルの演算の場合、ループ構造やメモリアクセスが複雑になるため、計算速度はメモリアクセスに制限され、高い演算性能が十分発揮できない。SX-4 では一演算の間に 2 要素のメモリアクセスが可能であったが、SX-7 ではこれが 1 要素に相対的に低減されたため特に問題が大きい。このため、高速化手法を組み込んでパフォーマンスの低下を避けた。
DO N=1,NMAX
DO I=1,NF(N)
WR=PHR(I,N)*WKR(I,N)-PHI(I,N)*WKI(I,N)
WI=PHR(I,N)*WKI(I,N)+PHI(I,N)*WKR(I,N)
TAR(I)=TAR(I)-WI*SYR(I,N)
TAI(I)=TAI(I)+WR*SYR(I,N)
TBR(I)=TBR(I)+WR*SYI(I,N)
TBI(I)=TBI(I)+WI*SYI(I,N)
END DOEND DOPARAMETER (V_LENGTH=256)
DO IOUT=0,INT(KMAX/V_LENGTH)*V_LENGTH-1,V_LENGTHDO N=1,NMAX
IF (IOUT+V_LENGTH.LE.NF(N)) THEN
*VDIR NODEPDO IIN=1,V_LENGTH
I=IIN+IOUT
WR=PHR(I,N)*WKR(I,N)-PHI(I,N)*WKI(I,N)
WI=PHR(I,N)*WKI(I,N)+PHI(I,N)*WKR(I,N)
TAR(I)=TAR(I)-WI*SYR(I,N)
TAI(I)=TAI(I)+WR*SYR(I,N)
TBR(I)=TBR(I)+WR*SYI(I,N)
TBI(I)=TBI(I)+WI*SYI(I,N)
END DOEND IF
END DOEND DO
DO N=1,NMAX*VDIR NODEP*VDIR SHORTLOOPDO I=INT(NF(N)/V_LENGTH)*V_LENGTH+1,NF(N)
WR=PHR(I,N)*WKR(I,N)-PHI(I,N)*WKI(I,N)
WI=PHR(I,N)*WKI(I,N)+PHI(I,N)*WKR(I,N)
TAR(I)=TAR(I)-WI*SYR(I,N)
TAI(I)=TAI(I)+WR*SYR(I,N)
TBR(I)=TBR(I)+WR*SYI(I,N)
TBI(I)=TBI(I)+WI*SYI(I,N)
C
END DOEND DO(a)変更前(b)変更後：最内側ループをベクトル長でループさせる

図 2.1 スパース配列要素を含む計算の高速化の例
３．高速化プログラムの応用例−超純水の原子・電子構造解析と超精密加工プロセスへの応用―
高速化プログラムの応用例として、超純水の原子・電子構造解析について紹介する。
本誌でも既報しているが、我々の研究グループでは、超純水の化学的機能を利用した超精密加工プロセスの研究・開発を行っている1)2)。これは、従来のプロセスのように化学薬品や特殊ガスを使用せず、水の化学的機能だけを用いて加工を行うものであり、環境負荷とコスト面に優れた新しい生産プロセスの実現が期待できる。このような加工プロセスの開発にあたり、我々は、スーパーコンピュータによる高速計算が可能な、量子力学の第一原理に基づく分子動力学シミュレーションプログラムを独自に作成し、加工物表面の化学反応解析を行ってきた。例えば、超純水のみによる電気化学的加工法4)の開発においては、OHのみの化学的作用により金属表面原子が加工されることや、陽極では加工できないSiやAlを陰極とすることにより加工が可能になることをシミュレーションで予測するなどの成果を挙げ、実用プロセスの開発を側面からサポートしてきた。
このように、水には多くの「隠された」化学的機能が潜んでいるものを考えられるが、これを発見し応用するためには、まず水の原子・電子構造や固体表面との相互作用など、水の基礎物性の解明が必要である。ところが、水は各分子が水素結合で繋がった「クラスター構造」を有していると考えられており、水分子 1 個の原子・電子構造を解明するだけでは、「水」の物性を解明することができず、水のクラスター構造を再現できるような多数の分子から成る大規模シミュレーションが必要である。現在
我々は、水のクラスター構造や固体表面との相互作用などの、大規模シミュレーションの実現を目指した基礎研究を行っており、その現状について紹介する。
本研究で用いた第一原理分子動力学シミュレーションプログラムは、32768（32×32×32）個の平面波基底とノルム保存型擬ポテンシャルを用いており、カットオフエネルギは 458.1eVである。また、時間ステップは 1.0 fsとした。スーパーセルは一辺9.312Åの立方体とし、27 個の水分子を等間隔に配置した（図 3.1(a)）。密度は、1.0g/cm3である。まず絶対温度 0 度で、全エネルギーが低下するように構造安定化を行い、その後温度を 50℃まで上昇させた後、徐冷して最安定構造を求めた。
結果を図 3.1(b)に示す。水分子密度の偏りが生じたことがわかる。これは、水素結合によりクラスター構造が生じたためと考え、水素結合の有無を判別し可視化するために、すべての酸素原子と水素原子間のボンドポピュレーションを計算した5)。図 3.2に結果を示す。ボンドポピュレーションの値は、3 つのグループに分かれていることがわかる。このうち、ボンドポピュレーションの値が 0.4 以上の酸素―水素間結合が:O(a)初期状態:H:電子密度分布(b)構造最適化後:水素結合(c)可視化した水素結合

図 3.1 27 個の水分子の第一原理分子動力学シミュレーション結果

同一水分子内のものであること、

1000800600
4000.1 から 0.2 あたりに分布している水素原子間結合によるものであることがわかった。

1086420度数

結合が隣接する水分子間の酸素―

これを可視化したものを図 3.1(c)に示す。この結果から、ボンドポ水分子内 O−H 結合

隣接水分子間O−H 結合0.0ピュレーションの値により水素結合の有無が判定できることがわか

H2O内のO−H結合
水素結合

0.10.20.30.40.5

ボンドポピュレーション
図 3.2 ボンドポピュレーションの分布

り、シミュレーションで得られたたクラスター構造であることがわかった。また、図 3.3 は、ボンドポピュレーションの値と酸素―水素原子間距離との関係を示したも8原子間距離Å原子間距離(Å)

水分子の構造は、水素結合を介し4水分子内 O−H 結合

320.00.10.20.30.40.5

ボンドポピュレーション

図 3.3 ボンドポピュレーションと酸素−水素原子間距離との関係

シミュレーション結果

g(r)隣接水分子間O−H 結合50でも水素結合の有無が判定できる
水素結合は、隣接する水分子間に61のであり、酸素―水素原子間距離ことがわかる。以上の結果から、

H2O内のO−H結合
水素結合

7SiOH

電子密度分布

実測値r(Å)

図 3.4 酸素原子の動径分布関数の計算結果と実測値5)
図 3.5 クラスター構造を有する水と相互作用する Si(001)水素終端化表面

おいて、原子間距離 1.5〜2.0Å、ボンドポピュレーション 0.1〜0.2 程度の範囲にあると結論した。これは、水素結合の原子間距離が 1.7〜2.0ÅとするX線回折による結果6)とも一致する。また、図 3.4 は、得られた原子配置から酸素原子の動径分布関数を求めた結果であり、これも実測値6)と一致していることがわかる。
以上の結果から、第一原理シミュレーションにより水のクラスター構造がある程度再現できることがわかった。しかし、「水」は水分子が水素結合で繋がった巨大分子であることから、現状の分子数では不十分であることも明白である。現在、さらに分子数を増やしたモデルに対するシミュレーションとプログラムの高速化に取り組んでいる。また、図 3.5 に示すように、クラスター構造を有する水が固体表面と相互作用した場合に生じる化学反応のシミュレーションも行っている。
４． オーダーN プログラム開発の現状
従来の第一原理分子動力学の計算法では、計算量がモデルサイズの３乗に比例（オーダーN3法；O(N3)）するため、モデルサイズを少し大きくするだけで、実用的な時間内での計算が急激に困難となる。このため、計算を本質的な部分だけに絞ってオーダーNで計算する方法（オーダーN法； O(N)法）が提案されており、その手法の実現とプログラム化が期待されている。
4.1 オーダーN アルゴリズムとスーパーコンピュータ
従来法とオーダーN法のアルゴリズムを図 3.1 に示す。オーダーN法のアルゴリズ
ムを述べる前に従来法のものについて簡単に述べる。第一原理分子動力学法は、いう（実際には、密までもなく量子力学のシュレーディンガー方程式に基づいている7)~9)。
度汎関数理論に基づいて一電子近似した、コーン・シャム方程式といわれる方程式が用いられる。
）

HΨi = Ei Ψi
これはH（ハミルトニアン）に対して、固有値（エネルギーレベル）Eiとそれに対応する固有関数（波動関数）Ψiを求める固有値問題であり、モデルに含まれる電子数に応じた複数の固有値・固有関数の計算が必要である。コンピュータ上では、実空間グリッド上で定義したΨiを用いて微分を差分化することによって、上式をエルミット行列の固有値問題に変換して取り扱っている10)。Ψiはベクトルの形態をなすことになるが、複数あるので全体としては行列のようになる。実空間グリッド数は、計算する 3次元空間を例えば 100×100×100 に区切れば 106にもなるため、Ψiは 106、Hは 1012もの要素を持つことになる。幸いなことにHはスパースであるので、固有値問題の解法に逐次近似法を用いれば、計算量はそれほどは多くならない。

初期値の入力（原子の種類、電子数）

初期値の入力（原子の種類、電子数）

電子密度の初期値を設定する

変分関数の初期値を設定する

Schrödinger 方程式

Schrödinger 方程式

電荷密度をもとにハミルトニアンを算出する

エネルギー汎関数の値を減少させる変分関数の修正方向を算出

固有値問題を解き、新しい電荷密度を求める

変分関数を修正する

電荷密度を修正する

No

変分関数が収束?
分子動力学

No電荷密度が収束?
分子動力学

Yes

原子に働く力に従って原子を動かし座標を更新終わり?

No

Yes

No終わり?

原子配置と電子状態が決定
各種物理量の計算

Yes・電子分布・エネルギー・平衡原子間距離・結合の強さ

原子配置と電子状態が決定
各種物理量の計算
・電子分布・エネルギー・平衡原子間距離・結合の強さ

(a)

Yes

原子に働く力に従って原子を動かし座標を更新

(b)

従来法

オーダーN 法図 4.1 計算アルゴリズム
従来法で最終的に計算の負荷となる部分は、スパース行列に密な行列を掛ける作業と、多数のΨiベクトルを直行化する作業の 2 点になる。
これまでのスーパーコンピュータの能力では計算時間の中で前者の占める割合が大きかったが、計算能力の向上によりモデルサイズが大きくとられるようになるにつれて、後者が多くを占めの作業はモデルサイズの３乗に比例して計算量が増えるため、今後のモデルサイズを大きく計算時間るようになってきた。こ従来の計算手法
計算時間∝N3
オーダN
第一原理計算手法
計算時間∝Nしたいという要求に対して、この項がボトルネックとなることが容易に分かる。

原子数N
図 4.2 オーダーN 法による計算時間の短縮

以上のようなことから、計算時間がモデルサイズの１乗にしか比例しないオーダー
N法の必要性が明らかである。これを実現するためにいくつかの方法が提案されているが、我々はMauri とGiulia が提案している方法11)12)に基づいたプログラムを開発している。その方法のキーポイントは、電子が存在する固有ベクトルΨiをすべて求める必要はなく、実はそれが張るベクトル空間のみが分かればよい、という性質と、方程式を固有値問題ではなくエネルギー最小化問題として解く方法を導入することにある。前者の性質から、固有ベクトル群Ψiの代わりに、ベクトル空間が一致する「変分関数」群φiを求めればよいことになる。したがってφi選択には任意性があることになり、φiを空間のほとんどの領域で０の値で、限られた一部にのみ値を持つようなものを選ぶことも可能である。さらに、この「ほぼ０」を「完全な０」とする近似を導入し、φiをスパースベクトルとすることによって、計算量の増加をオーダーNに抑えることが可能になる。本手法に、我々が開発した高精度実空間差分法を適用することによって実用的な計算精度を得ることができると期待できる10)13)。
4.2 開発の現状と計算例
平成 14 年度は手法の有用性（スパース行列を用いた場合にどの程度の誤差が生じるか）を確認するためのプログラムを作成し、本手法の適用性を確かめた。プログラ
ムに関しては、固有値問題解法から最小化問題解法に変更したことで、まったく違うものになった。高速化については、確認作業に最小限必要なものについて行った。実用化を目指した高速化は、平成 15 年度において行うことになっている。
開発した性能確認プログラムを用いて、長く伸びた原子鎖の計算を行った。図 3.3に結果の一例を示す。この図は、1152 原子からなる金の螺旋ナノチューブ (直径 12.7Å，長さ 127Å)の電子状態を計算した例である。局在変分関数を用いた精度低下が問

12.7Å

軸断面図127Å

図 4.3 オーダーN 法により金原子（球）の一次元鎖の電子状態（雲）を計算した例

題となることなく計算が可能であることが見て取れる。今後、計算の高速化や、今回の計算では無視した詳細な項を含めた高精度計算をめざしたプログラムの改良が必要である。
５．おわりに
本稿では、新しく東北大学情報シナジーセンターに導入された SX-７により、高い演算性能を十分発揮させるための要点について述べるとともに、高速・高精度第一原
理分子動力学シミュレーションプログラムの開発と応用について述べた。従来プログ
ラムの応用については、水に潜んでいる多くの化学的機能を発見し活用するために行ったシミュレーションについて述べた。また、開発については、オーダーN 法という従来とはまったく違った計算原理に基づく新しいシミュレーション手法について述べ、その可能性を示した。本プログラムのさらなる改良によって、シミュレーション可能な計算規模の飛躍的な発展が期待される。
本研究は文部省中核的研究拠点(COE)形成プログラムである大阪大学・超精密加工研究拠点における共同研究の一環として行われた。また、計算ソフトウェアの並列化・高速化については東北大学情報シナジーセンター殿に多大な協力をいただいた。

参考文献
1) 広瀬喜久治、後藤英和、稲垣耕司、小野倫也：「共同研究成果」スーパーコンピュータを用いたシミュレーションによる超精密加工プロセスの原子レベルでの解析
Ⅱ、東北大学シナジーセンター 大規模科学計算システム広報 SENAC、34、3(2001)pp.27-35.
2) 広瀬喜久治、後藤英和、稲垣耕司、小野倫也、塚本 茂：スーパーコンピュータを用いたシミュレーションによる超精密加工プロセス・ナノ電子デバイスの解析、
東北大学情報シナジーセンター 大規模科学計算システム広報 SENAC、35、3(2002)pp.59-66.
3) 左近彰一：並列コンピュータ AzusA の高速化技法 ―メモリアクセスの効率化―、
SENAC、36,2(2003).
4) 森 勇藏，広瀬喜久治，後藤英和：超純水のみによる電気化学的加工法の研究、電気加工学会誌、37,84 (2003) pp.1-8.
5) 広瀬喜久治、後藤英和、土屋八郎、森 勇藏、遠藤勝義、山内和人：材料表面現象の第一原理分子動力学シミュレーション−シリコン単結晶(001)表面の水素終
端化反応−、精密工学会誌、60,3 (1994) pp.402-406.
6) 例えば、上平恒：水の分子工学、(1998) 講談社 など
7) W.Kohn and L.J.Sham: Self-consistent Equations Including Exchange andCorrelation Effects, Phys. Rev., 140, 4A (1965) 1133.
8) N.Troullier and L.Martins: Efficient Pseudopotentials for Plane-wave
Calculations, Phy.Rev., B43, 3 (1991) 1993.

9) M.C.Pyne, M.P.Teter, D.C.Allan, T.A.Arias and J.D. Joannopoulos: Iterative
Minimization Techniques for ab initio Total-energy Calculations: Molecular
Dynamics and Conjugate Gradients, Rev. Mod. Phys., 64,4 (1992) 1045.
10) Tomoya Ono and Kikuji Hirose: Timesaving Double-Grid Method forReal-Space Electronic-Structure Calculations, Phys. Rev. Lett. Vol.82,
No.25 ,p.5016 (1999).
11) Francesco Mauri, Giulia Galli, and Roberto Car :Orbital formulation for electronic-structure calculations with linear system-size scaling, Phys. Rev. B47,(1993)9973.
12) Francesco Mauri and Giulia Galli: Electronic-structure calculations and molecular-dynamics simulations with linear system-size scaling, Phys. Rev. B50,(1994)4316.
13) K.Hirose and T.Ono: Direct minimization to generate electronic states withproper occupation numbers, Phys. Rev. B, Vol. 64, No.8, 085105.