In Silico

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

分子動力学とは？原子の運動を解く仕組みとシミュレーションの限界

2026/6/26 （更新: 2026/10/11）

🤖 AIを使用して執筆

  • 分子動力学
  • 計算物質科学
  • 材料モデリング
  • レナード–ジョーンズ
  • Velocity-Verlet
  • 速度スケーリング
  • LAMMPS
  • 結晶
  • 融点
  • 空孔
  • 周期境界
  • 多体効果
  • 密度汎関数理論
  • ニューラルネットワーク
  • 六方配向秩序
  • 平均二乗変位
  • フェムト秒
  • Anton
  • ニュートンの法則
  • 換算単位
  • 三角格子
  • ファンデルワールス力
  • 粒界
  • 自由表面
  • 過熱
目次
背景・問い・要点
背景

温度とは、原子の動きの激しさである。熱いと原子は激しく揺れ、ある限界を超えると並びが崩れて融ける。その様子は、本物を顕微鏡でのぞかなくても、原子に運動の法則を当てて時間を一歩ずつ進めれば、コンピュータの中で作り直せる。法則から動きを計算で作り直すこの方法を「分子動力学（MD）」と呼ぶ。

計算で作り直せるなら、実験室で本物を合成する前に、温めたらどう融けるか、力を加えたらどう変形するかを、計算機の中で覗ける。ただし計算が返すのは、原子どうしの力をどう置いたかに従った動きであり、その置き方が現実とずれれば、計算の結果も現実とずれる。本稿は、182個の原子を並べた2次元の結晶を実際に計算し、温めて秩序を崩す例を取る。

問い

原子に運動の法則を当てて時間を進めるだけの計算で、結晶の秩序が崩れる様子をどこまで再現でき、その結果は何に左右されるのか。

要点

原子ごとに力を求めて少しだけ動かすことを繰り返すと、温度とともに固体の秩序が崩れていく過程を、計算の中で再現できる。 短いランで見えるのは、ある温度で一気に融ける瞬間ではなく、なめらかに秩序が失われる過程で、融点そのものは別の計算が要る。計算は、原子どうしの力のモデルに入れた物理以上のことを返さない。単純な2体の力のモデルは金属の融点や熱膨張を大きく外す。時間を数フェムト秒の刻みで進めるので、ゆっくり進む現象は直接は追いにくい。

モデル・例示

力を求め、少し動かし、繰り返す

※ この節の数値は説明のための仮定で、測定値ではありません。

手順は3行で書ける。

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

ここでは、物理法則を小刻みに積分していく。ただし温度を上げる操作だけは法則の外から加える。具体的には、毎ステップ、全原子の速度を狙いの温度にそろえ直す。

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

  • 近づきすぎると、急激に強く反発する（電子の雲が重なるのを嫌う。距離の −12 乗、と置くのが定番）。
  • 少し離れると、弱く引き合う（ファンデルワールス力。距離の −6 乗）。

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

刻みの細かさが、追える時間の長さを決めることも数で確かめられる。1歩を1フェムト秒とすると、1ピコ秒（1兆分の1秒）を進めるには1000歩、1ナノ秒には100万歩、1マイクロ秒には10億歩が要る。1歩ごとに手順1の力の計算をやり直すので、進められる時間は、力を何回計算できるかで頭打ちになる。

実際に走らせて結晶を融かす

下のシミュレーションでは、原子どうしの力に、近いと強く反発し少し離れると弱く引き合うだけの単純な式（レナード–ジョーンズ・ポテンシャル）を使い、原子はこの力だけで動く。原子ごとに力を求め、ごく短い時間だけ動かすことを繰り返し、温度だけは全原子の速度をそろえ直して外から上げる。

仕組みを数式で書く（クリックで展開）

このシミュレーションの計算を式で書くとこうなる（単位はレナード–ジョーンズの換算単位＝ε=σ=1\varepsilon=\sigma=1。だから温度も 0.10.1〜1.61.6 という無次元の値で書ける）。

  1. 力を計算する（レナード–ジョーンズ）。 距離 rr だけ離れた2原子のあいだのエネルギー V(r)V(r) は、反発の 1212 乗と引力の 66 乗の差で書ける。原子にはたらく力は、そのエネルギーの傾き − dV/dr-\,dV/dr である。
V(r)=4ε[(σr)12−(σr)6]V(r) = 4\varepsilon\left[\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6}\right] F(r)=−dVdr=24εr[ 2(σr)12−(σr)6 ]F(r) = -\frac{dV}{dr} = \frac{24\varepsilon}{r}\left[\,2\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6}\,\right]

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

  1. 少しだけ動かす（Velocity–Verlet 積分）。 力から加速度 a=F/ma = F/m を出し、位置と速度を微小時間 dtdt だけ進める。この積分法は、位置を先に更新し、新しい位置で力を測り直してから速度を仕上げる。定番の積分法で、「Velocity-Verlet」と呼ばれる1。位置と速度が同じ時刻で揃うので扱いやすく、丸め誤差にも強い形として知られる1。
x(t+dt)=x(t)+v(t) dt+12 a(t) dt2x(t+dt) = x(t) + v(t)\,dt + \tfrac{1}{2}\,a(t)\,dt^{2} v(t+dt)=v(t)+12 [ a(t)+a(t+dt) ] dtv(t+dt) = v(t) + \tfrac{1}{2}\,\bigl[\,a(t) + a(t+dt)\,\bigr]\,dt
  1. 温度を上げる（速度スケーリング）。 温度とは、原子の運動エネルギーの平均そのものだ（NN 原子・2次元、kB=1k_B=1 の換算単位）。狙った温度 TtargetT_\text{target} に寄せるには、全原子の速度を一律に Ttarget/T\sqrt{T_\text{target}/T} 倍すればよい。この TtargetT_\text{target} を 0.10.1 から 1.61.6 へ少しずつ上げていくと、ある所で格子がほどけて融ける。
