利用四元数进行蛋白质原子坐标旋转变换

四元数简介

四元数（Quaternion）是一种扩展复数的数学工具，用于描述三维空间中的旋转操作。它比传统的旋转矩阵具有一些优点，比如避免**万向节锁（Gimbal Lock）**问题，并且在插值操作中更加稳定。四元数常用于计算机图形学、机器人学、航空航天等领域来表示物体的旋转。

一个四元数 q 由一个实数部分和一个虚数部分组成，可以表示为：

q=w+xi+yj+zk

其中，w 是实数部分，x、y、z 是虚数部分的系数。

四元数用于表示旋转时通常写成：

q=[w,x,y,z]

其中 w=cos⁡(θ/2)，(x,y,z)是旋转轴的单位向量乘以 sin(θ/2)，并且 θ 是旋转角度。

四元数旋转的优点

• 不容易出现奇点：相比欧拉角的旋转矩阵，四元数在空间旋转中不会出现万向节锁。
• 效率高：四元数的计算比旋转矩阵的乘法更高效。
• 插值方便：四元数适合用在旋转插值（如球形线性插值，Slerp）等场合。

四元数旋转公式

若我们有一个向量 vv，要通过四元数 qq 旋转它，过程如下：

1. 将向量 v 转换成一个四元数 vq，其形式为 vq=[0,vx,vy,vz]。
2. 计算 v′=q⋅vq⋅q−1，其中 q−1 是 q 的共轭，即 q−1=[w,−x,−y,−z]。
3. 结果 v′就是经过旋转后的向量。

Python代码：利用四元数旋转蛋白质原子坐标

import numpy as np
 
def quaternion_multiply(q1, q2):
    """计算两个四元数的乘积"""
    w1, x1, y1, z1 = q1
    w2, x2, y2, z2 = q2
    w = w1*w2 - x1*x2 - y1*y2 - z1*z2
    x = w1*x2 + x1*w2 + y1*z2 - z1*y2
    y = w1*y2 - x1*z2 + y1*w2 + z1*x2
    z = w1*z2 + x1*y2 - y1*x2 + z1*w2
    return np.array([w, x, y, z])
 
def quaternion_conjugate(q):
    """计算四元数的共轭"""
    w, x, y, z = q
    return np.array([w, -x, -y, -z])
 
def rotate_vector_by_quaternion(v, q):
    """通过四元数旋转向量 v"""
    v_q = np.array([0] + list(v))  # 将向量转换为四元数 [0, vx, vy, vz]
    q_conj = quaternion_conjugate(q)  # 计算 q 的共轭
    # 执行旋转: v' = q * v * q_conj
    v_rotated = quaternion_multiply(quaternion_multiply(q, v_q), q_conj)
    return v_rotated[1:]  # 只返回旋转后的向量部分
def rotate_coordinates_by_quaternion(coords, axis, angle):
    """
    使用四元数旋转一组原子坐标
    coords: (N, 3) numpy 数组，表示 N 个原子的坐标
    axis: (3,) numpy 数组，表示旋转轴
    angle: 旋转角度，单位为弧度
    """
    axis = axis / np.linalg.norm(axis)  # 归一化旋转轴
    w = np.cos(angle / 2)
    x, y, z = axis * np.sin(angle / 2)
    q = np.array([w, x, y, z])  # 构造四元数
    
    # 对每个原子坐标应用旋转
    rotated_coords = np.array([rotate_vector_by_quaternion(coord, q) for coord in coords])
    return rotated_coords
# 示例原子坐标 (单位：Å)
coords = np.array([[1.0, 0.0, 0.0],
                   [0.0, 1.0, 0.0],
                   [0.0, 0.0, 1.0]])
# 旋转轴（绕z轴旋转）
axis = np.array([0, 0, 1])
# 旋转角度（90度，单位：弧度）
angle = np.pi / 2
# 执行旋转
rotated_coords = rotate_coordinates_by_quaternion(coords, axis, angle)
print("旋转后的坐标：")
print(rotated_coords)

代码说明

• 四元数的乘法：quaternion_multiply 函数实现了两个四元数的乘法运算，计算两个四元数的乘积。
• 四元数的共轭：quaternion_conjugate 函数用于计算四元数的共轭，旋转时需要用到四元数的共轭来反向旋转。
• 旋转操作：rotate_vector_by_quaternion 使用四元数来旋转一个向量（原子坐标）。我们通过四元数和其共轭的乘法来完成旋转。
• 主函数：rotate_coordinates_by_quaternion 用于对多个原子坐标应用四元数旋转。

总结

• 四元数是用于描述旋转的工具，它比传统的旋转矩阵更加灵活且不易出现奇点问题。
• 四元数旋转的主要步骤是通过四元数的乘法与共轭来完成向量的旋转。
• 上述代码展示了如何用四元数对蛋白质的原子坐标进行旋转，可以扩展应用于三维分子结构的各种几何变换。