Task 2: 𝑅!! and 𝑅" # for N个frames for n_mol个分子 𝑅44 = 𝑅 𝑛 − 𝑅 0 计算 𝑅44 i，𝑡 𝑅5 ! = ' 8 ∑ #6' 0 𝑚 # 𝑟# − 𝑟-7 ! 计算𝑅9 ! i，𝑡 • 输入： MD轨迹文件，dr， 开始帧，结束帧，帧的间隔 • 输出： 1. 对于某一条高分子链（第i条）作图，体现𝑅!! 𝑡 和 𝑅" # 𝑡 随时间的变化 2. 对于一段时间（p个frames），对于所有的高分子链的𝑅!! 和 𝑅" # 作统计分布

首先，需要解析 MD 轨迹文件，获取每个分子的坐标信息。可以使用 Python 中的 MDAnalysis 库来实现：

import MDAnalysis as mda

u = mda.Universe("trajectory.xtc", "topology.pdb")

这里假设轨迹文件为 `trajectory.xtc`，拓扑文件为 `topology.pdb`。

然后，需要定义计算函数来计算 $R_{44}$ 和 $R_{5}$：

def compute_R44(u, i, t):
    """
    计算第 i 条高分子链在 t 时刻的 R44 值
    """
    sel = u.select_atoms(f"resid {i}")
    r_i = sel.center_of_geometry()
    r_0 = u.select_atoms("all").center_of_geometry()
    r_44 = r_i - r_0
    return r_44

def compute_R5(u, i, t):
    """
    计算第 i 条高分子链在 t 时刻的 R5 值
    """
    sel = u.select_atoms(f"resid {i}")
    r_i = sel.center_of_geometry()
    r_7 = u.select_atoms(f"around {dr} resid {i}").center_of_geometry()
    r_5 = r_i - r_7
    return r_5

其中，`i` 表示第几条高分子链，`t` 表示时间，`dr` 表示在计算 $R_5$ 时，选取周围分子的半径。

接下来，可以使用 Matplotlib 库来绘制 $R_{44}$ 和 $R_{5}$ 随时间的变化。

import numpy as np
import matplotlib.pyplot as plt

start_frame = 0
end_frame = 100
skip_frame = 10
p = 10  # 分布统计的帧数

# 计算 R44 和 R5
R44_data = []
R5_data = []
for ts in u.trajectory[start_frame:end_frame:skip_frame]:
    R44_i = compute_R44(u, i, ts.time)
    R5_i = compute_R5(u, i, ts.time)
    R44_data.append(R44_i)
    R5_data.append(R5_i)

# 绘制 R44 随时间的变化
R44_data = np.array(R44_data)
plt.plot(u.trajectory[start_frame:end_frame:skip_frame].time, R44_data[:, 2])
plt.xlabel("Time (ps)")
plt.ylabel("R44 (nm)")
plt.title("R44 vs. Time")
plt.show()

# 绘制 R5 随时间的变化
R5_data = np.array(R5_data)
plt.plot(u.trajectory[start_frame:end_frame:skip_frame].time, R5_data[:, 2])
plt.xlabel("Time (ps)")
plt.ylabel("R5 (nm)")
plt.title("R5 vs. Time")
plt.show()

# 绘制 R44 和 R5 的分布统计
R44_statistics = []
R5_statistics = []
for ts in u.trajectory[start_frame:end_frame:p]:
    R44_i = compute_R44(u, i, ts.time)
    R5_i = compute_R5(u, i, ts.time)
    R44_statistics.append(R44_i[2])
    R5_statistics.append(R5_i[2])

plt.hist(R44_statistics, bins=20, density=True, alpha=0.5, label="R44")
plt.hist(R5_statistics, bins=20, density=True, alpha=0.5, label="R5")
plt.xlabel("Distance (nm)")
plt.ylabel("Probability Density")
plt.title("Distribution of R44 and R5")
plt.legend()
plt.show()

这里用到了 `numpy` 库来处理数据，`start_frame`、`end_frame` 和 `skip_frame` 分别表示轨迹文件中的开始帧、结束帧和帧的间隔，`p` 表示分布统计的帧数。