计算化学公社

标题: 求助通过笛卡尔坐标，计算键角二面角的公式 [打印本页]

作者
Author: Laura_ 时间: 2020-8-27 16:40
标题: 求助通过笛卡尔坐标，计算键角二面角的公式
老师们好，请问通过原子的笛卡尔坐标（每个原子的绝对坐标，即XYZ）求键角与二面角的计算公式？感谢！键长的公式为sqrt((x1-x2)^2+(y1-y2)^2+(z1-z2)^2)
想请问一下键角与二面角的计算公式，谢谢各位老师！

作者
Author: granvia 时间: 2020-8-27 20:29
采用矢量运算就很简单。对于键角，设a=v1-v2, b=v3-v2，则矢量a与b的夹角就是键角，a^b=acos( a.b/|a|/|b| )
对于二面角，可令b=v4-v3, 再求a^b即可

作者
Author: highlight 时间: 2020-8-27 21:35
http://bbs.keinsci.com/thread-12861-1-1.html

作者
Author: 喵星大佬 时间: 2020-8-27 21:55
这不是最基本的解析几何问题嘛。。。。

作者
Author: sobereva 时间: 2020-8-28 04:27
注意你的贴图方式不对，仔细看置顶的新社员必读贴了解怎么正确贴图

作者
Author: sobereva 时间: 2020-8-28 04:28
Multiwfn源代码包的util.f90里的xyz2angle和xyz2dih子程序就是算键角和二面角的，直接输入原子的笛卡尔坐标就能用，一看便知

作者
Author: Daniel_Arndt 时间: 2020-8-28 11:40
我以前写过一段awk代码，就是矢量运算。纯粹的高中数学。

#The unit of the results given by function fangle is degree.
function fangle(a1,a2,a3, x1,y1,z1,x2,y2,z2,x3,y3,z3,cosangle)
{
x1 = Atom[a1,1];
y1 = Atom[a1,2];
z1 = Atom[a1,3];
x2 = Atom[a2,1];
y2 = Atom[a2,2];
z2 = Atom[a2,3];
x3 = Atom[a3,1];
y3 = Atom[a3,2];
z3 = Atom[a3,3];
cosangle = ((x1-x2)*(x3-x2)+(y1-y2)*(y3-y2)+(z1-z2)*(z3-z2))/sqrt(((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2))*((x3-x2)*(x3-x2)+(y3-y2)*(y3-y2)+(z3-z2)*(z3-z2)));
return 57.29577951*atan2(sqrt(1-cosangle*cosangle),cosangle);
}
#The unit of the results given by function fdihedral is degree.
function fdihedral(d1,d2,d3,d4, x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,norm1,x5,y5,z5,norm2,x6,y6,z6)
{
x1 = Atom[d1,1];
y1 = Atom[d1,2];
z1 = Atom[d1,3];
x2 = Atom[d2,1];
y2 = Atom[d2,2];
z2 = Atom[d2,3];
x3 = Atom[d3,1];
y3 = Atom[d3,2];
z3 = Atom[d3,3];
x4 = Atom[d4,1];
y4 = Atom[d4,2];
z4 = Atom[d4,3];
norm1 = sqrt((y1*(z2-z3)+y2*(z3-z1)+y3*(z1-z2))*(y1*(z2-z3)+y2*(z3-z1)+y3*(z1-z2))+(z1*(x2-x3)+z2*(x3-x1)+z3*(x1-x2))*(z1*(x2-x3)+z2*(x3-x1)+z3*(x1-x2))+(x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2))*(x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2)));
x5 = (y1*(z2-z3)+y2*(z3-z1)+y3*(z1-z2))/norm1;
y5 = (z1*(x2-x3)+z2*(x3-x1)+z3*(x1-x2))/norm1;
z5 = (x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2))/norm1;
norm2 = sqrt((x3-x2)*(x3-x2)+(y3-y2)*(y3-y2)+(z3-z2)*(z3-z2));
x6 = (y5*(z3-z2)-(y3-y2)*z5)/norm2;
y6 = (z5*(x3-x2)-(z3-z2)*x5)/norm2;
z6 = (x5*(y3-y2)-(x3-x2)*y5)/norm2;
return -57.29577951*atan2(x6*(y2*(z3-z4)+y3*(z4-z2)+y4*(z2-z3))+y6*(z2*(x3-x4)+z3*(x4-x2)+z4*(x2-x3))+z6*(x2*(y3-y4)+x3*(y4-y2)+x4*(y2-y3)),x5*(y2*(z3-z4)+y3*(z4-z2)+y4*(z2-z3))+y5*(z2*(x3-x4)+z3*(x4-x2)+z4*(x2-x3))+z5*(x2*(y3-y4)+x3*(y4-y2)+x4*(y2-y3)));
}

作者
Author: liuyuje714 时间: 2020-8-28 11:53
这都简单，注意二面角的取值范围，如果是-180~180，你应该考虑方向。我这有个小程序可以直接计算键长，键角，二面角，支持xyz/gjf/gro/pdb,。https://github.com/liuyujie714/T ... aster/cal_bondangle

作者
Author: granvia 时间: 2020-8-28 16:32

Daniel_Arndt 发表于 2020-8-28 11:40
我以前写过一段awk代码，就是矢量运算。纯粹的高中数学。

#The unit of the results given by function ...

这代码可进一步优化，很多相同的表达式重复运算了

作者
Author: Daniel_Arndt 时间: 2020-9-1 09:22

granvia 发表于 2020-8-28 16:32
这代码可进一步优化，很多相同的表达式重复运算了

确实可以，但当时就是临时写来用的，我也不想再折腾自己了

。

作者
Author: Laura_ 时间: 2020-9-1 15:05
十分感谢大家！以后贴图的方式会更加注意，谢谢各位老师

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