プログラムの高速化・並列化サービスの事例

＊ここに掲載するのは、クイーンズ大学ベルファストのMichael A. Lysaght博士らによるHECToRレポート「Implementation of established algorithms to extend HELIUM$WP5: Extending HELIUM to treat many-electron atoms$, Michael A. Lysaght and Laura R. Moore, Department of Applied Mathematics and Theoretical Physics, Queen’s University Belfast, March 31, 2011」を要約したものです。

概要

物質中の相関する電子運動の特性時間スケールはアト秒であり、フェムト秒レーザーは実時間の核運動観察を可能にし、アト秒光パルスは原子系の電子の相関運動をそのスケールで観測し制御する技術発展を約束するものです。最近の超高速測定は、原子や分子内の電子の相関運動の理解に障害になっているのは現状の理論であることを示しています。相関電子運動はその定義から量子多体プロセスであり、その理論は広く用いられている単一粒子描像を超えて、超高速単一電子プロセスを記述を可能にしなくてはなりません。現在、単一粒子描像を超えた最も洗練された方法がHELIUMコードの基礎になっています。HELIUMは2電子系の相関電子運動の記述に成功してきました。このプロジェクトの目的は、高強度短パルス光に曝された他電子原子内で生じる相関電子運動を正確にくじゅつできるように、HELIUMを拡張することです。このレポートは、HELIUMを新たなFortran90/MPIコード”RMT”（時間依存R行列法）へ拡張した実装の基礎となる数値計算手法とアルゴリズムを記述します。

過去10年間に渡り、短高強度光パルスによる多電子原子分子の詳細な応答を記述するために、時間依存シュレディンガー方程式を直接かつ正確に解くことが極めて重要となって来ました。この計算はアト秒科学[1,2]（1アト秒は10^-18秒）の成功に必須で、その目標は物質中の電子の相関運動をその運動時間において追跡・制御するという革新的なものです。アト秒科学の出現は、レーザー物理と電子衝突物理学に革命をもたらしています。アト秒時間スケール（ATS)の観測は、その進歩を牽引する研究室が存在する英国[3,4]、欧州[5,6,7]、北米[8,9,10,11]、中国[12]および日本[13,14]において現在進められています。アト秒光パルス（ALP）技術[2]が近い将来期待されるのは、バイオエナジェティクス[15]、分子エレクトロニクス[16]あるいは癌治療[17]など多様な分野への適用です。

この10年に多くのアト秒科学の技術的ブレイクスルーが生じましたが、理論の進歩としては、短および長波長の高強度フェムト秒（1フェムト秒は10^-15秒）レーザーパルスに対する原子分子の有効一電子：Single Active Electron（SAE）応答の定量的理解を可能にしたことがその一つです。長波長領域ではTi:sapphireレーザー技術の進歩[19]は高高調波生成（HHG)[20]の様なSAE現象の理解を進め、最近になって僅か100アト秒のみ持続する極紫外（XUV)光パルスを生成するに至りました[21]。

SAEベースの理論がこの10年来研究されてきた中で、ALP技術は現在原子分子内の複数電子の運動を励起してその超高速追跡を可能にしつつあります[5,7]。こうして、高強度短レーザーパルスに曝された多電子原子や分子内の電子の相関運動を扱うことのできる追跡手法の開発が極めて重要になります。現状で原子内電子の相関運動を記述する最も洗練されたコードはHELIUM [24,25]です。HELIUMはアト秒科学における先駆的コードです。これは重要な科学的発見（その後実験で確認された）[25]のみならず、多次元TDSEの正確な解の伝播の研究に用いられて来ました。このコードは特に、高強度レーザー光に曝された1，2電子原子やイオンを記述するように構成されています。これは当初からMPPアーキテクチャ向けに設計されたため、コードの成功には疑う余地がありません[24]（この領域では必須と認識されています）。このコードのカーネルは有限差分法（FD）と組み合わせた高次Arnoldiプロパゲータです。HELIUMはプロセッサー上で適切に分散させることで最近接のみで通信が可能になり、HECToR上の16,000コア以上で高い効率でスケールさせることが可能です。

