Python双原子合金化MD模拟代码：原子无聚集问题修正求助

核心问题分析 #

代码无法产生预期团簇结构的关键原因在于物理机制的多处错误：

• 温度单位混淆（使用摄氏温度而非开尔文）
• LJ力计算逻辑错误（直接用势能值而非势能梯度）
• 混合原子对的相互作用参数缺失
• 时间步长单位完全不符合MD模拟规范
• 周期性边界条件未正确应用于距离计算
• 未设置LJ相互作用截断距离
• 温度控制方式错误（直接累加温度值而非速度缩放）

修正后的完整代码 #

import numpy as np
import matplotlib.pyplot as plt

# 物理常数
Avogadro = 6.02214086e23
Boltzmann = 1.38064852e-23  # J/K
fs_to_s = 1e-15

# Pd/Au原子参数（kg）
pd_mass = 106.42 / Avogadro
au_mass = 196.97 / Avogadro

# 原子数量
n_pd = 50
n_au = 50
natoms = n_pd + n_au

# 模拟盒尺寸（m）- 初始值对应约1nm³的盒子，更符合MD尺度
box_size = np.array([3e-9, 3e-9, 3e-9])

# 生成初始位置：区分Pd和Au原子
def generate_initial_positions(n_pd, n_au, box_size):
    pd_pos = np.random.rand(n_pd, 3) * box_size
    au_pos = np.random.rand(n_au, 3) * box_size
    return np.vstack([pd_pos, au_pos])

# 生成初始速度：按原子质量分别遵循Maxwell-Boltzmann分布
def generate_initial_velocities(n_pd, n_au, temperature_K):
    # Pd原子速度
    pd_speed_std = np.sqrt(Boltzmann * temperature_K / pd_mass)
    pd_vel = np.random.normal(0, pd_speed_std, (n_pd, 3))
    # Au原子速度
    au_speed_std = np.sqrt(Boltzmann * temperature_K / au_mass)
    au_vel = np.random.normal(0, au_speed_std, (n_au, 3))
    # 移除整体动量
    all_vel = np.vstack([pd_vel, au_vel])
    total_momentum = np.sum(all_vel * np.tile([pd_mass]*n_pd + [au_mass]*n_au, (3,1)).T, axis=0)
    all_vel -= total_momentum / (n_pd*pd_mass + n_au*au_mass)
    return all_vel

# LJ势能参数（sigma: m, epsilon: J）
lj_params = {
    'Pd-Pd': {'sigma': 2.75e-10, 'epsilon': 1.02e-20},
    'Au-Au': {'sigma': 2.88e-10, 'epsilon': 1.0e-20}
}
# 洛伦兹-贝塞罗特混合规则计算Pd-Au参数
lj_params['Pd-Au'] = {
    'sigma': (lj_params['Pd-Pd']['sigma'] + lj_params['Au-Au']['sigma']) / 2,
    'epsilon': np.sqrt(lj_params['Pd-Pd']['epsilon'] * lj_params['Au-Au']['epsilon'])
}

# 计算原子对的LJ参数
def get_lj_params(atom_i, atom_j):
    if atom_i < n_pd and atom_j < n_pd:
        return lj_params['Pd-Pd']
    elif atom_i >= n_pd and atom_j >= n_pd:
        return lj_params['Au-Au']
    else:
        return lj_params['Pd-Au']

# 周期性边界条件下的最短距离计算
def pbc_distance(pos_i, pos_j, box_size):
    delta = pos_i - pos_j
    delta -= box_size * np.round(delta / box_size)
    return np.linalg.norm(delta), delta

# LJ力计算（含截断与移位）
def compute_forces(positions, box_size, cutoff):
    forces = np.zeros_like(positions)
    for i in range(natoms):
        for j in range(i+1, natoms):
            r, delta = pbc_distance(positions[i], positions[j], box_size)
            if r < cutoff:
                params = get_lj_params(i, j)
                sigma = params['sigma']
                epsilon = params['epsilon']
                # LJ势能导数（dU/dr）
                lj_deriv = 4 * epsilon * (-12 * (sigma**12)/(r**13) + 6 * (sigma**6)/(r**7))
                # 力是负梯度：F = -dU/dr * (delta/r)
                force = -lj_deriv * (delta / r)
                forces[i] += force
                forces[j] -= force
    return forces

