python - 在笛卡尔坐标中对齐两个 3D 对象

Question

我有两个与 .xyz 文件相同的分子的副本。这意味着每个原子都有 X、Y 和 Z 坐标。但是，您可以旋转分子并获得每个原子的不同坐标，尽管相对位置相同并且分子保持不变。我想使用三个原子作为参考点来对齐两个分子。然而，我正在努力完全对齐这两个分子。

首先，我通过翻译单个原子来对齐两个分子。然后，我使用旋转矩阵进行两次后续旋转，如别处所述。出于某种原因，我需要取两个向量的叉积的负数，并使用正弦而不是余弦来使两个结构完美对齐（经过大量试验和错误后我发现了这一点）。

对于第二次旋转，我将两个要对齐的向量投影在由旋转向量定义的平面上。这是必要的，因为我不想沿着两个向量的叉积旋转来对齐，因为这会使分子的其余部分不对齐。相反，我沿着两个已经对齐的向量旋转。该项目允许我找到两个向量之间的平面角度，从而找到必要的旋转。

但是，此代码没有正确对齐两个分子。

"group1[0]" 包含要在列表中对齐的三个原子的 XYZ 坐标。同样对于“group2 [0]”和结构 2。

#Point 1: align the functional groups to the origin

O1 = np.array(coords1[group1[0][0]])
O2 = np.array(coords2[group2[0][0]])

mat_2 = np.zeros((len(atoms2), 3))

for ind, c in enumerate(coords1):
    coords1[ind] = np.array(c) - O1

for ind, c in enumerate(coords2):
    coords2[ind] = np.array(c) - O2
    mat_2[ind] = coords2[ind]

#Point 2: align according to a first vector
v1 = np.array(coords1[group1[0][1]])#Since atom 1 is the origin, the coordinates is the vector already
v2 = np.array(coords2[group2[0][1]])#Since atom 1 is the origin, the coordinates is the vector already

v1 = v1/np.linalg.norm(v1)
v2 = v2/np.linalg.norm(v2)

#Let v be the axis of rotation
v = -np.cross(v1, v2)#why do I need a minus here?

if np.linalg.norm(v) != 0:
    a = np.arccos(np.dot(v1, v2)/(np.linalg.norm(v1)*np.linalg.norm(v2)))
    #c = np.dot(v1, v2)*np.cos(a)

    c = np.dot(v1, v2)*np.sin(a)#The internet says cos, but this works perfectly

    vx = np.array([[0, -v[2], v[1]], [v[2], 0, -v[0]], [-v[1], v[0], 0]])

    rot_mat = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) + vx + vx.dot(vx)*(1-c)/(1-c**2)
    mat_2 = np.array(mat_2)

    R_mat_rot = np.matmul(rot_mat, mat_2.T).T
else:
    exit(0) 

coords3 = R_mat_rot.copy()

#I get exactly what I want up until here

#Point 3: Rotate along atom2-atom1 (v1) to align the third atom
v = -v1.copy()

v2 = np.array(coords3[group2[0][2]]) - np.array(coords3[group2[0][0]]) #Since atom 1 is the origin, the coordinates is the vector already
v2 = v2/np.linalg.norm(v2)

v1 = np.array(coords1[group1[0][2]]) - np.array(coords1[group1[0][0]]) #Since atom 1 is the origin, the coordinates is the vector already
v1 = v1/np.linalg.norm(v1)

if np.linalg.norm(v) != 0:
    #consider v to be the vector normal to a plane
    #we want the projection of v1 and v2 unto that plane
    vp1 = np.cross(v, np.cross(v1, v)) - np.array(coords1[group1[0][0]])
    vp1 = vp1/np.linalg.norm(vp1)

    vp2 = np.cross(v, np.cross(v2, v)) - np.array(coords3[group2[0][0]])
    vp2 = vp2/np.linalg.norm(vp2)

    #we find the angle between those vectors on the plane
    a = np.arccos(np.dot(vp1, vp2))/(np.linalg.norm(vp1)*np.linalg.norm(vp2))

    #rotation of that amount
    c = np.dot(v1, v2)*np.cos(a)
    vx = np.array([[0, -v[2], v[1]], [v[2], 0, -v[0]], [-v[1], v[0], 0]])

    rot_mat = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) + vx + np.dot(vx, vx)*(1-c)/(1-c**2)

    R_mat_rot = np.matmul(rot_mat, coords3.T).T

coords4 = R_mat_rot.copy()#Final coordinates

标签： python