读写文件

输入

从拓扑和坐标文件读取信息创建Universe:

1
2
import MDAnalysis as mda
u = mda.Universe('topology.gro', 'trajectory.xtc')

拓扑文件对于一个“Universe”来说总是必需的，而坐标文件则是可选的。有些文件格式同时提供拓扑和坐标信息。MDAnalysis支持多种格式，这些格式会根据文件扩展名自动检测。例如，上述代码加载了一个 GROMACS XTC 轨迹。可以加载多个坐标文件，如下所述；以下代码加载了两个CHARMM/NAMD DCD文件并将它们连接起来：

1
u = mda.Universe('topology.psf', 'trajectory1.dcd', 'trajectory2.dcd')

某些格式可以使用特定于格式的关键字参数进行加载，例如 LAMMPS DATA 原子样式规范。

读取多条轨迹

一个Universe可以加载多个轨迹，这些轨迹会按照给出的顺序进行连接。但有一个例外情况，即对于XTC和TRR文件。如果将“continuous=True”标志传递给Universe对象，MDAnalysis将会尝试将它们连接起来，以使轨迹尽可能地在时间上保持连续性。这意味着不会有重复的时间帧，也不会出现时间倒流的情况。
例如，以下所示的轨迹被分为四个部分。该列表示时间。如您所见，有些段落会重叠。当 continuous=True 时，只有标有“+”的帧才会被读取。

1
2
3
4
part01:  ++++--
part02:      ++++++-
part03:            ++++++++
part04:                         ++++

然而，可能会出现时间上的间隔（即某些帧似乎缺失了）。最终，确保这些轨迹能够有意义地拼接起来的责任在于用户。
注：当使用XTC/TRR且设置continuous=True，你不能混合使用不同类型的轨迹。
更多信息见ChainReader API文档。

轨迹格式

如果未提供格式关键字，ChainReader将会根据每个文件的扩展名来尝试猜测其格式。您可以通过使用格式关键字来强制ChainReader对每个文件都采用相同的格式。您还可以通过传入一系列（文件名，格式）元组的方式来指定针对每个文件应使用的格式。

1
2
3
4
5
6
import MDAnalysis as mda
from MDAnalysis.tests.datafiles import PDB, GRO

u = mda.Universe(PDB, [(GRO, 'gro'), (PDB, 'pdb'), (GRO, 'gro')])
u.trajectory
## <ChainReader containing adk_oplsaa.gro, adk_oplsaa.pdb, adk_oplsaa.gro with 3 frames of 47681 atoms>

内存中轨迹

将轨迹读取到内存中

如果您的设备有足够的内存能够将整个轨迹加载到内存中，那么将轨迹转移到内存中就能极大地加快分析速度。这使得能够使用现有的MDAnalysis工具直接处理原始坐标。此外，它还允许用户对轨迹中的坐标进行修改（例如通过 AtomGroup.positions），而无需将整个状态写入文件。
最简单直接的方法是将参数in_memory=True传递给“Universe”，它会自动将轨迹存储到内存中：

1
2
3
from MDAnalysis.tests.datafiles import TPR, XTC

universe = mda.Universe(TPR, XTC, in_memory=True)

MDAnalysis使用MemoryReader类读取这种数据。

将轨迹迁移到内存

可以在任何时候使用transefer_to_memory()方法将轨迹迁移到内存中。这个操作需要一些时间(传入verbose=True，transfer_to_memory()会显示一个进度条)。然而，一旦完成，后续操作速度就会大大加快。

在内存中构建轨迹

MemoryReader还可以直接将轨迹生成为一个 numpy 数组的形式。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
from MDAnalysisTests.datafiles import PDB
from MDAnalysis.coordinates.memory import MemoryReader
import numpy as np

universe = mda.Universe(PDB)

universe.atoms.positions
## array([[ 52.017,  43.56 ,  31.555],
##        [ 51.188,  44.112,  31.722],
##        [ 51.551,  42.828,  31.039],
##        ...,
##        [105.342,  74.072,  40.988],
##        [ 57.684,  35.324,  14.804],
##        [ 62.961,  47.239,   3.753]], shape=(47681, 3), dtype=float32)

“load_new()”方法可用于将坐标加载到宇宙中，从而替换掉原有的坐标：

1
2
3
4
5
6
7
8
9
10
11
12
13
coordinates = np.random.rand(len(universe.atoms), 3)

universe.load_new(coordinates, format=MemoryReader);