T=12N∑i=1Nmi ∣vi∣2vi←vi TtargetTT = \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.6T=1.6）でも、最初の隣どうしの約7割はまだ保たれている。もっとも、この保持率は原子が元の隣から離れるのに時間がかかるぶん遅れて下がる量なので、短いランの終端に残る約7割だけでは、固体のままか液体に達したかは決められない。そこで結晶の向きのそろい具合を別に測ると、完全な結晶の1から終端では約0.14まで落ちており、同じ温度で計算を続けると原子は元の位置から離れ続け、そろい具合もほぼ0になる（液体に見られるふるまいだが、2次元でのこの量は融解の目安であって相の判定ではない）。

「温度とは原子の動きの激しさだ」という言葉が、ここでは計算の中で実際に起きている。融点そのものを求めるなら、別の計算が要る。銅の融解を分子動力学で調べた Lutsko らは、完全結晶と液相の自由エネルギーを解析して融点を独立に決めている2。

何の役に立つのか

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

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

ここで一つ注意がある。

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

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

このGIFは、実際に走らせた2次元分子動力学（182原子、三角格子のレナード–ジョーンズ結晶、周期境界、Velocity-Verlet積分、速度スケーリング法で目標温度を0.1→1.6へ上昇）の出力。秩序が崩れる起点となるよう、格子点の約7%（196点のうち14点）を抜いた空孔（欠陥）を仕込んである。この空孔濃度は現実の結晶の平衡空孔濃度よりはるかに高い。これは、自由表面を持たない周期境界の小さな系を、短いランのあいだに崩すための人為的な設定である。銅を対象とした分子動力学の計算では、完全な無欠陥結晶は過熱してもなかなか融けない一方、自由表面・粒界・ある大きさ以上の空隙といった欠陥を導入すると、融点を超えた温度でそこから融解が核生成すると報告されている2。本稿は周期境界で自由表面を持たないため、その代役として高濃度の空孔を置いた（同論文では孤立した空孔は融解を起こしておらず、同じ機構ではない）。このシミュレーションは決定論的（固定シード）で、コードとパラメータは管理下にあり再現可能である。

各時刻の格子の乱れ（初期位置からのずれと近傍の保持率）を測ると、温度の上昇とともに単調に秩序が失われる（終端 T=1.6T=1.6 で最近接ペアの保持率は約0.68）。急峻な転移点は、この短いランでは観測されない。結晶の向きのそろい具合は、各原子から最近接の原子へ向かう向きを6回対称で平均した量（六方配向秩序 ψ6）の系全体の大きさで測った（完全な三角格子で1、向きがばらばらなら0に近い）。こちらは昇温中に単調ではなく、目標温度1.06から1.30のあいだで0.55から0.22へと他の区間より急に下がるが、1シード・182原子の揺らぎを超えるかは判定していない。ランの終端から目標温度1.6のまま計算を続けると、換算時間2.5・5・10で平均二乗変位は0.97・1.42・4.22と増え続け、ψ6 は0.09・0.07・0.015へ下がる。

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

出典6件

  1. Swope ほか「A computer simulation method for the calculation of equilibrium constants…」J. Chem. Phys. 76, 1982. https://pubs.aip.org/aip/jcp/article/76/1/637/397800 https://archive.org/details/DTIC_ADA1030958 — 付録で速度形の Verlet 法を示し、丸め誤差への強さと、位置と速度が同時刻で揃う点を利点とする。 ↩ ↩2

  2. Lutsko ほか「Molecular-dynamics study of lattice-defect-nucleated melting…」Phys. Rev. B 40, 1989. https://journals.aps.org/prb/abstract/10.1103/PhysRevB.40.2841 — 無欠陥結晶は融点を200K超でも融けず、自由表面・粒界・13原子の空隙は融点直上で融解を核生成した。単一空孔は1400Kでも融けない。 ↩ ↩2

  3. Thompson ほか「LAMMPS - a flexible simulation tool…」Comput. Phys. Commun. 271, 2022. https://doi.org/10.1016/j.cpc.2021.108171 https://www.lammps.org/ — 概説論文。公式サイトは材料モデリングに重心を置く古典 MD コードで、金属・半導体向けのポテンシャルを持つと書く。 ↩

  4. Daw & Baskes「Embedded-atom method: Derivation and application…」Phys. Rev. B 29, 1984. https://doi.org/10.1103/PhysRevB.29.6443 https://doi.org/10.1098/rspa.2003.1190 — 多体効果を入れる EAM。物性値は Srinivasan & Baskes（Proc. R. Soc. A 460, 2004）表3。 ↩

  5. Behler & Parrinello「Generalized Neural-Network Representation of High-Dimensional Potential-Energy Surfaces」PRL 98, 2007. https://doi.org/10.1103/PhysRevLett.98.146401 — DFT のエネルギー面をネットワークで表し、力を DFT より数桁速く与える。 ↩

  6. Dror ほか「Exploring atomic resolution physiology…」J. Gen. Physiol. 135, 2010. https://doi.org/10.1085/jgp.200910373 — 1歩は数フェムト秒で1マイクロ秒に約10億歩、並列化と Anton でミリ秒級と報告。刻みを広げる生成モデル（Diez ほか 2026）の適用先は小分子とペプチド。 https://doi.org/10.1126/sciadv.aed2333 https://arxiv.org/abs/2510.07589 ↩ ↩2

この記事はAIを使用して執筆されました。内容には誤りが含まれる可能性があります。事実や数字の出典は脚注から確かめられます。