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

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

2026/6/26 （更新: 2026/9/28）シリーズ「シミュレーションで材料を試運転する」第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個の原子を並べ、少しずつ温度を上げていく）。色は速さを表す。青い原子はゆっくり（冷たい）、赤い原子は速い（熱い）。

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

最初、原子はきれいな格子に並び、青いまま、その場で小さく震えている。この状態が固体である。温度を上げていくと揺れが大きくなり、色が暖色に変わり、格子の像がしだいにぼやけていく。ただしこの短いランで見えるのは、ある温度を境に一気に融ける「瞬間」ではない。温度とともになめらかに秩序が失われていく過程である。実際、初期位置からのずれと近傍の保持率を測ると、秩序は単調に崩れていくだけで、飛びはない。ランの終端（ $T=1.6$ ）でも、最初の隣どうしの約7割はまだ保たれている。もっとも、この保持率は原子が元の隣から離れるのに時間がかかるぶん遅れて下がる量なので、短いランの終端に残る約7割だけでは、固体のままか液体に達したかは決められない。「温度とは原子の動きの激しさだ」という言葉が、ここでは計算の中で実際に起きている。融点そのものを求めるなら、別の計算が要る。² の著者たちは、完全結晶と液相の自由エネルギーを解析して融点を独立に決めている。

何の役に立つのか

これは、材料を作る前に“試運転”するということだ。実験室で本物を合成しなくても、計算機の中で「温めたらどの温度で融けるか」「力を加えたらどう変形・破壊するか」「別の物質がどう染み込んでいくか（拡散）」を覗ける。分子動力学は、新しい材料の候補を作る前にふるいにかける工程のうち、精密な段にあたる。こうした計算を走らせる公開コードの一つに LAMMPS があり、材料モデリングに重心を置いた古典分子動力学コードで、金属や半導体といった固体材料向けのポテンシャルを備える³。

シミュレーションは「入れた物理」しか返さない

ここで一つ注意がある。

このシミュレーションの心臓は、「原子どうしがどう力を及ぼし合うか」というモデルだ。今回は単純な模型（レナード–ジョーンズ）を使った。もしこのモデルが現実とずれていれば、どんなに精密に計算しても、出てくる“映画”は現実とずれる。実際、3次元の面心立方格子でこの単純な2体間の力のモデルが与える物性を実験値と並べると、融点は金属の実測範囲の上限の約2倍、熱膨張係数は下限の約8分の1で、体積弾性率も範囲の外に出る。原子が集団で及ぼす多体効果を加えて初めて、多くの面心立方金属で実験値の範囲に入る、と報告されている⁴。

シミュレーションは、入れた物理以上のことは教えない。だから本当の勝負は、「正しい力のモデルを使えているか」を確かめることにある。経験的な式の代わりに、量子力学の計算（密度汎関数理論）が与えるエネルギーと力をニューラルネットワークで表し、その計算より数桁速く求める方法も提案されている⁵。もう一つの現実的な制約は、時間の刻み幅そのものだ。原子の振動を追うには数フェムト秒の刻みが要り、1マイクロ秒を計算するだけでも10億回近い逐次ステップが要る⁶。原子スケールの計算が、ゆっくり進む現象を直接は追いにくい理由の一つだ。ただしこれは越えられない壁ではない。並列化の進展で1マイクロ秒を超えるランが汎用クラスタでも実用的になり、専用計算機 Anton ではミリ秒級に届いている⁶。

この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へ上げている。箱の大きさは固定（体積一定）なので、加熱しても格子は膨張しない。ただし空孔を抜いた分だけ系は完全結晶より約7%疎になっており、この密度低下も、空孔がもたらす局所的な乱れと同程度に秩序の喪失を早めている（同じコードで対照を取ると、空孔を入れずに箱だけ同じ密度まで広げた場合も、空孔を入れて密度を完全結晶に揃えた場合も、終端の保持率は約0.95から0.8前後へ下がる。数シードの目安）。各温度で平衡化させていないので、ここから融点そのものを読み取ることはできない。分子動力学・融解は計算物質科学の標準的な手法に基づく。

