4.8.1. 从轨迹生成密度 MDAnalysis.analysis.density

作者:

奥利弗·贝克斯坦

年:

2011

版权所有:

GNU公共许可证v3

该模块提供类和函数来生成和表示体积数据，特别是密度。

在 2.0.0 版本发生变更: 已弃用 density_from_Universe() ， density_from_PDB() ，以及 Bfactor2RMSF() 现已被移除。

4.8.1.1. 从MD轨迹生成密度

一个常见的用例是分析感兴趣的蛋白质周围的溶剂密度。密度的计算公式为 DensityAnalysis 在模拟单元格的固定坐标系中。因此，有必要相对于盒子坐标系对蛋白质进行定向和固定。在实践中，这意味着将蛋白质一帧接一帧地居中并叠加在参考结构上，并与蛋白质一起平移和旋转模拟的所有其他组件。这样，溶剂就会出现在蛋白质的参照系中。

输入轨迹必须

1. 已经集中在感兴趣的蛋白质上；

2. 所有跨越周期边界的分子都是完整的吗？ [1];

3. 重新映射溶剂分子，使它们最接近溶质(当使用三斜晶胞，如十二面体或截断的八面体时，这一点很重要) [1].

4. 有一个固定的参照系；例如，通过将蛋白质叠加在参照物结构上，这样人们就可以研究它周围的溶剂密度 [2].

要产生蛋白质周围水分子的密度(假设轨迹已经针对周期性边界伪影进行了适当的处理，并适当地叠加以提供固定的参考系) [3] ：：

from MDAnalysis.analysis.density import DensityAnalysis
u = Universe(TPR, XTC)
ow = u.select_atoms("name OW")
D = DensityAnalysis(ow, delta=1.0)
D.run()
D.results.density.convert_density('TIP4P')
D.results.density.export("water.dx", type="double")

所有水中氧的位置( AtomGroup ow )在带有间距的网格上绘制了组织图 德尔塔 =1？最初，密度是以 \(\text{{Å}}^{{-3}}\) 。与 Density.convert_density() 方法，改变测量单位。在该示例中，我们现在测量的是TIP4P水模型在环境条件下相对于文献值的密度(请参阅中的值 MDAnalysis.units.water 有关详细信息，请参见)。最后，密度被写为 OpenDX 可读入的兼容文件 VMD, Chimera, 或 PyMOL.

这个 Density 对象可以作为 DensityAnalysis.results.density 属性。特别是，密度的数据以NumPy数组的形式存储在 Density.grid ，它可以以任何方式处理。

4.8.1.2. 创造密度

这个 DensityAnalysis 类生成一个 Density 来自原子群。

class MDAnalysis.analysis.density.DensityAnalysis(atomgroup, delta=1.0, metadata=None, padding=2.0, gridcenter=None, xdim=None, ydim=None, zdim=None)

体积密度分析。

轨迹是逐帧读取的，其中的原子 atomgroup 在带有间距的3D网格上绘制了柱状图 delta 。

参数:

• atomgroup (AtomGroup or UpdatingAtomGroup) -- 被分析的原子组(如所有的水、氧原子)。这可以是一个 UpdatingAtomGroup 用于每个时间步长都会更改的选择。

• delta (float (optional)) -- 密度栅格的存储箱大小(在x、y、z中相同)。

• padding (float (optional)) -- 在？ngström中通过填充(在初始框大小的顶部)来增加直方图尺寸。设置用户定义的网格时，填充被忽略。

• gridcenter (numpy ndarray, float32 (optional)) -- 3元素NumPy数组，详细说明了？ngström中用户定义的网格框中心的x、y和z坐标。

• xdim (float (optional)) -- ？ngström中用户定义的x尺寸框边缘。

• ydim (float (optional)) -- 用户在？ngström中定义的y尺寸框边缘。

• zdim (float (optional)) -- ？ngström中用户定义的z尺寸框边缘。

results.density

A Density 实例包含单位的物理密度 \(Angstrom^{{-3}}\) 。

