LAMMPS 実装編

※注意 ここに書かれた内容は、LAMMPSの開発にともない古い情報となる可能性があります

はじめに 言語 LAMMPSはC++で書かれている C++が読めるエディタ(シンタックスハイライト、入力補完など) の使用を推奨 ソース 関連するファイルはほぼ全てsrcディレクトリ直下に存在 ただし、デフォルトのLAMMPSは最小構成なので、場合によりパッケージを 入れたコンパイルが必要になる(多体ポテンシャルのmanybodyなど)

Overview ファイルの単位 モジュール側から見たLAMMPS全体 ほぼ1コマンド1ファイルとみてよい (pairはstyleとcoeffを合わせて1つ) モジュール側から見たLAMMPS全体 各機能に対して系全体の情報をもつ”LAMMPSクラス”(のポインタ)が渡される 例えば原子の位置や速度、定義されたfixの情報などにこのポインタを通してアクセスできる fix_nve(速度ベルレ法の実装)やcompute_temp(温度の計算)などを見るとわかりやすい

C++クラスによる実装 各モジュールはC++のクラスで実装されている 例えばfixコマンドは”fix抽象クラス”の継承として、ポテンシャルは”pair抽象クラス” の継 承として実装 抽象クラス自体を定義したファイルも当然ある。例えばfix.hやpair.h LAMMPS側は”fixクラスの集合”, “pairクラスの集合”…のように把握している LAMMPS側からは適切なタイミングで対応する関数を呼ぶので、クラス側で は対応する仮想関数をオーバーロードしておけばよい 例えば、fix_setforce.cppではinit()で初期化の処理、post_force()で力の計算後の処理を 行っている 用語がわからない場合は、先にC++の「継承と多態性」を調べると楽かもしれない

LAMMPSのrunコマンドの中身 右はrunが実行されたときの擬似コード loop over N timesteps: ev_set() fix->initial_integrate() fix->post_integrate() nflag = neighbor->decide() if nflag: fix->pre_exchange() domain->pbc() domain->reset_box() comm->setup() neighbor->setup_bins() comm->exchange() comm->borders() fix->pre_neighbor() neighbor->build() else comm->forward_comm() force_clear() fix->pre_force() pair->compute() bond->compute() angle->compute() dihedral->compute() improper->compute() kspace->compute() comm->reverse_comm() fix->post_force() fix->final_integrate() fix->end_of_step() if any output on this step: output->write() 右はrunが実行されたときの擬似コード (LAMMPS Webページ マニュアルより) 基本的には力の計算(compute())を行っていて、 その前後でfixが担当する様々な計算を行っている 位置の更新もfixで行う(nveなど)

LAMMPSとマルチプロセス LAMMPSはマルチプロセスで動くことを前提にして設計されている LAMMPS内ではこの単位は”local”という用語で表現される 例えば、”プロセス内の原子数”はnlocalと書かれる local local local local local local

LAMMPSとマルチプロセス 1プロセスが持つ原子の情報は”領域内の原子”と”領域近傍の原子”の2種類 近傍原子がないとそれらとの相互作用が計算できないため これらはlocal, ghostとして管理されている。ghost原子はどこか別プロセスのlocal原子 (1プロセスで周期境界なら自分自身かもしれない) 通常、境界の向こう側の別プロセスとの処理をプログラマが気にする必要はない localとghostは各ステップで自動で構成してくれる 例えば、原子が移動して境界を飛び越えた場合、自動的にlocal原子から外してくれる ghost原子もカットオフ距離に応じて用意してくれる local ghost

ポテンシャルの実装 pairクラスを継承して作る pair_lj_cut.cppの実装が簡単で勉強しやすい 最低限やること: computeが呼ばれたときにatom->f (力の配列)を埋める 必要に応じてエネルギー、応力を計算する (大抵はev_tallyを呼んでやってもらう) pair_lj_cut.cppの実装が簡単で勉強しやすい 原子i, jについてループを回し、距離を計算し、力を計算し…というなじみのある実装 コードの書き方がC++というよりC89風 原子のリスト、各原子について近傍原子のリストなどが定義されていることがわかる localやghostの意味については前ページ参照 (rRESPAの実装のためコードが長くなっているが、computeのみ見れば良い)

ポテンシャルの実装: 発展編 (計算した)周囲の原子の値が必要なポテンシャル commクラスにより明示的に通信が行える 例: pair_eam.cpp EAM関数(MANYBODYパッケージ内にある) 力の計算には局所電子密度ρの値(とF(ρ)の微分)が必要だが… ghost原子のF(ρ)の微分値は自プロセスの計算だけではわからない →そこで、手動で隣接プロセスと通信を行いF(ρ)の微分値を交換する commクラスにより明示的に通信が行える commクラスの関数が呼ばれると、引数に渡されたクラスのpack_**_commを呼び出し、そこで作られた データを他プロセスに通信し、受け取ったデータをunpack_**_commに渡す forwardはlocal原子のデータを周囲のプロセスの対応するghostへ配る、 reverseは周囲のプロセスのghostにあるデータを対応する自らのlocal原子に集める 何を送るかは自分で決める(pack_**_commなどに書く) ghost local

ポテンシャルの実装: 確認編 (数学的にはエネルギーだけで良いとはいえ)実装にあたってはエネルギーと力と圧力を それぞれ別に計算するため、それらに齟齬があると物理的におかしくなる まずはNVEアンサンブルでエネルギーが一定かどうかを確認 エネルギーが発散していくor減衰していく→微分の計算を間違えた可能性 エネルギーがたまに飛ぶ、ずれる→カットオフが変な可能性(エネルギー曲線が折れている、不連続である) まれに原子が吹っ飛んでロストする→距離0でエネルギーが-∞のような不自然な関数形状になっている可能性 格子サイズを含めた構造最適化計算も行ってみる 圧力平衡点とエネルギー最小点がずれる→圧力の計算を一部抜かした可能性 収束はするが微妙に振動したりして遅い→カットオフに依存する可能性(1階あるいは2階微分が不連続など) 微分計算が重いので近似値を使い、「NVTアンサンブルしか行わないから細かい点は気にしない」という 思想のポテンシャルもある

拡張パッケージにする方法 USER-**ディレクトリを参考にするとよい 新しいディレクトリ内に必要なソースファイルを置く ソースファイルの内容が他とかぶらないようにする ファイル名、クラス名(+コンストラクタ・デストラクタ名)、インクルードガードの文字列、 LAMMPSで呼び出されるコマンド(pairなら.hファイル上部の#ifdef PAIR_CLASSと#elseの間) ディレクトリ内にInstall.sh(他のパッケージ参照)を置き、中にファイルを指定する →make yes-USER-PACKAGE-DIRという形で組み込めるようになる

高速化Tips 一般的なプログラムと同様にLAMMPSにプロファイラが使える 呼ばれる回数・時間の多い箇所が特定できる ボトルネックは特定箇所に集中することが多く、少ない労力で実行時間が改善できるかもしれない 例えばテーブルの使用、メモ化、関数の置換など とはいえ、(キャッシュヒットや分岐ハザードなど)凝った最適化をする前に、 以下の点を検討することをおすすめします コンパイラ・コンパイルオプションを変える カットオフを短くする タイムステップを大きくする