4.8.3. 水动力学分析 MDAnalysis.analysis.waterdynamics

作者:

亚历杭德罗·贝尔纳丁

年:

2014-2015年

版权所有:

GNU公共许可证v3

在 0.11.0 版本加入.

该模块提供分析水的动力学轨迹和水与其他分子的相互作用的功能。此模块中的功能包括：水定向松弛(WOR) [Yeh1999] ，氢键寿命(HBL) [Rapaport1983] ，角度分布(AD) [Grigera1996] ，均方位移(MSD) [Bródka1994] 和生存概率(SP) [Liu2004] 。

有关此类分析的更多信息，请参阅 [Araya-Secchi2014] (蛋白质腔中的水)和 [Milischuk2011] (纳米孔中的水)。

参考文献

[Yeh1999] (1,2,3,4,5,6)

Yu-ling Yeh and Chung-Yuan Mou. Orientational relaxation dynamics of liquid water studied by molecular dynamics simulation. The Journal of Physical Chemistry B, 103(18):3699–3705, 1999. doi:10.1021/jp984584r.

[Rapaport1983] (1,2,3)

D.C. Rapaport. Hydrogen bonds in water. Molecular Physics, 50(5):1151–1162, 1983. doi:10.1080/00268978300102931.

[Grigera1996] (1,2,3)

J. Raul Grigera, Susana G. Kalko, and Jorge Fischbarg. Wall−water interface. a molecular dynamics study. Langmuir, 12(1):154–158, 1996. doi:10.1021/la9408681.

[Bródka1994] (1,2,3)

Aleksander Bródka. Diffusion in restricted volume. Molecular Physics, 82(5):1075–1078, 1994. doi:10.1080/00268979400100764.

[Liu2004] (1,2,3)

Pu Liu, Edward Harder, and B. J. Berne. On the calculation of diffusion coefficients in confined fluids and interfaces with an application to the liquid−vapor interface of water. The Journal of Physical Chemistry B, 108(21):6595–6602, 2004. doi:10.1021/jp0375057.

[Milischuk2011] (1,2,3)

Anatoli A. Milischuk and Branka M. Ladanyi. Structure and dynamics of water confined in silica nanopores. The Journal of Chemical Physics, 135(17):174709, 2011. doi:10.1063/1.3657408.

4.8.3.1. 分析类的用法示例

4.8.3.1.1. HydrogenBondLifetimes

要分析氢键寿命，请使用 MDAnalysis.analysis.hydrogenbonds.hbond_analysis.HydrogenBondAnalysis.lifetime() 。

参见

MDAnalysis.analysis.hydrogenbonds.hbond_analysis

4.8.3.1.2. WaterOrientationalRelaxation

分析水的定向松弛(WOR) WaterOrientationalRelaxation 。在这种情况下，我们正在分析水分子旋转/改变方向的速度有多快。如果WOR非常稳定，我们可以假设水分子旋转/改变方向非常慢，另一方面，如果WOR衰减非常快，我们可以假设水分子旋转/改变方向非常快：

import MDAnalysis
from MDAnalysis.analysis.waterdynamics import WaterOrientationalRelaxation as WOR

u = MDAnalysis.Universe(pdb, trajectory)
select = "byres name OH2 and sphzone 6.0 protein and resid 42"
WOR_analysis = WOR(universe, select, 0, 1000, 20)
WOR_analysis.run()
time = 0
#now we print the data ready to plot. The first two columns are WOR_OH vs t plot,
#the second two columns are WOR_HH vs t graph and the third two columns are WOR_dip vs t graph
for WOR_OH, WOR_HH, WOR_dip in WOR_analysis.timeseries:
      print("{time} {WOR_OH} {time} {WOR_HH} {time} {WOR_dip}".format(time=time, WOR_OH=WOR_OH, WOR_HH=WOR_HH,WOR_dip=WOR_dip))
      time += 1

#now, if we want, we can plot our data
plt.figure(1,figsize=(18, 6))

#WOR OH
plt.subplot(131)
plt.xlabel('time')
plt.ylabel('WOR')
plt.title('WOR OH')
plt.plot(range(0,time),[column[0] for column in WOR_analysis.timeseries])

