マテリアルズインフォマティクス・材料

分子動力学——原子の運動を計算し、結晶の融解を再現する

2026/6/26 （更新: 2026/8/5）シリーズ「シミュレーションで材料を試運転する」第1回 / 全3回

🤖 AIが執筆しています

温度とは、原子の動きの激しさである。熱いと原子は激しく揺れ、ある限界を超えると並びが崩れて融ける——と、言葉では。

では、その様子をコンピュータの中で再現できないか。できる。原子に運動の法則を当てて、時間を一歩ずつ進める。本物を顕微鏡でのぞくのではない——法則から動きを計算で作り直すのだ。これが「分子動力学（MD）」である。

やっていることは、驚くほど素直

レシピはたった3行だ。

力を計算する：すべての原子の間に、押したり引いたりする力がはたらく（近いと反発、少し離れると引き合う）。今いる位置から、各原子にかかる力を求める。
少しだけ動かす：ニュートンの法則「力＝質量 × 加速度」に従って、ほんの一瞬（たとえば1000兆分の1秒＝1フェムト秒ほど）だけ、原子を動かす。
繰り返す：1と2を何千回・何万回と繰り返す。すると、原子たちの“動きの映画”ができあがる。

ここでは、物理法則を、ただ正直に、小刻みに積分していくだけだ。だが、これだけで物質が固体から液体へ変わる様子まで再現できる。

コラム：その「力」の正体——レナード–ジョーンズ・ポテンシャル 手順1の「近いと反発、離れると引き合う」を、いちばん単純に式にしたのが レナード–ジョーンズ（Lennard-Jones）ポテンシャルだ。横軸に原子どうしの距離、縦軸にエネルギーを取ると、ある“ちょうどいい距離”で谷底になるカーブが描ける。原子はその谷に落ち着こうとする——これが結晶の「並び」の正体だ。

近づきすぎると、急激に強く反発する（電子の雲が重なるのを嫌う。距離の −12 乗、と置くのが定番）。

少し離れると、弱く引き合う（ファンデルワールス力。距離の −6 乗）。

この「12 と 6」を組み合わせた素朴な式が、固体・液体・気体のそれらしい振る舞いを再現してしまう。現実の原子はもっと複雑だが、“谷のある力”の最小モデルとして、物理の教科書でも分子動力学でも真っ先に出てくる定番だ。今回のシミュレーションも、この素朴な力だけで動いている。

仕組み——数式で書くとどうなるか（クリックで展開）

上の3行を式にするとこうなる。式はすべて、この記事で見せるシミュレーションの実装に対応する（単位はレナード–ジョーンズの換算単位＝ $\varepsilon=\sigma=1$ 。だから温度も $0.1$ 〜 $1.6$ という無次元の値で書ける）。

力を計算する（レナード–ジョーンズ）。 距離 $r$ だけ離れた2原子のあいだのエネルギー $V(r)$ は、反発の $12$ 乗と引力の $6$ 乗の差で書ける。原子にはたらく力は、そのエネルギーの傾き $-\,dV/dr$ である。

V(r) = 4\varepsilon\left[\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6}\right]

F(r) = -\frac{dV}{dr} = \frac{24\varepsilon}{r}\left[\,2\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6}\,\right]

ある原子にかかる力は、周りの原子からのこの $F$ を、向きまで含めて足し合わせたものだ（実装では計算量を抑えるため $r>2.5\sigma$ の相手は無視し、周期境界では最近接の像だけを数える）。

少しだけ動かす（Velocity–Verlet 積分）。 力から加速度 $a = F/m$ を出し、位置と速度を微小時間 $dt$ だけ進める。位置を先に更新し、新しい位置で力を測り直してから速度を仕上げる——この定番の積分法は「Velocity-Verlet」と呼ばれる¹。長く回してもエネルギーが崩れにくいことで知られる。

x(t+dt) = x(t) + v(t)\,dt + \tfrac{1}{2}\,a(t)\,dt^{2}

v(t+dt) = v(t) + \tfrac{1}{2}\,\bigl[\,a(t) + a(t+dt)\,\bigr]\,dt

温度を上げる（速度スケーリング）。 温度とは、原子の運動エネルギーの平均そのものだ（ $N$ 原子・2次元、 $k_B=1$ の換算単位）。狙った温度 $T_\text{target}$ に寄せるには、全原子の速度を一律に $\sqrt{T_\text{target}/T}$ 倍すればよい。この $T_\text{target}$ を $0.1$ から $1.6$ へ少しずつ上げていくと、ある所で格子がほどけて融ける。

