vdw_potential_inter Subroutine

    subroutine vdw_potential_inter(vdw1, vdw2, V)
        !! Compute the dispersion repulsion energy for the whole system
        !! using a double loop algorithm

        use mod_io, only : fatal_error
        use mod_constants, only: eps_rp
        implicit none

        type(ommp_nonbonded_type), intent(in), target :: vdw1, vdw2
        !! Nonbonded data structure
        real(rp), intent(inout) :: V
        !! Potential, result will be added

        integer(ip) :: i, j, ineigh
        real(rp) :: eij, rij0, rij, ci(3), cj(3), vtmp
        type(ommp_topology_type), pointer :: top1, top2
        procedure(vdw_term), pointer :: vdw_func

        top1 => vdw1%top
        top2 => vdw2%top

        if(vdw1%vdwtype /= vdw2%vdwtype .or. &
           vdw1%radrule /= vdw2%radrule .or. &
           vdw1%epsrule /= vdw2%epsrule) then
            call fatal_error("Requested VdW potential between two incompatible &
                             &VdW groups.")
       end if

        select case(vdw1%vdwtype)
            case(OMMP_VDWTYPE_LJ)
                vdw_func => vdw_lennard_jones
            case(OMMP_VDWTYPE_BUF714)
                vdw_func => vdw_buffered_7_14
            case default
                vdw_func => vdw_buffered_7_14
                call fatal_error("Unexpected error in vdw_potential_inter")
        end select

        !$omp parallel do default(shared) schedule(dynamic) &
        !$omp private(i,j,ci,cj,ineigh,Rij,Rij0,Eij,vtmp) reduction(+:v) 
        do i=1, top1%mm_atoms
            if(abs(vdw1%vdw_f(i) - 1.0_rp) < eps_rp) then
                ci = top1%cmm(:,i)
            else
                ! Scale factors are used only for monovalent atoms, in that
                ! case the vdw center is displaced along the axis connecting
                ! the atom to its neighbour
                if(top1%conn(1)%ri(i+1) - top1%conn(1)%ri(i) /= 1) then
                    call fatal_error("Scale factors are only expected for &
                                     &monovalent atoms (top1)")
                end if
                ineigh = top1%conn(1)%ci(top1%conn(1)%ri(i))

                ci = top1%cmm(:,ineigh) + (top1%cmm(:,i) - top1%cmm(:,ineigh)) &
                     * vdw1%vdw_f(i)
            endif
                
            do j=1, top2%mm_atoms
                Rij0 = get_Rij0_inter(vdw1, vdw2, i, j)
                eij = get_eij_inter(vdw1, vdw2, i, j)
            
                if(abs(vdw2%vdw_f(j) - 1.0_rp) < eps_rp) then
                    cj = top2%cmm(:,j)
                else
                    ! Scale factors are used only for monovalent atoms, in that
                    ! case the vdw center is displaced along the axis connecting
                    ! the atom to its neighbour
                    if(top2%conn(1)%ri(j+1) - top2%conn(1)%ri(j) /= 1) then
                        call fatal_error("Scale factors are only expected &
                                         & for monovalent atoms (top2) atom")
                    end if
                    ineigh = top2%conn(1)%ci(top2%conn(1)%ri(j))

                    cj = top2%cmm(:,ineigh) + &
                         (top2%cmm(:,j) - top2%cmm(:,ineigh)) * vdw2%vdw_f(j)
                endif
                Rij = norm2(ci-cj)
                if(Rij < eps_rp) then
                    call fatal_error("Requesting inter non-bonded potential for two atoms &
                                     &placed in the same point, this could be &
                                     &an internal bug or a problem in your input &
                                     &file, please check.")
                end if
                
                vtmp = 0.0_rp
                call vdw_func(Rij, Rij0, Eij, vtmp)
                v = v + vtmp
            end do
        end do
    end subroutine vdw_potential_inter