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