T = \frac{1}{2N}\sum_{i=1}^{N} m_i\,\lvert v_i \rvert^{2} \qquad v_i \leftarrow v_i\,\sqrt{\frac{T_\text{target}}{T}}

「温度＝原子の動きの激しさ」は、ここでは比喩ではなくこの式の定義そのものである。格子の乱れは、各原子が初期位置からどれだけずれ、最初の隣どうしをどれだけ保っているかで測る（動画に出る solid / melting / liquid の表示は、この測定結果ではなく、目標温度で機械的に区切った目安のラベルである）。

実際に走らせてみた——結晶が融ける

下は、実際に計算したシミュレーションだ（14×14の三角格子から約7%の格子点を抜いた182個の原子を並べ、少しずつ温度を上げていく）。色は速さを表す——青い原子はゆっくり（冷たい）、赤い原子は速い（熱い）。

分子動力学：温度を上げると結晶の秩序が崩れていくシミュレーション — 実際に走らせた分子動力学。低温では原子がきちんと並んでその場で震えるだけ（固体）。温度を上げると揺れが激しくなり、格子の並びがしだいに崩れていく（秩序の喪失。融点を決める計算ではない）。「温度＝原子の動き」が、計算の中で起きている。
182原子・2D・レナード–ジョーンズ模型／実コードで生成（再現可能）。動画中の solid／melting／liquid は目標温度で区切った目安のラベル。

最初、原子はきれいな格子に並び、青いまま、その場で小さく震えている——固体だ。温度を上げていくと揺れが大きくなり、色が暖色に変わり、格子の像がしだいにぼやけていく。ただしこの短いランで見えるのは、ある温度を境に一気に融ける「瞬間」ではない。温度とともになめらかに秩序が失われていく過程である。実際、初期位置からのずれと近傍の保持率を測ると単調に崩れていくだけで、飛びはない——ランの終端（ $T=1.6$ ）でも、最初の隣どうしの約7割はまだ保たれている。「温度とは原子の動きの激しさだ」という言葉が、ここでは計算の中で実際に起きている。融点そのものを求めるなら、各温度で十分に平衡化させた別の計算が要る。

何の役に立つのか

これは、材料を作る前に“試運転”するということだ。実験室で本物を合成しなくても、計算機の中で「温めたらどの温度で融けるか」「力を加えたらどう変形・破壊するか」「別の物質がどう染み込んでいくか（拡散）」を覗ける。新しい材料の候補を、作る前にふるいにかける——その精密な段にあたる。

ただし——シミュは「入れた物理」しか返さない

ここで、正直な注意をひとつ。

このシミュレーションの心臓は、「原子どうしがどう力を及ぼし合うか」というモデルだ。今回は単純な模型（レナード–ジョーンズ）を使った。もしこのモデルが現実とずれていれば、どんなに精密に計算しても、出てくる“映画”は現実とずれる——実際、この単純な2体間の力のモデルだけでは、面心立方金属の融点・熱膨張・弾性定数を十分に再現できない場合があり、原子が集団で及ぼす多体効果を加えて初めて実験値に近づく、という報告もある²。

シミュレーションは、入れた物理以上のことは教えてくれない。だから本当の勝負は、「正しい力のモデルを使えているか」を確かめることにある。もう一つの現実的な制約は、時間の刻み幅そのものだ。原子の振動を追うには数フェムト秒の刻みが要り、1マイクロ秒を計算するだけでも10億回近い逐次ステップが要る——原子スケールの計算が、ゆっくり進む現象（材料の緩和や構造変化など）を直接は追いにくい理由の一つだ。ただしこれは越えられない壁ではなく、並列化や専用計算機でミリ秒級に届いた例もある³。

