四个点和三个点与面的输入的二面角的计算

1.计算咪唑环上三原子构成的面与坐标轴形成的面的二面角分布

使用二面角计算器计算一组数据以作参考

x0y

y0z

x0z

(1) 使用三点构成两个向量，计算法向量，利用向量的夹角计算二面角

验证二面角值的计算

N 31.01034171 62.77095584 53.16101503
N 32.60943575 63.81423837 54.26446791
C 30.93316986 64.10126884 52.83739902
用第一个分子的位置信息来验证代码的正确性。

与坐标轴形成的面的二面角	x0y	x0z	y0z
计算器参考值	139.9	81.38	51.25
本代码计算	40.06131242442094	81.38040749886316	128.74927894127075

观察发现，本代码的计算值有的与计算器一致，有的则是与计算器值互补。
这是因为，采用计算向量的cos，然后acos得到，计算得到的数值是0-180°，咪唑环上的法向量与坐标轴的法向量不在二面角的同侧时，就会出现计算结果互补的情况。

验证一帧的分布的计算

面y0z第一帧，二面角分布。通过使用EXCEL统计个数，验证计算相同角度个数的函数，没有问题。

验证多帧求和

最终结果
● x0y
● x0z
● y0z 与坐标轴三个面的6000帧的总体的分布情况

(2) 考虑到咪唑环处于同一个面，但是原子位置不同的时候，应当考虑为两种情况，有必要设法区分这两种状态，计算0-360°的二面角分布。

def dihedral_cos_360(p1, p2, p3, norm):
    V1 = vec_sub(p1, p2)
    V2 = vec_sub(p3, p2)
    norm1 = calc_normal_vector(V1, V2)
    norm2 = norm
    abs_norm1 = vec_length(norm1)
    abs_norm2 = vec_length(norm2)
    cos_theta = (norm1[0] * norm2[0] + norm1[1] * norm2[1] + norm1[2] * norm2[2])/(abs_norm1 * abs_norm2)
    theta_pi = math.acos(cos_theta)
    if ((V1[2] + V2[2]) > 0) :
        theta = theta_pi * (180 / math.pi)
    else:
        theta = -theta_pi * (180 / math.pi)
    return theta

简单的根据咪唑环上两个向量的和，当两个向量的和在Z轴的分量大于0，说明面的方向朝上。通过此方法，将二面角判定到-180 - 180°。其余计算方法仍然是按照cos的方法计算。

结果呈现

(3) 利用二面角新定义计算-180°到180°二面角

这里利用已知面的法向量，通过外积，找到了两个面相交的向量V2new，同时定义了一个能确定方向的V3new，按照三个向量构成二面角的计算公式，计算得到了二面角的分布函数。

def dihedral_tp_norm_tan(p1, p2, p3, norm): # 通过三个点的位置信息确定一个面，与已知法向量的面的二面角的夹角
    V1 = vec_sub(p1, p2)
    V2 = vec_sub(p3, p2)
    norm1 = calc_normal_vector(V1, V2)
    norm2 = norm
    V2new = calc_normal_vector(norm1, norm2)
    V3new = calc_normal_vector(V2new, norm2)
    V1new = V1
    m =calc_normal_vector(calc_normal_vector(V1new, V2new), calc_normal_vector(V2new, V3new))
    y = vec_dot(m , vec_scale(1/vec_length(V2new), V2new))
    x = vec_dot(calc_normal_vector(V1new, V2new), calc_normal_vector(V2new, V3new))
    DIHED = math.atan2(y,x) * 180/math.pi
    if (DIHED > 0):
        DIHED = 180 - DIHED
    else:
        DIHED = -180 - DIHED
    return DIHED

(4) 验证CONFIG文件（分子的初始位置）的分布

计算的结果显示，初始文件中就存在集中在90°附近。

2.计算阴离子C-S-S-C四个原子形成的两个面的二面角分布

验证二面角

