基于Lennard-Jones+FENE势的双链蒙特卡洛模拟代码修改求助

我帮你修改了代码，现在可以生成两条各20个单体的线性链，满足你提出的所有约束条件。下面是修改后的完整代码和关键改动说明：

关键改动点 #

• 调整粒子总数为40（2×20），为第二条链设置初始间距，避免初始重叠
• 固定两条链的首个单体（不参与MC移动），确保它们锚定在表面
• 更新键能计算逻辑，同时处理两条链的FENE相互作用
• 优化轨迹输出，用不同元素区分两条链，方便可视化
• 重构势能函数，将参数作为函数传入，提升代码可维护性和Numba兼容性

修改后的完整代码 #

import numpy as np
import numba as nb

@nb.jit()
def gen_chain(N_chain, chain_spacing=2.0):
    total_N = 2 * N_chain
    x = np.zeros(total_N)
    y = np.zeros(total_N)
    z = np.zeros(total_N)
    
    # 第一条链：起始于(0,0,0)，沿z轴初始排列
    z[:N_chain] = np.linspace(0, N_chain * 0.9, num=N_chain)
    # 第二条链：起始于(chain_spacing, 0, 0)，沿z轴排列，保持初始间距
    x[N_chain:] = chain_spacing
    z[N_chain:] = np.linspace(0, N_chain * 0.9, num=N_chain)
    
    return np.column_stack((x, y, z))

@nb.jit()
def lj(rij2, sigma, epsilon, rcutoff_sq):
    if rij2 >= rcutoff_sq:
        return 0.0
    sig_by_r6 = np.power(sigma**2 / rij2, 3)
    sig_by_r12 = np.power(sigma**2 / rij2, 6)
    lje = 4 * epsilon * (sig_by_r12 - sig_by_r6)
    return lje

@nb.jit()
def fene(rij2, K, R, r0):
    rij = np.sqrt(rij2)
    return (-0.5 * K * np.power(R, 2) * np.log(1 - ((rij - r0) / R)**2))

@nb.jit()
def total_energy(coord, N_chain, sigma, epsilon, rcutoff_sq, K, R, r0):
    total_N = 2 * N_chain
    e_nb = 0.0
    # 非键相互作用：覆盖所有粒子对（链内+链间）
    for i in range(total_N):
        for j in range(i):
            rij = coord[i] - coord[j]
            rij2 = np.dot(rij, rij)
            e_nb += lj(rij2, sigma, epsilon, rcutoff_sq)
    
    e_bond = 0.0
    # 计算第一条链的FENE键能（0-1, 1-2, ..., 18-19）
    for i in range(1, N_chain):
        rij = coord[i] - coord[i-1]
        rij2 = np.dot(rij, rij)
        e_bond += fene(rij2, K, R, r0)
    # 计算第二条链的FENE键能（20-21, ..., 38-39）
    for i in range(N_chain + 1, total_N):
        rij = coord[i] - coord[i-1]
        rij2 = np.dot(rij, rij)
        e_bond += fene(rij2, K, R, r0)
    
    return e_nb + e_bond

@nb.jit()
def move(coord, N_chain, max_delta):
    total_N = 2 * N_chain
    trial = np.ndarray.copy(coord)
    
    # 跳过两条链的首个单体（索引0和N_chain），它们固定在表面
    for i in range(total_N):
        if i == 0 or i == N_chain:
            continue
        
        while True:
            delta = (2 * np.random.rand(3) - 1) * max_delta
            trial[i] += delta
            # 确保z分量始终为正
            if trial[i, 2] > 0.0:
                break
            trial[i] -= delta
    
    return trial

@nb.jit()
def accept(delta_e, T):
    beta = 1.0 / T
    if delta_e < 0.0:
        return True
    random_number = np.random.rand(1)
    p_acc = np.exp(-beta * delta_e)
    return random_number < p_acc

if __name__ == "__main__":
    # FENE势能参数
    K = 40.0
    R = 0.3
    r0 = 0.7
    # L-J势能参数
    sigma = 0.5716
    epsilon = 1.0
    # MC模拟参数
    N_chain = 20  # 每条链的单体数
    chain_spacing = 2.0  # 两条链首个单体的初始间距
    total_N = 2 * N_chain
    rcutoff = 2.5 * sigma
    rcutoff_sq = rcutoff * rcutoff
    max_delta = 0.01
    n_steps = 100000
    T = 10
    
    # 生成初始双链结构
    coord = gen_chain(N_chain, chain_spacing)
    energy_current = total_energy(coord, N_chain, sigma, epsilon, rcutoff_sq, K, R, r0)
    
    # 输出轨迹文件
    traj = open('2GN_2x20_T_10.xyz', 'w')
    traj_txt = open('2GN_2x20_T_10.txt', 'w')
    
    for step in range(n_steps):
        if step % 1000 == 0:
            traj.write(f"{total_N}\n\n")
            for i in range(total_N):
                # 用不同元素区分两条链：第一条链用C，第二条用O
                element = "C" if i < N_chain else "O"
                traj.write(f"{element} {coord[i][0]:10.5f} {coord[i][1]:10.5f} {coord[i][2]:10.5f}\n")
                traj_txt.write(f"{coord[i][0]:10.5f} {coord[i][1]:10.5f} {coord[i][2]:10.5f}\n")
        
        # 执行MC移动尝试
        coord_trial = move(coord, N_chain, max_delta)
        energy_trial = total_energy(coord_trial, N_chain, sigma, epsilon, rcutoff_sq, K, R, r0)
        delta_e = energy_trial - energy_current
        
        if accept(delta_e, T):
            coord = coord_trial
            energy_current = energy_trial
        
        if step % 1000 == 0:
            print(f"Step {step}, Current Energy: {energy_current:.2f}")
    
    traj.close()
    traj_txt.close()

细节说明 #

1. 链生成逻辑：gen_chain函数现在支持指定每条链的单体数和链间距，第二条链的x坐标被设为间距值，确保初始状态下两条链不会重叠。
2. 固定单体处理：在move函数中明确跳过两条链的首个单体（索引0和20），确保它们始终锚定在表面，不会被移动。
3. 能量计算优化：分别遍历两条链的键对计算FENE键能，非键相互作用自动包含链间的LJ相互作用，符合模拟需求。
4. 可视化友好输出：用"C"和"O"分别标识两条链的粒子，你可以用VMD、PyMOL等软件直接打开XYZ轨迹文件，清晰观察两条链的坍缩过程。

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