使用 MDAnalysis 获得径向分布函数

问题描述

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

这可能吗？我在 Jupyter Notebook 中使用以下代码为 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？我不知道文档中是否明确说明了这一点，因此如果这是一个明显的问题，我深表歉意。

感谢您的任何建议！

解决方法

为了重用 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() 函数的开始/停止/步长参数，它允许您缩小要具体使用的帧的范围。

chemistry mdanalysis physics physics python