1

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

我想通过实施一种方法来解开包装盒边缘的离子的包装坐标来纠正这个问题。

本质上，我在想，对于我轨迹中的每一帧，MDAnalysis 会检查 ION 1 在第 1 帧中的位置。然后在第 2 帧中，它会再次检查相同的离子并将其与之前的位置进行比较。例如，如果它从 + 坐标变为 - 坐标，那么我会有一个添加 +1 的计数，这意味着它包裹了一次。如果它从 - 变为 +，我会让它减去 1。然后在所有帧结束时，我会有一个数字，可以帮助我确定如何执行我的分析。

但是我的编码技能并不乏味，我想知道我将如何实现这一点？我基本上已经倒计时了，但是帧之间的比较是我感到困惑的地方。我将如何进行这种比较？

提前致谢

4

2 回答 2

1

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

基本上我在想，对于我轨迹中的每一帧，MDAnalysis 会检查 ION 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 有一个名为“nojump”的解包实用程序，它可以解包跨越盒子边缘的原子，例如

gmx trjconv -f my_trajectory.xtc -s my_topology.gro -pbc nojump -o my_unwrapped_trajectory.xtc
于 2021-12-07T18:22:18.870 回答
0

正如 Lily 提到的，您可以编写自己的分析来执行此操作或使用 GROMACS。但是，Lily 的示例和“nojump”的 GROMACS 实现都无法解释 NPT 合奏下的框大小波动（假设您使用了 NPT）。冯布洛等人。几年前写过这个普遍存在的问题。据我所知，解决盒子大小波动的 nojump 展开的唯一实现是在LiPyphilic中（免责声明：我是 LiPyphilic 的作者）。

使用 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 Universe 进行进一步分析时，原子将在每一帧自动展开。

于 2022-01-19T15:13:03.830 回答