#WOR HH
plt.subplot(132)
plt.xlabel('time')
plt.ylabel('WOR')
plt.title('WOR HH')
plt.plot(range(0,time),[column[1] for column in WOR_analysis.timeseries])

#WOR dip
plt.subplot(133)
plt.xlabel('time')
plt.ylabel('WOR')
plt.title('WOR dip')
plt.plot(range(0,time),[column[2] for column in WOR_analysis.timeseries])

plt.show()

其中t0=0，tf=1000，dtmax=20。这样，我们创建了20个窗口时间步长(x轴上的20个值)，第一个窗口以1000个时间步长平均值(1000/1)创建，第二个窗口以500个时间步长平均值(1000/2)创建，第三个窗口以333个时间步长平均值(1000/3)创建，依此类推。

4.8.3.1.3. AngularDistribution

分析角度分布(AD) AngularDistribution 对于OH向量、HH向量和偶极向量。它返回带有向量方向偏好的线条直方图。输出图中的直线表示水分子没有择优取向。在这种情况下，我们正在分析水分子是否具有某种取向偏好，通过这种方式，我们可以看到水分子是否处于电场下，或者它们是否正在与某种东西(残基、蛋白质等)相互作用：

import MDAnalysis
from MDAnalysis.analysis.waterdynamics import AngularDistribution as AD

u = MDAnalysis.Universe(pdb, trajectory)
selection = "byres name OH2 and sphzone 6.0 (protein and (resid 42 or resid 26) )"
bins = 30
AD_analysis = AD(universe,selection,bins)
AD_analysis.run()
#now we print data ready to graph. The first two columns are P(cos(theta)) vs cos(theta) for OH vector ,
#the seconds two columns are P(cos(theta)) vs cos(theta) for HH vector and thirds two columns
#are P(cos(theta)) vs cos(theta) for dipole vector
for bin in range(bins):
      print("{AD_analysisOH} {AD_analysisHH} {AD_analysisDip}".format(AD_analysis.graph0=AD_analysis.graph[0][bin], AD_analysis.graph1=AD_analysis.graph[1][bin],AD_analysis.graph2=AD_analysis.graph[2][bin]))

#and if we want to graph our results
plt.figure(1,figsize=(18, 6))

#AD OH
plt.subplot(131)
plt.xlabel('cos theta')
plt.ylabel('P(cos theta)')
plt.title('PDF cos theta for OH')
plt.plot([float(column.split()[0]) for column in AD_analysis.graph[0][:-1]],[float(column.split()[1]) for column in AD_analysis.graph[0][:-1]])

#AD HH
plt.subplot(132)
plt.xlabel('cos theta')
plt.ylabel('P(cos theta)')
plt.title('PDF cos theta for HH')
plt.plot([float(column.split()[0]) for column in AD_analysis.graph[1][:-1]],[float(column.split()[1]) for column in AD_analysis.graph[1][:-1]])

#AD dip
plt.subplot(133)
plt.xlabel('cos theta')
plt.ylabel('P(cos theta)')
plt.title('PDF cos theta for dipole')
plt.plot([float(column.split()[0]) for column in AD_analysis.graph[2][:-1]],[float(column.split()[1]) for column in AD_analysis.graph[2][:-1]])

plt.show()

哪里 P(cos(theta)) 是角分布或角概率。

4.8.3.1.4. MeanSquareDisplacement

均方位移分析(MSD) MeanSquareDisplacement 对于水分子来说。在这种情况下，我们正在分析水分子在蛋白质内部沿XYZ方向(半径为11的圆柱形区域)移动的平均距离 [nm] ，Zmax 4.0 [nm] 和Zmin-8.0 [nm] )。强上升意味着水分子的快速运动，弱上升意味着颗粒的缓慢运动：

import MDAnalysis
from MDAnalysis.analysis.waterdynamics import MeanSquareDisplacement as MSD

u = MDAnalysis.Universe(pdb, trajectory)
select = "byres name OH2 and cyzone 11.0 4.0 -8.0 protein"
MSD_analysis = MSD(universe, select, 0, 1000, 20)
MSD_analysis.run()
#now we print data ready to graph. The graph
#represents MSD vs t
time = 0
for msd in MSD_analysis.timeseries:
      print("{time} {msd}".format(time=time, msd=msd))
      time += 1

