如何计算四个原子的二面角

首先给出Python实例，计算得到的结果可以在各种分子可视化软件上检查：

给定四个原子坐标，就可以返回二面角的数值。

import math

import numpy as np


def calculate_dihedral_angle(A, B, C, D):
    a = [B[0]-A[0], B[1]-A[1], B[2]-A[2]]
    b = [C[0]-B[0], C[1]-B[1], C[2]-B[2]]   # core vector 决定了旋转轴
    c = [D[0]-C[0], D[1]-C[1], D[2]-C[2]]
    n1, n2 = np.cross(a, b), np.cross(b, c)
    b_ = np.cross(n1, n2)  # 决定了两个原子的旋转方向
    if sum(b*b_) > 0:
        return math.acos(np.dot(n1, n2)/(np.linalg.norm(n1)*np.linalg.norm(n2))) * 180 / math.pi
    else:
        return -math.acos(np.dot(n1, n2)/(np.linalg.norm(n1)*np.linalg.norm(n2))) * 180 / math.pi

# 举例说明
A = [1.00, 2.92, 2.0]
B = [1.24, 2.942, 2.0]
C = [3.40, 2.15, 5.0]
D = [1.420, 1.9, 8.0]
calculate_dihedral_angle(A, B, C, D)

其中A，B，C，D是以列表形式给出的直角坐标 (例如： A = [x, y, z] )

1. 给出4个原子的坐标A、B、C、D（应该是笛卡尔坐标，内坐标的话也就不需要算了）

2. 注意4个原子的顺序（这涉及到二面角的大小，正负，可以参考各种分子可视化软件的选择二面角方法）

3. 求三个向量a = B-A, b = C-B, c = D-C（这一步需要注意减法的顺序，因为这决定了向量的方向，从而决定了平面法向量的方向）

4. 使用数学上向量的方法叉乘，n1=a×b，n2=b×c；（注意顺序，决定了法向量的方向，方向错了二面角会错）

5. 得出的两个法向量用余弦定理求角度 θ = arccos（n1·n2）/（|n1|*|n2|)

6. 由于之前我们严格定义了各种矢量的方向，所以得出的即四个原子的二面角（写程序的时候注意弧度与角度之间的转化）

实际应用中，二面角具有360°的旋转角度，在一些软件中是（-180°~180°），也有的是（0° ~360°），而此时我们用余弦计算角度的结果只能返回0~180，那如何实现正负的区分呢？

这里需要比较两个平面法向量外积与旋转轴的方向，如果同向，结果取正，如果反向，结果取负。

数学原理：

叉乘会根据给定的两个向量的方向和大小给出法向方向的向量（方向：右手定则）

A，B，C 就对应着AB，BC两个向量，叉乘给出的法向量方向是唯一的，而根据平面法向量和平面上的向量相互正交的性质列方程解方程组得到的法向量是有正负两个方向的（因为正负两个方向都满足正交的条件，算二面角的时候就会出现问题）

同理，B，C，D也会给出一个唯一方向的法向量，两个法向量之间会有（0，180°）的夹角，算出的cosθ=（n1·n2）/（|n1|*|n2|)也会有正负的分别，但根据余弦图像可得，每一个值都对应着一个在（0，180°）内确定的角度。

然而这还不是最终的结果，例如，对于一个计算得到的90°的角，他可能是旋转了90°，也可能是反方向转了270°，虽然两个“面”的夹角确实只会出现0°~180°，但在许多问题上并非如此。例如（判断化学键的旋转方向）下面以一个顺反异构的化学现象作为例子，对于分子ABCD和ABCD‘，如果按照之前的方法求n1和n2角度，那么都会得到90°的结果（两个平面也确实是相互垂直的）。但事实上如果以BC向量的方向为参照，它们旋转的方向是不同的。此时我们要额外增加一步判断旋转方向的步骤：

就是用n1和n2的外积与BC对比方向，如果同向，结果取正号，反向结果取负号（这一步是因为ABCD的顺序被固定了，BC的方向就是B指向C）此时，计算的结果就可以与分子可视化软件中的二面角对应上。