分析之后(请参阅 run() 方法)，则结果密度存储在 results.density 属性作为 Density 举个例子。注意：这将取代现在不推荐使用的 density 属性。

类型:

Density

density

的别名 results.density 。

自 2.0.0 版本弃用: 将在MDAnalysis 3.0.0中删除。请使用 results.density 取而代之的是。

类型:

Density

抛出:

• ValueError -- 如果原子组为空且未提供用户定义的网格，或者如果未提供或未正确提供用户定义的网格

• UserWarning -- 如果原子组为空并且提供了用户定义的网格

参见

pmda.density.DensityAnalysis

一个平行版本的 DensityAnalysis

备注

如果 gridcenter 和 x/y/zdim 不提供参数， DensityAnalysis 将尝试从“原子组”中的原子自动生成一个网格框(见示例)。

正常 AtomGroup 实例代表对原子的静态选择。如果你愿意的话 动态更改选择 (如“NAME OW AND约4.0(蛋白质而不是NAME H*)”，即在蛋白质重原子的4？范围内的水氧原子)然后创建 UpdatingAtomGroup (请参见示例)。

DensityAnalysis 将在以下情况下失败 AtomGroup 实例不包含任何原子选择，即使在 updating 设置为 True 。在这种情况下，可以提供用户定义的框限制以生成 Density 。尽管如此，用户仍有责任确保所提供的网格限制包含要在所有轨迹框架上选择的原子。

示例

一个常见的用例是分析感兴趣的蛋白质周围的溶剂密度。密度的计算公式为 DensityAnalysis 在模拟单元格的固定坐标系中。因此，有必要相对于盒子坐标系对蛋白质进行定向和固定。在实践中，这意味着将蛋白质一帧接一帧地居中并叠加在参考结构上，并与蛋白质一起平移和旋转模拟的所有其他组件。这样，溶剂就会出现在蛋白质的参照系中。

输入轨迹必须

1. 已经集中在感兴趣的蛋白质上；

2. 所有跨越周期边界的分子都是完整的吗？ [1];

3. 重新映射溶剂分子，使它们最接近溶质(当使用三斜晶胞，如十二面体或截断的八面体时，这一点很重要) [1];

4. 有一个固定的参照系；例如，通过将蛋白质叠加在参照物结构上，这样人们就可以研究它周围的溶剂密度 [2].

产生密度

要产生蛋白质周围水分子的密度(假设轨迹已经针对周期性边界伪影进行了适当的处理，并适当地叠加以提供固定的参考系) [3], 首先创建 DensityAnalysis 对象，然后使用 run() 方法：

from MDAnalysis.analysis import density
u = Universe(TPR, XTC)
ow = u.select_atoms("name OW")
D = density.DensityAnalysis(ow, delta=1.0)
D.run()
D.results.density.convert_density('TIP4P')

所有水氧化合物的位置都用有间距的网格作了组织图。 德尔塔 =1？并存储为 Density 属性中的对象 DensityAnalysis.results.density 。

与密度一起工作

A Density 包含文档中列出的大量方法和属性。在这里，我们使用 Density.convert_density() 将密度从反立方体转换为相对于标准条件下TIP4P水的体积密度的密度。(MDAnalysis在以下位置存储了许多文献价值 MDAnalysis.units.water 。)

您可以通过以下方式以3D NumPy数组的形式直接访问密度 Density.grid 。

默认情况下， Density 返回的对象包含以？为单位的物理密度 ^-3 。如果您有兴趣恢复潜在的 概率密度 ，只需除以总和：：

probability_density = D.results.density.grid / D.results.density.grid.sum()

同样，如果要恢复包含 原子计数直方图 ，只需乘以体积 dV 每个面元(或体素)的密度；在这种情况下，需要确保物理密度的测量单位为 ^-3 通过转换它：：

import numpy as np

# ensure that the density is A^{-3}
D.results.density.convert_density("A^{-3}")

dV = np.prod(D.results.density.delta)
atom_count_histogram = D.results.density.grid * dV

将密度写入文件