#Plot
plt.xlabel('time')
plt.ylabel('MSD')
plt.title('MSD')
plt.plot(range(0,time),MSD_analysis.timeseries)
plt.show()

4.8.3.1.5. SurvivalProbability

生存概率(SP)分析 SurvivalProbability 分子的。在这种情况下，我们正在分析水分子在半径为12.3的球体中停留多长时间，该球体以剩余物42和26的几何中心为中心。SP的缓慢衰减意味着水分子在区域中的持久时间较长，而快速衰减意味着较短的持久时间：：

import MDAnalysis
from MDAnalysis.analysis.waterdynamics import SurvivalProbability as SP
import matplotlib.pyplot as plt

universe = MDAnalysis.Universe(pdb, trajectory)
select = "byres name OH2 and sphzone 12.3 (resid 42 or resid 26) "
sp = SP(universe, select, verbose=True)
sp.run(start=0, stop=101, tau_max=20)
tau_timeseries = sp.tau_timeseries
sp_timeseries = sp.sp_timeseries

# print in console
for tau, sp in zip(tau_timeseries, sp_timeseries):
      print("{time} {sp}".format(time=tau, sp=sp))

# plot
plt.xlabel('Time')
plt.ylabel('SP')
plt.title('Survival Probability')
plt.plot(tau_timeseries, sp_timeseries)
plt.show()

需要注意的是， stop 上例中使用的关键字有一个 exclusive 行为，即这里使用的最终帧将是100而不是101。这一行为与 AnalysisBase 但目前不同于其他 MDAnalysis.analysis.waterdynamics 类，这些类都展示了 inclusive 最终帧选择的行为。

另一个例子适用于处理许多不同“残留物”的情况。在这里，我们计算钾离子在膜中每种脂类周围的SP，并对结果进行平均。在这个例子中，如果SP分析在不单独处理每种脂质的情况下运行，钾离子可能从一种脂类跳到另一种脂类，并且仍然被算作残留在指定区域。也就是说，将计算整个膜周围的钾离子的生存概率。

请注意，对于本例，建议使用 Universe(in_memory=True) 要确保模拟不会重新加载到内存中，请执行以下操作：

import MDAnalysis as mda
from MDAnalysis.analysis.waterdynamics import SurvivalProbability as SP
import numpy as np

u = mda.Universe("md.gro", "md100ns.xtc", in_memory=True)
lipids = u.select_atoms('resname LIPIDS')
joined_sp_timeseries = [[] for _ in range(20)]
for lipid in lipids.residues:
    print("Lipid ID: %d" % lipid.resid)

    select = "resname POTASSIUM and around 3.5 (resid %d and name O13 O14) " % lipid.resid
    sp = SP(u, select, verbose=True)
    sp.run(tau_max=20)

    # Raw SP points for each tau:
    for sps, new_sps in zip(joined_sp_timeseries, sp.sp_timeseries_data):
        sps.extend(new_sps)

# average all SP datapoints
sp_data = [np.mean(sp) for sp in joined_sp_timeseries]

for tau, sp in zip(range(1, tau_max + 1), sp_data):
    print("{time} {sp}".format(time=tau, sp=sp))

4.8.3.2. 输出

4.8.3.2.1. WaterOrientationalRelaxation

按窗口时间步长返回水定向松弛(WOR)数据，该数据存储在 WaterOrientationalRelaxation.timeseries ：：

results = [
    [ # time t0
        <WOR_OH>, <WOR_HH>, <WOR_dip>
    ],
    [ # time t1
        <WOR_OH>, <WOR_HH>, <WOR_dip>
    ],
    ...
 ]

4.8.3.2.2. AngularDistribution

每个向量返回角度分布(AD)数据，该数据存储在 AngularDistribution.graph 。事实上，AngularDistribution返回一个直方图：：

results = [
    [ # OH vector values
      # the values are order in this way: <x_axis  y_axis>
        <cos_theta0 ang_distr0>, <cos_theta1 ang_distr1>, ...
    ],
    [ # HH vector values
        <cos_theta0 ang_distr0>, <cos_theta1 ang_distr1>, ...
    ],
    [ # dip vector values
       <cos_theta0 ang_distr0>, <cos_theta1 ang_distr1>, ...
    ],
 ]

