用metropolis模拟高分子链

Question

我试图模拟聚合物链，但在函数total_energy中出现错误“包含多个元素的数组的真值不明确。请使用a.any()或a.all()”。在我看来，我在定义数组和列表之间的区别时犯了一个错误。谁能指出问题所在？

def gen_chain(N, R):
    x = np.linspace(0, (N - 1) * R * 0.8, num=N)
    y = np.zeros(N)
    z = np.zeros(N)
    return np.column_stack((x, y, z))

def lj(rij): # Lennard-Jones potential
    lje = 4.0 * epsilon* (sigma**12/rij**12 - sigma**6/r**6)

    return lje

def fene(rij):

    return (-0.5 * K * R0**2 * np.log(1 - ((rij - r0) / R)**2))

def total_energy(coord):
    # Non-bonded
    e_nb = 0
    for i in range(N):
        for j in range(i-1):
            ri = coord[i]
            rj = coord[j]
            rij = ri - rj
            if (rij < rcutoff): # HERE THE ERROR COMES
                e_nb += lj(rij)

    e_bond = 0
    for i in range(1, N):
        ri = coord[i]
        rj = coord[i-1]
        rij = ri - rj
        e_bond += fene(rij)
    return e_nb + e_bond

def move(coord):
    trial = np.ndarray.copy(coord)
    for i in range(N):
        delta = (2.0 * np.random.rand(3) - 1) * max_delta
        trial[i] += delta
    return trial

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

if __name__ == "__main__":

    # FENE parameters
    K = 40
    R0 = 0.3
    r0 = 0.7

    # LJ parameters
    sigma = r0 / 2**1 / 6  # 0.572
    epsilon = 1.0

    # MC parameters
    N = 50  # number of particles
    rcutoff = (2.5 * sigma)**2
    max_delta = 0.01
    n_steps = 10000
    T = 0.5

    coord = gen_chain(N, R0)
    energy_current = total_energy(coord)

    traj = open('traj_T_05.xyz', 'w')

    for step in range(n_steps):
        if step % 1000 == 0:
            traj.write(str(N) + '\n\n')
            for i in range(N):
                traj.write("C %10.5f %10.5f %10.5f\n" % (coord[i][0], coord[i][1], coord[i][2]))
            print(step, energy_current)
        coord_trial = move(coord)
        energy_trial = total_energy(coord_trial)
        delta_e = energy_trial - energy_current
        if accept(delta_e):
            coord = coord_trial
            energy_current = energy_trial

    traj.close()

Answer 1

您得到的特定错误是由于rij是一个数组，而rcutoff是一个数字。您可以计算rij的长度，并将其与截止长度进行比较：

if np.sqrt(np.dot(rij, rij)) < rcutoff:

或者使用define rcutoff_sqr = rcutoff**2并使用它来保存一些平方根计算：

if np.dot(rij, rij) < rcutoff_sqr:

请注意，在此之后，您的代码将无法正常工作，因为您必须修复更多的错误(未定义的变量，除以零，甚至更多)。