この効率性こそが、時間依存の原子分子系の実験的な対象に沿った、高強度短光パルスに曝された多電子原子を扱う数値計算手法に用いる理由です。特にHELIUMコードが基礎にするFD手法は、一価あるいは二価イオン化で多電子原子分子から放出された1，2電子を効率的かつ正確に記述するのに非常に適しています。この能力を用いて、HELIUMの設計者はFD手法がこうした状況を正確に記述することを実証してきました。これは、空間分割コンセプト[26]の利用を、原子分子過程で成功しているR行列法を用いて可能であることが示されていました[27]。空間分割コンセプトでは、電子の占有する配位空間は2つの分離された領域に分割され、FD法による積分で多電子原子の超高速イオン化を記述可能にするだけでなく、原子の多電子記述を核近傍の小さな領域のみに制限する有効な方法でもあります。この有限の領域においては、R行列法が多電子波動関数の正確かつ扱いやすい表現を提供することがよく知られています[27]。実際に、HELIUMコードが用いるArnoldiプロパゲータおよびR行列法の基底関数とFD手法の組み合わせは、効率的かつ正確に、高強度短パルス光に対する多電子原子系の全応答を記述し、これまでに既存の手法を超える能力を持つことが証明されてきました。

このために、2009年初頭にQueen’s University Belfast $QUB$,University College London $UCL$, the Open University and STFC CSEの4機関により5年間のEPSRCソフトウェア開発予算 $the UKRAMP Project$が実行され、電子-原子および電子-分子散乱と時間依存/非依存レーザー-原子分子相互作用の非経験的理論の構築が行われました。そのプロジェクトの主要な成果は、HELIUMコードで高強度短パルス光に対する多電子原子分子系の全応答を記述可能にすることです。この目的の第一の主要かつ重要なアルゴリズム開発が、HELIUMコードにR行列法を組み込み、高強度短パルス光に曝された一般的な多電子原子分子の一価イオン化応答を記述することでした。

·プロジェクトの目的

プロジェクトの目的は、高強度短パルス光に曝された一般的な多電子原子分子の一価イオン化応答を正確かつ効率的に計算するアプリケーションを開発することです。これは一価イオン化だけでなく多価イオン化を記述できる能力を持ちます。またこれは一電子水素原子に対する基底関数手法とFD手法をインターフェイスした最近の仕事[26]を大きく拡張します。このために以下の3つの作業項目が必要です（24人月以上かかります）：

1. 最初に、原子核近傍の有限領域（R行列法では「内殻領域」と呼ばれます）における多電子構造を記述することが必要です。このために、高精度かつ定評のあるシリアルR行列法コードRMATRIX2[28]を用います。RMATRX2は、HECToR[30]での利用が現在求められているPRMAT[29]のパーツです。RMATRX2は基底関数手法を用いて、角運動量のある範囲における内殻領域内の一般的な多電子原子の固有値と固有ベクトルを解き、全固有値間の近似双極子相互作用行列要素を持ちます。このデータは非時間依存のため、新たな時間依存多電子コードの初期値として用いることが出来ます。最初にスタンドアローンのコードとして、内殻領域における時間依存の一般的な多電子原子の正確な応答を記述するよう構築しました（以降、これをTD-RAを呼びます）。この作業には9人月の作業工数を要しました。
2. 次に、内殻領域を超えて、より広い範囲の外殻領域におけるイオン化した電子の運動を記述する必要があります。この電子はレーザー場と残った多電子イオンの場の中を移動します。FD手法はこうした領域での波動関数記述に適した数値計算法です。最終的に、HYDROとして知られる時間依存HELIUMコードの1電子バージョンを用いて新たなスタンドアローンコード（TD-OUTERと呼びます）を構築しました。残りのイオンが存在する中で電子運動を記述するために、HYDROは長距離多重極ポテンシャル項を含めるよう大きく更新されました。この項もRMATRX2で生成された初期値の一部となります。外部レーザー場は放出された電子のみならず、残りのイオンでも感じるため、TD-OUTERにもこのレーザー項を追加しました。この作業には6人月の作業工数を要しました。
3. 最後に、これら2つのコードを一つのコードRMT（時間依存R行列）に統合します。RMTの構築においてはArnoldiプロパゲータの実装に特に重点を置きました。これはMPPアーキテクチャの利点を活用します。この作業では、1電子手法プロトタイプに波動関数を時間発展させるアルゴリズムの拡張を中心的に行いました[26]。元の手法では、内殻と外殻領域内の波動関数の時間発展には、テイラー展開法が用いられていました。TD-RAとTD-OUTERの統合版ではArnoldiプロパゲータが実装され、内殻領域と外殻領域間の波動関数の受け渡しに用いられます。この作業の中でコードはMPI並列化され、HECToR Cray XT6へ移植されました。レーザー-原子相互作用に残りのイオンの複数状態を含めた、800nmおよびアト秒XUVレーザーパルスに対するネオン原子の応答がHECToRで計算可能になりました。この作業には9人月の作業工数を要しました。

