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											:: 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 ) 

th = acos( costh )

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

return

end subroutine IntraBend