0

我想生成两个彼此相距一定距离的 20 个单体的线性链。以下代码生成单个链。有人可以帮助我如何生成第二条链吗？两条链固定在一个表面上，即链的第一个单体是固定的，其余的单体在 xyz 方向上自由移动，但单体的 z 分量应该是正的。

像这样的东西： 在此处输入图像描述

import numpy as np
import numba as nb
#import pandas as pd

@nb.jit()
def gen_chain(N):

    x = np.zeros(N)
    y = np.zeros(N)
    z = np.linspace(0, (N)*0.9, num=N)

    return np.column_stack((x, y, z)), np.column_stack((x1, y1, z1))

    #coordinates = np.loadtxt('2GN_50_T_10.txt', skiprows=199950)
    #return coordinates

@nb.jit()
def lj(rij2):

    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):

    return (-0.5 * K * np.power(R, 2) * np.log(1 - ((np.sqrt(rij2) - r0) / R)**2))

@nb.jit()
def total_energy(coord):
    # Non-bonded energy.
    e_nb = 0.0
    for i in range(N):
        for j in range(i - 1):
            ri = coord[i]
            rj = coord[j]
            rij = ri - rj
            rij2 = np.dot(rij, rij)
            if (rij2 < rcutoff_sq):
                e_nb += lj(rij2)

    # Bonded FENE potential energy.
    e_bond = 0.0
    for i in range(1, N):
        ri = coord[i]
        rj = coord[i - 1] # Can be [i+1] ??
        rij = ri - rj
        rij2 = np.dot(rij, rij)
        e_bond += fene(rij2)
    return e_nb + e_bond

@nb.jit()
def move(coord):

    trial = np.ndarray.copy(coord)
    for i in range(1, N):
        while True:
            delta = (2 * np.random.rand(3) - 1) * max_delta

            trial[i] += delta
        #while True:
            if trial[i,2] > 0.0:
                break
            trial[i] -= delta



    return trial

@nb.jit()
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 potential parameters.
    K = 40.0
    R = 0.3
    r0 = 0.7
    # L-J potential parameters
    sigma = 0.5716
    epsilon = 1.0
    # MC parameters
    N = 20   # Numbers of monomers
    rcutoff = 2.5 * sigma
    rcutoff_sq = rcutoff * rcutoff
    max_delta = 0.01
    n_steps = 100000
    T = 10

    # MAIN PART OF THE CODE
    coord = gen_chain(N)
    energy_current = total_energy(coord)

    traj = open('2GN_20_T_10.xyz', 'w')
    traj_txt = open('2GN_20_T_10.txt', '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]))
                traj_txt.write("%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()

我除了粒子链坍缩成一个球体。

4

3 回答 3

0

您可以在下面找到代码的修改版本。为了让事情更简单，在这个版本中没有盒子和边界条件，只有自由空间中的一条链。链被初始化为粒子的线性序列，每个粒子距离R0下一个粒子的 80%，因为R0FENE 键的最大长度。该代码认为粒子i被结合i+1并且结合没有被破坏。这段代码只是一个概念证明。

#!/usr/bin/python

import numpy as np

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(rij2):
    sig_by_r6 = np.power(sigma/rij2, 3)
    sig_by_r12 = np.power(sig_by_r6, 2)
    lje = 4.0 * epsilon * (sig_by_r12 - sig_by_r6)
    return lje

def fene(rij2):
    return (-0.5 * K * R0**2 * np.log(1-(rij2/R0**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
            rij2 = np.dot(rij, rij)
            if (rij2 < rcutoff):
                e_nb += lj(rij2)
    # Bonded
    e_bond = 0
    for i in range(1, N):
        ri = coord[i]
        rj = coord[i-1]
        rij = ri - rj
        rij2 = np.dot(rij, rij)
        e_bond += fene(rij2)
    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 = 1.5

    # LJ parameters
    sigma = 1.0
    epsilon = 1.0

    # MC parameters
    N = 50 # number of particles
    rcutoff = 3.5
    max_delta = 0.01
    n_steps = 10000000
    T = 1.5

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

    traj = open('traj.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()

代码在每一步打印当前配置，您可以将其加载到 VMD 上并查看其行为方式。键最初不会在 VMD 上正确显示，您必须对粒子使用珠子表示并手动定义键或使用 VMD 中的脚本定义键。在任何情况下，您都不需要看到债券就可以注意到链条没有崩溃。

请记住，如果要以一定的密度模拟链，则需要小心生成正确的拓扑。我推荐EMC封装以在所需的热力学条件下有效地生成聚合物。这绝不是一个微不足道的问题，尤其是对于较大的连锁店。

顺便说一句，您的代码在 FENE 能量评估中出现错误。rij2已经平方了，你又平方了。

下面您可以看到总能量作为步数的函数对于T = 1.0、N = 20、rcutoff = 3.5以及 10000 步后的最后当前配置如何表现。

能量与步骤

10000步的当前配置

下面是N = 50, T = 1.5, max_delta = 0.01, K = 40, R = 1.5, rcutoff = 3.5, 和 1000 万步。这是当前的最后一个配置。

在此处输入图像描述

完整的“轨迹”，它不是真正的轨迹，因为这是 MC，你可以在这里找到（它小于 6 MB）。

于 2019-09-09T22:53:19.237 回答
0

您正在实施的 MC 的逻辑存在一些问题。要执行 MC，您需要尝试移动，评估新状态的能量，然后根据随机数接受/拒绝。在您的代码中，没有任何试图移动粒子的迹象。您需要移动一个（或多个），评估能量，然后更新您的坐标。顺便说一句，我想这不是你的全部代码。有许多参数未定义，例如您的 fene 电位中的“k”和“R0”

于 2019-09-09T09:39:28.230 回答
0

FENE 势模型键相互作用。您的代码是说截止内的所有粒子都由 FENE 弹簧粘合，并且键不是固定的，而是由截止定义的。当 r_cutoff = 3.0，大于 LJ 井的平衡距离时，您实际上是在考虑每个粒子都可能与许多其他粒子结合。您将 FENE 电位视为非键合电位。

对于键相互作用，您应该忽略截止值，只评估根据您的拓扑键合的实际对的能量，这意味着您首先需要定义一个拓扑。我建议在一个足够大以包含整个拉伸分子的盒子中生成一个 N 原子的线性分子，并考虑第 i 个原子与第 (i-1) 个原子键合，其中 i = 2，...， N. 这样，拓扑就得到了很好的定义和持久化。然后分别考虑这两种相互作用，非键合和键合，并将它们添加到最后。

像这样的东西，在伪代码中：

e_nb = 0
for particle i = 1 to N:
    for particle j = 1 to i-1:
        if (dist(i, j) < rcutoff):
            e_nb += lj(i, j)

e_bond = 0
for particle i = 2 to N:
    e_bond += fene(i, i-1)

e_tot = e_nb + e_bond
于 2019-09-09T10:32:19.713 回答