计算化学公社

标题: 关于提取二面角的python求助 [打印本页]

作者
Author: zhuyoucai 时间: 2023-6-27 21:13
标题: 关于提取二面角的python求助

各位老师好，本人目前刚开始学习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)

作者
Author: wzkchem5 时间: 2023-6-27 21:32
用z的正负判断。至于正的z对应正的二面角还是负的二面角，和一个已知能正确判断二面角符号的可视化软件对比一下就知道了

作者
Author: sobereva 时间: 2023-6-28 00:06
可以参考Multiwfn源代码包里util.f90里的xyz2dih_sign函数

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

作者
Author: zhuyoucai 时间: 2023-6-28 08:53

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

好的，谢谢老师

作者
Author: zhuyoucai 时间: 2023-6-28 08:53

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

好的，谢谢老师

作者
Author: zhuyoucai 时间: 2023-6-28 08:53

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

好的，谢谢老师

欢迎光临计算化学公社 (http://bbs.keinsci.com/)