轨迹

在MDAnalysis中，静态数据存于您的universe拓扑中，而动态数据则来自“Universe.trajectory”中的轨迹数据。这些数据通常是从轨迹文件中加载的，其中包括以下信息：

  • 原子坐标(Universe.atoms.positions)
  • 盒子尺寸(Universe.dimensions)
  • 速度和力(Universe.atoms.velocities)

尽管这些属性看起来是静态的，但它们实际上是动态的，其中包含的数据可能会发生变化。 为了保持内存效率，MDAnalysis 不会一次将轨迹的每一帧加载到内存中。相反，一个宇Universe有一个状态：它当前在轨迹中与之关联的特定时刻。当时刻发生变化时，上述属性中的数据也会相应地移动。

改变时刻的常见方法是通过索引来实现。Universe.trajectory可以被视为一系列时刻的列表，这是一种存储当前时间范围相关信息的数据结构。例如，您可以查询其长度。

1
2
3
4
5
import MDAnalysis as mda
from MDAnalysis.tests.datafiles import PSF, DCD
u = mda.Universe(PSF, DCD)
len(u.trajectory)
#98

首次从文件加载轨迹时，默认情况下，它设置为第一帧（索引为 0）。索引轨迹会返回该帧的时刻，并将Universe设置为指向该帧，直到下一次更改时刻。
许多任务涉及将函数应用于轨迹的每一帧。对于这些，您需要迭代帧，即使您不直接使用时间步长。这是因为迭代行为将Universe移动到下一帧，从而改变动态原子坐标。
如果您只想处理帧的子集，也可以对轨迹进行切片。

1
2
3
4
5
6
7
8
9
10
11
protein = u.select_atoms('protein')

for ts in u.trajectory[:20:4]:
    rad = protein.radius_of_gyration()
    print('frame={}: radgyr={}'.format(ts.frame, rad))

## frame=0: radgyr=16.669018368649777
## frame=4: radgyr=16.743960893217544
## frame=8: radgyr=16.78938645874581
## frame=12: radgyr=16.87231363208217
## frame=16: radgyr=17.003316543310998

请注意，在迭代轨迹后，帧始终设置回第一帧，即使循环在轨迹结束之前停止也是如此。

由于MDAnalysis将直接从读取的文件中提取轨迹数据，因此一旦更改帧，对原子坐标和长方体尺寸的更改将不会持续存在。使这些更改永久化的唯一方法是将轨迹加载到内存中，或者为每一帧写入一个新的轨迹到文件中。例如，要为每一帧设置立方框大小并将其写入文件：

有时您可能希望只转换部分轨迹，或者不写出文件。在这些情况下，MDAnalysis 支持在读取帧时对帧执行的“及时”转换。

切片轨迹

MDAnalysis轨迹可以被索引返回Timestep，或者被切片返回FrameIterator。

1
2
3
4
5
6
7
8
9
10
import MDAnalysis as mda

from MDAnalysis.tests.datafiles import PSF, DCD

u = mda.Universe(PSF, DCD)

u.trajectory[4], u.trajectory[:4], 

## (< Timestep 4 >,
##  <MDAnalysis.coordinates.base.FrameIteratorSliced at 0x213387af140>)

索引轨迹就是将Universe指向特定帧，更新动态数据，如Universe.atoms.positions。
通过切片轨迹创建的FrameIterator不是将Universe指向新帧，而是迭代切片轨迹会将轨迹倒回第一帧。

1
2
3
4
fiter = u.trajectory[10::10]
frames = [ts.frame for ts in fiter]
print(frames, u.trajectory.frame)
## [10, 20, 30, 40, 50, 60, 70, 80, 90] 0

您还可以通过布尔索引和多重索引来创建一个切片轨迹。布尔索引允许您仅选择满足特定条件的帧，方法是传递一个与原始轨迹长度相同的ndarray。只有布尔值为True的帧才会包含在生成的FrameIterator中。例如，要仅选择RMSD小于2 Å的轨迹帧：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
from MDAnalysis.analysis import rms
protein = u.select_atoms('protein')
rmsd = rms.RMSD(protein, protein).run()
bools = rmsd.results.rmsd.T[-1] < 2
print(bools)
## [ True  True  True  True  True  True  True  True  True  True  True  True
##   True  True False False False False False False False False False False
##  False False False False False False False False False False False False
##  False False False False False False False False False False False False
##  False False False False False False False False False False False False
##  False False False False False False False False False False False False
##  False False False False False False False False False False False False
##  False False False False False False False False False False False False
##  False False]