universe.atoms.positions
## array([[0.66662693, 0.938834  , 0.78954625],
##        [0.00282533, 0.22149895, 0.16824494],
##        [0.04847972, 0.7253196 , 0.03443228],
##        ...,
##        [0.7031928 , 0.84474987, 0.11842076],
##        [0.3746822 , 0.10468951, 0.01548544],
##        [0.4736413 , 0.12117212, 0.83662224]],
##       shape=(47681, 3), dtype=float32)

或者可以在创建Universe时直接加载：

1
2
3
4
5
6
7
8
9
10
universe2 = mda.Universe(PDB, coordinates, format=MemoryReader)
universe2.atoms.positions
## array([[0.66662693, 0.938834  , 0.78954625],
##        [0.00282533, 0.22149895, 0.16824494],
##        [0.04847972, 0.7253196 , 0.03443228],
##        ...,
##        [0.7031928 , 0.84474987, 0.11842076],
##        [0.3746822 , 0.10468951, 0.01548544],
##        [0.4736413 , 0.12117212, 0.83662224]],
##       shape=(47681, 3), dtype=float32)

内存轨迹的原子选择

要创建原子选择的轨迹，需要进行适当的单位转换。在使用Merge创建新Universe时通常会需要这样做，因为坐标不会自动加载进来。

输出

帧和轨迹

MDAnalysis 中的“Universe”对象可以通过“write()”方法以多种格式进行输出。例如，要将当前帧以 PDB 格式保存出来：

1
2
3
4
from MDAnalysis.tests.datafiles import PDB, TRR
u = mda.Universe(PDB, TRR)
ag = u.select_atoms("name CA")
ag.write("c-alpha.pdb")

传入"frame"关键字写出轨迹。 ag.write(‘c-alpha_all.xtc’, frames=‘all’)
切片或索引轨迹来选择写出的帧。

1
2
ag.write('c-alpha_skip2.trr', frames=u.trajectory[::2])
ag.write('c-alpha_some.dcd', frames=u.trajectory[[0,2,3]])

或者，可以逐帧使用 Writer() 对轨迹进行遍历。这需要您传入要写入的原子数量。

1
2
3
with mda.Writer('c-alpha.xyz', ag.n_atoms) as w:
    for ts in u.trajectory:
        w.write(ag)

您可以将关键字参数传递给某些格式编写器。例如，LAMMPS 的数据格式允许通过"lengthunit"和"timeunit"这两个关键字来指定输出单位。

序列化

MDAnalysis 能够对其大部分数据结构和轨迹格式进行处理。未被支持的属性可在 PR #2887 中找到。

1
2
3
4
5
6
7
8
import pickle

from MDAnalysis.tests.datafiles import PSF, DCD

psf = mda.Universe(PSF, DCD)

pickle.loads(pickle.dumps(psf))
## <Universe with 3341 atoms>

对于MDAnalysis.core.groups.AtomGroup类，在序列化过程中，它会与关联的MDAnalysis.core.universe.Universe类一起被序列化。这意味着在反序列化后，会创建一个新的MDAnalysis.core.universe.Universe对象，并将其附加到新的MDAnalysis.core.groups.AtomGroup对象上。如果Universe是与MDAnalysis.core.groups.AtomGroup一起进行序列化的，那么之后它们仍然会保持绑定关系：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
import pickle

from MDAnalysis.tests.datafiles import PSF, DCD

u = mda.Universe(PSF, DCD)

g = u.atoms

g_pickled = pickle.loads(pickle.dumps(g))

print("g_pickled.universe is u: ", u is g_pickled.universe)
## g_pickled.universe is u:  False

g_pickled, u_pickled = pickle.loads(pickle.dumps((g, u)))

print("g_pickled.universe is u_pickled: ",
       u_pickled is g_pickled.universe)

## g_pickled.universe is u_pickled:  True

如果多个MDAnalysis.core.groups.AtomGroups被绑定到同一个MDAnalysis.core.universe.Universe对象上，那么在进行序列化操作后，它们也会被绑定到同一个对象上：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
u = mda.Universe(PSF, DCD)
g = u.atoms[:2]
h = u.atoms[2:4]
g_pickled = pickle.loads(pickle.dumps(g))
h_pickled = pickle.loads(pickle.dumps(h))
print("g_pickled.universe is h_pickled.universe : ",
       g_pickled.universe is h_pickled.universe)
## g_pickled.universe is h_pickled.universe :  False

g_pickled, h_pickled = pickle.loads(pickle.dumps((g, h)))

print("g_pickled.universe is h_pickled.universe: ",
       g_pickled.universe is h_pickled.universe)
#g_pickled.universe is h_pickled.universe:  True

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