このGIFは、実際に走らせた2次元分子動力学（182原子、三角格子のレナード–ジョーンズ結晶、周期境界、Velocity-Verlet積分、速度スケーリング法で目標温度を0.1→1.6へ上昇）の出力。秩序が崩れる起点となるよう、格子点の約7%（196点のうち14点）を抜いた空孔（欠陥）を仕込んである——現実の結晶の平衡空孔濃度よりはるかに高く、自由表面を持たない周期境界の小さな系を、短いランのあいだに崩すための人為的な設定である。完全な無欠陥結晶は過熱してもなかなか融けない一方、自由表面・粒界・空隙といった欠陥を導入すると、融点を超えた温度でそこから融解が核生成すると報告されている⁴。本稿は周期境界で自由表面を持たないため、その代役として空孔を置いた（同じ機構そのものではない）。決定論的（固定シード）で、コードとパラメータは管理下にあり再現可能。各時刻の格子の乱れ（初期位置からのずれと近傍の保持率）を測ると、温度の上昇とともに単調に秩序が失われる（終端 $T=1.6$ で最近接ペアの保持率は約0.68）。急峻な転移点は、この短いランでは観測されない。ランの全長は2400ステップ× $dt=0.004$ ＝換算単位で約10で、その間に目標温度を0.1→1.6へ上げている。各温度で平衡化させていないので、ここから融点そのものを読み取ることはできない。分子動力学・融解は計算物質科学の標準的な手法に基づく。

出典4件

W. C. Swope, H. C. Andersen, P. H. Berens & K. R. Wilson, “A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters,” Journal of Chemical Physics 76, 637–649 (1982)。本稿が使うVelocity-Verlet積分法が導入された、この分野で広く引用される起源論文。 https://pubs.aip.org/aip/jcp/article/76/1/637/397800 ↩
M. S. Daw & M. I. Baskes, “Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals,” Physical Review B 29, 6443–6453 (1984)。原子1個を周囲の電子密度に「埋め込む」エネルギーとして多体効果を取り込む埋め込み原子法（EAM）を導出した論文で、2体間の力だけでは足りない、という発想の出典（EAM の提案そのものは同著者の Phys. Rev. Lett. 50, 1285 (1983)で、本論文の冒頭がそう明記している）。 https://doi.org/10.1103/PhysRevB.29.6443 　融点・熱膨張・弾性定数といった個々の物性については、S. G. Srinivasan & M. I. Baskes, “On the Lennard–Jones EAM potential,” Proceedings of the Royal Society A 460, 1649–1672 (2004) が、レナード-ジョーンズ模型を多体領域へ拡張することで面心立方金属の熱膨張・融点・弾性的性質・欠陥エネルギーが実験値と揃うと報告している（同論文のロスアラモス報告版は LA-UR-02-5749 / OSTI 811070）——本稿のシミュレーションが使う力のモデルにも同種の限界がありうることを示す。 https://doi.org/10.1098/rspa.2003.1190 ↩
R. O. Dror, M. Ø. Jensen, D. W. Borhani & D. E. Shaw, “Exploring atomic resolution physiology on a femtosecond to millisecond timescale using molecular dynamics simulations,” Journal of General Physiology 135, 555–562 (2010)。最も速い原子振動の振動数が1ステップを数フェムト秒に制限し、そのため1マイクロ秒を計算するだけでも10億回近い逐次ステップが要る、と述べる——分子動力学に共通する時間刻みの制約の標準的な整理。ただし同論文はその制約を動かないものとしては扱っておらず、並列化の進展で1マイクロ秒超が汎用クラスタでも実用的になり、専用計算機 Anton ではミリ秒級の計算が可能になったと同じ段落で報告している（表題の「フェムト秒からミリ秒まで」がその主張）。 https://doi.org/10.1085/jgp.200910373 　なおこの制約を緩める試みも続いており、近年の例としては生成モデルで刻み幅をナノ秒級まで引き上げる手法が提案されている（検証は小分子とペプチドで、結晶での実証ではない。“Transferable Generative Models Bridge Femtosecond to Nanosecond Time-Step Molecular Dynamics,” arXiv:2510.07589v1, 2025）。 https://arxiv.org/abs/2510.07589 ↩
J. F. Lutsko, D. Wolf, S. R. Phillpot & S. Yip, “Molecular-dynamics study of lattice-defect-nucleated melting in metals using an embedded-atom-method potential,” Physical Review B 40, 2841–2855 (1989)。完全な無欠陥結晶（約1000原子）は融点を200K以上超えても融けずに残る（過熱される）一方、格子欠陥を導入すると融点をわずかでも超えた温度から融解が核生成し、これが実際の融解を支配する機構だと報告する。ただし調べられたのは銅の(001)面上の粒界・自由表面・空隙の平面配列という広がりを持つ欠陥であって、本稿が仕込む孤立した空孔ではない——本稿の設定は、この知見に着想を得た代役である。 https://journals.aps.org/prb/abstract/10.1103/PhysRevB.40.2841 ↩

この記事はAIが執筆しています。内容には誤りが含まれる可能性があります。ご注意ください。