有轨迹文件用python按照半径计算聚苯乙烯体系的密度分布

按照半径计算聚苯乙烯体系的密度分布可以使用径向分布函数（Radial Distribution Function，RDF）来实现。以下是一个简单的Python示例代码，演示如何使用MDAnalysis库计算聚苯乙烯体系的RDF：

import numpy as np
import MDAnalysis as mda

# 读取轨迹文件和拓扑文件
u = mda.Universe('topology.pdb', 'trajectory.dcd')

# 获取聚苯乙烯分子
polymer = u.select_atoms('resname POL')

# 设置计算RDF的参数
n_bins = 100  # 网格数量
r_max = 10.0  # 最大半径
box = u.dimensions[:3]  # 模拟盒子大小

# 计算RDF
rdf, edges = np.histogramdd(mda.lib.distances.distance_array(polymer.positions, polymer.positions, box=box), bins=(n_bins, n_bins, n_bins), range=[(0, r_max), (0, r_max), (0, r_max)])
rdf = rdf / (4 * np.pi * edges[0][:-1]**2 * np.diff(edges[0])[0] * len(polymer) * u.trajectory.n_frames / box[0] / box[1] / box[2])

# 导出RDF数据
np.savetxt('rdf.txt', rdf)

在上面的示例代码中，我们首先使用MDAnalysis库读取了聚苯乙烯体系的轨迹文件和拓扑文件，并选择了聚苯乙烯分子。然后，我们设置了计算RDF的参数，包括网格数量、最大半径和模拟盒子大小。接下来，我们使用MDAnalysis库中的distance_array函数计算了聚合物分子之间的距离，并使用numpy库中的histogramdd函数计算了RDF。最后，我们对RDF进行了归一化，并将结果导出到文件中。

需要注意的是，上述示例代码仅适用于周期性边界条件下的聚苯乙烯体系，如果聚合物分子太大或者模拟盒子太小，可能需要使用更高级的计算方法来获得更准确的结果。