1

我一直在尝试使用一种特定离子的部分电荷在 mdanalysis 中进行计算。我已经尝试过（这只是我知道抛出错误的代码片段）：

Cl = u.select_atoms('resname CLA and prop z <= 79.14') 
Lz = 79.14                          #Determined from system set-up
Q_sum = 0
COM = 38.42979431152344             #Determined from VMD
file_object1 = open(fors, 'a')
print(dcd, file = file_object1)
for ts in u.trajectory[200:]:
    frame = u.trajectory.frame
    time = u.trajectory.time
    for coord in Cl.positions:
        q= Cl.total_charge(Cl.position[coord][2])
        coords = coord - (Lz/COM)
        q_prof = q * (coords + (Lz / 2)) / Lz
        Q_sum = Q_sum + q_prof
        print(q)

但我不断收到与此相关的错误。

当它通过循环以在 MD 分析中获得电荷时，我将如何选择这个特定的原子？在我将 q 设置为等于一个常数并且代码运行良好之前，所以我知道只有这一行会引发错误：

q = Cl.total_charge(Cl.position[coord][2])

谢谢您的帮助！

4

1 回答 1

0

我想通了：

def Q_code(dcd, topo):
Lz = u.dimensions[2]
Q_sum = 0
count = 0
CLAs = u.select_atoms('segid IONS or segid PROA or segid PROB or segid MEMB')
ini_frames = -200
n_frames = len(u.trajectory[ini_frames:])
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
with open('Q_av.txt', 'a') as f:
    print('The Q_av for {}   is   {}'.format(s, Q_av), file = f)
return Q_av
于 2021-12-16T03:23:43.820 回答