0

我一直在试图弄清楚这个非常依赖于收费的过程。但是，我发现，对于水，我得到的水分子电荷总和的值不等于 0。为了对此进行调查，我想比较 MDAnalysis 和我的拓扑文件使用的费用值，看看问题可能出在哪里。

在我的初始代码中，我可以使用以下方法来收取费用：

    for ts in u.trajectory[ini_frames:]:
    count += 1
    membrane = u.select_atoms('segid PROA or segid PROB or segid MEMB')
    COM = membrane.atoms.center_of_mass()[2]
    q_prof = CLAs.atoms.charges * (CLAs.positions[:,2] + (Lz/2 - COM))/Lz
    Q_instant = np.sum(q_prof)
    Q_sum += Q_instant
Q_av = Q_sum / n_frames

这对我有用

CLAs = u.select_atoms("segid HETA or segid PROA or segid PROB or segid MEMB or segid IONS")

因此，为了查看这种差异可能来自哪里，我尝试了：

def Q_code(dcd, topo):
Lz = u.dimensions[2]
Q_sum = 0
count = 0
CLAs = u.select_atoms('segid TIP3')
ini_frames = -20
n_frames = len(u.trajectory[ini_frames:])
for ts in u.trajectory[ini_frames:]:
    for i in CLAs:
        with open('Q_atom_temp.csv', 'a') as f:
            new_line = [s, i, i.charges, CV1_dict[s], CV2_dict[s]]
            writes = csv.writer(f)
            writes.writerow(new_line)
            f.close()
fields = ['Window', 'Atom', 'Charge', 'CV1', 'CV2'] with
 open('Q_atom_temp.csv', 'a') as f:
      writer = csv.writer(f)
      writer.writerow(fields) f.close()

但我得到了错误： 在此处输入图像描述 那么我将如何制作一个文件来显示哪个原子对应于哪个电荷？

指标不一样吗？谢谢你的帮助！

4

1 回答 1

0

事实证明，这是以相同的方式索引的。我的问题是我将事物应用于每个单独的原子而不是其分子。当我这样做时，我得到了一个有意义的答案。

因此，对于阅读本文的其他人，您可以使用相同的索引来遍历您的头寸和费用。

for i in range(len(CLAs.positions[:,2])):
        #print(CLAs.atoms.charges[i], CLAs.atoms.names[i], CLAs.positions[i,2])
        total_charge += CLAs.atoms.charges[i]
        count += 1
        position += CLAs.positions[i,2]
        if count == 5:
            count = 0
            av_position = (position / 5)
            q_prof += total_charge * (av_position + (Lz/2 - COM))/Lz
            position = 0
            total_charge = 0

该代码查看单个电荷，将其映射到同一离子的位置，然后在水分子完成后，获取分子的电荷并将其乘以相对于盒子的平均位置。

于 2022-01-07T02:40:00.200 回答