如何在MDAnalysis中利用轨迹片段计算RDF？尝试方法结果异常

我有一条超长分子动力学（MD）轨迹，想选取轨迹的特定片段（比如第100至1000帧）计算径向分布函数（RDF）。已经完成全轨迹的RDF计算，但处理片段轨迹时出了问题——我提取片段帧坐标创建新的MDAnalysis Universe，再调用InterRDF模块计算，得到的结果完全没有RDF的特征。

我的尝试代码如下：

frame_count = len(all_p.trajectory)

half_frame_count = int(frame_count/2)
frames_first_half = []
frames_second_half = []

for ts in all_p.trajectory[:half_frame_count]:
    frames_first_half.append(ts.positions.copy())
for ts in all_p.trajectory[half_frame_count:]:
    frames_second_half.append(ts.positions.copy())

u_fragment_1 = mda.Universe(PDB, coordinates=frames_first_half)
A_1 = u_fragment_1.select_atoms('*selection A*') 
B_1 = u_fragment_1.select_atoms('*selection B*') 

u_fragment_2 = mda.Universe(PDB, coordinates=frames_second_half)
A_2 = u_fragment_2.select_atoms('*selection A*') 
B_2 = u_fragment_2.select_atoms('*selection B*') 

随后执行RDF计算：

irdf = rdf.InterRDF(CA_A_1, CA_B_1,
                    nbins=2000,
                    range=(0.0,100.0),
                    exclusion_block=(1,1),
                    )
irdf.run()

rdf_x_axis = irdf.results.bins
rdf_y_axis = irdf.results.rdf
plt.plot(rdf_x_axis, rdf_y_axis)
plt.show()
plt.clf()

irdf = rdf.InterRDF(CA_A_2, CA_B_2,
                    nbins=2000,
                    range=(0.0,100.0),
                    exclusion_block=(1,1),
                    )
irdf.run()

rdf_x_axis = irdf.results.bins
rdf_y_axis = irdf.results.rdf
plt.plot(rdf_x_axis, rdf_y_axis)
plt.show()

得到的图像完全不符合RDF特征，请问正确的片段轨迹RDF计算方法是什么？需要手动逐帧计算距离吗？

你的问题出在创建新Universe时没有传递轨迹的盒尺寸（unit cell dimensions）——RDF计算依赖于周期性边界条件，而你只复制了原子坐标，丢失了每个帧的盒子信息，导致InterRDF无法正确统计密度分布，结果自然不对。

不需要手动逐帧计算距离，用MDAnalysis可以更高效地处理，以下是两种正确的方法：

方法1：直接遍历轨迹片段并传递盒信息 #

不用创建新Universe，直接在原轨迹的指定片段上运行InterRDF，通过start、stop参数限定帧范围：

# 计算前半段轨迹的RDF
irdf_first = rdf.InterRDF(A, B,  # A、B是原Universe中已选好的原子组
                          nbins=2000,
                          range=(0.0, 100.0),
                          exclusion_block=(1,1))
# 限定只处理前half_frame_count帧
irdf_first.run(start=0, stop=half_frame_count)

# 计算后半段轨迹的RDF
irdf_second = rdf.InterRDF(A, B,
                           nbins=2000,
                           range=(0.0, 100.0),
                           exclusion_block=(1,1))
irdf_second.run(start=half_frame_count, stop=frame_count)

# 绘图
plt.plot(irdf_first.results.bins, irdf_first.results.rdf, label="First half")
plt.plot(irdf_second.results.bins, irdf_second.results.rdf, label="Second half")
plt.legend()
plt.show()

这种方法最简洁，直接复用原Universe的原子选择和轨迹元数据，避免信息丢失。

方法2：正确创建包含盒信息的片段Universe #

如果你确实需要创建独立的片段Universe，必须同时保存每个帧的positions和dimensions：

frames_first_half = []
dims_first_half = []
for ts in all_p.trajectory[:half_frame_count]:
    frames_first_half.append(ts.positions.copy())
    dims_first_half.append(ts.dimensions.copy())

# 创建Universe时同时传入坐标和盒尺寸
u_fragment_1 = mda.Universe(PDB, coordinates=frames_first_half, dimensions=dims_first_half)
A_1 = u_fragment_1.select_atoms('*selection A*') 
B_1 = u_fragment_1.select_atoms('*selection B*') 

# 后续RDF计算即可正常运行
irdf = rdf.InterRDF(A_1, B_1,
                    nbins=2000,
                    range=(0.0,100.0),
                    exclusion_block=(1,1),
                    )
irdf.run()

这样每个帧的周期性边界信息被保留，InterRDF能正确计算粒子数密度和分布。

额外注意事项 #

• 避免使用过大的nbins值：2000 bins对应0.05Å的分辨率，对于0-100Å的范围来说过于精细，会导致曲线噪声大，建议根据需求调整（比如200 bins，0.5Å分辨率）。
• 确保原子选择语句*selection A*是正确的，错误的选择会导致计算对象不对，结果异常。

内容的提问来源于stack exchange，提问作者rauk72