·プロジェクトの主要な成果

このプロジェクトにおいてRMTコードが構築されました。RMTコードは、高強度レーザーパルスに曝された多電子原子の一価イオン化や原子イオンの多電子反応を正確に記述する非経験的計算コードです。

RMTコードの際立った特徴の一つは、FLASH$the future European XFEL facility$および米国LCLSのような第4世代の光源のXUV周波数からTi:sapphireレーザーの赤外周波数までの光に曝された原子に特化して設計されていることです。RMTコードは幅広い周波数に対応できるだけでなく、その時間依存特性から、短いものから長い光パルス幅に対する原子系の応答を記述することが可能です。さらに光場の顕わな表現を用いることから、レーザー光に対する原子の応答上の時間的なパルス形状変化の効果を見る実験などにも幅広く対応可能です。

一般的な多電子原子のアト秒スケールのダイナミクスを記述するRMTコードの能力は、英国の研究グループ全てに計り知れない利点を持ちます。高強度アト秒光パルスは、分子内の電荷移動の実時間観測と制御で期待されており、特に分子エレクトロニクス[15]と癌研究[17]に大きなインパクトをもたらすと予想されます。

HELIUMの多電子原子への拡張：RMT法

出発点は時間依存シュレディンガー方程式です。

R行列法[27]においては電子が存在する空間は2つの領域に分割されます。内殻領域$領域1$は核周囲にあり、多電子波動関数が構築され、原子-レーザーハミルトニアンが陽に計算されます。外殻領域$領域2$では1つあるいは多くて2つの電子が存在する様に選ばれ、その電子はレーザー場と共に残りの原子からの長距離の多重極相互作用の影響を受けます。通常のR行列理論では、時間は陽に含まれません。領域1では、時間依存波動関数は場の無いハミルトニアンのR行列固有状態で展開されます。またその波動関数はN+1個の電子の動径座標を変数とします。領域2では時間依存波動関数は等間隔の有限差分グリッドポイント上の値で定義されます。

·TD-OUTERの実装

TD-OUTERの設計:
空間分割の多くの利点は、原子の多電子表現を空間領域に限ることで導かれるもので、特に核近傍の領域1において必ず必要です。こうすることで領域2では電子数を1に制限され、残りのイオンの各状態のTDSEの次元は高々3になり、計算は簡便になります。
TD-OUTER領域では、波動関数を放出された電子の動径変数を持つ関数Fp$r,t$と、残りのイオンのターゲット状態チャネル$p$の関数の積の線形結合で表します。解くべき方程式はFp$r,t$以外の変数を積分したTDSEです。

TD-OUTERスタンドアローンコードの詳細：
これはHELIUMの1電子バージョンHYDROを用いて構築されました。HYDROは1電子波動関数に時間発展に正確かつ効率的なArnoldiプロパゲータとFD法を用います。HYDROはFortran90で記述され、HELIUMと同じMPI並列実装を用います。TD-OUTERは多電子原子存在下での放出された1電子を記述するため、ハミルトニアンの相互作用項の計算には多電子原子波動関数データが必要です。このデータはRMATRX2[28]を用いて計算します。

HYDROは1電子TDSEを全領域で解くため、TD-OUTER用に、内殻境界距離bから外殻境界距離Rまでの範囲でグリッドを切りなおします。

HYDROは水素原子のTDSEを解くため、残りのイオンの波動関数を扱えません。このため、TD-OUTERは、残りのイオンが形成する長距離多重極ポテンシャルから力を受ける「有効」1電子TDSEを解くように設計し直します。放出された電子は、その角運動量や残りのイオンの対称性、動径波動関数Fp$r,t$自身による様々なやり方でイオンと相互作用します。また、相互作用項を通じてチャネル間でもカップリングします。

·TD-RAの実装

領域1ではN+1電子波動関数は、R行列固有状態基底関数の線形結合で与えられます。この線形結合係数の時間依存を解きます。

このTDSEは、量子ダイナミクスの多くの領域で見られる非同次方程式です。こうした方程式の解法に高次の陽解法が適用されたのはごく最近のことです。このような高次のプロパゲータの一つに、テイラー展開による近似があります[32]。これは1電子の実装に対する既存の手法です[26]。TD-RAで採られた時間発展手法は、標準的なArnoldi法とよく似た手法です。ここでは領域2で計算される非同次項を除いた計算を行います。こうして標準的な同次Arnoldi法をR行列基底関数と組み合わせて実装します。

