python - 如何使用 MDAnalysis 分析一组原子的principal_axes 和 moment_of_inertia?

标签 python python-3.x physics chemistry mdanalysis

我正在尝试使用 MDAnalysis (MDAnalysis.__version__ == 0.17.0) API 函数 principal_axes()moment_of_inertia() 来计算 这些矩阵用于一组选定的原子,如doc中所述:

import MDAnalysis
from MDAnalysis.tests.datafiles import PSF, DCD
import numpy as np

u = MDAnalysis.Universe(PSF, DCD)

CA = u.select_atoms("protein and name CA")

I = np.matrix(CA.moment_of_inertia())
U = np.matrix(CA.principal_axes())
print("center of mass", CA.center_of_mass())
print("moment of inertia", I)
print("principal axes", U)
print("Lambda = U'IU", np.transpose(U)*I*U)

输出:

center of mass [ 0.06873595 -0.04605918 -0.24643682]
moment of inertia [[ 393842.2070687     -963.01376005   -6050.68541811]
 [   -963.01376005  474434.9289629    -3902.61617054]
 [  -6050.68541811   -3902.61617054  520207.91703069]]
principal axes [[-0.04680878 -0.08278738  0.99546732]
 [ 0.01813292 -0.9964659  -0.08201778]
 [-0.99873927 -0.01421157 -0.04814453]]
Lambda = U'IU [[ 519493.24344558   -4093.3268841    11620.96444297]
 [  -4093.3268841   473608.1536763     7491.56715845]
 [  11620.96444297    7491.56715845  395383.6559404 ]]

这看起来不对,原因之一是 U'IU 不是 doc 中提到的对角线。 : enter image description here

<小时/>

也许我需要将蛋白质与质心对齐以计算与之相关的转动惯量。

最佳答案

教程中的文档AtomGroup.principal_axes()原则上是正确的,但令人困惑的是 AtomGroup.principal_axes() 的返回值不是矩阵 U 而是其转置 U.T:

AtomGroup.principal_axes() 方法返回一个数组 [p1, p2, p3],其中主轴 p1, p2p3 是长度为 3 的数组;选择这种布局作为向量是为了方便(这样就可以使用p1, p2, p3 = ag.principal_axes()提取向量)。为了形成一个矩阵U,其中主轴是列向量,就像主轴的通常处理一样,必须进行转置。例如:

import MDAnalysis
from MDAnalysis.tests.datafiles import PSF, DCD
import numpy as np

u = MDAnalysis.Universe(PSF, DCD)

CA = u.select_atoms("protein and name CA")

I = CA.moment_of_inertia()
UT = CA.principal_axes()

# transpose the row-vector layout UT = [p1, p2, p3]
U = UT.T

# test that U diagonalizes I
Lambda = U.T.dot(I.dot(U))

print(Lambda)

# check that it is diagonal (to machine precision)
print(np.allclose(Lambda - np.diag(np.diagonal(Lambda)), 0))

矩阵Lambda应该是对角线的(最后一个print应该显示True):

[[ 5.20816990e+05 -6.56706349e-10 -2.83491351e-12]
[-6.62283524e-10  4.74131234e+05 -2.06979926e-11]
[-6.56687024e-12 -2.07159142e-11  3.93536829e+05]]
True

最后,如果您想“手动”计算:

values, evecs = np.linalg.eigh(I)
indices = np.argsort(values)
U = evecs[:, indices]

这给出了 U 的主轴作为列向量。

关于python - 如何使用 MDAnalysis 分析一组原子的principal_axes 和 moment_of_inertia?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/49239475/

相关文章:

python - 使用CoLab在Google云端硬盘上OpenCV读取和写入速度太慢

python - 如何在 PySpark 中将列从字符串转换为数组

python - 删除 pandas where 中的 "duplicate rows"并附加条件

python-3.x - 如何在 python 中使用 while 循环来添加列表中的值,直到它们超过最大值?

c++ - 在 Box2D 中制作漩涡

python - 给定 xy 触摸位置,如何用颜色填充 RGBA 图像的透明部分?

python - 创建一个图形,其中包含按 gridspec 布局排列的绘图、 slider 和其他小部件

algorithm - 比较海量数据的最佳算法

c# - 如何制作类似驱散的效果?

swift - 更改 ShowPhysics 轮廓颜色