＊ここに掲載するのは、クイーンズ大学ベルファストのJonathan S. Parker博士によるHECToRレポート「Optimization of the MPI parallel RMT code for HECToR and likely successors, Jonathan S.Parker, October 3, 2011」を要約したものです。
概要
RMT法（時間依存R行列）は、高強度短レーザーパルス内の多電子原子分子系の時間依存シュレディンガー方程式（TDSE)に対する新しいab initio法です。近年、この他にも様々な時間依存R行列法が提案されて来ましたが[1,2,3,4,5,6]、RMT法は、原子のコアよりも遠方の数電子の波動関数に対して有限差分（FD)法を用いることで桁違いの効率を示します。RMT法は外殻領域のFDモデルと多電子の内殻領域を融合します[7]。プロパゲータのユニタリティを保ちつつ基底関数モデルと隣接するFDモデルを融合する問題は、この分野において以前からの課題になっていました。RMTは2008年のNikolopoulos, Parker and Taylor [8]による報告をベースにしています。このレポートではこの方法が計算的に安定で高い効率を持つことを示します。
RMTの実装は現在、ワークステーション、中規模並列マシンおよびHECToR上で並列化され計算されています。RMTは内殻に対して並列R行列計算コード、外殻は有限差分を用いるHELIUMコードを使用します[9,10]。HELIUM[11,12]は、15年以上前のCray T3D設置以来大規模並列マシン上で煩雑に利用されてきました。HELIUMは強い場の内の2電子原子あるいはイオンの完全次元のシュレディンガー方程式を正確に解きます。これまで、HECToRの16,000コアやその姉妹機であるオークリッジのJaguarの60,000コア以上における実行など、大規模並列マシンにおいて高い効率を実証されています。IBM Blue Geneにおいても同様の精度で良好な効率を示しました。HELIUMはシングルプロセッサマシンや、8コアワークステーション、ワークステーションクラスタおよび大規模Beowulf型クラスタで定期的に利用されています。
HELIUMはレーザー励起された2電子原子の完全TDSEの解を生成するように設計されました。HELIUMの実験で判明したことですが、高い精度を得るには、大きな積分体積、言い換えれば1000ボーア半径の外殻領域が必要です。同じく重要なのが、放出された電子のエネルギースペクトル計算は、それらがコアから遠方に運ばれる際に電子波動関数を高精度に解かない限り困難であるという事実です。このため、RMT法がHELIUMと同等の精度が要求される場合には、有限差分法とプロパゲータアルゴリズムを検討すべきです。
最近の5年間に渡り、レーザー技術の革新により前例のない時間解像度を持つ物質-レーザー相互作用の実験研究が可能になりました。現在、800nm波長およびVacuum-Ultraviolet (VUV)波長の250アト秒パルスによる、高強度Ti:Sapphireレーザー光の数サイクルパルスの実験が可能です。特にアト秒パルス生成は、超高速電子プロセスの研究に新たなフロンティアへの道を開きました。最近の特筆すべき固体内アト秒スペクトル測定[13]として、エレクトロニクス分野の理論的速度限界で生じる固体状態プロセスが測定されました。アト秒パルスは最近になって、原子内電子トンネル現象[14]や内殻電子のオージェ崩壊[15]のリアルタイム測定を可能にしました。アト秒パルスはアルゴンの1電子イオン化のストロボスコピック測定研究[16]も可能にしました。
数サイクルのレーザーパルスの間に生じる相互作用は、極度に大きなパルス強度や短いパルス長を要因としますが、これまでの時間依存手法では信頼に足るモデリングが困難です。このような状況下で理論が意味のある役割を担い、特に予測を行うためには、多電子原子の構造とレーザーによる時間依存の多電子応答を正確に記述することが可能な、より洗練された計算手法が必要です。
R行列法は多電子原子や多電子分子とレーザー間の相互作用のモデル化に成功してきましたが、時間依存性を記述しません。HELIUMは時間依存の原子-レーザー相互作用のモデル化に成功していますが、2電子原子に限られています。RMT法はこの両方の制約から自由になります。R行列法を利用すれば、全電子が存在する核近傍の内殻領域の動的多電子記述が可能になります。これが重要なのは、高周波数極限（VUVやXUV）において内殻励起が相互作用の多くを占めると予想できるからです。HELIUMの有限差分法を利用すれば、二価イオン化において最大2つの電子が存在し、内殻領域よりも大きな空間を占める外殻領域の高精度記述が可能になります。RMTの成功にとって等しく重要なのが、HELIUMが持つ並列コンピュータ上での高い効率とスケーラビリティです。
まとめると、RMT開発の目標は、近年の実験の進歩に対して、既存の方法では不可能だった高い信頼性を持つ理論解析を可能にすることです。その対象には、アト秒時間スケールのイオン化現象の時間解像研究や、二価イオン化での電子放出時間差の研究、複雑な原子における内殻励起と崩壊、自由電子X線レーザーによるXUV極限での高強度場の原子-レーザー相互作用、および原子分子内の高調波生成などが含まれます。
目的と成果
作業項目1では、RMTの内殻領域（I.R.)と外殻領域（O.R.)の負荷バランス最適化のための開発とテストを行います。100から10,000コアを用いた場合の負荷バランスのための最適なパラメータ選択を議論します。テストはI.R.に100から1,000コア、O.R.に100から10,000コアを用いました。
I.R.とO.R.間の通信順序の最適化（これを以降はRed-Black最適化と呼びます）により、500I.R.コア以上では積分スピードが1.7倍に高速化しました。
成功した最適化は2つです：1）O.R.からI.Rへ計算を移動（こうして各時間ステップでの領域間通信が削減されます）、2）領域間通信に専用コアを一つ割り当てる。この最適化により大きく高速化しました。既存のMPI設計をコレクティブ通信で置き換える検討を行いましたが、Red-Blackスキームには適しませんでした。
作業項目2の当初の目標は、プロパゲータのアルゴリズムを改善し、空間グリッド間隔の最小極限において効率を改善することでした。当初はハミルトニアンの近似を累乗することが含まれていましたが、そのオーバーヘッドが効率の改善を上回ってしまいました。代わって開発したのは、最小二乗有限差分演算子をベースにしたもので、結果として速度が改善されました。この最小二乗有限差分法は既存の5点有限差分の精度を維持しつつ、有限差分グリッド上の高周波モードを抑制して安定性を改善することで、より大きな積分ステップと速度向上を可能にしました。最小二乗法を用いることで、5点有限差分に比べて1.8倍高速化しました。
作業項目1 MPI実装の最適化
こうした大きな領域を積分する時間依存R行列コードは、その高次元偏微分方程式を解くのに計算時間が多くを占めるため1万コア以上の並列処理が必要となりますが、純粋なR行列形式では直截的には困難です。
RMTコードをHECToR上で効率的に実行するために、この作業では、内殻領域(I.R.)と外殻領域(O.R.)間の通信を効率化して、複数コアに対する負荷バランスの改善に焦点を当てます。
負荷バランス化とは、実行中にHECToRのコアのアイドル時間を最小化することです。内殻と外殻領域に対するコア上の計算アルゴリズムは大きく異なるため、1ステップ先へ波動関数を更新する計算負荷も異なります。よってRMTコードを効率化するには、異なる条件下でのコード実行試験を通じて、両方の領域間の計算負荷の違いを理解する必要があります。内殻領域では時間の前進には行列-ベクトル積が含まれます。
この行列はブロック三重対角で、1ブロックは全角運動量の2つの値から成る行列要素からなります。I.R.はこれらブロックで並列化され、全角運動量はゼロからある最大値までの範囲で変化し、最小コア数はこの値+1を用います。この時各コアは行列-ベクトル積の一部分を担当します。
外殻領域では計算はグリッド上での領域分割を用います。各コアは波動関数の一部を扱い、それをNgridXNchannels次元の配列に格納します。ここでNgridはグリッド数、Nchannelsはチャネル数です。チャネル数は全角運動量と計算において残るターゲット状態数に依存します。
任意の全角運動量と残ターゲット状態数に対して、コア当たりの最適なグリッド数があれば、2領域のコアでの最良の計算バランスが可能です。全角運動量と残ターゲット状態数の値は対象とする問題に依存します。問題のタイプの指標を得るために、まずネオン原子を考えます。ネオンはエネルギーが似通った多くの残イオンターゲット状態が存在します。これらは8(21)つあり、角運動量は24(30)までで500(1500)チャネルを有し、I.R.ハミルトニアンは25,000(120,000)次元です。波動関数を高強度赤外レーザーの10サイクルに対して正確に更新するには、20,000-35,000ボーア半径程度に動径グリッドを伸ばす必要があります。場合によっては共鳴構造を解く必要もあります：こうした極端なケースでは、100以上の周期が必要となり、その結果動径グリッドは1から2オーダーさらに延長する必要があります。こうした場合は、I.R.では100-1000コア、O.R.では200-10,000コアが必要になると算段しました。
·200-10,000コア上の負荷バランス性
問題は、RMTの実際的な限界において計算のバランス化に要する最適なパラメータを選択することです。特に興味があるのは、10,000コア以上を用いたO.R.のケースです。これは高強度場による高調波生成研究での典型的なコア数です。
ここでは、内殻領域は100コアに固定して、O.R.コアを100-10,000コアで変化させてテストします。Red-Black最適化では、I.R.のコアは100コアまで変化させました。
O.R.コアが計算を素早く完了した場合にはアイドルし、I.R.コアの計算を待ちます。もしO.R.コアが時間当たりの計算が多すぎた場合は、I.R.コアO.R.コアの計算が終わるまでアイドルします。
I.R.を決めるパラメータは問題の物理に依存します。簡単に調節できるパラメータは、外殻領域のコア当たりのグリッド数です。これはI.R.とO.R.の実行時のスケール性から直接計算できます。
しかしながら実際には、10,000コアに近づくにつれてMPI通信オーバーヘッドが性能を劣化させます。
·I.R.とO.R.計算順序の最適化
RMTの波動関数更新では、時間ステップの最初に両方の領域において波動関数が確定していなくてはなりません。言い換えれば、各時間ステップの最初と最後で両方の領域は同期します。計算内容は、ハミルトニアンと波動関数の掛け算です。I.R.はO.R.と独立に計算できますが、O.R.は計算開始時にI.R.の波動関数を必要とします。そこでO.R.はこのデータを受け取るまでアイドルします。受信が完了したらO.R.は計算を開始し、I.R.はアイドルします。O.R.は一次の更新を行い、I.R.へ情報を送ります。ここでI.R.は一次の更新を開始し、O.R.は2次の更新をI.R.と並列で行います。高次ではより大きな割合が並列処理されます（即ち、より多くのI.R.が独立に動作します）。この処理はArnoldiプロパゲータの最高次数まで繰り返されます。最後にI.R.は最高次数まで完了し、2つの領域は同期します。
内殻と外殻領域が独立に動く際に、内殻コアを2つに分割して並列自由度を増やすことができます。ここで、それぞれをRed内殻コア、Black内殻コアと呼びますが、これらは独立にO.R.から情報を受け取ることが可能です。ある極限においてこの独立性は、片方のコアセットがアイドルする代わりに計算を開始しもう一方のコアセットがO.R.と同期することで、I.R.並列効率を改善します。これはより多くのコアを利用する能力を高め、積分全体の効率も改善します。
外殻領域のオーバーヘッドはコア当たりちょうど32グリッドを用いた場合に最小になります。
また、最大角運動量=23の場合、この最適化の結果、240内殻コアまでほぼリニアにスケールします。240コア以上では、Red-Blackを適用しないと性能向上は認められませんが、Red-Blackを用いれば（傾きは下降しますが）線形性は継続されます。580コアの場合、Red-Blockはそれを用いない場合に対し1.45倍高速化します。コア数が2倍以上ですが、45％の高速化にとどまります。
実際には外殻コアは1,000-10,000コアを用いるため、45％の性能向上においては内殻コアの増加分コストは安いものです。600I.R.コアに増やしたところこのケースでは改善は示しませんでした。
また、8次のプロパゲータのケースでは1.7倍になりました。通常はArnoldiプロパゲータの次数を倍にすると実行時間も倍になりますが、積分スピードの効率は倍以上になります。実際には、メモリー許容範囲で用いられる最高次数は通常12-16が標準です。
·I.R.とO.R.領域間通信の最適化
時間ステップ毎にハミルトニアン・波動関数乗算が数多く計算され、計算の度にI.R.コアとO.R.コア間で情報が好感されます。この処理の効率改善の方法を検討しました。
方程式の右辺はO.R.で計算され、I.R.へ送信されます。上述のように書くI.R.コアは行列・ベクトル積を行うため、O.R.データは全てのI.R.コアで必要になります。ここで考えられる最適化は、1対多のMPI集団通信を用いるようにすることです。しかしながらこの最適化は上述のRed-Black最適化と適合させることが出来ません。
領域間通信でもう一つ障害となるのは、I.R.で用いるsurface amplitudesのO.R.における計算です。各時間ステップでO.R.において、チャネル波動関数の動径微分の時間に依存しないsurface amplitudes への射影によるベクトルが繰り返し計算されて、I.R.へ送信されます。この計算はO.R.からI.R.へ移動させられました。これは最適負荷バランス状態での観測から、I.R.コアが計算に用いることのできるアイドルするCPUサイクルが余分にあったためです。
当初は領域間の通信サイズの削減を目論見ましたが、作業の結果、データパケットのサイズが大きくなってしまいました。しかしながら以下の最適化によりこの性能劣化が除去されました：一つのコアにデータ転送を任せ、他の時間が余っているコアにO.R.ハミルトニアンの更新を任せました。この結果O.R.コアは完全に同期しました。この最適化の結果、コードは高速化しました。
作業項目2 プロパゲータのアルゴリズムの改良
例えば核電荷Zの極限などの場合には、有限差分の外殻領域の積分は内殻領域の時間積分に比べて効率が極度に落ちます。直截的な改良も可能ですが、大きくイオン化されたモデルでは特に問題で、グリッド間隔を小さくする必要が出てきます。HELIUMを用いた典型的な問題では、Z=2の2電子原子ではδr=0.25auで十分です。一方ネオン、Z-10で電荷が8にイオン化された場合、δr=0.025auが必要です。
残念ながら、O.R.は1/δr2に比例する積分時間ステップδtが必要です。言い換えれば、δrを、δr=0.25auからδr=0.025auへ1/10にする場合には、δtを10,000倍に小さくしなくてはならず、プログラムは10,000倍遅くなります。この問題は、有限差分ハミルトニアンの最高エネルギー固有値が運動量の2乗に比例し、その最大値が10/δr2であることに起因します。この高いエネルギー固有値はグリッドに起因する最高のエネルギーを持つ平面波の存在によるものです。これは一般には物理的にあり得ない励起を意味します。これらはレーザーパルスで励起した物理的な平面波のエネルギーよりも何桁も大きな値です。
例えば、δr=0.25auでは、物理的な励起が10auのエネルギー以上になるのは稀です。しかしながらこの時、グリッド上の最高エネルギーは100auのオーダーです。これはδr=0.025auでは10,000auにもなります。ですが、プロパゲータはこうした非物理的なモードにも励起が存在するとして運動方程式を積分しなくてはなりません。一方、こうした偽励起は積分誤差を生じ、不安定になり計算はうまくいきません。HELIUMやRMTで用いられる陽的なプロパゲータでは、安定性が保てるのは、δtがグリッドに起因する平面波の最高エネルギー周期のオーダー、つまり1/δr2にスケールする場合です。
RMTコードの有限差分はO.R.部分です。O.R.は通常は原子のコアから20au以上の距離にあります。この距離においては、イオン化した電子と原子核間引力や、残りの電子との反発力は比較的に小さくなっています。ここでの有限差分ハミルトニアンの固有値は運動エネルギーKが支配します。つまりKの固有値スペクトルは二次微分演算子の固有値で支配されます。（求心力はrの逆2乗で小さくなり、スペクトルにはほぼ無視できます）。電場もスペクトルに影響を与えます。典型的な場の強さでは、自由場でのピーク固有値を、稀ですが5-20%程度シフトする可能性があります。極端な場合を除き、外殻領域ハミルトニアンの固有値スペクトルは、Kに現れる2次微分演算子の良い近似です。
δrの小さな極限では、積分ステップδtは、運動エネルギー演算子Kの2次微分有限差分演算子によるδ1/r2の依存性を持つことが期待されます。これは実際に、RMTコードの実行においても観察されます。内殻、外殻領域共に同じ時間ステップサイズδtを用いることから、Kの固有値スペクトル計算はRMTコードの実行効率に影響します。 そこで、Kの高い固有値の影響の緩和のための修正が作業項目2の目的となりました。Kのみに対して近似を適用する修正は、ステップサイズを大きくした場合の利得よりもオーバーヘッドが大きいことから採用しませんでした。その代わりに、最小2乗演算子を用いてKのピーク固有値を削減する手法を開発しました。この新しい手法により、僅かな追加コストでKのピーク固有値は4倍まで削減されました。実際にこの追加コストは、RMTの実行時間に全く影響を与えませんでした。
新しい手法は、δtを既存コードに比べて1.8倍大きくすることが可能でした。その積分は既存コードと同等の精度と安定性を持ちます。この手法は理想に近い解決策です。
2次微分計算で用いる多項式次数は、例えば9点ルールの場合、既存の方法は8次です。言い換えれば、9点各々を通るように8次多項式を用います。9点ルールでは、中心点でこの多項式の正確な2次微分が得られます。8次多項式は、任意に選ばれた9点を通ることのできる最小次数多項式です。最小二乗ルールは8以下の次数の多項式を用います。一般にこの多項式は9点を通りません。このため、exp(ikr)の線形結合として記述される関数の数値的な積分で生じる数値的な雑音に鈍感です。最小二乗プロセスは、低周波成分を大きく変えることなく、固有値スペクトルの高周波成分を劇的に除去します。
ピーク固有値は物理的な励起ではありません。レーザー励起された電子のエネルギーよりも10-100倍大きな値が示されます。これは有限差分の数値的アーチファクトです。これらを抑えることで、安定性と精度を確保しつつステップサイズを大きくとることが可能になります。
結果として、9点最小二乗演算子が、5点の非最小二乗ルールで可能なステップサイズよりも大きなステップサイズが可能でであることが示されました。
この手法の利点の一つは9点ルールの実行速度です。HECToR上でのテストでは、9点ルールは5点ルールと同じくらいの速度です。これは、有限差分計算のオーバーヘッドの多くが巨大な波動関数配列へのメモリーアクセスであることを意味しています。このメモリーフェッチのオーバーヘッドは9点と5点ルールで同じです。また、2次微分演算子は非物理的な固有値を計算しますが、その他の演算に比べれば小さな計算コストしか持ちません。このため、9点ルールのRMTは5点ルールと同程度の速度を持ち、δtが1.8倍に大きくなると、積分全体の速度も1.8倍になります。
謝辞
このプロジェクトは、NAG Ltd.が運営するHECToRの分散計算科学および工学(CSE)サービスの基に実行されました。英国の国立スーパーコンピューティング・サービスである、HECToR：英国リサーチ・カウンシル・ハイエンド計算サービスは、リサーチ・カウンシルを代行するEPSRCが管理しています。そのミッションは英国学術界の科学および工学の研究支援です。HECToRスーパーコンピューターは、UoE HPCx Ltd.およびNAG Ltd.のCSEサポートサービスにより管理運営されています。
文献
| [1] | van der Hart HW, Lysaght MA and Burke PG, Phys. Rev. A 76 043405 (2007). |
| [2] | van der Hart HW, Lysaght MA and Burke PG, Phys. Rev. A 77 065401 (2008). |
| [3] | Guan X, Noble CJ, Zatsarinny O, Bartschat K and Schneider BI, Phys. Rev. A 78 053402 (2008). |
| [4] | Lysaght MA, Burke PG and van der Hart HW, Phys. Rev. Lett. 101 253001 (2008). |
| [5] | Lysaght MA, van der Hart HW and Burke PG, Phys. Rev. A 79 053411 (2009). |
| [6] | Lysaght MA, Burke PG and van der Hart HW, Phys. Rev. Lett. 102 193001 (2009). |
| [7] | van der Hart HW & Burke PG (the adaptation);Burke VM, Noble CJ, Plummer M and Burke PG (RMATRIXII/RM95): both to be submitted to Comput. Phys. Comm. |
| [8] | Nikolopoulos LAA., Parker JS and Taylor KT, Phys. Rev. A 78 063420 (2008). |
| [9] | Lysaght MA, Moore LR, Nikolopoulos LAA, Parker JS, van der Hart HW and Taylor KT, 'Ab initio methods for few- and many-electron atomic systems in intense short-pulse laser light’ Quantum Dynamic Imaging: Theoretical and Numerical Methods, editors. A.D. Bandrauk and M. Ivanov, to appear (Springer) (2010) |
| [10] | Moore LR, Lysaght MA, Nikolopoulos LAA, Parker JS, van der Hart HW and Taylor KT, submitted to J. Mod. Opt. (2010) |
| [11] | Smyth ES, Parker JS and Taylor KT Comput Phys Comm 144 1 (1998). |
| [12] | Parker JS, Doherty BJS, Taylor KT et al, Phys. Rev. Lett. 96 133001 (2006). |
| [13] | Cavalieri AL et al Nature 449 1029 (2007). |
| [14] | Uiberacker M et al Nature 446 627 (2007). |
| [15] | Drescher M et al Nature 419 803 (2002). |
| [16] | Mauritsson J et al Phys Rev Lett 100 073003 (2008). |