密度可以是 exported to different formats 使用 Density.export() )多亏了这样一个事实 Density 是一个子类 gridData.core.Grid ，它提供了功能)。例如，要 write a DX file water.dx 可以使用VMD、PYMOL或Chimera：：

D.results.density.export("water.dx", type="double")

示例：整个模拟中的水密度

创建水密度的基本用途(只需使用水中的氧原子“OW”)：

D = DensityAnalysis(universe.select_atoms('name OW')).run()

示例：绑定位置中的水(更新选择)

如果您只对某个区域内的水感兴趣，例如，在结合位置附近，您可以使用一个选项，该选项通过使用 UpdatingAtomGroup ：：

near_waters = universe.select_atoms('name OW and around 5 (resid 156 157 305)',
              updating=True)
D_site = DensityAnalysis(near_waters).run()

示例：配体周围的小区域(手动框选择)

如果您对显式设置给定边大小和原点的网格框感兴趣，可以使用 gridcenter 和 xdim /ydim/`zdim`参数。例如，要绘制配体(在这种情况下，配体已被指定为残基名称“LIG”)5？内的水密度，请使用以配体的质心(COM)为中心的具有20？边的立方体网格：：

# Create a selection based on the ligand
ligand_selection = universe.select_atoms("resname LIG")

# Extract the COM of the ligand
ligand_COM = ligand_selection.center_of_mass()

# Create a density of waters on a cubic grid centered on the ligand COM
# In this case, we update the atom selection as shown above.
ligand_waters = universe.select_atoms('name OW and around 5 resname LIG',
                                      updating=True)
D_water = DensityAnalysis(ligand_waters,
                          delta=1.0,
                          gridcenter=ligand_COM,
                          xdim=20, ydim=20, zdim=20)

(应注意的是， padding 在分配用户定义的网格时不使用关键字)。

在 1.0.0 版本加入.

在 2.0.0 版本发生变更: _set_user_grid() 现在是一种 DensityAnalysis 。 Density 结果现在存储在 MDAnalysis.analysis.base.Results 实例。

static _set_user_grid(gridcenter, xdim, ydim, zdim, smin, smax)

用于将格网尺寸设置为用户定义的值的助手函数

参数:

• gridcenter (numpy ndarray, float32) -- 3包含栅格框中心的x、y和z坐标的元素ndarray

• xdim (float) -- X维度中的长方体边缘长度

• ydim (float) -- Y维度上的长方体边缘长度

• zdim (float) -- Y维度上的长方体边缘长度

• smin (numpy ndarray, float32) -- 输入选择的最小x，y，z坐标

• smax (numpy ndarray, float32) -- 输入选择的最大x，y，z坐标

返回:

• umin ( 麻木ndarray，浮动32 )--用户定义栅格的最小x、y、z坐标

• umax ( 麻木ndarray，浮动32 )--用户定义栅格的最大x、y、z坐标

在 2.0.0 版本发生变更: 现在是一种静态的方法 DensityAnalysis 。

run(start=None, stop=None, step=None, frames=None, verbose=None, *, progressbar_kwargs={})

执行计算

参数:

• start (int, optional) -- 分析的起始框架

• stop (int, optional) -- 停止分析框架

• step (int, optional) -- 要在每个分析的帧之间跳过的帧数

• frames (array_like, optional) -- 用于切片轨迹的整数或布尔数组； frames 只能使用 instead 的 start ， stop ，以及 step 。设置 both frames 并且至少有一项 start ， stop ， step 设置为非默认值将引发 ValueError 。。。版本已添加：：2.2.0

• verbose (bool, optional) -- 启用详细信息

• progressbar_kwargs (dict, optional) -- 带有有关进度条位置等的自定义参数的ProgressBar关键字；请参见 MDAnalysis.lib.log.ProgressBar 查看完整的列表。

在 2.2.0 版本发生变更: 添加了分析任意帧的功能，方法是在 frames 关键字参数。

在 2.5.0 版本发生变更: 增列 progressbar_kwargs 参数，允许修改tqdm进度条的描述、位置等

