mod_nonbonded.F90 Source File


This file depends on

sourcefile~~mod_nonbonded.f90~~EfferentGraph sourcefile~mod_nonbonded.f90 mod_nonbonded.F90 sourcefile~mod_memory.f90 mod_memory.F90 sourcefile~mod_nonbonded.f90->sourcefile~mod_memory.f90 sourcefile~mod_topology.f90 mod_topology.F90 sourcefile~mod_nonbonded.f90->sourcefile~mod_topology.f90 sourcefile~mod_neighbors_list.f90 mod_neighbors_list.F90 sourcefile~mod_nonbonded.f90->sourcefile~mod_neighbors_list.f90 sourcefile~mod_constants.f90 mod_constants.F90 sourcefile~mod_nonbonded.f90->sourcefile~mod_constants.f90 sourcefile~mod_adjacency_mat.f90 mod_adjacency_mat.F90 sourcefile~mod_nonbonded.f90->sourcefile~mod_adjacency_mat.f90 sourcefile~mod_io.f90 mod_io.F90 sourcefile~mod_nonbonded.f90->sourcefile~mod_io.f90 sourcefile~mod_profiling.f90 mod_profiling.F90 sourcefile~mod_nonbonded.f90->sourcefile~mod_profiling.f90 sourcefile~mod_jacobian_mat.f90 mod_jacobian_mat.F90 sourcefile~mod_nonbonded.f90->sourcefile~mod_jacobian_mat.f90 sourcefile~mod_memory.f90->sourcefile~mod_constants.f90 sourcefile~mod_memory.f90->sourcefile~mod_io.f90 sourcefile~mod_topology.f90->sourcefile~mod_memory.f90 sourcefile~mod_topology.f90->sourcefile~mod_constants.f90 sourcefile~mod_topology.f90->sourcefile~mod_adjacency_mat.f90 sourcefile~mod_topology.f90->sourcefile~mod_io.f90 sourcefile~mod_neighbors_list.f90->sourcefile~mod_memory.f90 sourcefile~mod_neighbors_list.f90->sourcefile~mod_constants.f90 sourcefile~mod_neighbors_list.f90->sourcefile~mod_adjacency_mat.f90 sourcefile~mod_neighbors_list.f90->sourcefile~mod_io.f90 sourcefile~mod_neighbors_list.f90->sourcefile~mod_profiling.f90 sourcefile~mod_adjacency_mat.f90->sourcefile~mod_memory.f90 sourcefile~mod_utils.f90 mod_utils.F90 sourcefile~mod_adjacency_mat.f90->sourcefile~mod_utils.f90 sourcefile~mod_io.f90->sourcefile~mod_constants.f90 sourcefile~mod_profiling.f90->sourcefile~mod_memory.f90 sourcefile~mod_profiling.f90->sourcefile~mod_constants.f90 sourcefile~mod_profiling.f90->sourcefile~mod_io.f90 sourcefile~mod_jacobian_mat.f90->sourcefile~mod_memory.f90 sourcefile~mod_jacobian_mat.f90->sourcefile~mod_constants.f90 sourcefile~mod_jacobian_mat.f90->sourcefile~mod_utils.f90 sourcefile~mod_utils.f90->sourcefile~mod_memory.f90 sourcefile~mod_utils.f90->sourcefile~mod_constants.f90

Files dependent on this one

sourcefile~~mod_nonbonded.f90~~AfferentGraph sourcefile~mod_nonbonded.f90 mod_nonbonded.F90 sourcefile~mod_mmpol.f90 mod_mmpol.F90 sourcefile~mod_mmpol.f90->sourcefile~mod_nonbonded.f90 sourcefile~mod_link_atom.f90 mod_link_atom.F90 sourcefile~mod_mmpol.f90->sourcefile~mod_link_atom.f90 sourcefile~mod_qm_helper.f90 mod_qm_helper.F90 sourcefile~mod_qm_helper.f90->sourcefile~mod_nonbonded.f90 sourcefile~mod_qm_helper.f90->sourcefile~mod_mmpol.f90 sourcefile~mod_qm_helper.f90->sourcefile~mod_link_atom.f90 sourcefile~mod_prm.f90 mod_prm.F90 sourcefile~mod_qm_helper.f90->sourcefile~mod_prm.f90 sourcefile~mod_link_atom.f90->sourcefile~mod_nonbonded.f90 sourcefile~mod_link_atom.f90->sourcefile~mod_prm.f90 sourcefile~mod_prm.f90->sourcefile~mod_nonbonded.f90 sourcefile~mod_iohdf5.f90 mod_iohdf5.F90 sourcefile~mod_iohdf5.f90->sourcefile~mod_nonbonded.f90 sourcefile~mod_iohdf5.f90->sourcefile~mod_mmpol.f90 sourcefile~mod_interface.f90 mod_interface.F90 sourcefile~mod_interface.f90->sourcefile~mod_nonbonded.f90 sourcefile~mod_interface.f90->sourcefile~mod_mmpol.f90 sourcefile~mod_interface.f90->sourcefile~mod_qm_helper.f90 sourcefile~mod_interface.f90->sourcefile~mod_link_atom.f90 sourcefile~mod_interface.f90->sourcefile~mod_prm.f90 sourcefile~mod_interface.f90->sourcefile~mod_iohdf5.f90 sourcefile~mod_inputloader.f90 mod_inputloader.F90 sourcefile~mod_interface.f90->sourcefile~mod_inputloader.f90 sourcefile~mod_geomgrad.f90 mod_geomgrad.F90 sourcefile~mod_interface.f90->sourcefile~mod_geomgrad.f90 sourcefile~mod_polarization.f90 mod_polarization.F90 sourcefile~mod_interface.f90->sourcefile~mod_polarization.f90 sourcefile~mod_inputloader.f90->sourcefile~mod_mmpol.f90 sourcefile~mod_inputloader.f90->sourcefile~mod_prm.f90 sourcefile~mod_c_interface.f90 mod_c_interface.F90 sourcefile~mod_c_interface.f90->sourcefile~mod_mmpol.f90 sourcefile~mod_c_interface.f90->sourcefile~mod_qm_helper.f90 sourcefile~mod_c_interface.f90->sourcefile~mod_interface.f90 sourcefile~mod_geomgrad.f90->sourcefile~mod_mmpol.f90 sourcefile~mod_geomgrad.f90->sourcefile~mod_polarization.f90 sourcefile~mod_polarization.f90->sourcefile~mod_mmpol.f90

Contents

Source Code


Source Code

