python实现四个氢原子的一维氢原子链的分子动力学模拟

好的，下面是一个简单的 Python 代码实现四个氢原子的一维氢原子链的分子动力学模拟：

import numpy as np

# 定义一些常数
m = 1.0  # 氢原子质量
k = 1.0  # 弹簧系数
T = 300.0  # 温度
dt = 0.001  # 时间步长
N = 4  # 原子数
L = 1.0  # 盒子长度

# 初始化位置和速度
x = np.linspace(0, L, N, endpoint=False)
v = np.random.randn(N) * np.sqrt(T / m)

# 初始化加速度
a = np.zeros(N)

# 定义一个计算势能的函数
def potential_energy(x):
    U = 0.0
    for i in range(N - 1):
        U += 0.5 * k * (x[i+1] - x[i])**2
    return U

# 计算初始势能和动能
U = potential_energy(x)
K = 0.5 * m * np.sum(v**2)

# 定义一个更新速度和位置的函数
def update(x, v, a, dt):
    x += v * dt + 0.5 * a * dt**2
    a_new = np.zeros(N)
    a_new[1:N-1] = k * (x[0:N-2] - 2*x[1:N-1] + x[2:N]) / m
    a = a_new
    v += 0.5 * (a + a_new) * dt
    return x, v, a

# 开始模拟
for i in range(100000):
    # 更新位置和速度
    x, v, a = update(x, v, a, dt)
    
    # 计算势能和动能
    U_new = potential_energy(x)
    K_new = 0.5 * m * np.sum(v**2)
    
    # 计算总能量
    E = U_new + K_new
    
    # 打印能量和位置信息
    if i % 100 == 0:
        print("Step {}: Energy = {:.4f}, Position = {}".format(i, E, x))

这个代码使用 Verlet 算法来更新位置和速度，同时计算势能和动能。它模拟了 100000 步，每 100 步打印一次能量和位置信息。你可以根据需要调整模拟参数和步数。