プログラムの高速化・並列化サービスの事例:チューニングレポート<要約>：線形スケール密度汎関数法Conquestの金属計算に対応した最適化

＊ここに掲載するのは、University College London $UCL$ のLianheng Tong博士によるHECToRレポート「Metal Conquest, Lianheng Tong, 2nd March 2011」を要約したものです。

[2017年2月掲載]

概要

今日興味を持たれている金属系を含む物性科学やナノテクノロジーにおける多くの問題の理解には、大規模システム上での計算が要求されます。絶縁体や半導体に対する高スケール性を持つ効率的な解析コードは大きく進歩してきましたが、金属系に対してはその努力は比較的に少ないものでした。絶縁体や半導体用の計算コードをスケーラブルにするアルゴリズムのほとんどは、効率的なハミルトニアン計算やコンパクトで柔軟な基底関数など、金属系にも適用可能ですが、金属系に対して同等のスケール性を達成するのは困難です。何故ならこれらの系の密度行列は一般に長距離相関を持つためで、短距離相関の絶縁体や半導体に対する通常の線形スケール手法はここでは適用できません。このプロジェクトの目的は、絶縁体や半導体向けに設計された線形スケール密度汎関数法コードConquest[4]を、効率的に金属を計算できるように修正することです。

Conquestがこの修正に適しているのは、SCaLAPACK v1.7を用いたハミルトニアンの直接対角化を行う金属向け密度行列ソルバが既に備わっているためです。さらにこれは、プロセッサ数に関して驚異的なスケール性を持つ成熟したオープンソースの線形スケールコードです。Conquestは、Bスプラインまたは疑原子軌道のいずれかを使用して、波動関数を拡張するための柔軟な関数系（サポート関数）を生成します。Bスプラインを用いた場合は、少数の関数で平面波の精度を再現できます。これは金属計算に重要です。何故なら対角化のコストは関数の数に依存するからです。さらにConquestは既に、HECToR上の4096プロセッサを用いて絶縁体のシミュレーションに成功しています（UKCPグラントでのプロジェクトe89で行われた）。

金属向けに効率的な計算を可能にするためにコードに適用すべき変更は以下の通りです：

既存の線形スケール手法ではなく、密度行列対角化へ並列線形代数ライブラリを適用します。これは既にConquestに備わっています。
効率的なブリルアンゾーン積分法 $Methfessel-Paxton method[14]$ を適用して、セルフコンシステント計算に必要な対角化の数を削減します。これが金属系で重要な理由は、部分的に占有されたバンドがフェルミエネルギーにおいて占有関数の非連続性を生じ、ブリルアンゾーンの数値積分精度、延いてはセルフコンシステンシーに悪影響を与えるためです。
大規模系におけるセルフコンシステント計算により効率的な手法を導入します。シミュレーション対象の系が大きくなると、charge-sloshingと呼ばれる電子密度の振動が生じて収束が困難になります。これは金属系ではよくあることです。この問題は、非経験的コードVASP[12]においてKerkerプレコンディショナとWave-dependent metricを用いて対処が成功しています。Conquestでも同じ手法を適用します。
金属系は一般に鋭く急激に変化するバンド構造を持つため、正確なブリルアンゾーン積分には多くのk点が必要です。また各k点では個別に対角化が必要です。さらに複数のk点での対角化を並列処理すれば効率的な計算が可能です。
金属の計算での性能プロファイルを採取し、最適な入力パラメータを決定します。

テストではバルクのアルミニウムを全体を通して使用します。アルミニウムは周期表に最初に出現する金属のため、相対論的あるいは擬ポテンシャルに対する補正項の考慮が不要です。同時にコア電子からの十分なスクリーニングが働くため擬ポテンシャルもそれほど複雑ではありません。コードの効率の実証用のテストのため、余計な複雑性によるオーバーヘッドは不要です。交換項の相関関数にはLDAを用います。これが適切な理由としては第1に、Gaudoin et al.[6, 7]により報告された通り、LDA近似はバルクのアルミニウムに適しています。第2には目的はコード性能の測定であい、バンド構造の精密性ではありません。アルミニウムの擬ポテンシャルは、Rappe et. al.[17, 18]の方法によるOPIUM[1]を用いて計算しました。この方法は一般に、標準的なTroullier-Martins method[19]よりもソフトな擬ポテンシャルを与えます。
アルミニウムでは2つの3sおよび一つの3p電子を価電子として、部分内殻補正（non-linear coreとも呼びます）や相対論補正は含まれません。基底関数はDZPの数値的なPAO基底を用います。この基底関数はSiesta[2]を用いて100meVのエネルギーシフトで生成しました。現状では基底関数当たり1サポート関数を用います。格子定数は、格子定数を3.50から4.40 Áまで変化させたAl4個のユニットセルの基底状態計算結果を内挿して、a₀=4.01155を得ました。また、実空間の各方向に32グリッド点、逆格子空間では25k点とれば（Monkhorst-Park[15]法を用いた）、smearing温度0.001Haで10^-6HAまでの精度に十分であることが判りました。また、結果をサイズでスケールして、32原子であれば各方向に実空間で64点、逆格子空間で13k点、108原子の場合は各方向に実空間で128点、逆格子空間で7kとしました。