s1
C -23.17885806 -21.76534081 15.50402625
C -26.95533703 -22.86314650 12.17431623
S -24.63272254 -22.89017219 15.37690433
S -25.50135682 -21.96109934 12.97211456
s2
C 7.90251561 22.78898166 -12.37276684
C 7.73501196 19.99085744 -9.29050914
S 7.12525099 21.19031621 -12.98969494
S 7.09221965 19.28633844 -10.90054852
s3
C 26.01163484 -19.58667385 6.09858539
C 27.44622182 -23.77730376 5.48816844
S 27.89430524 -19.79330731 6.03452357
S 28.60242134 -22.36978861 5.23343755

利用四点坐标余弦定理计算二面角

四点二面角	s1	s2	s3
二面角计算器	5.279	174.2	157.5
四点余弦定理法结果	5.526709306978917	176.34033176022595	155.3379774373199
新定义1	144.3939	27.07115	-74.363496

这个误差是因为计算器是采用了保留一位小数的值计算的。

利用二面角新定义计算

由于一个平面可以有几种定义方式(例如，通过它们中的向量或点，或它们的法向量)，二面角有几种等价的定义。
任何平面都可以由两个非共线的向量定义;求它们的外积并归一化得到平面的法向量。因此，二面角可以由四个成对非共线向量定义。
我们还可以定义三个非共线向量 ${\mathbf{b}}_1$、${\mathbf{b}}_2$和${\mathbf{b}}_3$的二面角(分别在图1中显示为红色、绿色和蓝色)。向量${\mathbf{b}}_1$和${\mathbf{b}}_2$定义了第一个平面，而${\mathbf{b}}_2$和${\mathbf{b}}_3$定义了第二个平面。二面角对应于一个外球面角(图1)，它是一个定义良好的有符号的量。

二面角的第二定义

参考网页1
$$\varphi =atan2\left(|{\mathbf{b}}_2|{\mathbf{b}}_1\cdot [{\mathbf{b}}_2\times {\mathbf{b}}_3],[{\mathbf{b}}_1\times {\mathbf{b}}_2]\cdot [{\mathbf{b}}_2\times {\mathbf{b}}_3]\right)$$
参考网页2
以上文字来自这两个参考页面，给出的文字信息一直，但是都没有图。
公式略有差别。
$$\phi =atan2\left(\left([{\mathbf{b}}_1\times {\mathbf{b}}_2]\times [{\mathbf{b}}_2\times {\mathbf{b}}_3]\right)\cdot \frac{{\mathbf{b}}_2}{|{\mathbf{b}}_2|},[{\mathbf{b}}_1\times {\mathbf{b}}_2]\cdot [{\mathbf{b}}_2\times {\mathbf{b}}_3]\right)$$
下文的计算，以公式2作为标准。

def calc_dihedral(p1, p2, p3, p4): # 通过四个点的位置信息计算二面角的夹角，四个点的位置命名为相邻的1234
    V1 = vec_sub(p1, p2)
    V2 = vec_sub(p3, p2)
    V3 = vec_sub(p4, p3)
    norm1 = calc_normal_vector(V1, V2)
    norm2 = calc_normal_vector(V2, V3)
    abs_norm1 = vec_length(norm1)
    abs_norm2 = vec_length(norm2)
    norm1 = vec_scale((1/abs_norm1) , norm1)
    norm2 = vec_scale((1/abs_norm2) , norm2)
    y = vec_dot(calc_normal_vector(norm1, norm2) , vec_scale((1/vec_length(V2)),V2))
    x = vec_dot(norm1, norm2)
    DIHED = math.atan2(y,x) * 180/math.pi
    if (DIHED > 0):
        DIHED = 180 - DIHED
    else:
        DIHED = -180 - DIHED
    return DIHED

结果分析

四点二面角	s1	s2	s3
二面角计算器	5.279	174.2	157.5
四点余弦定理法结果	5.526709306978917	176.34033176022595	155.3379774373199
新定义1	-174.5	-3.659668239775243	24.66202256268005

这个结果是与四点余弦定理结果一致的。
计算分布函数

与之前使用travis软件计算的结果对比

是非常接近的，在180°左右的误差可能是由于精度的原因导致的。这证明了这一新二面角定义，及其公式的正确性。