MDシミュレーション実行中に、以下のエラーメッセージが出て計算が止まってしまうことがあります。
このエラーは、「計算が数値的に破綻（発散）し、原子がシミュレーションボックスの彼方へ吹き飛んでしまった」状態になったことを指摘しています。本記事では、この現象が起きるメカニズムと、具体的な対処方法をいくつか紹介します。
■原因
多くの場合、原子同士が異常に近づきすぎたことが原因です。原子同士が近づきすぎると、反発力により、無限大に近い巨大な力が働きます。 この巨大な力を受けた次の瞬間のステップで、原子が異常な速度で弾き飛ばされ、座標が計算不能な値になってしまいます。
■対応策
構造最適化をする
MDシミュレーションを開始する際の初期構造が十分に緩和できておらず、原子に大きな力が働いている可能性があります。 構造最適化のステップ数を増やすか、収束判定条件（fmaxなど）をより厳しく設定し、十分に緩和した状態になるまで計算を行ってください。
timestepを小さくする
十分に構造を緩和しているにも関わらずエラーが出る場合、MDのタイムステップが大きすぎることが疑われます。timestepが大きいと、数値積分の精度が粗くなります。本来ならポテンシャルの壁（斥力）を感じて手前で止まるべき原子が、1ステップの移動距離が長すぎるために、相手の原子に近づき過ぎてしまうことがあります。すると次のステップで異常な反発力を受けて発散します。
軽い元素（H, Liなど）を含む系や高温でのシミュレーションでは原子の速度が大きくなりやすいため、特に注意が必要です。一般的には、timestep = 1.0 fsを用いますが、0.5 fs や 0.25 fsと小さくして試してみてください。
以下の図は、timestepを変えて1次元の2原子の運動をシミュレーションした際の軌跡です。 横軸は原子間距離、縦軸はポテンシャルエネルギーを表しています。timestepが小さい場合は安定した軌道を描いていますが、timestepが大きい場合は原子間距離が不自然に急増し、発散してしまっている様子が見て取れます。
本記事の図は、以下のコードを用いて作成しました。初期条件やtimestepの値を変更して、どのような条件下で発散が起きるのか、ぜひ実際に実験してみてください
import numpy as np
import matplotlib.pyplot as plt
class LennardJonesSimulator:
"""
A class to handle Lennard-Jones potential calculations.
"""
def __init__(self, epsilon=1.0, sigma=1.0):
"""
Initialize the simulator with potential parameters.
Args:
epsilon (float): Depth of the potential well.
sigma (float): Finite distance at which the potential is zero.
"""
self.epsilon = epsilon
self.sigma = sigma
def potential(self, r):
"""Calculates the Lennard-Jones potential V(r)."""
return 4 * self.epsilon * ((self.sigma / r)**12 - (self.sigma / r)**6)
def force(self, r):
"""Calculates the Lennard-Jones force F = -dV/dr."""
# Calculate force using class parameters
return 24 * self.epsilon * (2 * (self.sigma / r)**13 - (self.sigma / r)**7) / self.sigma
def simulate_verlet(dt, steps, potential, r=2.5, v=-5.0, m=1.0):
"""
Performs a 1D simulation using the Velocity Verlet integration method.
Args:
dt (float): Time step size.
steps (int): Maximum number of steps.
potential: Instance containing the force method.
r (float): Initial position.
v (float): Initial velocity.
m (float): Mass of the particle.
Returns:
list: Trajectory of the particle's position.
"""
# Initial conditions: Particle moving from right to left towards the wall
r_traj = [r]
for _ in range(steps):
# 1. Update velocity by half a step
f = potential.force(r)
v_half = v + 0.5 * f * dt / m
# 2. Update position by a full step
r_new = r + v_half * dt
# 3. Calculate force at the new position
f_new = potential.force(r_new)
# 4. Update velocity by the remaining half step
v_new = v_half + 0.5 * f_new * dt / m
# Update state variables
r = r_new
v = v_new
r_traj.append(r)
# Terminate if the particle reflects and moves far enough away
if r > 3.0:
break
return np.array(r_traj)
small_timestep = 0.01
large_timestep = 0.1
potential = LennardJonesSimulator()
r_success = simulate_verlet(dt=small_timestep, steps=50, potential=potential)
r_fail = simulate_verlet(dt=large_timestep, steps=50, potential=potential)
plt.figure(figsize=(10, 6))
r_grid = np.linspace(0.8, 10.0, 300)
plt.plot(r_grid, potential.potential(r_grid), 'k-', alpha=0.3, label='LJ Potential')
plt.plot(r_success, potential.potential(r_success), 'b-', label='Small timestep')
plt.plot(r_fail, potential.potential(r_fail), 'ro--', label='Large timestep')
plt.ylim(-1.5, 15)
plt.xlim(0.5, 10.0)
plt.axhline(0, color='gray', lw=0.5)
plt.xlabel('Distance r')
plt.ylabel('Potential Energy V(r)')
plt.title('Effect of Large Timestep')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()