KerkerプリコンディショナーとWave-dependent Metric

Pulay mixing[16]は、非経験的DFTコードに実装される標準的なコードで、基底状態のセルフコンシステント計算で用いられます。
charge sloshing[9]と呼ばれる現象は、セルフコンシステント・ループ内の入力密度の小さな振動が出力密度に大きな影響を与える場合を言い、金属系の計算でしばしば生じます。これは、セルフコンシステント・ループの非常に長い繰り返しをもたらし、金属系の高速計算に対する主な障害の1つです。charge sloshingの原因が、ハートレーポテンシャルの振動を支配する密度内の実空間の長距離性の変動（逆格子空間の短距離変化）に由来することに注意すれば、問題の解決は可能です。Mixing手法を少し修正して、密度内の長距離変動をフィルタリングすることで、charge sloshingの影響を抑えてより速い収束へ導くことが可能です。この前処理法の名称はその開発者Kerker[9]に由来し、既に非経験的電子構造コードで成功しています[11,12]。この他の方法として、mixing時の残差項に重みを追加することで、密度の短距離変動の重要度を増します。この方法をWave-dependent Metric法と呼ぶことにします。

·テスト結果

アルミニウム・バルクの対称性のため、charge sloshingは生じません。更新された密度は全て同じ対称性を持ち、格子方向は区別がなく、出力密度も同じで、ハートレーポテンシャルの対称性を不変に保ちます。32原子および108原子のユニットセルを用いた計算は、Kerker前処理やwave-dependent metricを使わなくとも12回のPulay mixingステップで収束します。

ここで対称性を破るために、バルクアルミニウムに原点座標の１原子を取り去った欠損を導入します。標準的な線形混合法では混合パラメータ0.1から0.8を用いても収束しませんが、Kerker前処理を追加すると収束性が大きく向上します。原点に欠損を持つ32原子から成るバルクアルミニウムを、単純な線形混合法を用いるとcharge sloshingが生じます。これは基底状態エネルギー計算における振動として観測されます。Kerker前処理の長距離相関パラメータq0を増やすと、charge sloshingの効果は減衰します。

smearing温度300K（ほぼ0.001Ha）と混合パラメータ0.5の単純線形混合を用いて計算すると、q0<0.4では収束せず、q0=1.0*2π/Bohr=1.188/Åが最適であることが今回のケースにおいては見出されました。これは文献[11]の結果と同等です。

同じケースでwavedependent metric methodをテストしましたが、Pulay mixing単独で（5ステップで）十分な収束に達しました。Kerker前処理やwavedependent metric methodの使用は、共に収束に影響を与えませんでした。

ステップ関数のMethfessel-Paxton近似

金属系の計算に関わるもう一つの大きな問題に、バンドが部分的にしか専有されていないことから生じるブリルアンゾーン上の不連続関数の積分があります。数値積分では、ブリルアンゾーン内に非常に細かい空間メッシュ $kメッシュ$ が必要になります。この問題の改善方法の一つが、占有関数のステップ関数を有限温度Tのフェルミディラック関数で置き換える方法です。こうすると被積分関数がブリルアンゾーン全体で微分可能となり、kメッシュの収束性を改善します。ただし得られた結果は、実際には有限温度下での系のわずかに異なる問題であることに注意すべきです。

MethfesselとPaxtonの方法[14]は、より洗練された占有関数の近似方法です。この方法は、デルタ関数を直交エルミート多項式で近似することから始めてステップ関数を近似する方法です。その後にデルタ関数を積分することでN次の近似S_N $x$ が求まります。０次の場合には、この方法はフェルミディラック様のsmearingに相当します。

Methfessel-Paxton smearingを用いて計算した基底状態エネルギーEは、もはや変分的ではありません。その代わりに基底状態の密度は自由エネルギーFを最小化したものです[11]。FとEの差は、大きなk_BTにおいても僅かです。これがMethfessel-Paxton法が金属系の計算に望ましい理由です。こうして、それほど大きなFとEの差を生じさせずに、比較的高いsmearingファクターk_BTを用いてより粗い逆格子グリッドを用いることが可能になります。もしF-Eが大きければ、計算された基底状態エネルギーが（ここで計算しようとしている）真のゼロ温度エネルギーから逸脱していることを意味しており、結果は信頼できません。

