关于提取二面角的python求助

zhuyoucai

各位老师好，本人目前刚开始学习python辅助计算化学，下面是我写的一个提取二面角的脚本。但是出现一个问题就是我提取的是原子1，2，3，4的二面角，最终输出的结果永远为正值，具体数值是正确的。所以想请教各位，该如何修改才能正确输出正负值呢，谢谢各位老师！

import math
import numpy as np
from numpy.linalg import norm
from numpy import dot, arccos, degrees, cross


def calculate_angle(coord1, coord2, coord3, coord4):
    x1, y1, z1 = coord1
    x2, y2, z2 = coord2
    x3, y3, z3 = coord3
    x4, y4, z4 = coord4
    a = [x2-x1,y2-y1,z2-z1]
    b = [x3-x2,y3-y2,z3-z2]
    c = [x4-x3,y4-y3,z4-z3]
    n1 = cross(a,b)
    n2 = cross(b,c)
    z = dot(n1,n2)/(norm(n1)*norm(n2))
    angle_rad = math.acos(z)
    angle_deg = math.degrees(angle_rad)
    return angle_deg

coord1 = [5,2,3]
coord2 = [2,3,2]
coord3 = [3,4,3]
coord4 = [4,1,4]

angle = calculate_angle(coord1, coord2, coord3, coord4)
print(angle)

wzkchem5

用z的正负判断。至于正的z对应正的二面角还是负的二面角，和一个已知能正确判断二面角符号的可视化软件对比一下就知道了

sobereva

可以参考Multiwfn源代码包里util.f90里的xyz2dih_sign函数

highlight

得不到负值就对了，你用的公式就不能输出负值
如果不明白，想想 arccoss 的值域
换一个定义式，使用 arctan 的那个
参考 https://stackoverflow.com/questi ... ordinates-in-python
或者我在 http://bbs.keinsci.com/thread-12861-1-1.html 下面的回复

zhuyoucai

wzkchem5 发表于 2023-6-27 21:32
用z的正负判断。至于正的z对应正的二面角还是负的二面角，和一个已知能正确判断二面角符号的可视化软件对比 ...

好的，谢谢老师

zhuyoucai

highlight 发表于 2023-6-28 08:43
得不到负值就对了，你用的公式就不能输出负值
如果不明白，想想 arccoss 的值域
换一个定义式，使用 arct ...

好的，谢谢老师

zhuyoucai

sobereva 发表于 2023-6-28 00:06
可以参考Multiwfn源代码包里util.f90里的xyz2dih_sign函数

好的，谢谢老师