計算兩個平面夾角的小程序

2021-Mar-3補充：后來在Multiwfn里加入了對動力學軌跡計算指定的兩個平面間夾角隨幀號變化的功能，見《計算分子動力學軌跡中兩個環平面間的距離和夾角》（http://www.shanxitv.org/590）。

---

今天思想家公社QQ群有人問怎么計算兩個平面的夾角。如果本身是二面角好辦直接量就行了，但如果不屬于這種情況就稍微麻煩，雖然有一些可視化程序如diamond可以量但終究為這點事安裝一個也麻煩。這里提供個自寫的不足掛齒的程序twoplane2angle來計算夾角。先輸入定義第一個平面的三個點的坐標，再輸入定義第二個平面的三個點的坐標，然后就會輸出夾角。

下載地址：

/usr/uploads/file/20150820/20150820222744_68251.rar

源代碼

program twoplane2angle
implicit real*8 (a-h,o-z)
write(*,*) "twoplane2angle: Get angle between two planes"
write(*,*) "Programmed by Sobereva (sobereva@sina.com), 2015-Aug-20"
write(*,"(a)") " Usage: Input three points for plane 1 and that for plane 2, then program outputs the angle between the two planes"
write(*,*)
do while(.true.)
 write(*,*) "Input XYZ coordinates of point 1 for plane 1, e.g. 0.2,0.4,-3.6"
 read(*,*) x1,y1,z1
 write(*,*) "Input XYZ coordinates of point 2 for plane 1, e.g. 0.2,0.4,-3.6"
 read(*,*) x2,y2,z2
 write(*,*) "Input XYZ coordinates of point 3 for plane 1, e.g. 0.2,0.4,-3.6"
 read(*,*) x3,y3,z3
 call pointABCD(x1,y1,z1,x2,y2,z2,x3,y3,z3,A1,B1,C1,D1)
 write(*,*) "Input XYZ coordinates of point 1 for plane 2, e.g. 0.2,0.4,-3.6"
 read(*,*) x1,y1,z1
 write(*,*) "Input XYZ coordinates of point 2 for plane 2, e.g. 0.2,0.4,-3.6"
 read(*,*) x2,y2,z2
 write(*,*) "Input XYZ coordinates of point 3 for plane 2, e.g. 0.2,0.4,-3.6"
 read(*,*) x3,y3,z3
 call pointABCD(x1,y1,z1,x2,y2,z2,x3,y3,z3,A2,B2,C2,D2)
 pt1norm=dsqrt(A1**2+B1**2+C1**2)
 pt2norm=dsqrt(A2**2+B2**2+C2**2)
 vecprod=A1*A2+B1*B2+C1*C2
 write(*,"(' Norm of normal vector 1:',f12.6)") pt1norm
 write(*,"(' Norm of normal vector 2:',f12.6)") pt2norm
 write(*,"(' Product of the two vectors:',f12.6)") vecprod
 ang=acos(vecprod/pt1norm/pt2norm)
 write(*,"(a,f12.4)") " The angle between the two planes is",ang/3.141592653589793238462D0*180
 write(*,*)
end do

end program

!!-------- Input three points, get ABCD of Ax+By+Cz+D=0
subroutine pointABCD(x1,y1,z1,x2,y2,z2,x3,y3,z3,A,B,C,D)
real*8 v1x,v1y,v1z,v2x,v2y,v2z,x1,y1,z1,x2,y2,z2,x3,y3,z3,A,B,C,D
v1x=x2-x1
v1y=y2-y1
v1z=z2-z1
v2x=x3-x1
v2y=y3-y1
v2z=z3-z1
! Solve determinant(Vector multiply) to get the normal vector (A,B,C):
!  i   j   k   //unit vector
! v1x v1y v1z
! v2x v2y v2z
A=v1y*v2z-v1z*v2y
B=-(v1x*v2z-v1z*v2x)
C=v1x*v2y-v1y*v2x
D=A*(-x1)+B*(-y1)+C*(-z1)
end subroutine