subroutine strtor_potential(bds, V)
use mod_constants
implicit none
type(ommp_bonded_type), intent(in) :: bds
! Bonded potential data structure
real(rp), intent(inout) :: V
real(rp) :: thet, costhet, dihef(3), dr(3), r(3), vst
integer(ip) :: i, j, k, ib1, ib2, ib3
if(.not. bds%use_strtor) return
!$omp parallel do default(shared) reduction(+:V) &
!$omp private(i,costhet,thet,j,dihef,ib1,ib2,ib3,r,dr,vst,k)
do i=1, bds%nstrtor
! Atoms that defines the dihedral angle
costhet = cos_torsion(bds%top, bds%strtorat(:,i))
thet = acos(costhet)
do j=1, 3
dihef(j) = 1.0 + cos(j*thet+bds%torsphase(j,bds%strtor_t(i)))
end do
ib1 = bds%strtor_b(1,i)
ib2 = bds%strtor_b(2,i)
ib3 = bds%strtor_b(3,i)
r(1) = norm2(bds%top%cmm(:, bds%bondat(1,ib1)) - &
bds%top%cmm(:, bds%bondat(2,ib1)))
r(2) = norm2(bds%top%cmm(:, bds%bondat(1,ib2)) - &
bds%top%cmm(:, bds%bondat(2,ib2)))
r(3) = norm2(bds%top%cmm(:, bds%bondat(1,ib3)) - &
bds%top%cmm(:, bds%bondat(2,ib3)))
dr(1) = r(1) - bds%l0bond(ib1)
dr(2) = r(2) - bds%l0bond(ib2)
dr(3) = r(3) - bds%l0bond(ib3)
do j=1,3
vst = 0.0
do k=1, 3
vst = vst + bds%strtork((j-1)*3+k,i) * dihef(k)
end do
V = V + vst * dr(j)
end do
end do
end subroutine strtor_potential