4.8.3.2.3. MeanSquareDisplacement

均方位移(MSD)数据在列表中返回，每个元素在其各自的窗口时间步长中表示一个MSD值。数据存储在 MeanSquareDisplacement.timeseries ：：

results = [
     #MSD values orders by window timestep
        <MSD_t0>, <MSD_t1>, ...
 ]

4.8.3.2.4. SurvivalProbability

生存概率(SP)计算两个列表：一个taus列表 (SurvivalProbability.tau_timeseries )和相应生存概率的列表 (SurvivalProbability.sp_timeseries )。

结果= [ tau1, tau2, ..., tau_n ] ， [ sp_tau1, sp_tau2, ..., sp_tau_n]

此外，还有一份名单 SurvivalProbability.sp_timeseries_data ，其中包含为每个tau计算的所有SP的列表。这可以用来计算SP的分布或时间相关性等。

4.8.3.3. 班级

class MDAnalysis.analysis.waterdynamics.WaterOrientationalRelaxation(universe, select, t0, tf, dtmax, nproc=1)

水位向松弛分析

评价叶玉玲、牟仲远提出的水位向松弛的作用 [Yeh1999] 。水取向松弛表示水分子旋转或改变方向的速度有多快。这是由以下公式给出的时间相关函数：

\[C_{\hat u}(\tau)=\lAngel\mathit{P}_2 [\mathbf{{\hat{{u}}}}(t_0)\cdot\mathbf{{\hat{{u}}}}(t_0+\tau)] \范围\]

哪里 \(P_2=(3x^2-1)/2\) 是二阶勒让德多项式，并且 \(\hat{{u}}\) 是沿HH、OH或偶极向量的单位向量。

参数:

• universe (Universe) -- 宇宙对象

• selection (str) -- 选水管柱 [‘BYRES NAME OH2’] 。

• t0 (int) -- 分析开始的帧

• tf (int) -- 分析结束处的框架

• dtmax (int) -- 最大DT大小， dtmax < tf 否则它就会崩溃。

在 0.11.0 版本加入.

在 1.0.0 版本发生变更: 变化 selection 关键字至 select

static lg2(x)

第二类勒让德多项式

run(**kwargs)

分析轨迹并制作时间序列

class MDAnalysis.analysis.waterdynamics.AngularDistribution(universe, select, bins=40, nproc=1, axis='z')

角分布函数分析

角分布函数(AD)被定义为余弦的分布概率 \(\theta\) 水分子的OH向量、HH向量或偶极向量与向量形成的角度 \(\hat n\) 平行于选定的轴(z为默认值)。余弦定义为 \(\cos \theta = \hat u \cdot \hat n\) ，在哪里 \(\hat u\) 是OH、HH或偶极向量。它创建直方图并返回列表列表，请参见 Output. AD也称为角概率(AP)。

参数:

• universe (Universe) -- 宇宙对象

• select (str) -- 用于计算其角度分布的选择字符串 [‘BYRES NAME OH2’]

• bins (int (optional)) -- 要通过以下方式创建直方图的仓数 numpy.histogram()

• axis ({'x', 'y', 'z'} (optional)) -- 用向量(HH、OH或偶极)创建角度并计算余弦角 ['z'] 。

在 0.11.0 版本加入.

在 1.0.0 版本发生变更: 变化 selection 关键字至 select

run(**kwargs)

求cos(Theta)角分布的函数

class MDAnalysis.analysis.waterdynamics.MeanSquareDisplacement(universe, select, t0, tf, dtmax, nproc=1)

均方位移分析

求均方位移的函数 (MSD) 。MSD给出了粒子移动的平均距离。MSD由以下人员提供：

\[\lAngel\增量r(T)^2\rAngel=2nDt\]

哪里 \(r(t)\) 是粒子在时间上的位置 \(t\) ， \(\Delta r(t)\) 是滞后后的位移吗？ \(t\) ， \(n\) 是维度，在本例中 \(n=3\) ， \(D\) 是扩散系数和 \(t\) 是时候了。

参数:

• universe (Universe) -- 宇宙对象

• select (str) -- 选水管柱 [‘BYRES NAME OH2’] 。

