十分感谢大家！以后贴图的方式会更加注意，谢谢各位老师

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

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

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

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

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

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

我以前写过一段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)));
}

Multiwfn源代码包的util.f90里的xyz2angle和xyz2dih子程序就是算键角和二面角的，直接输入原子的笛卡尔坐标就能用，一看便知

注意你的贴图方式不对，仔细看置顶的新社员必读贴了解怎么正确贴图

这不是最基本的解析几何问题嘛。。。。

http://bbs.keinsci.com/thread-12861-1-1.html

采用矢量运算就很简单。对于键角，设a=v1-v2, b=v3-v2，则矢量a与b的夹角就是键角，a^b=acos( a.b/|a|/|b| )
对于二面角，可令b=v4-v3, 再求a^b即可