1

我正在尝试使用 MDAnalysis ( MDAnalysis.__version__ == 0.17.0) API 函数principal_axes()并moment_of_inertia()为一组选定的原子计算这些矩阵，如文档中所述：

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不是文档中提到的对角线： 在此处输入图像描述

也许我需要将蛋白质与质心对齐，以计算相对于它的惯性矩。

4

2 回答 2

1

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

该AtomGroup.principal_axes()方法返回一个数组，[p1, p2, p3]其中主轴p1,是长度为 3 的数组；为方便起见，选择了这种作为行向量的布局（以便可以使用 提取向量）。为了形成一个矩阵，其中主轴是列向量，就像在主轴的通常处理中一样，必须进行转置。例如：p2p3p1, 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了主轴作为列向量。

于 2018-03-14T00:59:34.767 回答
1

问题是在他们说的文档中U'IU，但是 U 是返回值的转置CA.principal_axes()（参见源代码）：

    # Sort
    indices = np.argsort(e_val)[::-1]
    # Return transposed in more logical form. See Issue 33.
    return e_vec[:, indices].T

Matlab 确认：

>> I=[ 393842.2070687     -963.01376005   -6050.68541811 ;  -963.01376005  474434.9289629    -3902.61617054;  -6050.68541811   -3902.61617054  520207.91703069];
>> U=[-0.04680878 -0.08278738  0.99546732; 0.01813292 -0.9964659  -0.08201778;-0.99873927 -0.01421157 -0.04814453];
>> U*I*U'

ans =

   1.0e+05 *

    5.2082    0.0000   -0.0000
    0.0000    4.7413   -0.0000
   -0.0000   -0.0000    3.9354
于 2018-03-13T05:51:31.730 回答