MDAnalysis比较离子在弹道中的位置

Question

因此，我有一个系统，我需要能够确定我的离子的确切位置，并在这个离子的平均位置上运行一个方程。我发现我的离子位置是不一致的，因为一些离子绕过了周期边界，严重改变了那个窗口的位置。当离子在+40和-40之间移动时，我的平均值是+20。

我想通过实现一种方法来纠正这一点，在我的盒子边缘打开我包裹好的离子坐标。

本质上，我在想，对于我轨道上的每一帧，MDAnalysis会检查帧1中离子1的位置，然后在第2帧中，它会再次检查相同的离子，并将其与前一个位置进行比较。例如，如果它从+坐标转到-坐标，那么我就会有一个计数，它添加了+1，这意味着它包装了一次。如果从-到+，我会让它减去1。然后在所有帧的末尾，我会得到一个数字，可以帮助我识别我如何进行分析。

然而，我的编码技能并不那么出色，我想知道我将如何实现这一点？实际上，我已经把计数降低了，但是帧之间的比较让我感到困惑。我该怎么做这个比较？

提前感谢

Answer 1

回答这个问题有几种方法。首先，

本质上，我在想，对于我轨道上的每一帧，MDAnalysis会检查帧1中离子1的位置，然后在第2帧中，它会再次检查相同的离子，并将其与前一个位置进行比较。例如，如果它从+坐标转到-坐标，那么我就会有一个计数，它添加了+1，这意味着它包装了一次。如果从-到+，我会让它减去1。然后在所有帧的末尾，我会得到一个数字，可以帮助我识别我如何进行分析。

你可以编写自己的分析类。一种未经测试的方法是下面的原型--本教程将更多地介绍每个方法(_prepare、_conclude等)所做的工作。

from MDAnalysis.analysis.base import AnalysisBase
import numpy as np

class CountWrappings(AnalysisBase):

    def __init__(self, universe, select="name NA"):
        super().__init__(universe.universe.trajectory)
        # these are your selected ions
        self.atomgroup = universe.select_atoms(select)
        self.n_atoms = len(self.atomgroup)

    def _prepare(self):
        # self.results is a dictionary of results
        self.results.wrapping_per_frame = np.zeros((self.n_frames, self.n_atoms), dtype=bool)
        self._last_positions = self.atomgroup.positions
    
    def _single_frame(self):
        # does sign change for any element in 2D array?
        compare_signs = np.sign(self.atomgroup.positions) == np.sign(self._last_positions)
        sign_changes_any_axis = np.any(compare_signs, axis=1)
        # _frame_index is the relative index of the frame being currently analyzed
        self.results.wrapping_per_frame[self._frame_index] = sign_changes_any_axis
        self._last_positions = self.atomgroup.positions
    
    def _conclude(self):
        self.results.n_wraps = self.results.wrapping_per_frame.sum(axis=0)


n_wraps = CountWrappings(my_universe, select="name NA CL MG")
n_wraps.run()
print(n_wraps.results.wrapping_per_frame)
print(n_wraps.results.n_wraps)

然而，我不确定这是否符合你的实际目标：

我想通过实现一种方法来纠正这一点，在我的盒子边缘打开我包裹好的离子坐标。

你在计算离子的相对位置吗？您可以在每个离子和中心之间添加键，这样就可以使用AtomGroup.unwrap()函数。或者，您的数据是否与GROMACS兼容？GROMACS有一个名为“no跳转”的展开实用程序，它可以解开跨越方框边缘的原子。

gmx trjconv -f my_trajectory.xtc -s my_topology.gro -pbc nojump -o my_unwrapped_trajectory.xtc

Answer 2

正如莉莉提到的，您可以编写自己的分析来完成这一任务，或者使用GROMACS。然而，莉莉的例子和GROMACS实现的'nojump‘没有考虑到不扩散条约组合下的盒子大小波动(假设你已经使用了NPT)。等人几年前就写过这个普遍存在的问题。据我所知，只有在LiPyphilic (免责声明:我是LiPyphilic的作者)中才能实现能解释框大小波动的nojump解包装。

使用LiPyphilic，您可以像这样打开轨道：

import MDAnalysis as mda
from lipyphilic.transformations import nojump

u = mda.Universe(pdb, xtc)
ions = u.select_atoms('name NA CLA')

u.trajectory.add_transformations(
    nojump(
        ag=ions,
        nojump_x=True,
        nojump_y=True,
        nojump_z=True)
)

然后，当你用你的MDAnalysis宇宙做进一步的分析时，原子将在每个框架中自动展开。