intrabend.f90

来自「蒙特卡罗的一个程序分析 与大家分享 共同研究」· F90 代码 · 共 58 行

subroutine IntraBend( X, Y, Z, kbend, theta_eq, Uintra )

implicit none

! X, Y, Z contain the coordinates of the existing beads.

real, dimension(3), intent(in)					:: X, Y, Z

! kbend is the angle bending modulus and theta_eq is the equilibrium angle.

real, intent(in)								:: kbend, theta_eq

! Uintra is the intramolecular bending energy.

real, intent(out)								:: Uintra

! Local Variables

real, parameter									:: pi = 3.14159265359
real											:: costh, th
real, dimension(3)								:: V, W

V(1) = X(1) - X(2)
V(2) = Y(1) - Y(2)
V(3) = Z(1) - Z(2)

W(1) = X(3) - X(2)
W(2) = Y(3) - Y(2)
W(3) = Z(3) - Z(2)

V = V / sqrt( dot_product( V, V ) )
W = W / sqrt( dot_product( W, W ) )

costh = dot_product( V, W ) 

if( costh > 1.0 .OR. costh < -1.0 ) then
	
	if( costh > 1.0 ) th = 0.0
	if( costh < -1.0 ) th = pi

else

	th = acos( costh )

end if

Uintra = 0.5 * kbend * ( th - theta_eq ) * ( th - theta_eq )

return

end subroutine IntraBend