利用MDAnalysis获取径向分布函数

Question

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

这个是可能的吗？我在木星笔记本中使用以下代码为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

为了在MDAnalysis中重用分析工具，可以使用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参数，该参数允许您缩小要具体使用的框架。