• t0 (int) -- 分析开始的帧

• tf (int) -- 分析结束处的框架

• dtmax (int) -- 最大DT大小， dtmax < tf 否则它就会崩溃。

在 0.11.0 版本加入.

在 1.0.0 版本发生变更: 变化 selection 关键字至 select

run(**kwargs)

分析轨迹并制作时间序列

class MDAnalysis.analysis.waterdynamics.SurvivalProbability(universe, select, verbose=False)

生存概率(SP)给出了一组粒子停留在某个区域的概率。SP由以下人员提供：

\[P(\tau) = \langle \frac{ N(t, t + \tau )} { N(t) }\rangle\]

哪里 \(\tau\) 是时间步长， \(N(t)\) 一次粒子的数量 \(t\) ，以及 \(N(t, t+\tau)\) 是每一帧的粒子数 \(t\) 至 \(t + \tau\) 。尖括号表示所有时间起点的平均值， \(t\) 。看见 MDAnalysis.lib.correlations.autocorrelation() 以获取技术详细信息。

参数:

• universe (Universe) -- 宇宙对象

• select (str) -- 选择字符串；允许任何选择。通过此选择，您可以定义要分析的区域/区域，例如：“RENAME SOL AND AUBLE 5(RESID 10)”。看见 SP-examples 。

• verbose (Boolean, optional) -- 如果为True，则将进度和备注打印到控制台。

备注

目前 SurvivalProbability 是唯一上演的 MDAnalysis.analysis.waterdynamics 要支持 exclusive 行为(即类似于 AnalysisBase 发送到 stop 关键字传递给 SurvivalProbability.run() 。不同于其他 MDAnalysis.analysis.waterdynamics 最终帧定义是 inclusive 。

在 0.11.0 版本加入.

在 1.0.0 版本发生变更: 使用MDAnalysis.lib.correlations.py执行间歇性和自相关计算。变化 selection 关键字至 select 。删除了对不推荐使用的 t0 ， tf ，以及 dtmax 关键字。相反，这些应该传递给 SurvivalProbability.run() 作为 start ， stop ，以及 tau_max 关键字分别为。这个 stop 传递给的关键字 SurvivalProbability.run() 现在已经改变了行为，并将在 exclusive 方式(而不是以前的方式 inclusive 行为)，

在 2.7.0 版本发生变更: 更新文档以与离散自相关函数对齐。

run(tau_max=20, start=None, stop=None, step=None, residues=False, intermittency=0, verbose=False)

计算并返回生存概率(SP)时间序列

参数:

• start (int, optional) -- 要分析的第一帧的从零开始的索引，默认：无(第一帧)。

• stop (int, optional) -- 要分析的最后一帧的从零开始的索引(独占)，默认：无(最后一帧)。

• step (int, optional) -- 每跳一次 step -第1帧。这是兼容的，但独立于所使用的taus，最好考虑使用 step 等于 tau_max 以消除重叠部分。请注意 step 和 tau_max 时断时续地持续工作。默认：无(使用每帧)。

• tau_max (int, optional) -- 生存概率的计算范围为1<= tau <= tau_max 。

• residues (Boolean, optional) -- 如果为真，则将对残留物(.resds)而不是原子(.ids)进行分析。一个原子就足以将残基归类为距离内的残基。

• intermittency (int, optional) -- 原子可以离开但如果它在下一帧返回时被视为存在的连续帧的最大数目。间歇性的 0 相当于一个连续的生存概率，这不允许原子的离开和返回。例如，对于 intermittency=2 任何给定的原子可以离开感兴趣的区域多达两个连续的帧，但被视为在所有帧都存在。默认值为连续(0)。

• verbose (Boolean, optional) -- 将进度打印到控制台。

返回:

• tau_timeseries ( list )--tau从1到 tau_max 。保存在tau_TimeSeries字段中。

• sp_timeseries ( list )--每个值的生存概率 tau 。保存在字段sp_timeseries中。

• sp_timeseries_data ( list )--从中获取平均值的原始数据点(Sp_Timeseries)。可以提取时间相关性和分布。

在 1.0.0 版本发生变更: 要计算其他分析方法， stop 关键字现在是独占的，而不是包含的。