使用四元数计算两个分子之间的RMSD（附Python代码）
使用四元数计算两个分子之间的RMSD（附Python代码）
本文将简要介绍如何使用四元数方法计算两个分子之间RMSD，同时附上简单的示例Python代码。
1. 什么是RMSD
RMSD（Root Mean Square Deviation）是指均方根偏差，在化学中一般用于衡量一个分子结构相对于参照分子的原子偏离位置。RMSD的值越小，说明当前分子结构越接近参照的分子结构。RMSD的数学定义为[1]：
从上式也可以看出，当两个分子的几何结构完全一样时，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]：
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
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列指标可以组成一个矩阵
其中每个矩阵元表示一对横纵坐标之间的“成本函数”，我们通过更改行或列指标的“分配”，使得最后的“总成本”降到最低。从数学上讲我们是计算了这样一个函数：
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]。在读取原子坐标、对齐几何中心、对齐原子编号完成之后，接下来要计算平动-转动矩阵。从数学上讲，我们是要最小化这样一个函数：
原函数将简化为：
求极小值的函数可以写为：
需要优化的部分还剩下
这一项，通过第二节中纯四元数性质
代入可得：
# 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
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（参考了部分代码）