你也可以切片一个FrameIterator来创建一个新的FrameIterator。

及时转换

即时转换是一种在动态数据加载至内存时自动对其进行修改的功能（通常是指坐标数据）。每个时刻都会调用该转换函数，以将数据转换为所需的表示形式。转换函数还必须返回当前时间步信息，因为转换通常是按顺序串联在一起的。
MDAnalysis.transformation模块包含了一系列转换功能。例如，fit_rot_trans()可以对一个AtomGroup进行质量加权对齐，使其与参考对象对齐。

1
2
3
4
5
6
7
import MDAnalysis as mda
from MDAnalysis.tests.datafiles import TPR, XTC
from MDAnalysis import transformations as trans
u = mda.Universe(TPR, XTC)
protein = u.select_atoms('protein')
align_transform = trans.fit_rot_trans(protein, protein, weights=protein.masses)
u.trajectory.add_transformations(align_transform)

其他已实现的变换包括用于translate, rotate, fit一个AtomGroup到参考, 以及在晶胞内wrap或unwrap组。（请参阅 MDAnalysis 实时变换博客文章，其中包含了对这些拟合和包裹函数更完整的介绍。）
尽管只能调用一次 add_transformations() 函数，但可以将多个变换以列表形式传入，它们将按顺序执行。
还有一个变换教程，更详细地介绍了如何使用变换。下面给出了几个简单的示例。

示例

示例工作流如下:

  • 使所有分子完整（在周期性边界条件下解开它们），
  • 将蛋白质置于盒子的中心
  • 将水包回盒子里
1
2
3
4
5
6
7
u = mda.Universe(TPR, XTC)
protein = u.select_atoms('protein')
water = u.select_atoms('resname SOL')
workflow = [trans.unwrap(u.atoms),
            trans.center_in_box(protein, center='geometry'),
            trans.wrap(water, compound='residues')]
u.trajectory.add_transformations(*workflow)

如果你的转换不依赖于Universe的某些东西(如选择AtomGroup)，你可以在创建Universe时直接使用变换。下面的代码将坐标沿z轴向上怕平移1 Å:
u = mda.Universe(TPR, XTC, transformations=[trans.translate([0, 0, 1])])

自定义变换

从本质上讲，转换函数必须只将时刻作为其输入并返回时刻作为输出。

1
2
3
4
5
6
7
8
9
import numpy as np

def up_by_2(ts):
    """Translates atoms up by 2 angstrom"""
    ts.positions += np.array([0.0, 0.0, 0.2])
    return ts


u = mda.Universe(TPR, XTC, transformations=[up_by_2])

如果您的转换需要其他参数，则需要使用可以接受其他参数的包装函数包装核心转换。

1
2
3
4
5
6
7
8
9
10
11
def up_by_x(x):
    """Translates atoms up by x angstrom"""
    def wrapped(ts):
        """Handles the actual Timestep"""
        ts.positions += np.array([0.0, 0.0, float(x)])
        return ts
    return wrapped


## load Universe with transformations that move it up by 7 angstrom
u = mda.Universe(TPR, XTC, transformations=[up_by_x(5), up_by_x(2)])

或者，你可以使用functools.partial()来提交其他参数。

1
2
3
4
5
6
7
8
9
10
import functools

def up_by_x(ts, x):
    ts.positions += np.array([0.0, 0.0, float(x)])
    return ts


up_by_5 = functools.partial(up_by_x, x=5)

u = mda.Universe(TPR, XTC, transformations=[up_by_5])

及时转换函数可以应用于时刻的任何属性，而不仅仅是原子位置。例如，要为轨迹的每一帧提供一个框：

1
2
3
4
5
6
def set_box(ts):
    ts.dimensions = [10, 20, 30, 90, 90, 90]
    return ts


u = mda.Universe(TPR, XTC, transformations=[set_box])

翻译https://userguide.mdanalysis.org/stable/trajectories/trajectories.html