# Verlet积分
def integrate(positions, velocities, forces, masses, dt):
    # 更新位置
    positions += velocities * dt + 0.5 * forces / masses[:, np.newaxis] * dt**2
    # 应用周期性边界条件
    positions %= box_size
    # 计算新力
    new_forces = compute_forces(positions, box_size, lj_cutoff)
    # 更新速度
    velocities += 0.5 * (forces + new_forces) / masses[:, np.newaxis] * dt
    return positions, velocities, new_forces

# 速度缩放控制温度
def scale_velocities(velocities, masses, current_temp, target_temp):
    # 计算当前温度
    kinetic_energy = 0.5 * np.sum(masses[:, np.newaxis] * velocities**2)
    current_temp = (2 * kinetic_energy) / (3 * natoms * Boltzmann)
    # 缩放因子
    scale_factor = np.sqrt(target_temp / current_temp)
    velocities *= scale_factor
    return velocities

# 初始化模拟
temp_start_C = 100
temp_start_K = temp_start_C + 273.15
temp_end_C = 800
temp_end_K = temp_end_C + 273.15
heating_rate_K_per_ns = 100 / 1e9  # 100K/ns（合理的MD升温速率）
dt_fs = 2  # 2fs时间步长
dt = dt_fs * fs_to_s
lj_cutoff = 1.5e-9  # LJ截断距离（~5倍sigma）

# 初始化位置、速度、质量数组
positions = generate_initial_positions(n_pd, n_au, box_size)
velocities = generate_initial_velocities(n_pd, n_au, temp_start_K)
masses = np.array([pd_mass]*n_pd + [au_mass]*n_au)

# 模拟参数
steps_per_snapshot = 5000  # 每5000步输出一次快照（10ps）
total_steps = int((temp_end_K - temp_start_K) / (heating_rate_K_per_ns * dt * 1e9))

current_temp_K = temp_start_K
for step in range(total_steps):
    # 计算目标温度
    target_temp_K = temp_start_K + heating_rate_K_per_ns * step * dt * 1e9
    
    # 计算力
    forces = compute_forces(positions, box_size, lj_cutoff)
    # 积分
    positions, velocities, forces = integrate(positions, velocities, forces, masses, dt)
    # 速度缩放控温
    velocities = scale_velocities(velocities, masses, current_temp_K, target_temp_K)
    
    # 输出快照
    if step % steps_per_snapshot == 0:
        current_temp_C = target_temp_K - 273.15
        fig = plt.figure(figsize=(8,6))
        ax = fig.add_subplot(111, projection='3d')
        # 绘制Pd原子
        ax.scatter(positions[:n_pd,0], positions[:n_pd,1], positions[:n_pd,2], 
                   c='r', label='Pd', s=50)
        # 绘制Au原子
        ax.scatter(positions[n_pd:,0], positions[n_pd:,1], positions[n_pd:,2], 
                   c='b', label='Au', s=50)
        ax.set_xlim(0, box_size[0])
        ax.set_ylim(0, box_size[1])
        ax.set_zlim(0, box_size[2])
        ax.set_xlabel('X (m)')
        ax.set_ylabel('Y (m)')
        ax.set_zlabel('Z (m)')
        ax.set_title(f'Pd-Au合金化模拟 @ {current_temp_C:.1f}°C')
        ax.legend()
        plt.show()

    current_temp_K = target_temp_K

关键修正点详解 #

1. 温度单位修正：所有热力学计算统一使用开尔文温度，避免能量计算错误。
2. LJ力计算修复：基于势能梯度计算力，正确实现原子间的吸引/排斥作用。
3. 混合原子对参数：采用洛伦兹-贝塞罗特规则计算Pd-Au的LJ参数，符合真实原子间相互作用规律。
4. 时间步长调整：使用飞秒级时间步长，符合MD模拟的时间尺度要求。
5. 周期性边界条件：计算原子间距离时考虑镜像原子，确保跨边界原子对的相互作用正确。
6. LJ截断与移位：设置合理的截断距离，提升计算效率并避免势能不连续。
7. 温度控制：采用速度缩放法维持目标温度，确保系统按设定速率升温。
8. 初始条件优化：区分不同原子的速度分布，移除整体动量，保证模拟初始状态符合物理规律。

内容的提问来源于stack exchange，提问作者Leili_D