第3章：向量、矩阵、旋转与平移

各位同学，咱们今天不讲枯燥的数学定义。我直接说——你手里拿到的蛋白质坐标，本质上就是一堆向量。你要计算 RMSD，就是在比较两堆向量之间的差异。所以，这一章我们把最核心的线性代数工具过一遍。不用怕，都是高中和大学低年级的内容，我帮你重新组织一下，让它跟 RMSD 直接挂钩。

3.1 从原子坐标到向量空间

先看最基础的概念。一个原子的三维坐标 (x, y, z)，在数学上就是一个 三维向量。一整条蛋白质的 N 个原子，就是 N 个这样的向量。

RMSD 的本质是：

• 对每个原子，算它在参考结构和目标结构中的位置差。
• 这个差值也是个向量（差向量）。
• 把所有差向量的长度平方加起来，取平均，再开方。

用数学语言写出来：

对于第 i 个原子：
  差向量 d_i = p_i - q_i
  其中 p_i 是参考结构的坐标，q_i 是目标结构的坐标
  向量长度平方 ||d_i||² = (dx)² + (dy)² + (dz)²

RMSD = sqrt( (1/N) * Σ ||d_i||² )

看到了吗？整个 RMSD 就是在处理向量。所以，对向量有感觉，比死记公式重要得多。

3.2 平移——先让两个结构碰在一起

你想想看，两个结构可能放在空间中不同的位置。如果不先对齐，计算出来的 RMSD 值会非常大，但它反映的不是构型差异，而是整体位置偏移。这就没意义了。

解决办法很简单——把两个结构的质心移动到同一个点。一般做法是都移到原点。

步骤：

1. 计算每个结构的质心（所有原子坐标的平均值）。
2. 每个原子的坐标都减去质心坐标。
3. 这样两个结构就“悬浮”在原点附近，只比较形状差异。

这就是平移。在数学上，平移就是向量加法。

3.3 旋转——用正交矩阵让形状重合

平移对齐质心后，两个结构可能方向还不一致。这时候就需要旋转。

在三维空间中，旋转可以用一个 3×3 的正交矩阵 R 来表示。对一个向量 v 做旋转，就是矩阵乘法 v' = R · v。

旋转矩阵的特点是：

• 行列式为 +1（保定向旋转）。
• 逆矩阵等于转置矩阵：R⁻¹ = Rᵀ。
• 旋转不改变向量的长度：||R·v|| = ||v||。

在 RMSD 计算中，我们做的正是：找到一个旋转矩阵，让目标结构旋转后，与参考结构尽可能重合。这就是 Kabsch 算法 的核心。

绕 x 轴旋转 θ 角：
  R_x = [[1, 0,     0     ],
         [0, cosθ, -sinθ ],
         [0, sinθ,  cosθ  ]]

绕 y 轴旋转 θ 角：
  R_y = [[ cosθ, 0, sinθ ],
         [ 0,    1, 0    ],
         [-sinθ, 0, cosθ ]]

绕 z 轴旋转 θ 角：
  R_z = [[ cosθ, -sinθ, 0 ],
         [ sinθ,  cosθ, 0 ],
         [ 0,     0,    1 ]]

3.4 矩阵——RMSD 优化中的核心运算

好，现在把向量和矩阵串起来。RMSD 的最优化问题，可以抽象成这么几步：

1. 用矩阵 P 表示参考结构（N×3 矩阵）。
2. 用矩阵 Q 表示目标结构（N×3 矩阵）。
3. 找到旋转矩阵 R 和平移向量 t，使得 min ||P - (R·Q + t)||。

这个 ||·|| 就是 Frobenius 范数，说白了就是把所有对应原子距离的平方加起来再开方。而求解这个问题的标准方法，就是我们后面会详细讲的 Kabsch 算法——它本质上是做奇异值分解（SVD）。

我整理了一个对比表，帮你理清概念：

3.5 整个计算流程一图看明白

我把上面讲的串成一张流程图。你跟着箭头走，就是完整的 RMSD 计算逻辑：

3.6 实战小代码——验证你的理解

光说不练假把式。我习惯在正式写 Kabsch 之前，先试一段快速验证代码，确认自己没搞错向量和矩阵的基本操作。

import numpy as np

# 构造两个简单的测试结构（3个原子）
P = np.array([[0, 0, 0], [1, 0, 0], [0, 1, 0]], dtype=float)
Q = np.array([[1, 1, 1], [2, 1, 1], [1, 2, 1]], dtype=float)

# 1. 计算质心
centroid_P = np.mean(P, axis=0)
centroid_Q = np.mean(Q, axis=0)
print(f"质心 P: {centroid_P}, 质心 Q: {centroid_Q}")

# 2. 质心对齐
P_centered = P - centroid_P
Q_centered = Q - centroid_Q
print(f"对齐后 P 第一原子: {P_centered[0]}")
print(f"对齐后 Q 第一原子: {Q_centered[0]}")

# 3. 不做旋转，直接算 RMSD（你猜值是多少？）
diff = P_centered - Q_centered
rmsd_naive = np.sqrt(np.mean(np.sum(diff**2, axis=1)))
print(f"未旋转的 RMSD: {rmsd_naive:.4f}")

# 4. 构建协方差矩阵
H = P_centered.T @ Q_centered
print(f"协方差矩阵 H:\n{H}")

# 5. 做 SVD
U, S, Vt = np.linalg.svd(H)
R = Vt.T @ U.T
print(f"最优旋转矩阵 R:\n{R}")

# 6. 应用旋转后计算 RMSD
Q_rotated = (R @ Q_centered.T).T
diff_rotated = P_centered - Q_rotated
rmsd_final = np.sqrt(np.mean(np.sum(diff_rotated**2, axis=1)))
print(f"旋转对齐后的 RMSD: {rmsd_final:.4f}")

你运行一下这段代码，就能直观感受到：平移把结构拉到原点，旋转让形状最佳重合，最后计算向量差的长度平方，就是 RMSD。每一步都离不开向量和矩阵。