urey_potential Subroutine

    subroutine urey_potential(bds, V)
        !! Compute the Urey-Bradley potential.  
        !! This is basically a virtual bond, with its stretching harminic 
        !! potential that connect two otherwise un-connected bonds. The same
        !! potential formula used for normal stretching is used.

        use mod_constants, only : eps_rp

        implicit none

        type(ommp_bonded_type), intent(in) :: bds
        ! Bonded potential data structure
        real(rp), intent(inout) :: V
        !! Urey-Bradley potential, result will be added to V

        integer :: i
        logical(lp) :: use_cubic, use_quartic
        real(rp) :: dr(3), l, dl, dl2
        
        if(.not. bds%use_urey) return

        use_cubic = (abs(bds%urey_cubic) > eps_rp)
        use_quartic = (abs(bds%urey_quartic) > eps_rp)

        if(.not. use_cubic .and. .not. use_quartic) then
            ! This is just a regular harmonic potential
            !$omp parallel do default(shared) reduction(+:V) &
            !$omp private(i,dr,l,dl)
            do i=1, bds%nurey
                dr = bds%top%cmm(:,bds%ureyat(1,i)) - &
                     bds%top%cmm(:,bds%ureyat(2,i))
                l = sqrt(dot_product(dr, dr))
                dl = l - bds%l0urey(i)
                V = V + bds%kurey(i) * dl * dl
            end do
        else
            !$omp parallel do default(shared) reduction(+:V) &
            !$omp private(i,dr,l,dl,dl2)
            do i=1, bds%nurey
                dr = bds%top%cmm(:,bds%ureyat(1,i)) - &
                     bds%top%cmm(:,bds%ureyat(2,i))
                l = sqrt(dot_product(dr, dr))
                dl = l - bds%l0urey(i)
                dl2 = dl * dl

                V = V + bds%kurey(i)*dl2 * (1.0_rp + bds%urey_cubic*dl + &
                                            bds%urey_quartic*dl2)
            end do
        end if
    end subroutine urey_potential