用mdp分析法求径向分布函数

Question

我正在GROMOS54a7中运行一个简单的苯模拟。 我想用MDAnalysis 1.0.0计算每个苯分子质心的RDF

这可能吗？我在Jupyter笔记本中使用以下代码为C分子g_cc（r）创建了rdf：

import MDAnalysis
import numpy as np
%matplotlib inline
import MDAnalysis.analysis.rdf as mda
import matplotlib.pyplot as plt

u = MDAnalysis.Universe("739-c6h6-MolDynamics.tpr","739-c6h6-MolDynamics_good-pbc.xtc")
s1 = u.select_atoms("resid 0 and type CAro")
s2 = u.select_atoms("not (resid 0) and type CAro")
rdf = mda.InterRDF(s1, s2)
rdf.run()

我想取每个苯分子（每个苯分子在我的模拟中都是一个残基），计算它的COM，然后在上面运行一个脚本。有可能做那样的事吗

关于RDF的一个普遍问题：我上面使用的方法是否使用我轨迹的每一帧构造RDF？我不知道文档中是否明确说明了这一点，因此，如果这是一个明显的问题，我深表歉意

谢谢你的建议

Answer 1

为了在MDA分析中重用分析工具，可以使用CG基团作为天然原子

下面是一个模拟MDAnalysis组的快速修复程序，它提供了一个新的positions属性。新的positions提供几何体的中心，而不是实际位置。我还覆盖了len，表示CG元素只使用了一个珠子

import MDAnalysis as mda
import numpy as np
import MDAnalysis.analysis.rdf
import matplotlib.pyplot as plt

u = mda.Universe('sys_solv.pdb','prod.dcd')

class CG:
    def __init__(self, ag):
        self.ag = ag
        self.universe = self.ag.universe
        self.trajectory = self.ag.universe.trajectory

    @property
    def positions(self):
        return np.array([self.ag.center_of_geometry()])

    def __len__(self):
        return 1

cg_selection = u.select_atoms('resid 1')
cg_atom = CG(cg_selection.atoms)
waters = u.select_atoms('name O')

rdf = MDAnalysis.analysis.rdf.InterRDF(cg_atom, waters)
rdf.run() 
plt.plot(rdf.bins, rdf.rdf)

验证：我为CG珠选择了一个原子，它复制了原始RDF

MDAnalysis使用整个轨迹。您可以在文档中为.run（）函数指定start/stop/step参数，该函数允许您缩小要专门使用的帧的范围