もっとも基本的な分子動力学ワークフロー
`examples/lj/in.lj.nve` サンプルを使って, LAMMPS2001 が分子動力学計算をどのように準備・ループ・出力するかを, データファイル・主要サブルーチン・数式の観点から追跡します.
まず押さえておきたい用語
以下のサンプルでは分子動力学 (Molecular Dynamics; MD) の基礎概念が多数登場します. 段階を追って読むためのキーワードを先に整理しておきましょう.
- 分子動力学 (MD)
- 原子・分子を質点として扱い, ニュートンの運動方程式 \( m \frac{\mathrm{d}^2 \mathbf{r}}{\mathrm{d} t^2} = \mathbf{F} \) を数値積分する手法.
- Lennard-Jones (LJ) ポテンシャル
- 希ガスなどを近似するペアポテンシャル. 反発と引力を \(\left(\sigma / r\right)^{12}\) と \(\left(\sigma / r\right)^6\) の項で表します. LAMMPS の
units ljは質量・長さ・エネルギーを無次元化した単位系です. - NVE アンサンブル
- 粒子数 (N)・体積 (V)・内部エネルギー (E) が一定の扱い. 外部の熱浴とは接続せず, 総エネルギー保存を重視します.
- 温度制御 (Thermostat)
- 速度スケーリングなどで温度を一定に保つ操作. サンプルでは単純な速度リスケールを NVE に重ねて目標温度を維持します.
- 近接リスト (Neighbor List)
- 一定距離内に入る原子ペアをまとめたテーブル. 力計算を全組み合わせではなく近傍だけに限定し, 計算量を削減します.
対象とするサンプルと目的
本稿では LAMMPS2001 の examples/lj/in.lj.nve を題材に, 典型的な非結合 Lennard-Jones 系の NVE 計算がどのように進むかをステップごとに解説します. 前節で整理した用語を念頭に置きながら, 以下の観点を押さえることで Fortran 実装の読み解きが加速します.
- どのデータファイルを読み込み, どの配列に配置するか (
read_data.f,setup.f) - どのループで空間分割・近接リスト・力計算が実行されるか (
integrate.f,neighbor.f,force_many.f) - どの数式が使われ, MathJax で表したときに何を意味するか
- どの出力が得られ, 何をデバッグ指標にできるか (
thermo.f, ログファイル)
入力スクリプトとデータファイル
入力スクリプト (examples/lj/in.lj.nve)
最小限の NVE 計算を構成するコマンドはつぎのとおりです.
# LJ system with NVE ensemble
units lj
maximum cutoff 2.5
neighbor 0.3 1 1 5 1
read data data.lj
create temp uniform 1.0 87287
coulomb style none
nonbond style lj/cutoff 2.5 1
nonbond coeff 1 1 1.0 1.0 1.12246
...timestep 0.005
temp control rescale 1.0 1.0 20 0.02 1.0
thermo flag 50
thermo style 1
run 500input.f が上から順に解釈し, 必要なサブルーチンを呼び出します. たとえば read data は read_data, nonbond style は force_many.f 内で使われるフラグ類を設定します.
データファイル (examples/lj/data.lj)
read_data.f はセクション見出しを読み取ってから x, type, mass などの配列に値を格納します. 最初の部分は以下のとおりです.
LAMMPS Description
2048 atoms
3 atom types
-7.52829 7.52829 xlo xhi
-7.52829 7.52829 ylo yhi
-7.52829 7.52829 zlo zhi
Masses
1 1
2 1
3 1
Atoms
1 0 3 0 -7.52829 -7.52829 -7.52829
... 座標は reduced LJ 単位系で, 境界条件は 3 次元周期. read_data は scan_data を通じて各セクションの行数を先読みし, のちほど setup が原子配列サイズと領域情報を確定します.
プログラム側の全体フロー
プログラムは lammps.f の program lammps から始まり, 入力データの読み込み→シミュレーション設定→時間発展→出力の順に流れます. 大まかな関係を流れ図で整理すると次のようになります.
入力スクリプト ─┐
├─▶ input.f ──┐
データファイル ─┘ │
▼
read_data / setup
│
▼
初期配列・条件差し込み
│
▼
integrate (時間ループ)
│
┌──── 力計算 ────┴─── 近接リスト更新 ──┐
│ │
▼ ▼
nve_v / nve_x neighbor.f
│ │
└──────────── thermo.f 出力 ────────────┘
ステップ別の処理
- 初期化段階
startとinitializeが MPI 環境, 乱数シード, データ構造を初期化します.setupはドメイン分割の指標や送受バッファを確定し,globalモジュールの配列サイズを決めます. - 入力解析
input.fがスクリプトを 1 行ずつ読み取り, キーワードごとに専用サブルーチン（例：read_data,set_style_nonbond,timestep）を呼び出します. ここで力場係数やシミュレーション条件がglobalの配列に格納されます. - 時間発展ループ
integrate.fのdo itime = 1, nstepsがメインループです. 各ステップで位置・速度更新, 近接リスト検査, 力計算, 温度操作, 統計出力の順に処理が進みます. - 終了処理
finish.fが累積統計を出力し, MPI リソースを解放します. ログに Loop time などのプロファイルがまとめられるのはこの段階です.
時間ループ内部のチェックポイント
時間ループでは, integrate が力計算前後で速度を 2 回更新し, 途中で check_neighbor が原子の移動量を監視します. 要求があれば neighbor に処理が渡り, 新しい近接リストが作成されてから force_many などの力計算サブルーチンが再開されます.
続くセクションでは, このループの内部構造と近接リスト生成アルゴリズムをコード断片とともに詳しく掘り下げます.
主要ループと生成されるリスト
時間ステップループ (integrate.f)
integrate では以下の順序で計算が行われます.
- 半ステップ速度更新 —
nve_vが前ステップで求めた力f(:, i)を用いて速度を半ステップだけ更新します. ここで「半分」にとどめるのは, 後半に新しい力で仕上げの更新を行うためです. - 位置更新 —
nve_xが更新済み速度を使って位置x(:, i)を \(\Delta t\) だけ進めます. 周期境界条件ではボックス外に出た粒子をpbcが折り返します. - 近接リストの判定・再構築 —
check_neighborが原子の移動距離を計測し, スキン距離を超えた場合だけneighborを呼び出してリストを再構築します. 不要な再計算を避ける省力化ポイントです. - 通信とゴースト原子更新 — 並列計算では
communicate/exchange/bordersが各プロセス間で原子をやり取りし, 力計算に必要なゴースト原子情報をそろえます. シリアル実行でも同じ関数が呼ばれますが通信量はゼロです. - 力計算 —
zero_force_virialで配列を初期化後,force_many.fなど力場ごとのサブルーチンが近接リストを走査してポテンシャルと力を求めます. - 温度調整 — 本サンプルでは 20 ステップごとに
temp_rescaleが現在温度をチェックし, 必要に応じて速度の大きさだけをスケールします. - 後半の速度更新 — 新しい力を使って再び
nve_vが呼ばれ, 速度をもう半ステップ前進. ここで初めて\(\mathbf{v}(t+\Delta t)\)が決まり, 次のステップに進みます.
近接リスト (neighbor.f)
原子が互いに与える力は距離が離れると急速に小さくなるため, 一定距離内の組み合わせだけを調べる「近接リスト」を作ってから力計算を行います. 処理の流れは次の通りです.
- 移動量の評価 —
check_neighborが前回リスト作成時の座標xholdと現在座標を比較し, スキン距離neighbor 0.3 ...で設定したtriggersqを超える原子があるか判定します. - 位置情報の保管 — リストを更新する場合, 現在座標を
xholdにコピーして次回判定に備えます. - 原子ペアの探索 —
neighborは指定された方式に応じてneighbor_nsq（全対比較）かneighbor_bin_*（セル法）を呼び出します. 本サンプルではneighstyle = 0のためneighbor_nsqが使われます. - Newton 条件と特殊結合の扱い — 作用反作用を 1 回だけ計算する「Newton 法」が有効なら, ゴースト原子とのペアを片側だけに登録します. また
special配列を使って 1-2・1-3 などの結合近接に重みをかけたり除外したりします.
neighbor_nsq の核となる二重ループは以下のような擬似コードで表せます.
npnt = 0
do i = 1, nlocal
nnfirst(i) = npnt + 1
do j = i + 1, nlocal + nghost
if (newton rule でスキップする条件) cycle
rsq = |x(i) - x(j)|^2
if (rsq <= cutneighsq(itype, jtype)) then
if (special pair) then
npnt = npnt + 1; nlist(npnt) = scaled index
else
npnt = npnt + 1; nlist(npnt) = j
endif
endif
enddo
nnlast(i) = npnt
enddoこのループで作成されたテーブルは以下の配列に格納され, force_many.f の力計算ルーチンで逐次参照されます.
nnfirst(i)/nnlast(i)— 原子iの隣接ペアがnlistのどこに並んでいるかを示す開始・終了インデックス.nlist(npnt)— 条件を満たした隣接原子の ID または結合種別情報.cutneighsq— 原子タイプごとの「力を計算する最大距離」の二乗値. スキン距離だけ余裕を持たせた値が設定され, 頻繁なリスト再構築を避けます.xhold— 最後にリストを作ったときの座標を保存し, 次回の移動量チェックに使用.
スキン幅を大きくするとリスト更新頻度が下がる一方で nlist の要素数が増えます. バランスを確認するにはログファイル末尾の「Number of reneighborings」や「Total # of neighbors」を指標にすると良いでしょう.
力計算の数式と実装
非結合相互作用は force_many.f → force_lj 系ルーチンで計算されます. Lennard-Jones ポテンシャルは以下の式です.
$$ U_{ij}(r) = 4 \varepsilon_{ij} \left[ \left( \frac{\sigma_{ij}}{r} \right)^{12} - \left( \frac{\sigma_{ij}}{r} \right)^6 \right], $$ $$ \mathbf{F}_{ij} = -\frac{\mathrm{d} U_{ij}}{\mathrm{d} r}\, \hat{\mathbf{r}}_{ij} = 24 \varepsilon_{ij} \left[ 2 \left( \frac{\sigma_{ij}}{r} \right)^{12} - \left( \frac{\sigma_{ij}}{r} \right)^6 \right] \frac{\hat{\mathbf{r}}_{ij}}{r}. $$
force_many.f 内では平方距離 rsq を使って係数を事前計算し, Newton の作用反作用を利用して f(:, i) と f(:, j) を同時に更新します. カットオフ半径は nonbond style lj/cutoff 2.5 で与えられ, rsq <= cutsq(itype, jtype) の条件で力が計算されています.
速度ベルレ法
NVE 積分器は速度ベルレ法 (velocity Verlet) を採用しています. MathJax で書くと以下の手順です.
$$ \mathbf{v}\left(t + \frac{\Delta t}{2}\right) = \mathbf{v}(t) + \frac{\Delta t}{2m} \mathbf{F}(t), $$ $$ \mathbf{r}(t + \Delta t) = \mathbf{r}(t) + \Delta t\, \mathbf{v}\left(t + \frac{\Delta t}{2}\right), $$ $$ \mathbf{F}(t + \Delta t) = -\nabla U\left(\mathbf{r}(t + \Delta t)\right), $$ $$ \mathbf{v}(t + \Delta t) = \mathbf{v}\left(t + \frac{\Delta t}{2}\right) + \frac{\Delta t}{2m} \mathbf{F}(t + \Delta t). $$
実装上は dthalf = 0.5 * dt を共通化し, nve_v と nve_x が上記の半ステップ更新を担当します.
温度制御と補助操作
NVE アンサンブルでは本来エネルギーが保存されますが, 初期速度の揺らぎや有限精度の影響で温度が目標値から外れる場合があります. サンプルでは temp control rescale 1.0 1.0 20 0.02 1.0 を指定し, 20 ステップごとに単純な速度リスケールをかけています.
ensemble.f 内の temp_rescale は瞬時温度 \(T\) と目標温度 \(T_0\) からスケール係数 \(s\) を次の式で求め, すべての速度ベクトルに乗じます.
$$ s = \sqrt{ \frac{T_0}{T} }, \qquad \mathbf{v}_i \leftarrow s \, \mathbf{v}_i. $$
リスケールが実行されるタイミングでは以下のチェックが走ります：
- 「最小ステップ間隔」パラメータが満たされているか（上記コマンドの第 3 引数）.
- 温度偏差が許容範囲（第 4 引数）を超えているか.
- 速度を調整する原子グループ（ここでは全原子）に合わせて再分布するか.
この操作はエネルギー保存を厳密には崩しますが, 初心者が短時間で目標温度に合わせたい場合に便利です. より高度な熱浴（Langevin, Nose-Hoover など）を試す場合は入力スクリプトの fix コマンドを選び替えることで変更できます.
出力と検証
thermo flag 50 により 50 ステップごとに thermo.f がエネルギーや圧力を表示します. thermo style 1 の出力例は以下です.
Step Temp E_nbond E_bond E_long E_total Pressure Volume
0 1.000000 1.328286 0.000000 0.000000 2.828286 7.505132 3413.3357
50 0.997381 1.731314 0.000000 0.000000 3.227385 9.520529 3413.3357
... 計算終了後, ログファイル examples/lj/log.in.lj.nve. * にはリスト再構築回数や通信時間がまとまっており, neighbor_nsq が生成した隣接ペア数（Total # of neighbors）も確認できます.
同じ設定で自分のビルドを走らせたい場合は, ソースをコンパイルしたディレクトリで以下を実行します.
lmp_serial < examples/lj/in.lj.nve（実行バイナリ名はビルド設定に依存します. MPI 版の場合は mpirun -np 4 lmp_mpi < ... のように呼び出してください. ）
押さえておきたいポイント
- データ読み込み → 近接リスト → 力計算 → 積分 → 出力の順に Fortran サブルーチンが連携している.
- 隣接リストの構造 (
nnfirst/nnlast/nlist) を理解すると力計算ループの読み解きが容易になる. - 速度ベルレ法の式を追いながら
nve_v/nve_xを読むと, 温度制御などの追加操作の位置づけが明確になる. - ログ出力に含まれるペア数・通信統計はデバッグだけでなく負荷見積もりにも役立つ.