没有拓扑文件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），对于所有的高分子链的𝑅!! 和 𝑅" # 作统计分布

如果没有拓扑文件，可以使用 `MDAnalysis` 库中的 `Universe.from_trajectory()` 方法来生成 Universe 对象。该方法会自动从轨迹文件中读取拓扑信息。

具体实现可以参考以下代码：

import MDAnalysis as mda

u = mda.Universe("trajectory.xtc")

这里假设轨迹文件为 `trajectory.xtc`。

接下来，计算 $R_{44}$ 和 $R_{5}$ 的方法与之前相同，可以参考前面的代码。

绘制 $R_{44}$ 和 $R_{5}$ 随时间的变化和分布统计的方法也与之前相同，可以参考前面的代码。需要注意的是，在计算 $R_{5}$ 值时，需要使用 `u.select_atoms()` 方法来选取周围的分子，具体方法可以参考以下代码：

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} center ({r_i[0]}, {r_i[1]}, {r_i[2]})").center_of_geometry()
    r_5 = r_i - r_7
    return r_5

这里使用 `u.select_atoms(f"around {dr} center ({r_i[0]}, {r_i[1]}, {r_i[2]})")` 来选取距离 $r_i$ 不超过 $dr$ 的周围分子。

完整代码如下：

import numpy as np
import matplotlib.pyplot as plt

u = mda.Universe("trajectory.xtc")

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} center ({r_i[0]}, {r_i[1]}, {r_i[2]})").center_of_geometry()
    r_5 = r_i - r_7
    return r_5

start_frame = 0
end_frame = 100
skip_frame = 10
p = 10  # 分布统计的帧数
dr = 0.8  # 选取周围分子的半径

# 计算 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()

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