本内容为西安交通大学 MATL501202 Introduction of computational materials science: PartI 课程课程作业内容，仅供参考

简介

Velocity Verlet算法是在分子动力学中使用较为普遍的一种算法，其特点是通过半时刻变量来抵抗系统扰动，提升模拟精度。同时，模拟过程符合能量守恒。Velocity Verlet算法的基本原理如下图：

可逆性证明

要证明可逆性，即证明该算法可以在已知$t+\delta t$时刻状态时叠加$-\delta t$进行反推到$t$时刻状态。

位置更新算法可逆性证明

$$ \begin{gathered} r(t)=r(t+\delta t)+v(t+\delta t)(-\delta t)+\frac{1}{2} a(t+\delta t)(-\delta t)^2 \\ \Rightarrow r(t)=r(t+\delta t)-v(t+\delta t) \delta t+\frac{1}{2} a(t+\delta t)(\delta t)^2 \\ \Rightarrow r(t)=r(t+\delta t)-\delta t\left[v(t+\delta t) \delta t-\frac{1}{2} a(t+\delta t) \cdot \frac{1}{2} \delta t\right] \\ \Rightarrow r(t)=r(t+\delta t)-\delta t \cdot v\left(t+\frac{1}{2} \delta t\right) \end{gathered} $$

速度更新算法可逆性证明

在上一部分已经将位置矢量进行了逆向推演，此时即可知道$t$时刻的$F$和$a$，由此：

$$ v(t)=v(t+\delta t)-\frac{1}{2} a(t+\delta t) \delta t-\frac{1}{2} a(t) \delta t $$

$$ \Rightarrow v(t)=v(t+\delta t)-\frac{1}{2}[a(t)+a(t+\delta t)] \delta t $$

由此，对比正逆向算法：

代码实现

构建力场

此处，我们举例$CO_2$分子的 Velocity-Verlet 算法模拟，使用 Lenard-Jones Pontential 构建力场：

$$ \begin{aligned} &\text { Calculation of L-J Potential: } V(r)=4 \epsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^6\right]\\ &\Rightarrow F(r)=-\frac{d V(r)}{d r}=24 \epsilon\left[2\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^6\right] \frac{1}{r} \end{aligned} $$

# Lennard-Jones 势能计算
def lennard_jones_force(position1, position2, epsilon=1.0, sigma=1.0):
    r_vec = position1 - position2
    r = np.linalg.norm(r_vec)
    if r == 0:
        return np.zeros_like(r_vec)
    force_magnitude = 48 * epsilon * ((sigma**12 / r**13) - 0.5 * (sigma**6 / r**7))
    force = force_magnitude * r_vec / r
    return force

Velocity-Verlet核心代码

    # 更新位置和速度（Velocity-Verlet）
    positions += velocities * dt + 0.5 * forces / mass * dt**2
    new_forces = np.zeros_like(forces)
    for i in range(num_particles):
        for j in range(i + 1, num_particles):
            force = lennard_jones_force(positions[i], positions[j])
            new_forces[i] += force
            new_forces[j] -= force
    velocities += 0.5 * (forces + new_forces) / mass * dt

完整代码

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from mpl_toolkits.mplot3d import Axes3D

# Lennard-Jones 势能计算
def lennard_jones_force(position1, position2, epsilon=1.0, sigma=1.0):
    r_vec = position1 - position2
    r = np.linalg.norm(r_vec)
    if r == 0:
        return np.zeros_like(r_vec)
    force_magnitude = 48 * epsilon * ((sigma**12 / r**13) - 0.5 * (sigma**6 / r**7))
    force = force_magnitude * r_vec / r
    return force

# 初始化条件
num_particles = 10  # 粒子数量
box_size = 10.0  # 模拟盒子的大小
mass = 44.01  # CO₂ 分子的质量近似
dt = 0.01  # 时间步长
num_steps = 1000  # 模拟步数

# 初始化粒子位置和速度
positions = np.random.rand(num_particles, 3) * box_size
velocities = (np.random.rand(num_particles, 3) - 0.5) * 2

# 设置图形
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.set_xlim(0, box_size)
ax.set_ylim(0, box_size)
ax.set_zlim(0, box_size)
scat = ax.scatter(positions[:, 0], positions[:, 1], positions[:, 2])

# 更新函数
def update(frame):
    global positions, velocities
    
    # 计算所有粒子的力
    forces = np.zeros_like(positions)
    for i in range(num_particles):
        for j in range(i + 1, num_particles):
            force = lennard_jones_force(positions[i], positions[j])
            forces[i] += force
            forces[j] -= force

    # 更新位置和速度（Velocity-Verlet）
    positions += velocities * dt + 0.5 * forces / mass * dt**2
    new_forces = np.zeros_like(forces)
    for i in range(num_particles):
        for j in range(i + 1, num_particles):
            force = lennard_jones_force(positions[i], positions[j])
            new_forces[i] += force
            new_forces[j] -= force
    velocities += 0.5 * (forces + new_forces) / mass * dt

    # 周期性边界条件
    positions = positions % box_size

    # 更新散点图
    scat._offsets3d = (positions[:, 0], positions[:, 1], positions[:, 2])

# 创建动画
ani = FuncAnimation(fig, update, frames=range(num_steps), interval=50)
ax.set_xlabel('X Position')
ax.set_ylabel('Y Position')
ax.set_zlabel('Z Position')
ax.set_title('3D Real-time Simulation of CO₂ Molecule Motion')
plt.show()