本文将简要介绍如何使用四元数方法计算两个分子之间RMSD，同时附上简单的示例Python代码。
1. 什么是RMSD
RMSD（Root Mean Square Deviation）是指均方根偏差，在化学中一般用于衡量一个分子结构相对于参照分子的原子偏离位置。RMSD的值越小，说明当前分子结构越接近参照的分子结构。RMSD的数学定义为[1]：
其中是原子的数量。用于比较的分子可以分别看成包含个矢量的集合，，因此RMSD可具体定义为：
从上式也可以看出，当两个分子的几何结构完全一样时，RMSD的值是0。
在量子化学中，xyz文件是一种比较通用的记录分子几何结构的文件格式，其内容如下：
1 原子数量
2 标题
3 原子1 x1 y1 z1
4 原子2 x2 y2 z2
5 原子3 x3 y3 z3
…………
上面每行开头的数字为行数的辅助标记，不在xyz文件中出现。我们的目标是使用四元数方法，写出一个可以计算A、B两个分子之间RMSD值的Python脚本rmsd.py，即在给出两个坐标文件a.xyz和b.xyz后，输入如下命令：
$ ./rmsd.py a.xyz b.xyz
$ rmsd: x.xxxxxx
即可得到A、B分子之间的RMSD值。
2. 基本思路
RMSD的计算公式很简单，主要难点在于怎样将两个分子放在尽可能”相近“的位置上计算。换言之，RMSD会随着两个分子的相对位置变化而变化，我们需要找到RMSD最小的时候对应的相对位置。举个例子，假设我们有两个水分子1.xyz和2.xyz：
1.xyz
3
H2O-1
H -0.78397589 0.44324751 0.00000000
O 0.00000000 -0.11081188 0.00000000
H 0.78397589 0.44324751 0.00000000
2.xyz
3
H2O-2
H -0.81664155 0.46171616 0.00000000
O 0.00000000 -0.11542904 0.00000000
H 0.81664155 0.46171616 0.00000000
上面第二个水分子的键长做了微小调整，根据公式2可得出1.xyz和2.xyz的RMSD值为0.307549。假如我们对第二个水分子做一些平动和转动（为保持直观感受，假设两个分子都固定在xy平面），那么RMSD的值就会发生改变：
除了平动和转动会影响RMSD，原子之间的编号顺序也会产生影响，比如下图：
假设四个灰色原子是同样类型的原子，如果将2、3号原子的编号交换然后和未交换编号的分子做比较，RMSD就会从0变为1.088944。
由此我们可以看出，在计算两个分子RMSD值之前，还至少需要四个步骤：确认两个分子的原子类型和数量相等、优化同类原子的编号顺序、优化分子的平动和优化分子的转动。
3. 什么是四元数
我们首先看看四元数是什么。四元数和复数类似，可以表示成如下形式[2]：
其中i、j、k满足如下乘法表：
式(3)中的称为四元数的实部，实部为零的四元数称为纯四元数。四元数也可以写成矢量的形式，比如
四元数的一些基本运算：
对于纯四元数和，计算可以有一些简化：
注意四元数的乘法不满足交换律，但是满足结合律和分配律。此外，四元数的乘法运算也可以写成矩阵的形式[3]：
上式中矩阵的下表L和R分别表示左乘和右乘。两个矩阵可以表示为[4]：
上面两个矩阵在我们写程序的时候会特别有用。纯四元数可以表示三维空间中的一个矢量，单位纯四元数（也就是）可以用于表示三维空间中的旋转：
是四元数代表的三维转动矩阵。
4. 开始计算RMSD
我们通过xyz文件读取原子坐标的信息，所以可以先写一个简单的xyz文件解析器：
# xyz file parser
def parse_xyz(_xyz):
# check number of atoms
atom_number = 0
# line status
line_number = 0
# buffer
atom_name_buffer = []
atom_coord_buffer = []
with open(_xyz, "r") as xyz_file:
for line in xyz_file:
line_number += 1
# the first line: number of atoms
if (line_number == 1):
if (len(line.split()) == 1):
atom_number = int(line.split()[0])
else:
print("parse xyz file failed: line ", line_number)
exit()
# the second line: ignore title
if (line_number == 2):
pass
# other lines: atom names and coordinates
if (line_number > 2):
if (len(line.split()) == 4):
atom_name_buffer.append(str(line.split()[0]))
line_buffer = [0.0, 0.0, 0.0]
line_buffer[0] = float(line.split()[1])
line_buffer[1] = float(line.split()[2])
line_buffer[2] = float(line.split()[3])
atom_coord_buffer.append(line_buffer)
else:
print("parse xyz file failed: line ", line_number)
exit()
# check
if (not (len(atom_name_buffer) == (len(atom_coord_buffer)) == atom_number)):
print("wrong atom number")
exit()
# change to numpy array after check
atom_name = numpy.array(atom_name_buffer)
atom_coord = numpy.array(atom_coord_buffer)
return atom_name, atom_coord
这个函数的功能是解析形参中的xyz文件，比较第一行原子数目、原子符号数目和坐标数目是否一致，然后将原子符号放入一维矩阵atom_name中、将原子坐标放入的矩阵atom_coord中并返回。RMSD的计算需要比较两个分子，因此在函数调用时应该解析两个xyz文件：
def main():
# 1. arguments check
if (len(sys.argv) != 3):
print("wrong number of arguments")
print("usage: ./rmsd.py [a.xyz] [b.xyz]")
exit()
# 2. parse xyz files
a_xyz_file = sys.argv[1]
b_xyz_file = sys.argv[2]
a_atom, a_coord = parse_xyz(a_xyz_file)
b_atom, b_coord = parse_xyz(b_xyz_file)
if __name__ == "__main__":
main()
为了给原子编号优化和转动矩阵计算做准备，我们首先把两个分子的几何中心全都移动到坐标原点（后面会解释为什么），并且添加函数调用：
def centroid_translation(_atom_coord):
centroid = _atom_coord.mean(axis=0)
_atom_coord -= centroid
return _atom_coord
def main():
# 1. arguments check
......
# 2. parse xyz files
......
# 3. move centroid
a_coord = centroid_translation(a_coord)
b_coord = centroid_translation(b_coord)
现在我们的两个分子的几何中心已经重合，并且都在坐标原点。接下来我们要进行第一个优化步骤，尽可能对齐两个分子的原子编号，也就是纠正第2节中图2的那种编号错位。对齐原子编号可以使用匈牙利算法（Hungarian algorithm），匈牙利算法所解决的问题可以抽象为如下数学模型[5]：假设M个行指标和N列指标可以组成一个矩阵
其中每个矩阵元表示一对横纵坐标之间的“成本函数”，我们通过更改行或列指标的“分配”，使得最后的“总成本”降到最低。从数学上讲我们是计算了这样一个函数：
其中矩阵是指，当行指标和列指标相匹配时，其余矩阵元为零的矩阵。在化学分子的RMSD计算中，横纵指标可以理解为两个分子的原子编号顺序，“成本函数”是两个分子之间对应编号的原子的距离函数，优化“总成本函数”意味着我们要找到一种两个分子间原子编号的对应关系，在这种对应关系下，对应编号原子的距离之和最小。在Python中我们不需要单独实现匈牙利算法，可以借用scipy中的函数进行计算：
def reorder(_a_atom, _a_coord, _b_atom, _b_coord):
# result reorder index container
reorder_indx = numpy.zeros(_b_atom.shape, dtype=int)
# reordered by atom type
unique_atom = numpy.unique(_a_atom)
for atom in unique_atom:
a_atom_indx, = numpy.where(_a_atom == atom)
b_atom_indx, = numpy.where(_b_atom == atom)
# Hungarian assignment, cost matrix is distance matrix
a_atom_coord = _a_coord[a_atom_indx]
b_atom_coord = _b_coord[b_atom_indx]
distance_cost = scipy.spatial.distance.cdist(a_atom_coord, b_atom_coord)
a_indx, b_reorder_indx = scipy.optimize.linear_sum_assignment(
distance_cost)
# unique atom reorder index
reorder_indx[a_atom_indx] = b_atom_indx[b_reorder_indx]
# reorder operation
_b_atom = _b_atom[reorder_indx]
_b_coord = _b_coord[reorder_indx]
return _b_atom, _b_coord
上面代码通过scipy.spatial.distance.cdist函数计算“成本矩阵”，使用scipy.optimize.linear_sum_assignment函数进行指标分配。此外，在上面的计算中，我们是在同类型原子之间进行编号优化，这也很好理解，比如对于甲烷分子，把C原子和H原子进行编号交换是不合理的。
接下来就到了四元数参与的部分了[3]。在读取原子坐标、对齐几何中心、对齐原子编号完成之后，接下来要计算平动-转动矩阵。从数学上讲，我们是要最小化这样一个函数：
其中是转动矩阵，是平动矩阵，和分别表示两个分子中第个原子的坐标矢量。一般来说，为了同时描述刚体在三维空间中的平动和转动，需要使用二元四元数（dual number quaternion）来表示，也就是将四元数中的每一个实数都替换为二元数。但对于化学分子的RMSD计算来说并不需要这么复杂。我们对进行变分来求上述函数的极值（实际上就是对的最小值）：
如果对于两个矢量集合和 ，我们都将几何中心和移动到坐标原点：
原函数将简化为：
所以在RMSD计算中，我们将分子的几何中心移动到原点，就已经完成了对平动矩阵的优化。接下来主要讨论转动矩阵的优化问题，下面我们去掉波浪线分别使用和来代表几何中心在坐标原点时原子坐标的矢量集合（也就是移动后的原子坐标），作用在上的旋转可以写为：
求极小值的函数可以写为：
左右两侧同乘并展开等号右侧，可得：
考虑到单位四元数，且纯四元数，上式前两项可以进行简化：
需要优化的部分还剩下这一项，通过第二节中纯四元数性质代入可得：
上式使用了四元数乘积的结合律，脚标表示其内部结果仅第一个分量不为零。其中可以写为矩阵形式：
考虑到是纯四元数，因此将上式代入原函数，可得：
要取得的最小值，我们就要求二次型的最大值，至此我们已经将RMSD的转动矩阵优化问题变成了本征值问题，接下来我们只需求出矩阵的最大本征值即可。为了方便numpy索引纯四元数矢量的三个非零分量，下面的代码将改写为。首先我们要写出原子坐标生成对应的和矩阵的函数：
# generate A_L matrix
def A_L_matrix(q1, q2, q3, q4=0):
A_L = numpy.array([
[q4, -q3, q2, q1],
[q3, q4, -q1, q2],
[-q2, q1, q4, q3],
[-q1, -q2, -q3, q4],
])
return A_L
# generate A_R matrix
def A_R_matrix(q1, q2, q3, q4=0):
A_R = numpy.array([
[q4, q3, -q2, q1],
[-q3, q4, q1, q2],
[q2, -q1, q4, q3],
[-q1, -q2, -q3, q4],
])
return A_R
得到矩阵之后就可以求最大本征值对应的本征矢量，该本征矢量所表示的的四元数对应的转动矩阵就是我们优化的转动矩阵，在对原子坐标进行旋转操作之后即可计算两个分子之间的RMSD，代码如下：
def quaternion_rotation(_a_coord, _b_coord):
# generate A matrix
N = _a_coord.shape[0]
A = numpy.zeros([4, 4], dtype=float)
for k in range(N):
A_L = A_L_matrix(_a_coord[k, 0], _a_coord[k, 1], _a_coord[k, 2])
A_R = A_R_matrix(_b_coord[k, 0], _b_coord[k, 1], _b_coord[k, 2])
A_Lt_dot_A_R = numpy.dot(A_L.transpose(), A_R)
A += A_Lt_dot_A_R
# generate rotation matrix
eigen = numpy.linalg.eigh(A)
Q = eigen[1][:, eigen[0].argmax()]
A_Rt_rot = A_R_matrix(*Q).transpose()
A_L_rot = A_L_matrix(*Q)
rot = numpy.dot(A_Rt_rot, A_L_rot)[:3, :3]
# rotate coordinates: rotate _a_coord because U(q)x_k from the equation
_a_coord = numpy.dot(_a_coord, rot)
return _a_coord
最后一步，计算坐标变换之后的RMSD，所有步骤就全都完成了：
def rmsd(_a_coord, _b_coord):
diff = numpy.array(_a_coord) - numpy.array(_b_coord)
N = numpy.array(_a_coord).shape[0]
result = numpy.sqrt((diff * diff).sum() / N)
return result
上面提供了计算RMSD的一种方法，适用于分子结构较为简单、RMSD较小的情况。对于复杂分子、结构差别较大或含有手性中心等特殊情况，则需要额外引入其他计算方法[6]。
参考文献
[1] https://en.wikipedia.org/wiki/Root-mean-square_deviation_of_atomic_positions
[2] https://en.wikipedia.org/wiki/Quaternion
[3] Coutsias E. A.; Seok C.; Dill K. A. Using Quaternions to Calculate RMSD. J. Comput. Chem. 2004, 25, 1849-1857.
[4] Walker M. W.; Shao L.; Volz R. A. Estimating 3-D location parameters using dual number quaternions. CVGIP: Image Understanding. 1991, 54, 358-367.
[5] https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.linear_sum_assignment.html
[6] https://github.com/charnley/rmsd（参考了部分代码）