出典6件

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)。本稿が使う積分法を、この論文が付録 “Appendix: Velocity Form of the Verlet Algorithm” で書き下している。後年これが Velocity-Verlet と呼ばれるようになった、広く引用される出典。原典が速度形式の利点として挙げるのは、丸め誤差に強い総和形式の精度を保つことと、位置と速度が同時刻で揃うため確率的衝突を差し込みやすいことの2つである。出版社版は購読制だが、同一原稿の ONR 技術報告版が公開されている（ https://archive.org/details/DTIC_ADA1030958 ）。 https://pubs.aip.org/aip/jcp/article/76/1/637/397800 ↩ ↩²
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原子）は、この模型ポテンシャルについて別途決めた融点（1171±30K）を200K以上超えても融けずに残り（過熱）、約1450Kで力学的不安定により自発的に融ける。一方、銅の(001)面上の粒界・自由表面・空隙の平面配列を導入すると融点をわずかでも超えた温度から融解が核生成し、これが実際の融解を支配する機構だと報告する。ただし同論文は単一空孔と5原子空隙（濃度約0.5%）も走らせており、5原子空隙は600Kでの平衡化中に崩れ始めて昇温とともに動きやすい単一空孔5個へ解離したが、どちらも融点を約230K上回る1400Kでも融解を起こさなかった。空隙のうち融解を起こしたのは13原子のものだけで、著者らは、効くのは欠陥の広がりではなく局所的な乱れの量だと結論している。本稿の孤立した空孔（約7%）はこの試験範囲の外にあり、原典の実測ではなくその原理に着想を得た代役である。 https://journals.aps.org/prb/abstract/10.1103/PhysRevB.40.2841 ↩ ↩²
A. P. Thompson, H. M. Aktulga, R. Berger ほか, “LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales,” Computer Physics Communications 271, 108171 (2022)。LAMMPS の概説論文。公式サイトは LAMMPS を材料モデリングに重心を置いた古典分子動力学コードと説明し、金属・半導体などの固体材料向けポテンシャルを持ち、GPLv2 のオープンソースとして配布されると書く。 https://doi.org/10.1016/j.cpc.2021.108171 　https://www.lammps.org/ ↩
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）。本文の数値は同論文の表3（3次元の面心立方格子、純粋なレナード-ジョーンズ模型の列）による。融点 0.082 に対し面心立方金属は 0.024–0.041、体積膨張係数 0.34 に対し 2.64–8.67、体積弾性率 51.7 に対し 14–28（いずれも換算単位）。この比較は、本稿のシミュレーションが使う力のモデルにも同種の限界がありうることを示す。 https://doi.org/10.1098/rspa.2003.1190 ↩
J. Behler & M. Parrinello, “Generalized Neural-Network Representation of High-Dimensional Potential-Energy Surfaces,” Physical Review Letters 98, 146401 (2007)。密度汎関数理論（DFT）のような計算の重い方法では大きな系を長くシミュレーションできないとして、DFT のポテンシャルエネルギー面をニューラルネットワークで表し、任意の大きさの系で全原子の位置からエネルギーと力を DFT より数桁速く与える方法を提案した論文（要旨による。精度はシリコン結晶で経験的ポテンシャルおよび DFT と比べて示している）。 https://doi.org/10.1103/PhysRevLett.98.146401 ↩
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 　なおこの制約を緩める試みも続いており、近年の例としては生成モデルで刻み幅をナノ秒級まで引き上げる手法が提案されている（検証は小分子とペプチドで、結晶での実証ではない。J. V. Diez, M. Schreiner & S. Olsson, “Transferable generative models bridge femtosecond to nanosecond time-step molecular dynamics,” Science Advances 12, eaed2333 (2026)。プレプリントは arXiv:2510.07589）。 https://doi.org/10.1126/sciadv.aed2333 　https://arxiv.org/abs/2510.07589 ↩ ↩²

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