module mod_nonbonded

    use mod_memory, only: rp, ip, lp
    use mod_neighbor_list, only: ommp_neigh_list
    use mod_constants, only: OMMP_STR_CHAR_MAX
    use mod_adjacency_mat, only: yale_sparse
    use mod_topology, only: ommp_topology_type
    use mod_constants, only: OMMP_VDWTYPE_LJ, & 
                             OMMP_VDWTYPE_BUF714, &
                             OMMP_RADRULE_ARITHMETIC, &
                             OMMP_RADRULE_CUBIC, &
                             OMMP_RADTYPE_RMIN, & 
                             OMMP_EPSRULE_GEOMETRIC, &
                             OMMP_EPSRULE_HHG

    implicit none
    private

    integer(ip), parameter :: pair_allocation_chunk = 16

    type ommp_nonbonded_type
        !! Derived type for storing the information relative to the calculation
        !! of non-bonding interactions
        type(ommp_topology_type), pointer :: top
        !! Data structure for topology
        logical(lp) :: use_nl
        !! Flag for using neighbors list
        type(ommp_neigh_list) :: nl
        !! Neighbor list struture
        real(rp), allocatable, dimension(:) :: vdw_r
        !! VdW radii for the atoms of the system
        real(rp), allocatable, dimension(:) :: vdw_e
        !! Vdw energies for the atoms of fthe system
        real(rp), allocatable, dimension(:) :: vdw_f
        !! Scale factor for displacing the interaction center
    
        integer(ip) :: npair = 0
        !! Number of nonbonded pair parameters
        logical(lp), allocatable :: vdw_pair_mask_a(:,:), vdw_pair_mask_b(:,:)
        !! Mask to apply a pair parameter to a pair of atoms dimension is
        !! Natoms x Nparams
        real(rp), allocatable :: vdw_pair_r(:)
        !! Radii for the VdW atom pairs
        real(rp), allocatable :: vdw_pair_e(:)
        !! VdW energies for atom pairs
    
        real(rp), dimension(4) :: vdw_screening = [0.0, 0.0, 1.0, 1.0]
        !! Screening factors for 1,2 1,3 and 1,4 neighbours.
        !! Default vaules from tinker manual.
        real(rp) :: radf = 1.0
        !! Scal factor for atomic radii/diameters (1.0 is used for diameters,
        !! 2.0 for radii)
        integer(ip) :: radrule
        !! Flag for the radius rule to be used
        integer(ip) :: radtype
        !! Flag for the radius type
        integer(ip) :: vdwtype
        !! Flag for the vdW type
        integer(ip) :: epsrule
        !! Flag for eps rule !!TODO
    end type ommp_nonbonded_type
    
    abstract interface
    subroutine vdw_term(Rij, Rij0, Eij, V)
        use mod_memory, only: rp
        real(rp), intent(in) :: Rij
        real(rp), intent(in) :: Rij0
        real(rp), intent(in) :: Eij
        real(rp), intent(inout) :: V
    end subroutine
    end interface
   
    abstract interface
    subroutine vdw_gterm(Rij, Rij0, Eij, Rijgrad)
        use mod_memory, only: rp
        real(rp), intent(in) :: Rij
        real(rp), intent(in) :: Rij0
        real(rp), intent(in) :: Eij
        real(rp), intent(out) :: Rijgrad
    end subroutine
    end interface
   
    public :: ommp_nonbonded_type
    public :: vdw_init, vdw_terminate, vdw_set_pair, vdw_remove_potential
    public :: vdw_set_cutoff
    public :: vdw_potential, vdw_geomgrad
    public :: vdw_potential_inter, vdw_geomgrad_inter
    public :: vdw_potential_inter_restricted, vdw_geomgrad_inter_restricted

    contains

    subroutine vdw_init(vdw, top, vdw_type, radius_rule, radius_size, &
                        radius_type, epsrule, cutoff)
        !! Initialize the non-bonded object allocating the parameters vectors
        
        use mod_memory, only: mallocate
        use mod_io, only: fatal_error
        use mod_neighbor_list, only: nl_init
        use mod_constants, only: OMMP_DEFAULT_NL_SUB

        implicit none

        type(ommp_nonbonded_type), intent(inout) :: vdw
        type(ommp_topology_type), intent(in), target :: top
        character(len=*) :: vdw_type, radius_rule, radius_size, radius_type, &
                            epsrule
        real(rp) :: cutoff

        select case(trim(vdw_type))
            case("lennard-jones")
                vdw%vdwtype = OMMP_VDWTYPE_LJ
                continue
            case("buckingham")
                call fatal_error("VdW type buckingham is not implemented")
            case("buffered-14-7")
                vdw%vdwtype = OMMP_VDWTYPE_BUF714
                continue
            case("mm3-hbond")
                call fatal_error("VdW type mm3-hbond is not implemented")
            case("gaussian")
                call fatal_error("VdW type gaussian is not implemented")
            case default
                call fatal_error("VdW type specified is not understood")
        end select
        
        select case(trim(radius_rule))
            case("arithmetic")
                vdw%radrule = OMMP_RADRULE_ARITHMETIC
                continue
            case("geometric")
                call fatal_error("radiusrule geometric is not implemented")
            case("cubic-mean")
                vdw%radrule = OMMP_RADRULE_CUBIC
                continue
            case default
                call fatal_error("radiusrule specified is not understood")
        end select
        
        select case(trim(radius_size))
            case("radius")
                vdw%radf = 2.0
                continue
            case("diameter")
                vdw%radf = 1.0
                continue
            case default
                call fatal_error("radiussize specified is not understood")
        end select
        
        select case(trim(radius_type))
            case("sigma")
                call fatal_error("radiustype sigma is not implemented")
            case("r-min")
                vdw%radtype = OMMP_RADTYPE_RMIN
                continue
            case default
                call fatal_error("radiustype specified is not understood")
        end select
        
        select case(trim(epsrule))
            case("geometric")
                vdw%epsrule = OMMP_EPSRULE_GEOMETRIC
                continue
            case("arithmetic")
                call fatal_error("epsilonrule arithmetic is not implemented")
            case("harmonic")
                call fatal_error("epsilonrule harmonic is not implemented")
            case("w-h")
                call fatal_error("epsilonrule w-h is not implemented")
            case("hhg")
                vdw%epsrule = OMMP_EPSRULE_HHG
                continue
            case default
                call fatal_error("epsilonrule specified is not understood")
        end select

        vdw%top => top
        call mallocate('vdw_init [vdw_r]', top%mm_atoms, vdw%vdw_r)
        call mallocate('vdw_init [vdw_e]', top%mm_atoms, vdw%vdw_e)
        call mallocate('vdw_init [vdw_f]', top%mm_atoms, vdw%vdw_f)
        
        call mallocate('vdw_init [vdw_pair_mask_a]', &
                       top%mm_atoms, pair_allocation_chunk,vdw%vdw_pair_mask_a)
        call mallocate('vdw_init [vdw_pair_mask_b]', & 
                       top%mm_atoms, pair_allocation_chunk,vdw%vdw_pair_mask_b)
        call mallocate('vdw_init [vdw_pair_r]', pair_allocation_chunk, &
                       vdw%vdw_pair_r)
        call mallocate('vdw_init [vdw_pair_e]', pair_allocation_chunk, &
                       vdw%vdw_pair_e)

        vdw%vdw_f = 1.0_rp

        if(cutoff > 0.0) then
            vdw%use_nl = .true.
            if(vdw%use_nl) call nl_init(vdw%nl, top%cmm, cutoff, OMMP_DEFAULT_NL_SUB)
        else
            vdw%use_nl = .false.
        end if
    end subroutine vdw_init

    subroutine vdw_set_cutoff(vdw, cutoff, subdivision)
        use mod_neighbor_list, only: nl_terminate, nl_init
        implicit none

        type(ommp_nonbonded_type), intent(inout) :: vdw
        real(rp), intent(in) :: cutoff
        integer(ip), intent(in) :: subdivision

        if(vdw%use_nl) then
            ! This is re-initialization, old stuff should be cleaned.
            call nl_terminate(vdw%nl)
        end if
        if(cutoff > 0.0) then
            vdw%use_nl = .true.
            call nl_init(vdw%nl, vdw%top%cmm, cutoff, subdivision)
        else
            vdw%use_nl = .false.
        end if 
    end subroutine

    subroutine vdw_terminate(vdw)
        use mod_memory, only: mfree
        use mod_adjacency_mat, only: free_yale_sparse
        use mod_neighbor_list, only: nl_terminate

        implicit none
        
        type(ommp_nonbonded_type), intent(inout) :: vdw
        
        call mfree('vdw_terminate [vdw_r]', vdw%vdw_r)
        call mfree('vdw_terminate [vdw_e]', vdw%vdw_e)
        call mfree('vdw_terminate [vdw_f]', vdw%vdw_f)
        call mfree('vdw_terminate [vdw_pair_r]', vdw%vdw_pair_r)
        call mfree('vdw_terminate [vdw_pair_e]', vdw%vdw_pair_e)
        call mfree('vdw_terminate [vdw_pair_mask_b]', vdw%vdw_pair_mask_a)
        call mfree('vdw_terminate [vdw_pair_mask_b]', vdw%vdw_pair_mask_b)
        if(vdw%use_nl) call nl_terminate(vdw%nl)

    end subroutine

    subroutine vdw_remove_potential(vdw, i)
        !! Remove the VdW interaction from the specified atom
        !! the atom will not interact anymore with any other atom
        implicit none

        type(ommp_nonbonded_type), intent(inout) :: vdw
        integer(ip), intent(in) :: i

        vdw%vdw_e(i) = 0.0
        
        ! TODO do it also for pair interactions
    end subroutine

    subroutine vdw_set_pair(vdw, mask_a, mask_b, r, e)
        use mod_io, only: ommp_message, fatal_error
        use mod_constants, only: OMMP_VERBOSE_LOW
        use mod_memory, only: mallocate, mfree

        implicit none

        type(ommp_nonbonded_type), intent(inout) :: vdw
        !! Nonbonded data structure
        logical(lp), intent(in) :: mask_a(vdw%top%mm_atoms)
        logical(lp), intent(in) :: mask_b(vdw%top%mm_atoms)
        real(rp), intent(in) :: r
        !! Equilibrium distance for the pair
        real(rp), intent(in) :: e 
        !! Depth of the potential 

        integer(ip) :: oldsize, newsize
        logical(lp), allocatable :: ltmp(:,:)
        real(rp), allocatable :: rtmp(:)
        character(len=OMMP_STR_CHAR_MAX) :: msg

        oldsize = size(vdw%vdw_pair_r)
        if(vdw%npair >= oldsize) then
            newsize = oldsize + pair_allocation_chunk
            call mallocate('vdw_set_pair [rtmp]', oldsize, rtmp)
            call mallocate('vdw_set_pair [ltmp]',vdw%top%mm_atoms, oldsize, ltmp)
            
            rtmp = vdw%vdw_pair_r
            call mfree('vdw_set_pair [vdw%vdw_pair_r]', vdw%vdw_pair_r)
            call mallocate('vdw_set_pair [vdw%vdw_pair_r]', newsize, vdw%vdw_pair_r)
            vdw%vdw_pair_r(1:oldsize) = rtmp
            
            rtmp = vdw%vdw_pair_e
            call mfree('vdw_set_pair [vdw%vdw_pair_e]', vdw%vdw_pair_e)
            call mallocate('vdw_set_pair [vdw%vdw_pair_e]', newsize, vdw%vdw_pair_e)
            vdw%vdw_pair_e(1:oldsize) = rtmp

            ltmp = vdw%vdw_pair_mask_a
            call mfree('vdw_set_pair [vdw_pair_mask_a]', vdw%vdw_pair_mask_a)
            call mallocate('vdw_set_pair [vdw_pair_mask_a]', &
                           vdw%top%mm_atoms, newsize, vdw%vdw_pair_mask_a)
            vdw%vdw_pair_mask_a(:,1:oldsize) = ltmp
            
            ltmp = vdw%vdw_pair_mask_b
            call mfree('vdw_set_pair [vdw_pair_mask_b]', vdw%vdw_pair_mask_b)
            call mallocate('vdw_set_pair [vdw_pair_mask_b]', &
                           vdw%top%mm_atoms, newsize, vdw%vdw_pair_mask_b)
            vdw%vdw_pair_mask_b(:,1:oldsize) = ltmp

            call mfree('vdw_set_pair [rtmp]', rtmp)
            call mfree('vdw_set_pair [ltmp]', ltmp)
        end if
        vdw%npair = vdw%npair + 1
        vdw%vdw_pair_mask_a(:,vdw%npair) = mask_a
        vdw%vdw_pair_mask_b(:,vdw%npair) = mask_b
        vdw%vdw_pair_r(vdw%npair) = r
        vdw%vdw_pair_e(vdw%npair) = e

    end subroutine vdw_set_pair

    subroutine vdw_lennard_jones(Rij, Rij0, Eij, V)
        implicit none

        real(rp), intent(in) :: Rij
        real(rp), intent(in) :: Rij0
        real(rp), intent(in) :: Eij
        real(rp), intent(inout) :: V

        real(rp) :: sigma_ov_r

        sigma_ov_r = Rij0 / Rij
        V = V + Eij*(sigma_ov_r**12 - 2_rp*sigma_ov_r**6)

    end subroutine vdw_lennard_jones
    
    subroutine vdw_lennard_jones_Rijgrad(Rij, Rij0, Eij, Rijgrad)
        implicit none

        real(rp), intent(in) :: Rij
        real(rp), intent(in) :: Rij0
        real(rp), intent(in) :: Eij
        real(rp), intent(out) :: Rijgrad

        real(rp) :: sigma_ov_r

        sigma_ov_r = Rij0 / Rij
        Rijgrad = -12.0 * Eij * (sigma_ov_r ** 12 - sigma_ov_r ** 6) / Rij

    end subroutine vdw_lennard_jones_Rijgrad
    
    subroutine vdw_buffered_7_14(Rij, Rij0, Eij, V)
        !! Compute the dispersion-repulsion energy using the buffered 7-14 
        !! potential. Details can be found in ref: 10.1021/jp027815
        
        implicit none

        real(rp), intent(in) :: Rij
        real(rp), intent(in) :: Rij0
        real(rp), intent(in) :: Eij
        real(rp), intent(inout) :: V

        real(rp) :: rho
        real(rp), parameter :: delta = 0.07, gam = 0.12
        integer(ip), parameter :: n = 14, m = 7

        rho = Rij / Rij0
        V = V + Eij*((1+delta)/(rho+delta))**(n-m) * ((1+gam)/(rho**m+gam)-2)

    end subroutine vdw_buffered_7_14
    
    subroutine vdw_buffered_7_14_Rijgrad(Rij, Rij0, Eij, Rijgrad)
        !! Compute the gradient of vdw energy (using the buffered 7-14 
        !! potential ref: 10.1021/jp027815) with respect to the distance 
        !! between the two atoms (Rij).
        
        implicit none

        real(rp), intent(in) :: Rij
        real(rp), intent(in) :: Rij0
        real(rp), intent(in) :: Eij
        real(rp), intent(out) :: Rijgrad

        real(rp) :: rho, rhom
        real(rp), parameter :: delta = 0.07, gam = 0.12
        integer(ip), parameter :: n = 14, m = 7

        rho = Rij / Rij0
        rhom = rho**m
        Rijgrad = Eij * ((delta+1)/(delta+rho))**(n-m) * &
                  (rho*(gam+rhom)*(m-n)*(1-2*rhom-gam) - &
                  rhom*m*(delta+rho)*(gam+1)) / &
                  (rho*(delta+rho)*(gam+rhom)**2)
        Rijgrad = Rijgrad / Rij0

    end subroutine vdw_buffered_7_14_Rijgrad

    pure function get_Rij0(vdw, i, j) result(Rij0)
        use mod_constants, only: eps_rp
        
        implicit none

        type(ommp_nonbonded_type), intent(in), target :: vdw
        !! Nonbonded data structure
        integer(ip), intent(in) :: i, j
        !! Indices of interacting atoms
        real(rp) :: Rij0
        
        if(vdw%vdw_r(i) < eps_rp .and. vdw%vdw_r(j) < eps_rp) then
            Rij0 = 0.0
            return
        end if

        select case(vdw%radrule) 
            case(OMMP_RADRULE_ARITHMETIC)
                Rij0 = vdw%radf*(vdw%vdw_r(i) + vdw%vdw_r(j))/2
            case(OMMP_RADRULE_CUBIC)
                Rij0 = (vdw%vdw_r(i)**3 + vdw%vdw_r(j)**3) / &
                       (vdw%vdw_r(i)**2 + vdw%vdw_r(j)**2)
            case default
                Rij0 = 0.0
        end select
    end function
    
    pure function get_eij(vdw, i, j) result(eij)
        use mod_constants, only: eps_rp
        
        implicit none

        type(ommp_nonbonded_type), intent(in) :: vdw
        !! Nonbonded data structure
        integer(ip), intent(in) :: i, j
        !! Indices of interacting atoms
        real(rp) :: eij

        if(vdw%vdw_e(i) < eps_rp .and. vdw%vdw_e(j) < eps_rp) then
            eij = 0.0
            return
        end if
        
        select case(vdw%epsrule) 
            case(OMMP_EPSRULE_GEOMETRIC)
                eij = sqrt(vdw%vdw_e(i)*vdw%vdw_e(j))
            case(OMMP_EPSRULE_HHG)
                eij = (4*vdw%vdw_e(i)*vdw%vdw_e(j)) / &
                      (vdw%vdw_e(i)**0.5 + vdw%vdw_e(j)**0.5)**2
            case default
                eij = 0.0
        end select
    end function

    pure function get_Rij0_inter(vdw1, vdw2, i, j) result(Rij0)
        use mod_constants, only: eps_rp
        
        implicit none

        type(ommp_nonbonded_type), intent(in), target :: vdw1, vdw2
        !! Nonbonded data structure
        integer(ip), intent(in) :: i, j
        !! Indices of interacting atoms
        real(rp) :: Rij0

        if(abs(vdw1%radrule - vdw2%radrule) > eps_rp .or. &
           abs(vdw1%radf - vdw2%radf) > eps_rp) then
           Rij0 = 0.0
           return
        end if
        
        if(vdw1%vdw_r(i) < eps_rp .and. vdw2%vdw_r(j) < eps_rp) then
            Rij0 = 0.0
            return
        end if

        select case(vdw1%radrule) 
            case(OMMP_RADRULE_ARITHMETIC)
                Rij0 = vdw1%radf*(vdw1%vdw_r(i) + vdw2%vdw_r(j))/2
            case(OMMP_RADRULE_CUBIC)
                Rij0 = (vdw1%vdw_r(i)**3 + vdw2%vdw_r(j)**3) / &
                       (vdw1%vdw_r(i)**2 + vdw2%vdw_r(j)**2)
            case default
                Rij0 = 0.0
        end select
    end function
    
    pure function get_eij_inter(vdw1, vdw2, i, j) result(eij)
        use mod_constants, only: eps_rp
        
        implicit none

        type(ommp_nonbonded_type), intent(in) :: vdw1, vdw2
        !! Nonbonded data structure
        integer(ip), intent(in) :: i, j
        !! Indices of interacting atoms
        real(rp) :: eij

        if(vdw1%epsrule /= vdw2%epsrule) then
            eij = 0.0
            return
        end if
        
        if(vdw1%vdw_e(i) < eps_rp .and. vdw2%vdw_e(j) < eps_rp) then
            eij = 0.0
            return
        end if

        select case(vdw1%epsrule) 
            case(OMMP_EPSRULE_GEOMETRIC)
                eij = sqrt(vdw1%vdw_e(i)*vdw2%vdw_e(j))
            case(OMMP_EPSRULE_HHG)
                eij = (4*vdw1%vdw_e(i)*vdw2%vdw_e(j)) / &
                      (vdw1%vdw_e(i)**0.5 + vdw2%vdw_e(j)**0.5)**2
            case default
                eij = 0.0
        end select
    end function

    subroutine vdw_potential(vdw, V)
        !! Compute the dispersion repulsion energy for the whole system
        !! using a double loop algorithm

        use mod_io, only : fatal_error
        use mod_profiling, only: time_push, time_pull
        use mod_constants, only: eps_rp
        use mod_memory, only: mallocate, mfree
        use mod_neighbor_list, only: get_ith_nl
        implicit none

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

        integer(ip) :: i, j, jc, l, ipair, ineigh, nthreads, ithread, nn
        real(rp) :: eij, rij0, rij, ci(3), cj(3), s, vtmp
        type(ommp_topology_type), pointer :: top
        procedure(vdw_term), pointer :: vdw_func

        integer(ip), allocatable :: nl_neigh(:,:)
        real(rp), allocatable :: nl_r(:,:)

        integer :: omp_get_num_threads, omp_get_thread_num
        
        call time_push()
        !$omp parallel
        nthreads = omp_get_num_threads()
        !$omp end parallel

        top => vdw%top
        select case(vdw%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")
        end select

        if(vdw%use_nl) then
            call mallocate('vdw_potential [rneigh]', top%mm_atoms, nthreads, nl_r)
            call mallocate('vdw_potential [nl_neigh]', top%mm_atoms, nthreads, nl_neigh)
        end if

        !$omp parallel do default(shared) reduction(+:v)  schedule(dynamic) &
        !$omp private(i,j,jc,ineigh,ithread,nn,s,ci,cj,ipair,l,Eij,Rij0,Rij,vtmp)
        do i=1, top%mm_atoms
            ithread = omp_get_thread_num() + 1
            if(abs(vdw%vdw_f(i) - 1.0_rp) < eps_rp) then
                ci = top%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(top%conn(1)%ri(i+1) - top%conn(1)%ri(i) /= 1) then
                    call fatal_error("Scale factors are only expected for &
                                     &monovalent atoms")
                end if
                ineigh = top%conn(1)%ci(top%conn(1)%ri(i))

                ci = top%cmm(:,ineigh) + (top%cmm(:,i) - top%cmm(:,ineigh)) &
                     * vdw%vdw_f(i)
            endif
            
            ! If neighbor list are enabled get the one for the current
            if(vdw%use_nl) call get_ith_nl(vdw%nl, i, top%cmm, &
                                           nl_neigh(:,ithread), &
                                           nl_r(:,ithread), nn)

            do jc=1, top%mm_atoms
                ! If the two atoms aren't neighbors, just skip the loop
                if(vdw%use_nl) then
                    if(jc > nn) exit !! All neighbors done!
                    j = nl_neigh(jc,ithread)
                    if(j <= i) cycle
                else
                    ! Skip all iteration with j <= i
                    if(jc > i) then
                        j = jc
                    else
                        cycle
                    end if
                end if
                ! Compute the screening factor for this pair
                s = 1.0_rp
                do ineigh=1,4
                    ! Look if j is at distance ineigh from i
                    if(any(top%conn(ineigh)%ci(top%conn(ineigh)%ri(i): &
                                               top%conn(ineigh)%ri(i+1)-1) == j)) then
                       
                        s = vdw%vdw_screening(ineigh)
                        ! Exit the loop
                        exit 
                    end if
                end do
                
                if(s > eps_rp) then
                    ipair = -1
                    do l=1, vdw%npair
                        if((vdw%vdw_pair_mask_a(i,l) .and. vdw%vdw_pair_mask_b(j,l)) .or. &
                           (vdw%vdw_pair_mask_a(j,l) .and. vdw%vdw_pair_mask_b(i,l))) then
                            ipair = l
                            exit
                        end if
                    end do

                    if(ipair > 0) then
                        Rij0 = vdw%vdw_pair_r(ipair)
                        eij = vdw%vdw_pair_e(ipair)
                    else
                        Rij0 = get_Rij0(vdw, i, j)
                        eij = get_eij(vdw, i, j)
                    end if
                
                    if(abs(vdw%vdw_f(j) - 1.0_rp) < eps_rp) then
                        cj = top%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(top%conn(1)%ri(j+1) - top%conn(1)%ri(j) /= 1) then
                            call fatal_error("Scale factors are only expected &
                                             & for monovalent atoms")
                        end if
                        ineigh = top%conn(1)%ci(top%conn(1)%ri(j))

                        cj = top%cmm(:,ineigh) + &
                             (top%cmm(:,j) - top%cmm(:,ineigh)) * vdw%vdw_f(j)
                    endif
                    Rij = norm2(ci-cj)
                    if(Rij < eps_rp) then
                        call fatal_error("Requesting 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*s
                end if
            end do
        end do
        
        if(vdw%use_nl) then
            call mfree('vdw_potential [rneigh]', nl_r)
            call mfree('vdw_potential [nl_neigh]', nl_neigh)
        end if
        call time_pull('VdW potential calculation')
    end subroutine vdw_potential
    
    subroutine vdw_geomgrad(vdw, grad)
        !! Compute the dispersion repulsion geometric gradients for the whole system
        !! using a double loop algorithm

        use mod_io, only : fatal_error
        use mod_constants, only: eps_rp
        use mod_jacobian_mat, only: Rij_jacobian
        use mod_profiling, only: time_push, time_pull
        use mod_memory, only: mallocate, mfree
        use mod_neighbor_list, only: get_ith_nl
        implicit none

        type(ommp_nonbonded_type), intent(in), target :: vdw
        !! Nonbonded data structure
        real(rp), intent(inout) :: grad(3,vdw%top%mm_atoms)
        !! Gradients, result will be added

        integer(ip) :: i, j, l, ipair, ineigh, ineigh_i, ineigh_j, jc, &
                       nn, ithread, nthreads
        real(rp) :: eij, rij0, rij, ci(3), cj(3), s, J_i(3), J_j(3), Rijg, &
                    f_i, f_j
        logical :: skip
        type(ommp_topology_type), pointer :: top
        procedure(vdw_gterm), pointer :: vdw_gfunc
        
        integer(ip), allocatable :: nl_neigh(:,:)
        real(rp), allocatable :: nl_r(:,:)
        
        integer :: omp_get_num_threads, omp_get_thread_num
        
        call time_push()
        !$omp parallel
        nthreads = omp_get_num_threads()
        !$omp end parallel

        top => vdw%top
        select case(vdw%vdwtype)
            case(OMMP_VDWTYPE_LJ)
                vdw_gfunc => vdw_lennard_jones_Rijgrad
            case(OMMP_VDWTYPE_BUF714)
                vdw_gfunc => vdw_buffered_7_14_Rijgrad
            case default
                vdw_gfunc => vdw_buffered_7_14_Rijgrad
                call fatal_error("Unexpected error in vdw_geomgrad")
        end select
        
        if(vdw%use_nl) then
            call mallocate('vdw_geomgrad [rneigh]', top%mm_atoms, nthreads, nl_r)
            call mallocate('vdw_geomgrad [nl_neigh]', top%mm_atoms, nthreads, nl_neigh)
        end if

        !$omp parallel do default(shared) schedule(dynamic) &
        !$omp private(i,j,ci,cj,ineigh,ineigh_i,ineigh_j,f_i,f_j,s,ipair,l) &
        !$omp private(Eij,Rij0,Rijg,Rij,J_i,J_j,skip,jc,nn,ithread)
        do i=1, top%mm_atoms
            ithread = omp_get_thread_num() + 1
            if(abs(vdw%vdw_f(i) - 1.0) < eps_rp) then
                ci = top%cmm(:,i)
                ineigh_i = 0 ! This is needed later for force projection
                f_i = 1.0
            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(top%conn(1)%ri(i+1) - top%conn(1)%ri(i) /= 1) then
                    call fatal_error("Scale factors are only expected for &
                                     &monovalent atoms")
                end if
                ineigh_i = top%conn(1)%ci(top%conn(1)%ri(i))
                f_i = vdw%vdw_f(i)

                ci = top%cmm(:,ineigh_i) + (top%cmm(:,i) - &
                                            top%cmm(:,ineigh_i)) * f_i
            endif
                
            if(vdw%use_nl) call get_ith_nl(vdw%nl, i, top%cmm, &
                                           nl_neigh(:,ithread), &
                                           nl_r(:,ithread), nn)

            do jc=1, top%mm_atoms
                ! If the two atoms aren't neighbors, just skip the loop
                if(vdw%use_nl) then
                    if(jc > nn) exit !! All neighbors done!
                    j = nl_neigh(jc,ithread)
                    if(j <= i) cycle
                else
                    ! Skip all iteration with j <= i
                    if(jc > i) then
                        j = jc
                    else
                        cycle
                    end if
                end if
                ! Compute the screening factor for this pair
                s = 1.0_rp
                do ineigh=1,4
                    ! Look if j is at distance ineigh from i
                    if(any(top%conn(ineigh)%ci(top%conn(ineigh)%ri(i): &
                                               top%conn(ineigh)%ri(i+1)-1) == j)) then
                       
                        s = vdw%vdw_screening(ineigh)
                        ! Exit the loop
                        exit 
                    end if
                end do
                
                if(s > eps_rp) then
                    if(abs(vdw%vdw_f(j) - 1.0) < eps_rp) then
                        cj = top%cmm(:,j)
                        ineigh_j = 0 ! This is needed later for force projection
                        f_j = 1.0
                    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(top%conn(1)%ri(j+1) - top%conn(1)%ri(j) /= 1) then
                            call fatal_error("Scale factors are only expected &
                                             & for monovalent atoms")
                        end if
                        ineigh_j = top%conn(1)%ci(top%conn(1)%ri(j))
                        f_j = vdw%vdw_f(j)

                        cj = top%cmm(:,ineigh_j) + &
                             (top%cmm(:,j) - top%cmm(:,ineigh_j)) * f_j
                    endif
                    
                    ! if all atoms in the interaction are frozen 
                    ! just skip to next iteration
                    if(top%use_frozen) then
                        skip = .true.
                        skip = skip .and. top%frozen(i)
                        skip = skip .and. top%frozen(j)
                        if(ineigh_i > 0) skip = skip .and. top%frozen(ineigh_i)
                        if(ineigh_j > 0) skip = skip .and. top%frozen(ineigh_j)
                        if(skip) cycle
                    end if

                    ipair = -1
                    do l=1, vdw%npair
                        if((vdw%vdw_pair_mask_a(i,l) .and. vdw%vdw_pair_mask_b(j,l)) .or. &
                           (vdw%vdw_pair_mask_a(j,l) .and. vdw%vdw_pair_mask_b(i,l))) then
                            ipair = l
                            exit
                        end if
                    end do

                    if(ipair > 0) then
                        Rij0 = vdw%vdw_pair_r(ipair)
                        eij = vdw%vdw_pair_e(ipair)
                    else 
                        Rij0 = get_Rij0(vdw, i, j)
                        eij = get_eij(vdw, i, j)
                    end if

                    call Rij_jacobian(ci, cj, Rij, J_i, J_j)
                    if(Rij < eps_rp) then
                        call fatal_error("Requesting non-bonded gradients 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
                    call vdw_gfunc(Rij, Rij0, Eij, Rijg)

                    Rijg = Rijg * s

                    if(ineigh_i == 0) then
                        if(.not. (top%use_frozen .and. top%frozen(i))) then
                            !$omp atomic update
                            grad(1,i) =  grad(1,i) + J_i(1) * Rijg
                            !$omp atomic update
                            grad(2,i) =  grad(2,i) + J_i(2) * Rijg
                            !$omp atomic update
                            grad(3,i) =  grad(3,i) + J_i(3) * Rijg
                        end if
                    else
                        ! If the center is displaced, the forces should be 
                        ! projected onto the two atoms that determine the
                        ! position of the center
                        if(.not. (top%use_frozen .and. top%frozen(i))) then
                            !$omp atomic update
                            grad(1,i) = grad(1,i) + J_i(1) * Rijg * f_i
                            !$omp atomic update
                            grad(2,i) = grad(2,i) + J_i(2) * Rijg * f_i
                            !$omp atomic update
                            grad(3,i) = grad(3,i) + J_i(3) * Rijg * f_i
                        end if
                        if(.not. (top%use_frozen .and. top%frozen(ineigh_i))) then
                            !$omp atomic update
                            grad(1,ineigh_i) = grad(1,ineigh_i) + J_i(1) * Rijg * (1-f_i)
                            !$omp atomic update
                            grad(2,ineigh_i) = grad(2,ineigh_i) + J_i(2) * Rijg * (1-f_i)
                            !$omp atomic update
                            grad(3,ineigh_i) = grad(3,ineigh_i) + J_i(3) * Rijg * (1-f_i)
                        end if
                    end if

                    if(ineigh_j == 0) then
                        if(.not. (top%use_frozen .and. top%frozen(j))) then
                            !$omp atomic update
                            grad(1,j) =  grad(1,j) + J_j(1) * Rijg
                            !$omp atomic update
                            grad(2,j) =  grad(2,j) + J_j(2) * Rijg
                            !$omp atomic update
                            grad(3,j) =  grad(3,j) + J_j(3) * Rijg
                        end if
                    else
                        ! If the center is displaced, the forces should be 
                        ! projected onto the two atoms that determine the
                        ! position of the center
                        if(.not. (top%use_frozen .and. top%frozen(j))) then
                            !$omp atomic update
                            grad(1,j) = grad(1,j) + J_j(1) * Rijg * f_j
                            !$omp atomic update
                            grad(2,j) = grad(2,j) + J_j(2) * Rijg * f_j
                            !$omp atomic update
                            grad(3,j) = grad(3,j) + J_j(3) * Rijg * f_j
                        end if
                        if(.not. (top%use_frozen .and. top%frozen(ineigh_j))) then
                            !$omp atomic update
                            grad(1,ineigh_j) = grad(1,ineigh_j) + J_j(1) * Rijg * (1-f_j)
                            !$omp atomic update
                            grad(2,ineigh_j) = grad(2,ineigh_j) + J_j(2) * Rijg * (1-f_j)
                            !$omp atomic update
                            grad(3,ineigh_j) = grad(3,ineigh_j) + J_j(3) * Rijg * (1-f_j)
                        end if
                    endif
                end if
            end do
        end do
        call time_pull('VdW gradients calculation')
    end 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
    
    subroutine vdw_geomgrad_inter(vdw1, vdw2, grad1, grad2)
        !! 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
        use mod_jacobian_mat, only: Rij_jacobian
        implicit none

        type(ommp_nonbonded_type), intent(in), target :: vdw1, vdw2
        !! Nonbonded data structure
        real(rp), intent(inout) :: grad1(3,vdw1%top%mm_atoms), &
                                   grad2(3,vdw2%top%mm_atoms)
        !! Potential, result will be added

        integer(ip) :: i, j, ineigh_i, ineigh_j
        real(rp) :: eij, rij0, rij, ci(3), cj(3), Rijg, f_i, f_j, &
                    J_i(3), J_j(3)
        logical :: skip
        type(ommp_topology_type), pointer :: top1, top2
        procedure(vdw_gterm), pointer :: vdw_grad

        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_grad => vdw_lennard_jones_Rijgrad
            case(OMMP_VDWTYPE_BUF714)
                vdw_grad => vdw_buffered_7_14_Rijgrad
            case default
                vdw_grad => vdw_buffered_7_14_Rijgrad
                call fatal_error("Unexpected error in vdw_geomgrad_inter")
        end select
        
        !$omp parallel do default(shared) schedule(dynamic) &
        !$omp private(i,j,ci,cj,f_i,f_j,ineigh_i,ineigh_j,Eij,Rij0,Rij,Rijg,J_i,J_j,skip)
        do i=1, top1%mm_atoms
            if(abs(vdw1%vdw_f(i) - 1.0) < eps_rp) then
                ci = top1%cmm(:,i)
                ineigh_i = 0
                f_i = 1.0
            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")
                end if
                ineigh_i = top1%conn(1)%ci(top1%conn(1)%ri(i))
                f_i = vdw1%vdw_f(i)
                ci = top1%cmm(:,ineigh_i) + (top1%cmm(:,i) - top1%cmm(:,ineigh_i)) &
                     * f_i
            endif
                
            do j=1, top2%mm_atoms
                if(abs(vdw2%vdw_f(j) - 1.0) < eps_rp) then
                    cj = top2%cmm(:,j)
                    ineigh_j = 0
                    f_j = 1.0
                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")
                    end if
                    ineigh_j = top2%conn(1)%ci(top2%conn(1)%ri(j))
                    f_j = vdw2%vdw_f(j)

                    cj = top2%cmm(:,ineigh_j) + &
                         (top2%cmm(:,j) - top2%cmm(:,ineigh_j)) * f_j
                endif
                
                skip = top1%use_frozen .and. top2%use_frozen
                if(skip .and. top1%use_frozen) then
                    skip = skip .and. top1%frozen(i)
                    if(ineigh_i > 0) skip = skip .and. top1%frozen(ineigh_i)
                end if
                if(skip .and. top2%use_frozen) then
                    skip = skip .and. top2%frozen(j)
                    if(ineigh_j > 0) skip = skip .and. top2%frozen(ineigh_j)
                end if
                if(skip) cycle
                
                Rij0 = get_Rij0_inter(vdw1, vdw2, i, j)
                eij = get_eij_inter(vdw1, vdw2, i, j)
            
                call Rij_jacobian(ci, cj, Rij, J_i, J_j)
                if(Rij < eps_rp) then
                    call fatal_error("Requesting inter non-bonded gradients 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
                call vdw_grad(Rij, Rij0, Eij, Rijg)

                !$omp critical
                if(ineigh_i == 0) then
                    if(.not. (top1%use_frozen .and. top1%frozen(i))) &
                        grad1(:,i) =  grad1(:,i) + J_i * Rijg
                else
                    ! If the center is displaced, the forces should be 
                    ! projected onto the two atoms that determine the
                    ! position of the center
                    if(.not. (top1%use_frozen .and. top1%frozen(i))) &
                        grad1(:,i) = grad1(:,i) + J_i * Rijg * f_i
                    if(.not. (top1%use_frozen .and. top1%frozen(ineigh_i))) &
                        grad1(:,ineigh_i) = grad1(:,ineigh_i) + J_i * Rijg * (1-f_i)
                end if

                if(ineigh_j == 0) then
                    if(.not. (top2%use_frozen .and. top2%frozen(j))) &
                        grad2(:,j) =  grad2(:,j) + J_j * Rijg
                else
                    ! If the center is displaced, the forces should be 
                    ! projected onto the two atoms that determine the
                    ! position of the center
                    if(.not. (top2%use_frozen .and. top2%frozen(j))) &
                        grad2(:,j) = grad2(:,j) + J_j * Rijg * f_j
                    if(.not. (top2%use_frozen .and. top2%frozen(ineigh_j))) &
                        grad2(:,ineigh_j) = grad2(:,ineigh_j) + J_j * Rijg * (1-f_j)
                endif
                !$omp end critical
            end do
        end do
    end subroutine
    
    subroutine vdw_potential_inter_restricted(vdw1, vdw2, pairs, s, n, V)
        !! Compute the dispersion repulsion energy between two systems vdw1
        !! and vdw2 accounting only for the pairs pairs(1,i)--pairs(2,i) and scaling 
        !! each interaction by s(i).

        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
        integer(ip), intent(in) :: n
        !! number of pairs for which the interaction should be computed
        integer(ip), intent(in) :: pairs(2,n)
        !! pairs of atoms for which the interaction should be computed
        real(rp), intent(in) :: s(:)
        !! scaling factors for each interaction
        real(rp), intent(inout) :: V
        !! Potential, result will be added

        integer(ip) :: i, j, ineigh, ip
        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_restricted")
        end select

        !$omp parallel do default(shared) schedule(dynamic) &
        !$omp private(i,j,ip,ci,cj,ineigh,Rij,Rij0,Eij,vtmp) reduction(+:v) 
        do ip=1, n
            i = pairs(1,ip)
            j = pairs(2,ip)

            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")
                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
                
            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")
                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 * s(ip)
        end do
    end subroutine vdw_potential_inter_restricted
    
    subroutine vdw_geomgrad_inter_restricted(vdw1, vdw2, pairs, s, n, &
                                             grad1, grad2)
        !! 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
        use mod_jacobian_mat, only: Rij_jacobian
        implicit none

        type(ommp_nonbonded_type), intent(in), target :: vdw1, vdw2
        !! Nonbonded data structure
        integer(ip), intent(in) :: n
        !! number of pairs for which the interaction should be computed
        integer(ip), intent(in) :: pairs(2,n)
        !! pairs of atoms for which the interaction should be computed
        real(rp), intent(in) :: s(:)
        !! scaling factors for each interaction
        real(rp), intent(inout) :: grad1(3,vdw1%top%mm_atoms), &
                                   grad2(3,vdw2%top%mm_atoms)
        !! Potential, result will be added

        integer(ip) :: i, j, ineigh_i, ineigh_j, ip
        real(rp) :: eij, rij0, rij, ci(3), cj(3), Rijg, f_i, f_j, &
                    J_i(3), J_j(3)
        logical :: skip
        type(ommp_topology_type), pointer :: top1, top2
        procedure(vdw_gterm), pointer :: vdw_grad

        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_grad => vdw_lennard_jones_Rijgrad
            case(OMMP_VDWTYPE_BUF714)
                vdw_grad => vdw_buffered_7_14_Rijgrad
            case default
                vdw_grad => vdw_buffered_7_14_Rijgrad
                call fatal_error("Unexpected error in vdw_geomgrad_inter_restricted")
        end select
        
        do ip=1, n
            i = pairs(1,ip)
            j = pairs(2,ip)

            if(abs(vdw1%vdw_f(i) - 1.0) < eps_rp) then
                ci = top1%cmm(:,i)
                ineigh_i = 0
            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")
                end if
                ineigh_i = top1%conn(1)%ci(top1%conn(1)%ri(i))
                f_i = vdw1%vdw_f(i)
                ci = top1%cmm(:,ineigh_i) + (top1%cmm(:,i) - top1%cmm(:,ineigh_i)) &
                     * f_i
            endif
                
            if(abs(vdw2%vdw_f(j) - 1.0) < eps_rp) then
                cj = top2%cmm(:,j)
                ineigh_j = 0
            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")
                end if
                ineigh_j = top2%conn(1)%ci(top2%conn(1)%ri(j))
                f_j = vdw2%vdw_f(j)

                cj = top2%cmm(:,ineigh_j) + &
                        (top2%cmm(:,j) - top2%cmm(:,ineigh_j)) * f_j
            endif
                
            skip = top1%use_frozen .and. top2%use_frozen
            if(skip .and. top1%use_frozen) then
                skip = skip .and. top1%frozen(i)
                if(ineigh_i > 0) skip = skip .and. top1%frozen(ineigh_i)
            end if
            if(skip .and. top2%use_frozen) then
                skip = skip .and. top2%frozen(j)
                if(ineigh_j > 0) skip = skip .and. top2%frozen(ineigh_j)
            end if
            if(skip) cycle
                
            Rij0 = get_Rij0_inter(vdw1, vdw2, i, j)
            eij = get_eij_inter(vdw1, vdw2, i, j) * s(ip)
            
            call Rij_jacobian(ci, cj, Rij, J_i, J_j)
            if(Rij < eps_rp) then
                call fatal_error("Requesting inter non-bonded gradients 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
            call vdw_grad(Rij, Rij0, Eij, Rijg)

            if(ineigh_i == 0) then
                if(.not. (top1%use_frozen .and. top1%frozen(i))) &
                    grad1(:,i) =  grad1(:,i) + J_i * Rijg
            else
                ! If the center is displaced, the forces should be 
                ! projected onto the two atoms that determine the
                ! position of the center
                if(.not. (top1%use_frozen .and. top1%frozen(i))) &
                    grad1(:,i) = grad1(:,i) + J_i * Rijg * f_i
                if(.not. (top1%use_frozen .and. top1%frozen(ineigh_i))) &
                    grad1(:,ineigh_i) = grad1(:,ineigh_i) + J_i * Rijg * (1-f_i)
            end if

            if(ineigh_j == 0) then
                if(.not. (top2%use_frozen .and. top2%frozen(j))) &
                    grad2(:,j) =  grad2(:,j) + J_j * Rijg
            else
                ! If the center is displaced, the forces should be 
                ! projected onto the two atoms that determine the
                ! position of the center
                if(.not. (top2%use_frozen .and. top2%frozen(j))) &
                    grad2(:,j) = grad2(:,j) + J_j * Rijg * f_j
                if(.not. (top2%use_frozen .and. top2%frozen(ineigh_j))) &
                    grad2(:,ineigh_j) = grad2(:,ineigh_j) + J_j * Rijg * (1-f_j)
            endif
        end do
    end subroutine
    
end module mod_nonbonded