しかしながら、Methfessel-Paxton法もまた数値的なアーチファクトを生じます。これは、フェルミエネルギーE_Fが複数存在してしまうことです。

それでも、ここではフェルミエネルギーを正しい電子数を与える最小エネルギーであると定義し、最小エネルギーであることを保証するために、エネルギーの下限と上限を注意深く決めつつ最低のフェルミエネルギーを決定します。

フェルミディラックおよびMethfessel-Paxton smearing法を用いて、ユニットセルに4原子のバルクアルミニウム計算を比較しました。フェルミsmearingは、smearingファクターk_BTを増やすと異なる結果を生じます。この時、自由エネルギーと基底状態エネルギー（ここではHarrisエネルギーと呼びます）の差も大きく増大します。これはそのsmearingファクターが物理的な温度に相当し、k_BTが増加すると温度ゼロの領域から離れていくためです。一方、Methfessel-Paxton法ではこれは近似パラメータでしかありません。ここではk_BTの増加はエネルギーに小さな影響しか与えず、結果はエルミート多項式の次数を上げると改善します。つまり、Methfessel-Paxton近似では比較的大きなsmearingファクターを用いて、より少数のｋ点を用いる事が可能であることを意味します。
具体的には、k_BT=0.1Ha=2.72eV（smearing 温度31565.51Kに相当）では各逆格子方向に10点を用いて収束しました。全エネルギーは、小さなsmearingと大きなk点を用いた基底状態エネルギーの10^-5以内です。

ScaLAPACK性能プロファイリング

プロファイリングは、CrayPATおよびCray LibSci10.5.0を用いてHECToR XT5h上で行いました。ユニットセルに32原子を配置したバルクアルミニウムのテストを、プロセッサグリッド1×4のBLASを用いて、異なるScaLAPACKブロック次元数に振って計算しました。

最適なブロック次元は26×26で、比較対象のConquestのデフォルト104×104に比べて約1割改善していました。ボトルネックとして、対角化で利用されるScaLAPACKルーチンpzhegvx内のMPI_recv呼出しに、最大の負荷分散がある事が判りました。しかしながら大きなブロックサイズの場合は、MPI_recv呼び出しの負荷分散は、呼出しが減って部分的にオフになり、MPI_bcastが主要なボトルネックとなります。
つまり、ScaLAPACKブロックサイズの選択がその効率を決めています。32原子の場合はその最適値が26×26です。

K点並列

Conquestの実装は、系が周期境界条件を持つとは見做さず、系は実空間内の任意の場所に存在するものとして扱われます[5]。ユーザ入力セルはFundamental Simulation Cell $FSC$ として定義され、全ての格子方向に繰り返されて、FSC内の原子と相互作用する全ての原子が計算対象となります。全ての量子力学演算子は、サポート関数（これは行列表現を持つ基底関数です）を用いて行列として表現されます。サポート関数については文献[4,5]を参照してください。

Conquestの行列保存形では、FSCの原子はプロセッサで分割されて、各プロセッサが担当するFSC原子はプライマリーセットと呼ばれます。また柔軟なデータ転送のために、空間はパーティションに分割されて、各パーティション内の原子は特定のパーティションに属するものとして扱われます。異なる行列（演算子）が与えられた場合、原子間の相互作用範囲が異なり、与えられた行列（範囲）型のハロ原子は、プライマリーセットの任意の原子の相互作用範囲内にある任意の原子であると定義されます。

ハロ原子のセットは、プライマリーセット原子の近傍原子セットの集合（近傍リスト）に相当します。与えられた相互作用範囲のハロパーティションとは、少なくとも一つのハロ原子を含むパーティションです。次に、ハロは与えられた範囲に相当する全ハロパーティションの集合と定義されます。最後に、プロセッサのカバリングセットが、最も大きい（最も長い相互作用範囲を有する）ハロを含む最小斜方晶セルとして定義されます。与えられた原子を有する全ての項が保存されるわけではなく、実際に各プライマリーセット原子の近傍リストにあるものだけが行ごとにフォーマットされて保存されます。

Conquestは電子密度を解く際に、一般化固有値問題を周期境界条件で解きます。また、固有値問題は一度に一つのk点での方程式を解きます。更新されたハミルトニアン行列と重なり積分行列は、ScaLAPACKブロックフォーマットにでBLASプロセッサグリッドによる全てのプロセッサへ分配されます。次にScaLAPACKフォーマットからConquestのフォーマット[3]へ変換され、セルフコンシステント計算へ渡されます。

この異なるk点の固有値ベクトル計算は互いに独立です。異なるｋに対してプロセッサのサブグループを割当てる際には、ScaLAPACK計算に最適化された最適なパラメータを選択することが可能です。与えられたサイズの行列に対して、最適なプロセッサ数が存在しますが、多すぎると効率の悪い通信を引き継ぐことを意味します。したがって、k点並列計算は、理論的には、各ScaLAPACKサブルーチンコールに対して最適なプロセッサ数のグループを選ぶことで、より多くのプロセッサをより効果的に利用することを可能にします。金属計算に必要とされるk点は数千オーダーであり、これは特にHECToRなどのHPCシステム上で利用できる実際の自由度数です。