最初のステージで、RMATRX2が生成する行列要素を用いて非時間依存の多電子ハミルトニアンを構築します。次に多電子波動関数の時間発展用の標準的なArnoldi法を実装します。その後、時間依存ハミルトニアンを用いて、Arnoldi法により多電子波動関数の時間発展をさせます。

·TD-RAとTD-OUTERの統合：RMTプログラムコード

以下の3機能を実装します：
A:領域1の多電子波動関数の時間発展は、領域2（との境界）において波動関数の時間微分が得られれば開始可能なため、この非同次項の時間依存境界での計算を実装します。
B:領域2の1電子波動関数の時間発展は、領域1のFDポイント上の波動関数が得られれば開始可能なため、領域1内の波動関数を陽に求めます。
C:領域1のArnoldiアルゴリズムを、非同次項の時間発展ができるように拡張します。

RMTコードの並列化

並列RMTコードは、BLAS,LAPACKおよびMPIを用いて実装されました。並列プロセスは各時間ステップで、１）領域1のTD-RAグループと、２）領域2のTD-OUTERグループに割り当てられます。2つのグループ間のデータ交換は、各グループのマスターコア間の同期1対1通信を通して双方向に発生します。このやり方はスケール性が良く、計算負荷は適切にバランスします。計算結果は周期的にディスクに保存されます。

領域2での並列化：
TD-OUTERグループはグリッドポイントをコアに割り当てて並列化します。各コアはチャネル数×グリッド数の波動関数を扱います。ハミルトニアン内の微分演算は5点FDを用いるため、コアは再隣接通信のみが発生します。

領域1での並列化：
TD-RAの並列化には、計算負荷の高い2つのルーチンに焦点が当てられました。それは、非同次項のArnoldi法とFDポイント上の波動関数計算で呼ばれる行列ベクトル積サブルーチンです。

TD-RAとTD-OUTER間のMPI通信：
通信はTD-RAのマスターコア（MCR1)とTD-OUTERのマスターコア（MCR2)間で行われます。この通信はMPIのSEND/RECV 同期通信をベースにし、以下の2種類のシーケンスで行われます。

·TD-RAからTD-OUTERへの通信
各時間ステップの最初に、MCR1は境界内のFDポイント上の動径波動関数から成る配列をMCR2へ送信します。MCR2はこの配列を受信し、TD-OUTERで時間ステップを１つ更新します。つまり1つの時間ステップでは、1つのSEND/RECV呼び出しが発生します。

·TD-OUTERからTD-RAへの通信
各時間ステップ内で、MCR2からMCR1へのベクトルUjの送受信は複数回発生します（通常は12-14回）。

負荷バランス処理：
TD-RAのアイドル時間を最小にするには、MCR2は、MCR1が前回の計算をしたらすぐに次の計算が開始できるように同期してベクトルUjを送信する必要があります。RMTにはそれを調節するパラメータとして、TD-RAグループのコア数とTD-OUTERの各コアに割り当てられるFDポイント数があります。もしこれらを十分大きくとるとTD-RA内のアイドル時間が増加します。逆に小さくするとTD-OUTER内のアイドル時間が増加します。最適なサイズはCrayPatツールを用いた計測実験により選択が可能です。

RMTコードのHECToR上の性能

RMTおよびRMATRX2はHECToRへ2011年1月に移植されました。コンパイルはPGIで行いました。

最初のベンチマークでは、TD-OUTERの最大コア数を384としました。またTD-RAの最大コア数は48としました。この時のハミルトニアン次元数は約35000×35000です。

·弱スケール性能

RMTコードはTD-OUTERに割り当てられたコア数に対して線形にスケールすることが示されました。FDポイント数を4倍に増やすと時間当たりの繰り返し数は約4倍に増加しており、明らかな弱スケール性が示されました。

·強スケール性能

グリッドポイント数を固定した場合に、TD-OUTERに割り当てるコア数に対してRMTは線形にスケールしました。ここで、外殻サイズは7200a.u.、FDグリッド間隔は0.2a.u.に固定しました。全グリッド数は36000です。TD-RAに24コアを割り当てた場合はアイドル時間のオーバヘッドがありましたが、48コアを割り当てた場合には理想的な結果になりました。

·負荷バランス

CrayPatツールでの分析から、24$48$コアをTD-RAに割り当てた場合の最適なFDポイント数は144$89$であることが判りました。

