subroutine angtor_potential(bds, V)
implicit none
type(ommp_bonded_type), intent(in) :: bds
! Bonded potential data structure
real(rp), intent(inout) :: V
real(rp) :: thet, costhet, dihef(3), delta_a(2), vat, l1, l2, &
dr1(3), dr2(3), angle1, angle2
integer(ip) :: i, j, k, ia1, ia2
if(.not. bds%use_angtor) return
!$omp parallel do default(shared) reduction(+:V) &
!$omp private(i,costhet,thet,j,dihef,ia1,ia2,dr1,dr2,l1,l2,angle1,angle2,delta_a,vat,k)
do i=1, bds%nangtor
! Atoms that defines the dihedral angle
costhet = cos_torsion(bds%top, bds%angtorat(:,i))
thet = acos(costhet)
do j=1, 3
dihef(j) = 1.0 + cos(j*thet+bds%torsphase(j,bds%angtor_t(i)))
end do
ia1 = bds%angtor_a(1,i)
ia2 = bds%angtor_a(2,i)
dr1 = bds%top%cmm(:, bds%angleat(1,ia1)) - &
bds%top%cmm(:, bds%angleat(2,ia1))
dr2 = bds%top%cmm(:, bds%angleat(3,ia1)) - &
bds%top%cmm(:, bds%angleat(2,ia1))
l1 = norm2(dr1)
l2 = norm2(dr2)
angle1 = acos(dot_product(dr1, dr2)/(l1*l2))
dr1 = bds%top%cmm(:, bds%angleat(1,ia2)) - &
bds%top%cmm(:, bds%angleat(2,ia2))
dr2 = bds%top%cmm(:, bds%angleat(3,ia2)) - &
bds%top%cmm(:, bds%angleat(2,ia2))
l1 = norm2(dr1)
l2 = norm2(dr2)
angle2 = acos(dot_product(dr1, dr2)/(l1*l2))
delta_a(1) = angle1 - bds%eqangle(bds%angtor_a(1,i))
delta_a(2) = angle2 - bds%eqangle(bds%angtor_a(2,i))
do j=1,2
vat = 0.0
do k=1, 3
vat = vat + bds%angtork((j-1)*3+k,i) * dihef(k)
end do
V = V + vat * delta_a(j)
end do
end do
end subroutine angtor_potential