MDAnalysis:从同一个宇宙的读者中保存的轨迹是不正确的

Question

我有三个轨道复制(xtc)膜蛋白在模拟生理环境(水，离子，膜.)在MDAnalysis‘(2.2.0)宇宙中。我想保存另外三个只包含蛋白质轨迹(蛋白质原子)的xtcs，每一个原始xtc轨迹一个。当我试图遍历包含在宇宙中的三个MDAnalysis阅读器中的每一个时，第一个保存的轨迹似乎是正确的，但是其他两个在所有帧中都有相同的坐标。开始的，完整的轨迹是正确的。如果我的出发点必然是三个读者的宇宙，我如何正确和有效地做到这一点？

代码：

import MDAnalysis as mda
u = mda.Universe("11159_dyn_117.pdb", "11156_trj_117.xtc", "11157_trj_117.xtc", "11158_trj_117.xtc")
protein = u.select_atoms("protein")
protein.write("protein.pdb")

for num, reader in enumerate(u.trajectory.readers, 1):
    with mda.Writer(f"{num}.xtc", protein.n_atoms) as w:
        for ts in reader.trajectory:
            w.write(protein.atoms)

# Then check the generated individual trajectories by loading them in 
# Universes and checking the positions array. I checked them in PyMOL.

可下载的文件：https://submission.gpcrmd.org/dynadb/dynamics/id/117/ (模型文件和弹道文件)

Answer 1

您可以使用写一个轨迹方法直接从AtomGroup中提取AtomGroup.write(name, frames=trajectory_iterator)。使用私有的链式弹道属性访问ChainReader._start_frames中的开始/停止框架(未记录)。

import MDAnalysis as mda

# example data
from MDAnalysisTests import datafiles as data

# create a chained trajectory and select some atoms
u = mda.Universe(data.PSF, [data.DCD, data.DCD])
protein = u.select_atoms("protein")

# get start/stop frames: 
# array([  0,  98, 196]) for this example
sf = u.trajectory._start_frames

# write each subtrajectory of the chained trajectory
# to a new file in a different format (only containing
# the atoms of the selected AtomGroup)
for i, (start, stop) in enumerate(zip(sf[:-1], sf[1:])):
    protein.atoms.write(f"protein_{i}.xtc", frames=u.trajectory[start:stop])

这将产生轨迹protein_0.xtc和protein_1.xtc。如果要加载它们，请不要忘记创建包含所选内容的最小拓扑的文件。

protein.write("protein.gro")

这样你就可以用

p1 = mda.Universe("protein.gro", "protein_1.xtc")
p2 = mda.Universe("protein.gro", "protein_2.xtc")

Notes

• 我使用的是来自ChainReader._start_frames的开始/停止框架，而不是原始的轨迹长度，因为如果ChainReader 检测重叠帧在使用continuous=True时，它们会被适当地更新。
• 不使用frames=u.trajectory[start:stop]，也可以使用frames=slice(start, stop)。
• 不是写protein.atoms.write()，而是通常写protein.write()，这是等价的和更短的，但我想说明的是，去原子总是正确的，特别是如果一个人想要写一个完整的宇宙，在这种情况下，我们会使用u.atoms.write()。
• 不使用方便的方法，AtomGroup.write()也可以使用显式轨迹写入，在这里，您可以打开一个用mda.Writer("protein_1.xtc"，protein.n_atoms)作为W: ts in u.trajectorystart:stop: W.write(蛋白质)进行书写的轨迹。 它为每一步提供了更多的控制，但更冗长。