intratorsion.f90

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

subroutine IntraTorsion( X, Y, Z, model, C, Uintra	)

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

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

! model is the type of torsion potential being used.
! model = 1 for the OPLS model, Uintra = SUM( C(k) * cos(phi)**k ) for k=0..5
! model = 2 for the de Pablo model, 
! Uintra = C(O) + C(1)*(1+cosphi) + C(2)*(1-cos2phi) + C(3)*(1+cos3phi)

integer, intent(in)								:: model

! C contains the constants in the torsion potential.
! for model = 1, C has 6 entries, 0..5
! for model = 2, C has 4 entries, 0..3

real, dimension(0:5), intent(in)				:: C

! Uintra is the torsional energy.

real, intent(out)								:: UINTRA 

! Local Variables

integer											:: j
real											:: cosphi
real, dimension(3)								:: S, T, U, V, W

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

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

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

S(1) = U(2) * V(3) - U(3) * V(2)
S(2) = U(3) * V(1) - U(1) * V(3)
S(3) = U(1) * V(2) - U(2) * V(1)

T(1) = W(2) * V(3) - W(3) * V(2)
T(2) = W(3) * V(1) - W(1) * V(3)
T(3) = W(1) * V(2) - W(2) * V(1)

if( model == 1 ) then

	S = - S / sqrt( dot_product( S, S ) )
	T = T / sqrt( dot_product( T, T ) )

	cosphi = dot_product( S, T )

	Uintra = C(0)

	do j = 1, 5

		Uintra = Uintra + C(j) * cosphi ** j

	end do

else if( model == 2 ) then

!	S = S / sqrt( dot_product( S, S ) )
	S = S / sqrt( dot_product( S, S ) )
	T = T / sqrt( dot_product( T, T ) )

	cosphi = dot_product( S, T )

	Uintra = C(0) + 0.5 * &
			 ( C(1) * ( 1.0 + cosphi ) + &
			   2.0 * C(2) * (1.0 - cosphi * cosphi ) + &
			   C(3) * ( 1.0 + 4.0 * cosphi * cosphi * cosphi - 3.0 * cosphi ) )

end if

return

end subroutine IntraTorsion