torsion_potential Subroutine

    subroutine torsion_potential(bds, V)
        !! Compute torsion potential
        use mod_constants, only: pi, eps_rp

        implicit none

        type(ommp_bonded_type), intent(in) :: bds
        ! Bonded potential data structure
        real(rp), intent(inout) :: V
        !! torsion potential, result will be added to V
        real(rp) :: thet, costhet
        integer(ip) :: i, j
        
        if(.not. bds%use_torsion) return

        !$omp parallel do default(shared) &
        !$omp private(i,costhet,thet,j) reduction(+:V)
        do i=1, bds%ntorsion
            ! Atoms that defines the dihedral angle
            costhet = cos_torsion(bds%top, bds%torsionat(:,i))
            
            if(costhet + 1.0 <= eps_rp) then
                thet = pi
            else if(abs(costhet - 1.0) <= eps_rp) then
                thet = 0.0
            else
                thet = acos(costhet)
            end if

            do j=1, 6
                if(bds%torsn(j,i) < 1) exit
                V = V + bds%torsamp(j,i) * (1+cos(real(bds%torsn(j,i))*thet &
                                            - bds%torsphase(j,i)))
            end do
        end do

    end subroutine torsion_potential