4.8.1.3. 密度对象

密度创建函数的主要输出是 Density 实例，该实例派生自 gridData.core.Grid 。一个 Density 本质上是具有原点和长度的3D数组。

class MDAnalysis.analysis.density.Density(*args, **kwargs)

基类：Grid

类，表示正则笛卡尔网格上的密度。

参数:

• grid (array_like) -- 直方图或密度，通常为 numpy.ndarray

• edges (list) -- 数组列表，沿轴线的下仓边和上仓边

• parameters (dict) -- 类参数词典；与一起保存 Density.save() 。以下关键字对班级有意义。这些值的含义如下所示： IsDensity - False ：GRID是带有计数的直方图 [默认设置] - True ：应用密度 Density.make_density() 将其设置为 ``True` 。

• units (dict) -- 一张有钥匙的字典- 长度 ：栅格边的物理单位(埃或纳米) [阿格斯特罗姆] - 密度 ：密度单位，如果 isDensity=True 或 None 否则，密度的缺省值为“Angstrom^{-3}”(表示 \(\text{{Å}}^{{-3}}\) )。

• metadata (dict) -- 与密度关联的任意值的用户定义字典；类不接触 Density.metadata 但将其存储在 Density.save()

grid

计数或密度

类型:

阵列

edges

中每个单元格的边界 grid 沿所有轴(相当于 numpy.histogramdd() 返回)。

类型:

一维阵列列表

delta

每个维度的单元格大小。

类型:

阵列

origin

的坐标 居中 索引处的单元格的 grid[0, 0, 0, ..., 0] ，它被认为是左下角的前角。

类型:

阵列

units

长度和密度的单位；使用方法更改单位 convert_length() 或 convert_density() 。

类型:

DICT

备注

数据 (Density.grid )可以作为标准的NumPy数组进行操作。可以使用将更改保存到文件 Density.save() 方法。网格可以使用 Density.load() 方法，或者向构造函数提供文件名。

该属性 Density.metadata 包含可用于批注数据的用户定义词典。它还与一起保存 Density.save() 。

这个 Density.export() 方法始终导出3D对象(以这样的方式编写以便于读取 VMD, Chimera, 和 PyMOL) ，其余的应该适用于任何维度的数组。请注意 PyMOL 只理解DX数据类型为“DOUBLE”的DX文件(请参见 known issues when writing OpenDX files 并发布 MDAnalysis/GridDataFormats#35 有关详细信息，请参见)。使用关键字 type="double" 对于该方法， Density.export() ，用户可以确保DX文件以适合于 PyMOL.

如果输入直方图由每个单元格的计数组成，则 Density.make_density() 方法将栅格转换为物理密度。要获得概率密度，请将其除以 Density.grid.sum() 或使用 density=True 马上就到 histogramdd() 。

用户 应该 设置 参数 关键字(有关构造函数，请参阅文档)；尤其是，如果数据已经是密度，则必须设置 isDensity=True 因为没有可靠的方法来检测数据是否代表计数或密度。为了方便起见，如果从文件中读取数据，并且用户尚未设置 isDensity 然后假设数据实际上是一个密度。

参见

gridData.core.Grid

示例

典型用法：

1. 从直方图(即网格上的计数)：

   h,edges = numpy.histogramdd(...)
   D = Density(h, edges, parameters={'isDensity': False}, units={'length': 'A'})
   D.make_density()
   

2. 来自保存的密度文件(例如，OpenDX格式)，其中长度以埃为单位，密度以1/A**3：：

   D = Density("density.dx")
   

3. 根据保存的密度文件(例如，OpenDX格式)，其中长度以埃为单位，密度相对于环境条件下的水密度进行测量：

   D = Density("density.dx", units={'density': 'water'})
   

4. 从已保存的 直方图 (不太常见，但为了证明 参数 关键字)，其中长度以nm：：为单位

   D = Density("counts.dx", parameters={'isDensity': False}, units={'length': 'nm'})
   D.make_density()
   D.convert_length('Angstrom^{-3}')
   D.convert_density('water')
   

   在最后一步之后， D 将包含栅格上的密度(以ngström为单位测量)，密度值本身相对于水的密度进行测量。

Density 可以对对象进行代数操作(加、减、乘等)但是有一些 没有健全的检查 确保单元、元数据等兼容！

备注

建议从直方图构造网格对象，提供适当的长度单位，并使用 Density.make_density() 以获得密度。这确保了长度单位和密度单位彼此对应。

centers()

作为迭代器返回所有网格单元格的中心坐标。

参见

numpy.ndindex()

check_compatible(other)

检查是否 other 可以在算术运算中使用。

other 在以下情况下兼容

1. other 是一个标量

2. other 是在相同的边上定义的格网

为了使 other 兼容，在与此网格相同的网格上重新采样，使用 resample() 。

参数:

other (Grid or float or int) -- 另一个用于标准算术运算的对象，使用此 Grid

抛出:

TypeError -- 如果不兼容

参见

resample()

convert_density(unit='Angstrom')

将密度转换为给出的物理单位 unit 。

参数:

unit (str (optional)) -- 密度应转换为的目标单位。 unit 可以是以下类型之一：=单位的===============================================================名称说明===============================================================angstrom^{-3}颗粒/A 3 nm^{{-3}} particles/nm 3标准条件下SPC水的SPC密度TIP3P。看见 MDAnalysis.units.water TIP4P...看见 MDAnalysis.units.water 标准条件下实际水的水密度(0.997 g/cm**3)摩尔摩尔/l===============================================================

抛出:

• RuntimeError -- 如果密度一开始就没有与之关联的单位(即不是密度)，则不能进行转换。

• ValueError -- 对于未知 unit 。

备注

1. 仅当已存在与密度关联的长度单位时，此方法才有效；否则引发 RuntimeError

2. 转换总是返回到统一状态，因此可以对多个转换进行舍入和浮点处理。

convert_length(unit='Angstrom')

将网格对象转换为新的 unit 。

参数:

unit (str (optional)) -- 栅格应转换为的单位：“埃”、“纳米”之一

备注

这会更改边，但不会更改密度；如果栅格对象是根据密度构建的，则用户有责任提供适当的单位。建议从直方图和长度单位开始，并使用 make_density() 。

default_format = 'DX'

使用导出的默认格式 export() 。

export(filename, file_format=None, type=None, typequote='"')

使用给定格式将密度导出到文件。

也可以从文件名的后缀推导出格式，尽管 file_format 关键字优先。

的默认格式 export() 是‘DX’。使用‘dx’进行可视化。

实施的格式：

DX

OpenDX

泡菜

泡菜(使用 Grid.load() 恢复)； Grid.save() 比这更简单 export(format='python') 。

参数:

• filename (str) -- 输出文件的名称

• file_format ({'dx', 'pickle', None} (optional)) -- 输出文件格式，默认为“DX”

• type (str (optional)) -- 对于DX，设置输出DX数组类型，例如“DOUBLE”或“FLOAT”。默认情况下 (None )时，DX类型由网格数组的NumPy Dtype确定(这通常会导致“DOUBLE”)。。。添加的版本：：0.4.0

• typequote (str (optional)) -- 对于DX，设置用于引用类型字符串的字符；默认情况下，这是一个双引号字符‘“’。来自NAMD-GridForce(MDFF后端)的自定义解析器不需要引号，并且可以使用TypeQUOTE=‘’来安抚它们。.版本已添加：：0.5.0

property interpolated

数据网格(x，y，z)上的B-Spline函数。

这个 interpolated() 函数允许用户获取以下坐标的任意值的数据值：

interpolated([x1,x2,...],[y1,y2,...],[z1,z2,...]) -> F[x1,y1,z1],F[x2,y2,z2],...

插补顺序在中设置 Grid.interpolation_spline_order 。

插值函数只需计算一次即可缓存，以获得更好的性能。什么时候都行 interpolation_spline_order 是被修改的， Grid.interpolated() 是重新计算的。

未知数据的值在中设置 Grid.interpolation_cval (待办事项：也要重新计算何时 interpolation_cval 值已更改。)

示例

重采样的用法示例：

XX, YY, ZZ = numpy.mgrid[40:75:0.5, 96:150:0.5, 20:50:0.5]
FF = interpolated(XX, YY, ZZ)

备注

使用样条线函数对值进行插补。样条线可能会生成通常不会出现在数据中的值。例如，密度是非负的，但三次样条线插值可以生成负值，特别是在0和高值之间的边界处。

在内部，该函数使用 scipy.ndimage.map_coordinates() 使用 mode="constant" 从而通过用相同的常量值填充边缘之外的所有值来确定内插格网之外的内插值，该值由 interpolation_cval 参数，如果未设置该参数，则默认为内插栅格中的最小值。

在 0.6.0 版本发生变更: 格网外的插补现在是使用 mode="constant" 而不是 mode="nearest" ，在栅格之外进行插补时消除了挤出体积。

property interpolation_spline_order

数据的B-样条插值阶数。

3=立方体；还支持4和5

仅选择可接受的值 scipy.ndimage.spline_filter() 好了！

参见

interpolated

load(filename, file_format=None)

从加载保存的格线和边 filename

这个 load() 方法调用类的构造函数方法，并根据加载的数据完全重置所有值。

make_density()

将网格(直方图，单元格中的计数)转换为密度(计数/体积)。

这种方法不可逆转地改变了网格。

要获得概率密度，请手动除以 grid.sum() 。

如果这已经是密度，则发出警告且不执行任何操作，因此调用 make_density 多吃几次也没什么害处。

resample(edges)

将数据重采样到带边的新格网 边 。

此方法创建一个新网格，并将当前网格中的数据重采样为由 edges 。插补的顺序由以下设置 Grid.interpolation_spline_order ：更改值 在此之前 呼叫 resample() 。

参数:

edges (tuple of arrays or Grid) -- 新网格的边缘或 Grid 实例，它提供 Grid.edges

返回:

一个新的 Grid 在新的格网单元上插入数据

返回类型:

Grid

示例

提供 edges (三个数组的元组，指示每个网格单元格的边界)：

g = grid.resample(edges)

作为一种便利，一个也可以提供另一个 Grid 作为此方法的参数：：

g = grid.resample(othergrid)

边缘被取自 Grid.edges 。

resample_factor(factor)

重新采样到新的规则栅格。

参数:

factor (float) -- 格网单元格的数量按比例调整 factor 在每个维度，即， factor * N_i 沿每个维度i的单元格必须为正，并且不能导致沿维度的单元格少于2个。

返回:

内插网格 --重新采样的数据表示在 Grid 使用新的网格单元大小。

返回类型:

Grid

参见

resample

在 0.6.0 版本发生变更: 以前的实现不会改变重新采样的栅格边的范围。其结果是，网格边缘的值将稳步向内爬行。新的实现为每次重采样重新计算网格边的范围。

save(filename)

将栅格对象保存到 filename 并添加“.ickle”扩展名。

在内部，这调用 Grid.export(filename, format="python") 。可以使用以下命令从保存的数据重新生成网格：

g = Grid(filename="grid.pickle")

备注

Pickle格式取决于Python版本，因此不能保证用Python2.7保存的网格也能用Python3.5读取。在可移植性方面，OpenDX格式是更好的选择。

脚注

[1] (1,2,3,4)

使分子变得完整可以通过 MDAnalysis.core.groups.AtomGroup.wrap() 的 Universe.atoms (使用 compound="fragments" )。或中的PBC包装转换 MDAnalysis.transformations.wrap 。

[2] (1,2)

可以用以下形式进行叠加 MDAnalysis.analysis.align.AlignTraj 或中的拟合变换 MDAnalysis.transformations.fit 。

[3] (1,2)

请注意，示例中的轨迹 (XTC )是 not 正确地制作成整体并适合参考结构；省略了这些步骤，以清楚地显示实际密度计算所需的步骤。