·RMTによる実計算例

800 nm Ti:Sapphire レーザー場の存在下における、ネオン原子とXUVアト秒光パルスの相互作用を計算しました。実験においては予期しない2p電子と2s電子の放出時間差が見出されており、これが複雑な多電子効果であると推測されています。アト秒レーザーストリーキング法（attosecond streaking）を用いた実験ではこの時間差が21アト秒と見積もられていました[7]。これまでの理論ではこの実験値を再現することが出来ていません。

RMTは一回の計算で2pと2s電子のイオン化計算が可能な、世界で数少ないコードの一つであることが実証されました。

この結果の詳細は、XXVII International Conference on Photonic Electronic and Atomic Collisions $ICPEAC$ 2011 conferenceおよび論文（予定）で報告しています。

謝辞

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

文献

[1] P. Corkum and F. Krausz, Nat. Phys. 3, 381 $2007$
[2] F. Krausz and M. Ivanov, Rev. Mod. Phys. 81, 163 $2009$
[3] S. Baker et al Science, 312, 424 $2006$
[4] L. E. Chipperfield et al, Phys. Rev. Lett. 102, 063003 $2009$
[5] M. Uiberacker et al, Nature 446, 627 $2007$
[6] E. Goulielmakis et al, Nature 466, 739 $2010$
[7] M. Schultze et al, Science 328, 1658 $2010$
[8] N. Dudovich et al, Nature Physics 2, 781 $2006$
[9] H. J. W?rner et al, Nature 446, 604 $2010$
[10] K. P. Singh et al, Phys. Rev. Lett. 104, 023001 $2010$
[11] H.Wang et al, Phys. Rev. Lett. 105, 143002 $2010$
[12] Y. Zeng et al, Phys. Rev. Lett.103, 043904 $2010$
[13] Y. Nabekawa et al, Phys. Rev. Lett.102, 213904 $2009$
[14] E. Takahashi et al, Phys. Rev. Lett.104, 233901 $2010$
[15] J. Vura-Weis et al, Science 328, 1547 $2010$
[16] M. Stockman et al, Nature Photonics 1, 539 $2007$
[17] H. B. Gray and J. Halpern, Proc. Nat. Ac. Sci. 102, 3533 $2005$
[18] P. Corkum, Phys. Rev. Lett. 71, 1994 $1993$
[19] T. Brabec and F. Krausz, Rev. Mod. Phys 72, 545 $2000$
[20] M. Lewenstein et al, Phys. Rev. A 49, 2117 $1994$
[21] G. Sansone et al, Science 314, 443 $2006$
[22] W.E. Arnoldi, Quart. Appl. Math. 9, 17 $1951$
[23] T.J. Park, J.C. Light, J. Chem. Phys. 85, 5870 $1986$
[24] E.S. Smyth, J.S. Parker, K.T. Taylor, Comput. Phys. Commun. 114, 1 $1998$
[25] J.S. Parker et al, Phys. Rev. Lett. 96, 133001 $2006$
[26] L.A.A. Nikolopoulos, J.S. Parker, K.T. Taylor, Phys. Rev. A 78, 063420 $2008$
[27] P.G. Burke, K.A. Berrington, Atomic and Molecular Processes: An R-matrix Approach, $IOP Publishing,Bristol$ $1993$
[28] P.G. Burke, V.M. Burke, K.Dunseath, J. Phys. B: At. Mol. Opt. Phys. 27, 5341 $1994$
[29] A. G. Sunderland et al, Comp. Phys. Commun. 145, 311 $2002$
[30] A. G. Sunderland, C. J. Noble and M. Plummer, dCSE Technical Report $2010$
[31] B. Skaflestad, W.M.Wright, Appl. Numer. Math. 59, 783 $2009$
[32] M. Ndong, H. Tal-Ezer, R. Kosloff. C.P. Koch, J. Chem. Phys. 130, 124108 $2009$
[33] R. Sidje, ACM Trans. on Math. Soft. 24, 130 $1998$
[34] Y. Saad, SIAM J. Numer. Anal. 29, 209 $1992$
[35] J.S. Parker et al, J. Phys. B: At. Mol. Opt. Phys. 33, L239 $2000$
[36] D.H. Glass, P.G. Burke, J. Phys. B: At. Mol. Opt. Phys. 32, 407 $1999$
[37] P.G. Burke and K. T. Taylor, J. Phys. B: At. Mol. Opt. Phys. 8, 2620 $1975$