テストは、HECToR XT5h $およびCray LibSci 10.5.0$ と、ローカルのAMD Opteron 2214 $2 Core 2.2GHz$ のSunワークステーション $およびACML/LAPACK、ScaLAPACKはローカルでコンパイル$ で行いました。ユニットセルに32原子のバルクアルミニウムを用い、13×13×13k点メッシュ、0.001Haのフェルミ-ディラックsmearingで計算しました。全て4ノードを用いて非セルフ紺紙ステント計算です。
既存の実装に比べて、同等のリソースを用いてk点並列により大きな改善が見られました（ScaLAPCKブロック26×26を用いた場合に、HECToRでは4並列で1.27倍、ワークステーションでは4並列で2.64倍に高速化しました）。

謝辞

このプロジェクトは、nAG Ltd.が運営するHECToRの分散計算科学および工学 $CSE$ サービスの基に実行されました。英国の国立スーパーコンピューティング・サービスである、HECToR：英国リサーチ・カウンシル・ハイエンド計算サービスは、リサーチ・カウンシルを代行するEPSRCが管理しています。そのミッションは英国学術界の科学および工学の研究支援です。HECToRスーパーコンピューターは、UoE HPCx Ltd.およびnAG Ltd.のCSEサポートサービスにより管理運営されています。

文献

[1]	http://opium.sourceforge.net/.
[2]	http://www.icmab.es/siesta/
[3]	D. R. Bowler. Implementation of diagonalisation within conquest. Technical report, University College London, Oct. 2008.
[4]	D. R. Bowler, I. J. Bush, and M. J. Gillan. Practical methods for ab initio calculations on thousands of atoms. Int. J. Quantum Chem., 77 $5$ :831-842, 2000.
[5]	D. R. Bowler, T. Miyazaki, and M. J. Gillan. Parallel sparse matrix multiplication for linear scaling electronic structure calculations. Comput. Phys. Commun., 137 $2$ :255-273, June 2001.
[6]	R. Gaudoin and W. M. C. Foulkes. Ab initio calculations of bulk moduli and comparison with experiment. Phys. Rev. B, 66 $5$ :052104, Aug 2002.
[7]	R. Gaudoin, W. M. C. Foulkes, and G. Rajagopal. Ab initio calculations of the cohesive energy and the bulk modulus of aluminium. J. Phys.: Condens. Matter, 14 $38$ :8787-8793, 2002.
[8]	P. Hohenberg and W. Kohn. Inhomogeneous electron gas. Phys. Rev., 136 $3B$ :B864, 1964.
[9]	G. P. Kerker. Efficient iteration scheme for self-consistent pseudopotential calculations. Phys. Rev. B, 23 $6$ :3082-3084, Mar 1981.
[10]	W. Kohn and L. J. Sham. Self-consistent equations including exchange and correlation effects. Phys. Rev., 140 $4A$ :A1133, 1965.
[11]	G. Kresse and J. Furthmüller. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. Computational Materials Science, 6 $1$ :15-50, July 1996.
[12]	G. Kresse and J. Furthmüller. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B, 54 $16$ :11169-11186, Oct 1996.
[13]	R. M. Martin. Electronic Structure: Basic Theory and Practical Methods. Cambridge University Press, 2004.
[14]	M. Methfessel and A. T. Paxton. High-precision sampling for brillouin-zone integration in metals. Phys. Rev. B, 40 $6$ :3616-3621, Aug 1989.
[15]	H. J. Monkhorst and J. D. Pack. Special points for brillouin-zone integrations. Phys. Rev. B, 13 $12$ :5188-5192, Jun 1976.
[16]	P. Pulay. Convergence acceleration of iterative sequences. the case of scf iteration. Chem. Phys. Lett., 73 $2$ :393 - 398, 1980.
[17]	A. M. Rappe, K. M. Rabe, E. Kaxiras, and J. D. Joannopoulos. Optimized pseudopotentials. Phys. Rev. B, 41 $2$ :1227-1230, Jan 1990.
[18]	A. M. Rappe, K. M. Rabe, E. Kaxiras, and J. D. Joannopoulos. Erratum: Optimized pseudopotentials [phys. rev. b 41, 1227 $1990$ ]. Phys. Rev. B, 44 $23$ :13175-13176, Dec 1991.
[19]	N. Troullier and J. L. Martins. Efficient pseudopotentials for plane-wave calculations. Phys. Rev. B, 43 $3$ :1993-2006, Jan 1991.