mod_electrostatics.f90 Source File


This file depends on

sourcefile~~mod_electrostatics.f90~~EfferentGraph sourcefile~mod_electrostatics.f90 mod_electrostatics.f90 sourcefile~mod_io.f90 mod_io.f90 sourcefile~mod_electrostatics.f90->sourcefile~mod_io.f90 sourcefile~mod_profiling.f90 mod_profiling.f90 sourcefile~mod_electrostatics.f90->sourcefile~mod_profiling.f90 sourcefile~mod_adjacency_mat.f90 mod_adjacency_mat.f90 sourcefile~mod_electrostatics.f90->sourcefile~mod_adjacency_mat.f90 sourcefile~mod_constants.f90 mod_constants.f90 sourcefile~mod_electrostatics.f90->sourcefile~mod_constants.f90 sourcefile~mod_memory.f90 mod_memory.f90 sourcefile~mod_electrostatics.f90->sourcefile~mod_memory.f90 sourcefile~mod_topology.f90 mod_topology.f90 sourcefile~mod_electrostatics.f90->sourcefile~mod_topology.f90 sourcefile~mod_io.f90->sourcefile~mod_constants.f90 sourcefile~mod_profiling.f90->sourcefile~mod_io.f90 sourcefile~mod_profiling.f90->sourcefile~mod_constants.f90 sourcefile~mod_profiling.f90->sourcefile~mod_memory.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_memory.f90->sourcefile~mod_io.f90 sourcefile~mod_memory.f90->sourcefile~mod_constants.f90 sourcefile~mod_topology.f90->sourcefile~mod_io.f90 sourcefile~mod_topology.f90->sourcefile~mod_adjacency_mat.f90 sourcefile~mod_topology.f90->sourcefile~mod_constants.f90 sourcefile~mod_topology.f90->sourcefile~mod_memory.f90 sourcefile~mod_utils.f90->sourcefile~mod_constants.f90 sourcefile~mod_utils.f90->sourcefile~mod_memory.f90

Files dependent on this one

sourcefile~~mod_electrostatics.f90~~AfferentGraph sourcefile~mod_electrostatics.f90 mod_electrostatics.f90 sourcefile~mod_geomgrad.f90 mod_geomgrad.f90 sourcefile~mod_geomgrad.f90->sourcefile~mod_electrostatics.f90 sourcefile~mod_mmpol.f90 mod_mmpol.f90 sourcefile~mod_geomgrad.f90->sourcefile~mod_mmpol.f90 sourcefile~mod_polarization.f90 mod_polarization.f90 sourcefile~mod_geomgrad.f90->sourcefile~mod_polarization.f90 sourcefile~mod_c_interface.f90 mod_c_interface.f90 sourcefile~mod_c_interface.f90->sourcefile~mod_electrostatics.f90 sourcefile~mod_c_interface.f90->sourcefile~mod_mmpol.f90 sourcefile~mod_qm_helper.f90 mod_qm_helper.f90 sourcefile~mod_c_interface.f90->sourcefile~mod_qm_helper.f90 sourcefile~mod_interface.f90 mod_interface.f90 sourcefile~mod_c_interface.f90->sourcefile~mod_interface.f90 sourcefile~mod_solvers.f90 mod_solvers.f90 sourcefile~mod_solvers.f90->sourcefile~mod_electrostatics.f90 sourcefile~rotate_multipoles.f90 rotate_multipoles.f90 sourcefile~rotate_multipoles.f90->sourcefile~mod_electrostatics.f90 sourcefile~mod_link_atom.f90 mod_link_atom.f90 sourcefile~mod_link_atom.f90->sourcefile~mod_electrostatics.f90 sourcefile~mod_prm.f90 mod_prm.f90 sourcefile~mod_link_atom.f90->sourcefile~mod_prm.f90 sourcefile~mod_mmpol.f90->sourcefile~mod_electrostatics.f90 sourcefile~mod_mmpol.f90->sourcefile~mod_link_atom.f90 sourcefile~mod_prm.f90->sourcefile~mod_electrostatics.f90 sourcefile~mod_qm_helper.f90->sourcefile~mod_electrostatics.f90 sourcefile~mod_qm_helper.f90->sourcefile~mod_link_atom.f90 sourcefile~mod_qm_helper.f90->sourcefile~mod_mmpol.f90 sourcefile~mod_qm_helper.f90->sourcefile~mod_prm.f90 sourcefile~mod_interface.f90->sourcefile~mod_electrostatics.f90 sourcefile~mod_interface.f90->sourcefile~mod_geomgrad.f90 sourcefile~mod_interface.f90->sourcefile~mod_link_atom.f90 sourcefile~mod_interface.f90->sourcefile~mod_mmpol.f90 sourcefile~mod_interface.f90->sourcefile~mod_prm.f90 sourcefile~mod_interface.f90->sourcefile~mod_qm_helper.f90 sourcefile~mod_interface.f90->sourcefile~mod_polarization.f90 sourcefile~mod_inputloader.f90 mod_inputloader.f90 sourcefile~mod_interface.f90->sourcefile~mod_inputloader.f90 sourcefile~mod_polarization.f90->sourcefile~mod_electrostatics.f90 sourcefile~mod_polarization.f90->sourcefile~mod_solvers.f90 sourcefile~mod_polarization.f90->sourcefile~mod_mmpol.f90 sourcefile~mod_inputloader.f90->sourcefile~mod_electrostatics.f90 sourcefile~mod_inputloader.f90->sourcefile~mod_mmpol.f90 sourcefile~mod_inputloader.f90->sourcefile~mod_prm.f90 sourcefile~test_si_geomgrad_num.f90 test_SI_geomgrad_num.f90 sourcefile~test_si_geomgrad_num.f90->sourcefile~mod_mmpol.f90 sourcefile~test_si_geomgrad_num.f90->sourcefile~mod_interface.f90 sourcefile~test_si_geomgrad.f90 test_SI_geomgrad.f90 sourcefile~test_si_geomgrad.f90->sourcefile~mod_interface.f90 sourcefile~test_si_geomgrad.f90->sourcefile~test_si_geomgrad_num.f90 sourcefile~test_si_init.f90 test_SI_init.f90 sourcefile~test_si_init.f90->sourcefile~mod_interface.f90 sourcefile~test_si_potential.f90 test_SI_potential.f90 sourcefile~test_si_potential.f90->sourcefile~mod_interface.f90

Contents


Source Code

#include "f_cart_components.h"
module mod_electrostatics
    use mod_io, only: fatal_error, ommp_message
    use mod_constants, only: OMMP_VERBOSE_DEBUG
    use mod_profiling, only: time_push, time_pull
    use mod_memory, only: ip, rp, lp
    use mod_adjacency_mat, only: yale_sparse
    use mod_topology, only: ommp_topology_type

    !! TODO Check the signs in electrostatic elemental functions
    !! TODO [OPT] Use Laplace equation to simplify the calculations:
    !! TODO  1. Egrd(_zz_) = -(Egrd(_xx_) + Egrd(_yy_))
    !! TODO  2. EHes(_zzz_) = -(EHes(_xxz_) + EHes(_yyz_))
    !! TODO  3. EHes(_zzx_) = -(EHes(_xxx_) + EHes(_yyx_))
    !! TODO  4. EHes(_zzy_) = -(EHes(_xxy_) + EHes(_yyy_))
    !! TODO [OPT] Fundamental electrostatic functions should be pure/elemental
    !! TODO [BUG] Handling of flags gg
    implicit none 
    private

    type ommp_electrostatics_type
        integer(ip) :: def_solver
        !! Solver to be used by default for this eel object.
        integer(ip) :: def_matv
        !! Matrix-vector method to be used by default for this eel object.

        type(ommp_topology_type), pointer :: top
        !! Data structure containing all the topological informations
        integer(ip) :: pol_atoms 
        !! number of polarizable atoms
        logical(lp) :: amoeba
        !! True if AMOEBA FF is used
        integer(ip) :: ld_cart, ld_cder
        !! size of the cartesian multipolar distribution (i.e., (l+1)*(l+2)*(l+3)/6)
        !! this is 1 for AMBER (charges only), 10 for AMOEBA (up to quadrupoles). 
        !! this is also the size of the array that contains the electrostatic properties
        !! of the sources at the sources. ld_cder is the leading size of the derivative of
        !! such a distribution, which is 3 for AMBER and 19 for AMOEBA.
        integer(ip) :: n_ipd 
        !! number of induced point dipoles distributions 
        !! this is 1 for AMBER and 2 for AMOEBA
        
        real(rp) :: mscale(4)
        !! factors for charge-charge (or multipole-multipole) interactions

        real(rp) :: pscale(4)
        !! factors for chrage-ipd (or multipole-ipd) interactions.
        !! in AMOEBA, this is used to define the polarization field, i.e., the right-hand
        !! side to the polarization equations, and depends on the connectivity.

        real(rp) :: pscale_intra(4)
        !! Only used for AMOEBA, same as pscale but for atoms that belong to the 
        !! same polarization group
        
        real(rp) :: dscale(4)
        !! factors for multipoles-ipd interactions used to compute the direct field,
        !! which is used to define the polarization energy. these factors depend on 
        !! the polarization group "connectivity" (AMOEBA only)

        real(rp) :: uscale(4)
        !! factor for ipd-ipd interactions. these depend on the connectivity (AMBER)
        !! or on the polarization group " connectivity (AMOEBA)

        ! allocatable arrays which describe the polarizable system
        
        real(rp), allocatable :: thole(:)
        !! array to store the thole factors for computing damping functions
        
        real(rp) :: thole_scale
        !! Scale factor for thole damping (only used by non-AMOEBA FF); all
        !! the element of thole(:) are multiplied by thole_scale ** 0.5

        real(rp), allocatable :: cpol(:,:)
        !! Coordinates of polarizable atoms (3:pol_atoms)
        
        real(rp), allocatable :: q(:,:)
        !! Mutlipolar distribution (ld_cart:mm_atoms)
        !! For AMOEBA this is the rotated distribution.
        !! The order for the stored multipoles is
        !! q, px, py, pz, Qxx, Qxy, Qyy, Qxz, Qyx, Qzz.

        real(rp), allocatable :: q0(:,:)
        !! Unrotated utlipolar distribution (ld_cart:mm_atoms)
        !! (AMOEBA only)
        
        real(rp), allocatable :: pol(:)
        !! Polarizabilities for each polarizable atom
        
        integer(ip), allocatable :: mm_polar(:)
        !! for each mm atom: 0 if is not polarizable, index in 
        !! polarizable atom list otherwise

        integer(ip), allocatable :: polar_mm(:)
        !! positions of a polarizable atom in the mm atoms list

        integer(ip), allocatable :: C_polar_mm(:)
        !! [[polar_mm]] with 0-based C-indexing, only allocated at need.
        
        integer(ip), allocatable :: mmat_polgrp(:)
        !! Polarizability group index for each MM site

        type(yale_sparse) :: polgrp_mmat
        !! For each polarization group index, list all the MM atoms included.
        !! It basically is a sparse boolean matrix of dimension 
        !! N_polgroups x N_mmatoms

        type(yale_sparse), allocatable :: pg_conn(:)
        !! Adjacency and connectivity matytrices between polarizability groups.
        !! Two groups are said to be adjacent if they are connected by a chemical 
        !! bond. The 1st element is the identity matrix for code simplicity.
        
        ! parameters for the definition of the rotation matrices for the multipoles:
        integer(ip), allocatable :: mol_frame(:)
        !! definition of the molecular frame
        !! convention: 0 ... do not rotate
        !!             1 ... z-then-x
        !!             2 ... bisector
        !!             3 ... z-only
        !!             4 ... z-bisector
        !!             5 ... 3-fold

        integer(ip), allocatable :: ix(:), iy(:), iz(:)
        !! neighboring atoms used to define the axes of the molecular frame
        
        !- Intermediate data allocate here -!
        logical(lp) :: M2M_done = .false.
        !! flag to set when M2M electrostatic quantities are computed.
        logical(lp) :: M2Mgg_done = .false.
        !! flag to set when M2M electrostatic quantities for geometrical
        !! gradients are computed.
        real(rp), allocatable :: V_M2M(:)
        !! potential of MM permanent multipoles at MM sites; 
        real(rp), allocatable :: E_M2M(:,:)
        !! electric_field of MM permanent multipoles at MM sites; 
        real(rp), allocatable :: Egrd_M2M(:,:)
        !! electric_field gradient of MM permanent multipoles at MM sites;
        real(rp), allocatable :: EHes_M2M(:,:) 
        !! electric field Hessian of MM permanent multipoles at MM sites;

        logical(lp) :: M2D_done = .false.
        !! Flag to set when M2D electrostatics have been computed.
        logical(lp) :: M2Dgg_done = .false.
        !! Flag to set when M2D electrostatics for geometrical gradients 
        !! have been computed.
        real(rp), allocatable :: V_M2D(:,:)
        ! electrostatic potential of MM permanent multipoles at POL sites; unused.
        real(rp), allocatable :: E_M2D(:,:,:) ! TODO the third dimension is used?
        !! electric field of MM permanent multipoles at POL sites; 
        real(rp), allocatable :: Egrd_M2D(:,:,:)
        !! electric field of MM permanent multipoles at POL sites;
        real(rp), allocatable :: EHes_M2D(:,:,:)
        ! electric field Hessian of MM permanent multipoles at POL sites; unused.

        logical(lp) :: D2Mgg_done = .false.
        real(rp), allocatable :: V_D2M(:)
        real(rp), allocatable :: E_D2M(:,:)
        real(rp), allocatable :: Egrd_D2M(:,:)
        real(rp), allocatable :: EHes_D2M(:,:)

        logical(lp) :: D2Dgg_done = .false.
        real(rp), allocatable :: V_D2D(:,:)
        real(rp), allocatable :: E_D2D(:,:,:)
        real(rp), allocatable :: Egrd_D2D(:,:,:)
        real(rp), allocatable :: EHes_D2D(:,:,:)
        
        logical(lp) :: ipd_done = .false.
        !! Flag to set when IPD have been computed.
        logical(lp) :: ipd_use_guess = .false.
        !! Flag to set when current value of IPD can be
        !! used as guess for next solution of LS.
        real(rp), allocatable :: ipd(:,:,:)
        !! induced point dipoles (3:pol_atoms:ipd) 
    
        real(rp), allocatable :: TMat(:,:)
        !! Interaction tensor, only allocated for the methods that explicitly 
        !! requires it.
        
        logical(lp) :: screening_list_done = .false.
        !! Flag to check if screening list have already been prepared
        type(yale_sparse), allocatable :: list_S_S, list_P_P, list_S_P_P, list_S_P_D
        !! Sparse matrix containg the scale factors for the scaled elements
        !! of electrostatic interactions (all the element that are not 
        !! present in the sparse matrix have a scaling factor 1.0).
        logical(lp), dimension(:), allocatable :: todo_S_S, todo_P_P, todo_S_P_P, todo_S_P_D
        !! Logical array of the same dimension of column-index vector; true if 
        !! the scaling factor is zero, false otherwise
        real(rp), dimension(:), allocatable :: scalef_S_S, scalef_P_P, &
                                               scalef_S_P_P, scalef_S_P_D
        !! Array of the same dimension of column-index vector; contains the 
        !! value of the scaling factors different from 1.0
    end type ommp_electrostatics_type

    public :: ommp_electrostatics_type
    public :: electrostatics_init, electrostatics_terminate
    public :: set_def_solver, set_def_matv
    public :: thole_init, remove_null_pol, set_screening_parameters
    public :: screening_rules, make_screening_lists
    public :: damped_coulomb_kernel, field_extD2D
    public :: energy_MM_MM, energy_MM_pol
    public :: prepare_fixedelec, prepare_polelec
    public :: q_elec_prop, coulomb_kernel
    public :: potential_M2E, potential_D2E
    public :: field_M2E, field_D2E

    contains

    subroutine electrostatics_init(eel_obj, amoeba, pol_atoms, top_obj)
        use mod_memory, only: mallocate
        use mod_constants, only: OMMP_MATV_DEFAULT, OMMP_SOLVER_DEFAULT

        implicit none 

        logical(lp), intent(in) :: amoeba
        integer(ip), intent(in) :: pol_atoms
        type(ommp_topology_type), intent(in), target :: top_obj
        type(ommp_electrostatics_type), intent(inout) :: eel_obj

        integer(ip) :: mm_atoms

        mm_atoms = top_obj%mm_atoms
        eel_obj%amoeba = amoeba
        eel_obj%pol_atoms = pol_atoms
        eel_obj%top => top_obj
        eel_obj%def_solver = OMMP_SOLVER_DEFAULT
        eel_obj%def_matv = OMMP_MATV_DEFAULT

        if(amoeba) then
            eel_obj%ld_cart = 10_ip
            eel_obj%ld_cder = 19_ip
            eel_obj%n_ipd = 2_ip
        else
            eel_obj%ld_cart = 1_ip
            eel_obj%ld_cder = 3_ip
            eel_obj%n_ipd = 1_ip
        endif

        call mallocate('electrostatics_init [q]', eel_obj%ld_cart, &
                       mm_atoms, eel_obj%q)
        call mallocate('electrostatics_init [pol]', eel_obj%pol_atoms, &
                       eel_obj%pol)
        call mallocate('electrostatics_init [cpol]', 3_ip, eel_obj%pol_atoms, &
                       eel_obj%cpol)
        call mallocate('electrostatics_init [polar_mm]', eel_obj%pol_atoms, &
                       eel_obj%polar_mm)
        call mallocate('electrostatics_init [mm_polar]', mm_atoms, &
                       eel_obj%mm_polar)
        call mallocate('electrostatics_init [thole]', mm_atoms, &
                       eel_obj%thole)
        
        call mallocate('electrostatics_init [idp]', 3_ip, eel_obj%pol_atoms, &
                       eel_obj%n_ipd, eel_obj%ipd) 
        eel_obj%ipd_done = .false.
        eel_obj%ipd = 0.0_rp

        if (eel_obj%amoeba) then
            ! Extra quantities that should be allocated only
            ! for AMOEBA
            call mallocate('electrostatics_init [q0]', eel_obj%ld_cart, &
                           mm_atoms, eel_obj%q0)
            
            call mallocate('electrostatics_init [mmat_polgrp]', & 
                           mm_atoms, eel_obj%mmat_polgrp)

            call mallocate('electrostatics_init [mol_frame]', &
                           mm_atoms, eel_obj%mol_frame)
            call mallocate('electrostatics_init [ix]', mm_atoms, eel_obj%ix)
            call mallocate('electrostatics_init [iy]', mm_atoms, eel_obj%iy)
            call mallocate('electrostatics_init [iz]', mm_atoms, eel_obj%iz)
        end if

    end subroutine electrostatics_init

    subroutine electrostatics_terminate(eel_obj)
        use mod_memory, only: mfree
        use mod_adjacency_mat, only: matfree

        implicit none

        type(ommp_electrostatics_type), intent(inout) :: eel_obj
        integer(ip) :: i

        call mfree('electrostatics_terminate [q]', eel_obj%q)
        call mfree('electrostatics_terminate [pol]', eel_obj%pol)
        call mfree('electrostatics_terminate [cpol]', eel_obj%cpol)
        call mfree('electrostatics_terminate [polar_mm]', eel_obj%polar_mm)
        call mfree('electrostatics_terminate [mm_polar]', eel_obj%mm_polar)
        call mfree('electrostatics_terminate [thole]', eel_obj%thole)
        call mfree('electrostatics_terminate [idp]', eel_obj%ipd) 

        if (eel_obj%amoeba) then
            call mfree('electrostatics_terminate [q0]', eel_obj%q0)
            call mfree('electrostatics_terminate [mmat_polgrp]', &
                       eel_obj%mmat_polgrp)
            if(allocated(eel_obj%pg_conn)) then
                do i=1, size(eel_obj%pg_conn)
                    call matfree(eel_obj%pg_conn(i))
                end do
                deallocate(eel_obj%pg_conn)
            end if
            
            call matfree(eel_obj%polgrp_mmat)
            call mfree('electrostatics_terminate [mol_frame]', &
                       eel_obj%mol_frame)
            call mfree('electrostatics_terminate [ix]', eel_obj%ix)
            call mfree('electrostatics_terminate [iy]', eel_obj%iy)
            call mfree('electrostatics_terminate [iz]', eel_obj%iz)
        end if
        
        call mfree('electrostatics_terminate [E_M2D]', eel_obj%E_M2D)
        call mfree('electrostatics_terminate [V_M2M]', eel_obj%V_M2M)
        call mfree('electrostatics_terminate [E_M2M]', eel_obj%E_M2M)
        call mfree('electrostatics_terminate [Egrd_M2M]', eel_obj%Egrd_M2M)

        if(allocated(eel_obj%todo_S_S)) deallocate(eel_obj%todo_S_S)
        if(allocated(eel_obj%todo_P_P)) deallocate(eel_obj%todo_P_P)
        call mfree('electrostatics_terminate [scalef_S_S]', eel_obj%scalef_S_S)
        call mfree('electrostatics_terminate [scalef_P_P]', eel_obj%scalef_P_P)
        if(allocated(eel_obj%list_S_S)) then
            call matfree(eel_obj%list_S_S)
            deallocate(eel_obj%list_S_S)
        end if
        if(allocated(eel_obj%list_P_P)) then
            call matfree(eel_obj%list_P_P)
            deallocate(eel_obj%list_P_P)
        end if

    end subroutine electrostatics_terminate

    subroutine set_def_solver(eel_obj, solver)
        use mod_constants, only: OMMP_SOLVER_CG, OMMP_SOLVER_INVERSION, OMMP_SOLVER_DIIS
        implicit none
        
        type(ommp_electrostatics_type), intent(inout) :: eel_obj
        integer(ip), intent(in) :: solver

        if(solver /= OMMP_SOLVER_CG .and. &
           solver /= OMMP_SOLVER_INVERSION .and. &
           solver /= OMMP_SOLVER_DIIS) &
            call fatal_error("Unrecognized setting for default solver method")
        eel_obj%def_solver = solver
    end subroutine

    subroutine set_def_matv(eel_obj, matv)
        use mod_constants, only: OMMP_MATV_INCORE, OMMP_MATV_DIRECT
        implicit none
        
        type(ommp_electrostatics_type), intent(inout) :: eel_obj
        integer(ip), intent(in) :: matv

        if(matv /= OMMP_MATV_INCORE .and. &
           matv /= OMMP_MATV_DIRECT) &
            call fatal_error("Unrecognized setting for default matrix-vector method")
        eel_obj%def_matv = matv
    end subroutine
    
    subroutine set_screening_parameters(eel_obj, m, p, d, u, i)
        !! Subroutine to initialize the screening parameters
       
        implicit none

        type(ommp_electrostatics_type), intent(inout) :: eel_obj
        real(rp), intent(in) :: m(4), p(4), d(4), u(4)
        real(rp), optional, intent(in) :: i(4)
        
        eel_obj%mscale = m
        eel_obj%pscale = p
        eel_obj%dscale = d
        eel_obj%uscale = u
        
        if(present(i)) then
            if(eel_obj%amoeba) then
                eel_obj%pscale_intra = i
            else
                call fatal_error("Scale factors for atoms of the same group &
                                 &cannot be set outside AMOEBA FF")
            end if
        else
            if(eel_obj%amoeba) &
                call fatal_error("Scale factors for atoms of the same group &
                                 &should be defined in AMOEBA FF")
        end if
        
    end subroutine set_screening_parameters

    subroutine remove_null_pol(eel)
        !! Check which polarizabilities are close enough to 0 to be 
        !! just excluded from the calculation, and remove them.

        use mod_constants, only: OMMP_STR_CHAR_MAX, OMMP_VERBOSE_LOW, eps_rp
        use mod_memory, only: mallocate, mfree
        implicit none

        type(ommp_electrostatics_type), intent(inout) :: eel
        integer(ip), allocatable :: idx(:), polar_mm(:)
        integer(ip) :: i, nidx
        real(rp), allocatable :: tmp(:)
        character(len=OMMP_STR_CHAR_MAX) :: msg

        if(eel%pol_atoms > 0 .and. allocated(eel%pol)) then
            call mallocate('remove_null_pol [idx]', eel%pol_atoms, idx)
            nidx = 0
            do i=1, eel%pol_atoms
                if(abs(eel%pol(i)) > eps_rp) then
                    nidx = nidx + 1
                    idx(nidx) = i
                end if
            end do

            if(nidx < eel%pol_atoms) then
                write(msg, '(A,I0,A,I0)') "Removing ", eel%pol_atoms - nidx, &
                                          " polarizabilities out of ", eel%pol_atoms
                call ommp_message(msg, OMMP_VERBOSE_LOW)
                call mallocate('remove_null_pol [tmp]', eel%pol_atoms, tmp)
                ! Polarizabilities
                tmp = eel%pol
                call mfree('remove_null_pol [eel%pol]', eel%pol)
                call mallocate('remove_null_pol [eel%pol]', nidx, eel%pol)
                do i=1, nidx
                    eel%pol(i) = tmp(idx(i))
                end do
                call mfree('remove_null_pol [tmp]', tmp)
                    
                call mallocate('remove_null_pol [polar_mm]', eel%pol_atoms, polar_mm)
                polar_mm = eel%polar_mm
                call mfree('remove_null_pol [eel%polar_mm]', eel%polar_mm)
                call mallocate('remove_null_pol [eel%polar_mm]', nidx, eel%polar_mm)
                call mfree('remove_null_pol [eel%cpol]', eel%cpol)
                call mallocate('remove_null_pol [eel%cpol]', 3, nidx, eel%cpol)
                eel%mm_polar = 0
                do i=1, nidx
                    eel%polar_mm(i) = polar_mm(idx(i))
                    eel%mm_polar(eel%polar_mm(i)) = i
                    eel%cpol(:,i) = eel%top%cmm(:,eel%polar_mm(i))
                end do

                call mfree('remove_null_pol [polar_mm]', polar_mm)
                
                eel%pol_atoms = nidx
                call mfree('remove_null_pol [eel%ipd]', eel%ipd)
                call mallocate('electrostatics_init [idp]', 3_ip, eel%pol_atoms, &
                               eel%n_ipd, eel%ipd) 
                eel%ipd_done = .false.
                eel%ipd_use_guess = .false.
                eel%ipd = 0.0_rp
                
                call ommp_message("Removing null polarizable sites done", OMMP_VERBOSE_LOW)
            end if
            
            call mfree('remove_null_pol [idx]', idx)
        end if

    end subroutine

    subroutine make_screening_lists(eel)
        use mod_memory, only: mallocate, mfree
        use mod_adjacency_mat, only: compress_list, compress_data
        use mod_constants, only: eps_rp
        
        implicit none

        type(ommp_electrostatics_type), intent(inout) :: eel
        
        integer(ip) :: i, ineigh, ij, j, ns, ns_guess, n, npol, &
                       ipp, jp, ns_guess_grp, pg_i, igrp, grp, ib, ie
        logical :: to_do, to_scale
        real(rp) :: scalf
        integer(ip), allocatable :: itmp(:,:), ntmp(:)
        real(rp), allocatable :: rtmp(:,:)

        if(eel%screening_list_done) return

        n = eel%top%mm_atoms
        npol = eel%pol_atoms

        ! First estimate the maximum number of elements to scale
        ns_guess = 0
        !$omp parallel do default(shared) private(ns, ineigh, i)
        do i=1, n
            ns = 0
            do ineigh=1, 4
                ns = ns + eel%top%conn(ineigh)%ri(i+1) - eel%top%conn(ineigh)%ri(i)
            end do
            !$omp atomic
            ns_guess = max(ns, ns_guess)
        end do
        
        if(eel%amoeba) then
            ! If amoeba is used, some scaling (D-field) are performed using 
            ! polarization group, so a second guess should be done and the 
            ! maximum is used in memory allocation.
            ns_guess_grp = 0
            
            !$omp parallel do default(shared) & 
            !$omp private(i, pg_i, ineigh, igrp, grp, ns)
            do i=1, n
                pg_i = eel%mmat_polgrp(i)
                ns = 0
                do ineigh=1, 4
                    do igrp=eel%pg_conn(ineigh)%ri(pg_i), &
                            eel%pg_conn(ineigh)%ri(pg_i+1)-1

                        grp = eel%pg_conn(ineigh)%ci(igrp)
                        ns = ns + eel%polgrp_mmat%ri(grp+1)-eel%polgrp_mmat%ri(grp)
                    end do
                end do
                !$omp atomic
                ns_guess_grp = max(ns, ns_guess_grp)
            end do
            
            ns_guess = max(ns_guess_grp, ns_guess)
        end if

        allocate(eel%list_S_S)
        allocate(eel%list_P_P)
        allocate(eel%list_S_P_P)
        if(eel%amoeba) allocate(eel%list_S_P_D)

        !allocate(ltmp(ns_guess))
        call mallocate('make_screening_list [rtmp]', ns_guess, n, rtmp)
        call mallocate('make_screening_list [itmp]', ns_guess, n, itmp)
        call mallocate('make_screening_list [ntmp]', n, ntmp)


        ! Build S S list
        !$omp parallel do schedule(dynamic) default(shared) &
        !$omp private(i,ineigh,ij,j, scalf)
        do i=1, n
            ntmp(i) = 0
            do ineigh=1, 4
                do ij=eel%top%conn(ineigh)%ri(i), eel%top%conn(ineigh)%ri(i+1)-1
                    j = eel%top%conn(ineigh)%ci(ij)
                    scalf = screening_rules(eel, i, 'S', j, 'S', '-')
                    if(abs(scalf-1.0) > eps_rp) then
                        ntmp(i) = ntmp(i) + 1
                        itmp(ntmp(i),i) = j
                        rtmp(ntmp(i),i) = scalf
                    end if
                end do
            end do
        end do

        ! Compress the list
        call compress_list(n, itmp, ntmp, eel%list_S_S)
        call compress_data(eel%list_S_S, rtmp, eel%scalef_S_S)
        call mallocate('make_screening_list [todo_S_S]', &
                       size(eel%scalef_S_S), eel%todo_S_S)
        !$omp parallel do
        do i=1, size(eel%scalef_S_S)
            eel%todo_S_S(i) = (abs(eel%scalef_S_S(i)) > eps_rp)
        end do

        ! Build P P list
        ntmp = 0
        !$omp parallel do schedule(dynamic) default(shared) & 
        !$omp private(ipp,i,ineigh,ij,j,jp,scalf)
        do ipp=1, npol
            i = eel%polar_mm(ipp)
            do ineigh=1, 4
                do ij=eel%top%conn(ineigh)%ri(i), eel%top%conn(ineigh)%ri(i+1)-1
                    j = eel%top%conn(ineigh)%ci(ij)
                    jp = eel%mm_polar(j)
                    if(jp > 0) then
                        scalf = screening_rules(eel, ipp, 'P', jp, 'P', '-')
                        if(abs(scalf-1.0) > eps_rp) then
                            ntmp(ipp) = ntmp(ipp) + 1
                            itmp(ntmp(ipp),ipp) = jp
                            rtmp(ntmp(ipp),ipp) = scalf
                        end if
                    end if
                end do
            end do
        end do

        ! Compress the list
        call compress_list(npol, itmp, ntmp(1:npol), eel%list_P_P)
        call compress_data(eel%list_P_P, rtmp, eel%scalef_P_P)
        call mallocate('make_screening_list [todo_P_P]', &
                       size(eel%scalef_P_P), eel%todo_P_P)
        !$omp parallel do
        do i=1, size(eel%scalef_P_P)
            eel%todo_P_P(i) = (abs(eel%scalef_P_P(i)) > eps_rp)
        end do
        
        ! Build S P lists
        !$omp parallel do schedule(dynamic) default(shared) &
        !$omp private(i,ineigh,ij,j,jp,scalf)
        do i=1, n
            ntmp(i) = 0
            do ineigh=1, 4
                do ij=eel%top%conn(ineigh)%ri(i), eel%top%conn(ineigh)%ri(i+1)-1
                    j = eel%top%conn(ineigh)%ci(ij)
                    jp = eel%mm_polar(j)
                    if(jp > 0) then
                        ! Needed for either amoeba or amber
                        scalf = screening_rules(eel, i, 'S', jp, 'P', 'P')
                        if(abs(scalf-1.0) > eps_rp) then
                            ntmp(i) = ntmp(i) + 1
                            itmp(ntmp(i),i) = jp
                            rtmp(ntmp(i),i) = scalf
                        end if
                    end if
                end do
            end do
        end do

        ! Compress the list
        call compress_list(n, itmp, ntmp, eel%list_S_P_P)
        call compress_data(eel%list_S_P_P, rtmp, eel%scalef_S_P_P)
        call mallocate('make_screening_list [todo_S_P_P]', &
                       size(eel%scalef_S_P_P), eel%todo_S_P_P)
        !$omp parallel do
        do i=1, size(eel%scalef_S_P_P)
            eel%todo_S_P_P(i) = (abs(eel%scalef_S_P_P(i)) > eps_rp)
        end do

        if(eel%amoeba) then
            !$omp parallel do schedule(dynamic) default(shared) &
            !$omp private(i,ineigh,ij,j,jp,scalf,pg_i,igrp,grp)
            do i=1, n
                ntmp(i) = 0
                pg_i = eel%mmat_polgrp(i)
                
                do ineigh=1, 4
                    do igrp=eel%pg_conn(ineigh)%ri(pg_i), &
                            eel%pg_conn(ineigh)%ri(pg_i+1)-1
                        
                        grp = eel%pg_conn(ineigh)%ci(igrp)

                        do j=eel%polgrp_mmat%ri(grp), &
                             eel%polgrp_mmat%ri(grp+1)-1
                            jp = eel%mm_polar(j)
                            if(jp > 0) then
                                scalf = screening_rules(eel, i, 'S', jp, 'P', 'D')
                                if(abs(scalf-1.0) > eps_rp) then
                                    ntmp(i) = ntmp(i) + 1
                                    itmp(ntmp(i),i) = jp
                                    rtmp(ntmp(i),i) = scalf
                                end if
                            end if
                        end do
                    end do
                end do
            end do

            ! Compress the list
            call compress_list(n, itmp, ntmp, eel%list_S_P_D)
            call compress_data(eel%list_S_P_D, rtmp, eel%scalef_S_P_D)
            call mallocate('make_screening_list [todo_S_P_D]', &
                           size(eel%scalef_S_P_D), eel%todo_S_P_D)
            !$omp parallel do
            do i=1, size(eel%scalef_S_P_D)
                eel%todo_S_P_D(i) = (abs(eel%scalef_S_P_D(i)) > eps_rp)
            end do
        end if
        
        eel%screening_list_done = .true.

        call mfree('make_screening_list [rtmp]', rtmp)
        call mfree('make_screening_list [itmp]', itmp)
        call mfree('make_screening_list [ntmp]', ntmp)
    end subroutine

    subroutine thole_init(eel)
        ! This routine compute the thole factors and stores
        ! them in a vector. TODO add reference
        ! TODO in AMOEBA should be read from FF
       
        use mod_constants, only: OMMP_VERBOSE_LOW, eps_rp
        implicit none

        type(ommp_electrostatics_type), intent(inout) :: eel
        
        integer(ip) :: i, j
        
        eel%thole = 0.0_rp
        
        do i = 1, eel%pol_atoms
            j = eel%polar_mm(i)
            eel%thole(j) = eel%pol(i) ** (1.0_rp/6.0_rp)
        end do
        
        if(.not. eel%amoeba) then
            if(eel%thole_scale > eps_rp) then
                eel%thole = eel%thole * sqrt(eel%thole_scale)
            else
                call fatal_error("Scale factor for Thole damping should be &
                                 &greater than 0.0 when non-AMOEBA FF are used")
            end if
        else
            if(eel%thole_scale > eps_rp) then
                call ommp_message("Scale factor set for thole damping is &
                    &ignored because AMOEBA FF is used", OMMP_VERBOSE_LOW)
            end if
        end if
    end subroutine thole_init

    subroutine energy_MM_MM(eel, ene)
        !! This function computes the interaction energy of 
        !! static electric multipoles
        use mod_memory, only: mallocate, mfree

        implicit none

        type(ommp_electrostatics_type), intent(inout) :: eel
        !! Electrostatics data structure
        real(rp), intent(inout) :: ene
        !! Energy (results will be added)
        real(rp) :: eMM

        integer(ip) :: i

        call prepare_fixedelec(eel)
        eMM = 0.0

        if(eel%amoeba) then
            eMM = eMM + dot_product(eel%q(1,:), eel%V_M2M)
            do i=1, 3
                eMM = eMM - dot_product(eel%q(i+1,:), eel%E_M2M(i,:))
            end do

            do i=1,6
                if(i == _xx_ .or. i == _yy_ .or. i == _zz_) then
                    ! diagonal elements
                    eMM = eMM + dot_product(eel%q(i+4,:), eel%Egrd_M2M(i,:))
                else
                    ! off-diagonal elements (are stored once, but 
                    ! should be counted twice
                    eMM = eMM + 2.0 * dot_product(eel%q(i+4,:), eel%Egrd_M2M(i,:))
                end if
            end do
        else
            eMM = eMM + dot_product(eel%q(1,:), eel%V_M2M)
        end if

        ! Since potential is computed using all the sites each 
        ! interaction is counted twice
        ene = ene + 0.5_rp * eMM
    
    end subroutine energy_MM_MM
    
    subroutine energy_MM_pol(eel, ene)
        !! This function computes the interaction energy of 
        !! static electric multipoles
        use mod_memory, only: mallocate, mfree

        implicit none

        type(ommp_electrostatics_type), intent(inout) :: eel
        !! Electrostatics data structure
        real(rp), intent(inout) :: ene
        !! Energy (results will be added)
        real(rp) :: eMM

        integer(ip) :: i

        call prepare_polelec(eel)
        if(.not. eel%ipd_done) call fatal_error("IPD should be computed before&
                                                & computing MM-Pol energy.")
        eMM = 0.0
        
        if(eel%amoeba) then
            do i=1, 3
                eMM = eMM - dot_product(eel%ipd(i,:,_amoeba_D_), &
                                        eel%E_M2D(i,:,_amoeba_P_))
            end do
        else
            do i=1, 3
                eMM = eMM - dot_product(eel%ipd(i,:,1), eel%E_M2D(i,:,1))
            end do
        end if

        ene = ene + 0.5_rp * eMM

    end subroutine energy_MM_pol

    subroutine coulomb_kernel(dr, maxder, res)
        !! This function compute the coulomb kernel for the distance vector dr and 
        !! its derivative up to the value required by maxder.
        use mod_memory, only: ip, rp
        use mod_io, only: fatal_error
        use mod_constants, only: eps_rp

        implicit none

        real(rp), intent(in) :: dr(3)
        !! Distance vector from atom i to atom j
        integer(ip), intent(in) :: maxder
        !! Maximum derivative to be computed
        real(rp), intent(out), dimension(maxder+1) :: res
        !! Results vector
        
        integer(ip) :: ii
        real(rp) :: norm2_r, inv_norm_sq 

        if(maxder < 0) then
            ! Strange request, maybe a warning should be printed
            return 
        end if
        
        norm2_r = sqrt(dr(1)*dr(1) + dr(2)*dr(2) + dr(3)*dr(3))
        if(norm2_r < eps_rp) then
            call fatal_error("Requesting Coulomb kernel 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
        res(1) = 1.0_rp / norm2_r
        inv_norm_sq = res(1) * res(1)

        do ii = 1, maxder
            res(ii+1) = res(ii) * inv_norm_sq
        end do
        
    end subroutine coulomb_kernel

    subroutine damped_coulomb_kernel(eel, i, j, maxder, res, dr)
        !! This subroutine computes the damped coulomb kernel between two atoms.
        !! Note that this only makes sense between two MM atoms, as it is only used
        !! to compute the field that induces the point dipoles!

        use mod_memory, only: ip, rp
        use mod_constants, only: eps_rp
        
        implicit none

        type(ommp_electrostatics_type), intent(in) :: eel
        !! Electrostatics data structure
        integer(ip), intent(in) :: i, j
        !! Index of atoms (in MM atom list) for which the kernel should be
        !! computed
        integer(ip), intent(in) :: maxder
        !! Maximum derivative to be computed
        real(rp), intent(out), dimension(maxder+1) :: res
        !! Results vector
        real(rp), intent(out), dimension(3) :: dr
        !! Distance vector between i and j
        
        real(rp) :: s, u, u3, u4, fexp, eexp
        
        s = eel%thole(i) * eel%thole(j)
        
        ! Compute undamped kernels
        dr = eel%top%cmm(:,j) - eel%top%cmm(:,i)
        call coulomb_kernel(dr, maxder, res)

        if(abs(s) < eps_rp) then
            ! either thole(i) or thole(j) are zero, so the damped kernel
            ! is just the same as the regular one. Job is done.
            return
        end if

        u = 1.0_rp / (res(1) * s)
        u3 = u**3

        if(eel%amoeba .and. u3 * 0.39_rp < 50.0_rp) then
            ! 1.0 / (res(1) * s)**3 * 0.39_rp < 50_rp TODO Remove 0.39, this is a
            ! FF constants, ask someone why this condition is here.
            ! Here basically we multiply the standard interactions kernel for
            ! dumping coefficients. The equations implemented here correspond to 
            ! eq. (5) of 10.1021/jp027815
            fexp = -0.39_rp * u3
            eexp = exp(fexp)
            if(maxder >= 1) res(2) = res(2) * (1.0_rp - eexp)
            if(maxder >= 2) res(3) = res(3) * (1.0_rp - (1.0_rp - fexp) * eexp)
            if(maxder >= 3) res(4) = res(4) * &
                            (1.0_rp - (1.0_rp - fexp + 0.6_rp * fexp * fexp) * eexp)
            if(maxder >= 4) res(5) = res(5) * &
                            (1.0_rp - (1.0_rp - fexp + 18.0_rp/35.0_rp * fexp *  fexp - &
                            9.0_rp/35.0_rp * fexp * fexp * fexp) * eexp)
            if(maxder >= 5) then
                call fatal_error("Damped Coulomb kernel (AMOEBA) only supports up to the 5th derivative")
            end if
        else if(.not. eel%amoeba .and. res(1) > 1_rp/s) then
            ! TODO Again it is not clear to me why condition res(1) > 1_rp/s is here.
            u4 = u3*u
            if(maxder >= 1) res(2) = res(2) * (4.0_rp * u3 - 3.0_rp * u4)
            if(maxder >= 2) res(3) = res(3) * u4
            if(maxder >= 3) res(4) = res(4) * 0.2_rp * u4
            if(maxder >= 4) then
                call fatal_error("Damped Coulomb kernel (WANG) only supports up to the 4th derivative")
            end if
        end if
    end subroutine damped_coulomb_kernel

    subroutine q_elec_prop(q, dr, kernel, &
                           do_V, V, do_E, E, do_grdE, grdE, do_HE, HE)
        !! TODO
        !! Computes the electric potential of a charge \(q\) at position
        !! \(\mathbf{dr}\) from the charge itself. Pre-computed kernel should
        !! be provided as input. 
        !! The result is added to \(V\).
        !! $$ V_q = q \frac{f(r)}{||\mathbf{dr}||} $$

        implicit none

        real(rp), intent(in) :: q
        !! Charge
        real(rp), intent(in) :: dr(3)
        !! Distance vector
        real(rp), intent(in) :: kernel(:)
        !! Array of coulomb kernel (either damped or undamped)
        logical, intent(in) :: do_V, do_E, do_grdE, do_HE
        !! Flags to enable/disable calculation of different electrostatic 
        !! properties
        real(rp), intent(inout) :: V, E(3), grdE(6), HE(10)
        !! Electric potential
        
        if(do_V) then
            V =  V + kernel(1) * q
        end if
        
        if(do_E) then
            E = E + q * kernel(2) * dr
        end if

        if(do_grdE) then
            ! xx
            grdE(1) = grdE(1) + q * 3.0_rp * kernel(3) * dr(1) * dr(1) - q * kernel(2)
            ! xy
            grdE(2) = grdE(2) + q * 3.0_rp * kernel(3) * dr(1) * dr(2)
            ! yy
            grdE(3) = grdE(3) + q * 3.0_rp * kernel(3) * dr(2) * dr(2) - q * kernel(2)
            ! xz
            grdE(4) = grdE(4) + q * 3.0_rp * kernel(3) * dr(1) * dr(3)
            ! yz
            grdE(5) = grdE(5) + q * 3.0_rp * kernel(3) * dr(2) * dr(3)
            ! zz
            grdE(6) = grdE(6) + q * 3.0_rp * kernel(3) * dr(3) * dr(3) - q * kernel(2)
        end if

        if(do_HE) then
            ! xxx
            HE(1) = HE(1) + (15.0_rp * kernel(4) * dr(1) * dr(1) * dr(1) - &
                             9.0_rp * kernel(3) * dr(1)) * q
            ! xxy
            HE(2) = HE(2) + (15.0_rp * kernel(4) * dr(1) * dr(1) * dr(2) - &
                             3.0_rp * kernel(3) * dr(2)) * q
            ! xxz
            HE(3) = HE(3) + (15.0_rp * kernel(4) * dr(1) * dr(1) * dr(3) - &
                             3.0_rp * kernel(3) * dr(3)) * q
            ! xyy
            HE(4) = HE(4) + (15.0_rp * kernel(4) * dr(2) * dr(2) * dr(1) - & 
                             3.0_rp * kernel(3) * dr(1)) * q 
            ! xyz
            HE(5) = HE(5) + 15.0_rp * kernel(4) * dr(1) * dr(2) * dr(3) * q
            
            ! xzz
            HE(6) = HE(6) + (15.0_rp * kernel(4) * dr(3) * dr(3) * dr(1) - &
                             3.0_rp * kernel(3) * dr(1)) * q
            ! yyy
            HE(7) = HE(7) + (15.0_rp * kernel(4) * dr(2) * dr(2) * dr(2) - & 
                             9.0_rp * kernel(3) * dr(2)) * q 
            ! yyz
            HE(8) = HE(8) + (15.0_rp * kernel(4) * dr(2) * dr(2) * dr(3) - &
                             3.0_rp * kernel(3) * dr(3)) * q
            ! yzz
            HE(9) = HE(9) + (15.0_rp * kernel(4) * dr(3) * dr(3) * dr(2) - &
                             3.0_rp * kernel(3) * dr(2)) * q
            ! zzz
            HE(10) = HE(10) + (15.0_rp * kernel(4) * dr(3) * dr(3) * dr(3) - &
                               9.0_rp * kernel(3) * dr(3)) * q
        end if

    end subroutine q_elec_prop
    
    subroutine mu_elec_prop(mu, dr, kernel, &
                            do_V, V, do_E, E, do_grdE, grdE, do_HE, HE)
        
        implicit none

        real(rp), intent(in) :: mu(3)
        !! point dipole
        real(rp), intent(in) :: dr(3)
        !! Distance vector
        real(rp), intent(in) :: kernel(:)
        !! Array of coulomb kernel (either damped or undamped)
        logical, intent(in) :: do_V, do_E, do_grdE, do_HE
        !! Flags to enable/disable calculation of different electrostatic 
        !! properties
        real(rp), intent(inout) :: V, E(3), grdE(6), HE(10)
        !! Electric potential
        
        real(rp) :: mu_dot_dr
        
        mu_dot_dr = dot_product(mu, dr)

        if(do_V) then
            V = V + mu_dot_dr * kernel(2)
        end if
        
        if(do_E) then
            E = E + 3.0_rp * mu_dot_dr * dr * kernel(3) - mu * kernel(2)
        end if

        if(do_grdE) then
            ! xx
            grdE(1) = grdE(1) + mu_dot_dr * kernel(4) * 15.0 * dr(1) * dr(1) - &
                              (mu_dot_dr + 2.0 * mu(1) * dr(1)) * 3.0 * kernel(3)
            ! xy
            grdE(2) = grdE(2) + mu_dot_dr * kernel(4) * 15.0 * dr(1) * dr(2) - &
                              (mu(1)*dr(2) + mu(2)*dr(1)) * 3.0 * kernel(3)
            ! yy
            grdE(3) = grdE(3) + mu_dot_dr * kernel(4) * 15.0 * dr(2) * dr(2) - &
                              (mu_dot_dr + 2.0 * mu(2) * dr(2)) * 3.0 * kernel(3)
            ! xz
            grdE(4) = grdE(4) + mu_dot_dr * kernel(4) * 15.0 * dr(1) * dr(3) - &
                              (mu(1)*dr(3) + mu(3)*dr(1)) * 3.0 * kernel(3)
            ! yz
            grdE(5) = grdE(5) + mu_dot_dr * kernel(4) * 15.0 * dr(2) * dr(3) - &
                              (mu(2)*dr(3) + mu(3)*dr(2)) * 3.0 * kernel(3)
            ! zz
            grdE(6) = grdE(6) + mu_dot_dr * kernel(4) * 15.0 * dr(3) * dr(3) - &
                              (mu_dot_dr + 2.0 * mu(3) * dr(3)) * 3.0 * kernel(3)
        end if

        if(do_HE) then
            ! xxx
            HE(1) = HE(1) + 105.0_rp * mu_dot_dr * kernel(5) * dr(1)*dr(1)*dr(1) &
                          - 45.0_rp * kernel(4) * dr(1) * (mu(1)*dr(1) + mu_dot_dr) &
                          + 9.0_rp * kernel(3) * mu(1)
            ! xxy
            HE(2) = HE(2) + 105.0_rp * kernel(5) * mu_dot_dr * dr(1)*dr(1)*dr(2) &
                          - 15.0_rp * kernel(4) * (mu(2)*dr(1)*dr(1) + &
                                                   2.0_rp *mu(1)*dr(1)*dr(2) + &
                                                   mu_dot_dr*dr(2)) &
                          + 3.0_rp * kernel(3) * mu(2)

            ! xxz
            HE(3) = HE(3) + 105.0_rp * kernel(5) * mu_dot_dr * dr(1)*dr(1)*dr(3) &
                          - 15.0_rp * kernel(4) * (mu(3)*dr(1)*dr(1) + &
                                                   2.0_rp *mu(1)*dr(1)*dr(3) + &
                                                   mu_dot_dr*dr(3)) &
                          + 3.0_rp * kernel(3) * mu(3)
            ! xyy
            HE(4) = HE(4) + 105.0_rp * kernel(5) * mu_dot_dr * dr(2)*dr(2)*dr(1) &
                          - 15.0_rp * kernel(4) * (mu(1)*dr(2)*dr(2) + &
                                                   2.0_rp *mu(2)*dr(2)*dr(1) + &
                                                   mu_dot_dr*dr(1)) &
                          + 3.0_rp * kernel(3) * mu(1)
            ! xyz
            HE(5) = HE(5) + 105.0_rp * mu_dot_dr * kernel(5) * dr(1)*dr(2)*dr(3) &
                          - 15.0_rp * kernel(4) * (mu(1)*dr(2)*dr(3) + &
                                                   dr(1)*mu(2)*dr(3) + &
                                                   dr(1)*dr(2)*mu(3))
            ! xzz
            HE(6) = HE(6) + 105.0_rp * kernel(5) * mu_dot_dr * dr(3)*dr(3)*dr(1) &
                          - 15.0_rp * kernel(4) * (mu(1)*dr(3)*dr(3) + &
                                                   2.0_rp *mu(3)*dr(3)*dr(1) + &
                                                   mu_dot_dr*dr(1)) &
                          + 3.0_rp * kernel(3) * mu(1)
            ! yyy
            HE(7) = HE(7) + 105.0_rp * mu_dot_dr * kernel(5) * dr(2)*dr(2)*dr(2) &
                          - 45.0_rp * kernel(4) * dr(2) * (mu(2)*dr(2) + mu_dot_dr) &
                          + 9.0_rp * kernel(3) * mu(2)
            ! yyz
            HE(8) = HE(8) + 105.0_rp * kernel(5) * mu_dot_dr * dr(2)*dr(2)*dr(3) &
                          - 15.0_rp * kernel(4) * (mu(3)*dr(2)*dr(2) + &
                                                   2.0_rp *mu(2)*dr(2)*dr(3) + &
                                                   mu_dot_dr*dr(3)) &
                          + 3.0_rp * kernel(3) * mu(3)
            ! yzz
            HE(9) = HE(9) + 105.0_rp * kernel(5) * mu_dot_dr * dr(3)*dr(3)*dr(2) &
                          - 15.0_rp * kernel(4) * (mu(2)*dr(3)*dr(3) + &
                                                   2.0_rp *mu(3)*dr(3)*dr(2) + &
                                                   mu_dot_dr*dr(2)) &
                          + 3.0_rp * kernel(3) * mu(2)
            ! zzz
            HE(10) = HE(10) + 105.0_rp * mu_dot_dr * kernel(5) * dr(3)*dr(3)*dr(3) &
                            - 45.0_rp * kernel(4) * dr(3) * (mu(3)*dr(3) + mu_dot_dr) &
                            + 9.0_rp * kernel(3) * mu(3)
        end if
    end subroutine mu_elec_prop
    
    subroutine quad_elec_prop(quad, dr, kernel, &
                              do_V, V, do_E, E, do_grdE, grdE, do_HE, HE)
        
        implicit none

        real(rp), intent(in) :: quad(6)
        !! point quadrupole stored as (xx, xy, yy, xz, yz, zz)
        real(rp), intent(in) :: dr(3)
        !! Distance vector
        real(rp), intent(in) :: kernel(:)
        !! Array of coulomb kernel (either damped or undamped)
        logical, intent(in) :: do_V, do_E, do_grdE, do_HE
        !! Flags to enable/disable calculation of different electrostatic 
        !! properties
        real(rp), intent(inout) :: V, E(3), grdE(6), HE(10)
        !! Electric potential
        
        real(rp) :: quadxr(3), quadxr_dot_r

        quadxr(1) = quad(1)*dr(1) + quad(2)*dr(2) + quad(4)*dr(3)
        quadxr(2) = quad(2)*dr(1) + quad(3)*dr(2) + quad(5)*dr(3)
        quadxr(3) = quad(4)*dr(1) + quad(5)*dr(2) + quad(6)*dr(3)

        quadxr_dot_r = dot_product(quadxr, dr)

        if(do_V) then
            V = V + 3.0_rp * quadxr_dot_r * kernel(3)
        end if
        
        if(do_E) then
            E = E + 15.0_rp * quadxr_dot_r * dr * kernel(4) &
                - 6.0_rp * quadxr * kernel(3)
        end if

        if(do_grdE) then
            ! xx
            grdE(1) = grdE(1) + 105.0 * kernel(5) * quadxr_dot_r * dr(1)*dr(1) + &
                              6.0 * kernel(3) * quad(1) - &
                              15.0*kernel(4)*(quadxr_dot_r+4.0*quadxr(1)*dr(1))
            ! xy
            grdE(2) = grdE(2) + 105.0 * kernel(5) * quadxr_dot_r * dr(1)*dr(2) + &
                              6.0 * kernel(3) * quad(2) - &
                              30.0*kernel(4)*(quadxr(1)*dr(2)+quadxr(2)*dr(1))
            ! yy
            grdE(3) = grdE(3) + 105.0 * kernel(5) * quadxr_dot_r * dr(2)*dr(2) + &
                              6.0 * kernel(3) * quad(3) - &
                              15.0*kernel(4)*(quadxr_dot_r+4.0*quadxr(2)*dr(2))
            ! xz
            grdE(4) = grdE(4) + 105.0 * kernel(5) * quadxr_dot_r * dr(1)*dr(3) + &
                              6.0 * kernel(3) * quad(4) - &
                              30.0*kernel(4)*(quadxr(1)*dr(3)+quadxr(3)*dr(1))
            ! yz
            grdE(5) = grdE(5) + 105.0 * kernel(5) * quadxr_dot_r * dr(2)*dr(3) + &
                              6.0 * kernel(3) * quad(5) - &
                              30.0*kernel(4)*(quadxr(2)*dr(3)+quadxr(3)*dr(2))
            ! zz
            grdE(6) = grdE(6) + 105.0 * kernel(5) * quadxr_dot_r * dr(3)*dr(3) + &
                              6.0 * kernel(3) * quad(6) - &
                              15.0*kernel(4)*(quadxr_dot_r+4.0*quadxr(3)*dr(3))
        end if

        if(do_HE) then
            ! 3
            HE(_xxx_) = HE(_xxx_) + 945*kernel(6)*dr(_x_)**3 * quadxr_dot_r &
                          - 315*kernel(5)*dr(_x_)*(2*quadxr(_x_)*dr(_x_) + quadxr_dot_r) &
                          + 90*kernel(4)*(quad(_xx_)*dr(_x_) + quadxr(_x_))   
            
            HE(_yyy_) = HE(_yyy_) + 945*kernel(6)*dr(_y_)**3 * quadxr_dot_r &
                          - 315*kernel(5)*dr(_y_)*(2*quadxr(_y_)*dr(_y_) + quadxr_dot_r) &
                          + 90*kernel(4)*(quad(_yy_)*dr(_y_) + quadxr(_y_))   
            
            HE(_zzz_) = HE(_zzz_) + 945*kernel(6)*dr(_z_)**3 * quadxr_dot_r &
                          - 315*kernel(5)*dr(_z_)*(2*quadxr(_z_)*dr(_z_) + quadxr_dot_r) &
                          + 90*kernel(4)*(quad(_zz_)*dr(_z_) + quadxr(_z_))   
           
            ! 2 + 1
            HE(_xxy_) = HE(_xxy_) + 945*kernel(6)*dr(_x_)**2*dr(_y_)*quadxr_dot_r &
                          - 105*kernel(5)*(4*quadxr(_x_)*dr(_x_)*dr(_y_) + &
                                           2*quadxr(_y_)*dr(_x_)*dr(_x_) + &
                                           dr(_y_)*quadxr_dot_r) &
                          + 30*kernel(4)*(quad(_xx_)*dr(_y_) + 2*quad(_xy_)*dr(_x_) + quadxr(_y_)) 
            
            HE(_xxz_) = HE(_xxz_) + 945*kernel(6)*dr(_x_)**2*dr(_z_)*quadxr_dot_r &
                          - 105*kernel(5)*(4*quadxr(_x_)*dr(_x_)*dr(_z_) + &
                                           2*quadxr(_z_)*dr(_x_)*dr(_x_) + &
                                           dr(_z_)*quadxr_dot_r) &
                          + 30*kernel(4)*(quad(_xx_)*dr(_z_) + 2*quad(_xz_)*dr(_x_) + quadxr(_z_)) 

            HE(_yyx_) = HE(_yyx_) + 945*kernel(6)*dr(_y_)**2*dr(_x_)*quadxr_dot_r &
                          - 105*kernel(5)*(4*quadxr(_y_)*dr(_y_)*dr(_x_) + &
                                           2*quadxr(_x_)*dr(_y_)*dr(_y_) + &
                                           dr(_x_)*quadxr_dot_r) &
                          + 30*kernel(4)*(quad(_yy_)*dr(_x_) + 2*quad(_xy_)*dr(_y_) + quadxr(_x_)) 
            
            HE(_yyz_) = HE(_yyz_) + 945*kernel(6)*dr(_y_)**2*dr(_z_)*quadxr_dot_r &
                          - 105*kernel(5)*(4*quadxr(_y_)*dr(_y_)*dr(_z_) + &
                                           2*quadxr(_z_)*dr(_y_)*dr(_y_) + &
                                           dr(_z_)*quadxr_dot_r) &
                          + 30*kernel(4)*(quad(_yy_)*dr(_z_) + 2*quad(_zy_)*dr(_y_) + quadxr(_z_)) 

            HE(_zzx_) = HE(_zzx_) + 945*kernel(6)*dr(_z_)**2*dr(_x_)*quadxr_dot_r &
                          - 105*kernel(5)*(4*quadxr(_z_)*dr(_z_)*dr(_x_) + &
                                           2*quadxr(_x_)*dr(_z_)*dr(_z_) + &
                                           dr(_x_)*quadxr_dot_r) &
                          + 30*kernel(4)*(quad(_zz_)*dr(_x_) + 2*quad(_zx_)*dr(_z_) + quadxr(_x_)) 
            
            HE(_zzy_) = HE(_zzy_) + 945*kernel(6)*dr(_z_)**2*dr(_y_)*quadxr_dot_r &
                          - 105*kernel(5)*(4*quadxr(_z_)*dr(_z_)*dr(_y_) + &
                                           2*quadxr(_y_)*dr(_z_)*dr(_z_) + &
                                           dr(_y_)*quadxr_dot_r) &
                          + 30*kernel(4)*(quad(_zz_)*dr(_y_) + 2*quad(_zy_)*dr(_z_) + quadxr(_y_)) 
            
            ! 1 + 1 + 1
            HE(_xyz_) = HE(_xyz_) + 945*kernel(6)*dr(_x_)*dr(_y_)*dr(_z_)*quadxr_dot_r &
                          - 210*kernel(5)*(quadxr(_x_)*dr(_y_)*dr(_z_) + &
                                           quadxr(_y_)*dr(_x_)*dr(_z_) + &
                                           quadxr(_z_)*dr(_x_)*dr(_y_)) &
                          + 30*kernel(4)*(quad(_xy_)*dr(_z_) + quad(_xz_)*dr(_y_) + quad(_yz_)*dr(_x_))  
        end if
    end subroutine quad_elec_prop

    subroutine prepare_fixedelec(eel, arg_dogg) 
        !! This function allocate and populate array of electrostatic 
        !! properties of static multipoles at static multipoles sites.
        !! It should be called blindly before any calculation that requires
        !! V_M2M etc.
   
        use mod_memory, only: mallocate

        implicit none

        type(ommp_electrostatics_type), intent(inout) :: eel
        logical, optional, intent(in) :: arg_dogg

        integer(ip) :: mm_atoms
        logical :: do_gg

        mm_atoms = eel%top%mm_atoms

        !! TODO Improve logic
        do_gg = .false.
        if(present(arg_dogg)) then
            if(arg_dogg) do_gg = .true.
        end if

        if(.not. do_gg .and. eel%M2M_done) return
        if(do_gg .and. eel%M2M_done .and. eel%M2Mgg_done) return

        if(eel%amoeba) then
            if(.not. allocated(eel%V_M2M)) then
                call mallocate('prepare_fixedelec [V_M2M]', mm_atoms, eel%V_M2M)
            end if

            if(.not. allocated(eel%E_M2M)) then
                call mallocate('prepare_fixedelec [E_M2M]', 3, mm_atoms, eel%E_M2M)
            end if

            if(.not. allocated(eel%Egrd_M2M)) then
                call mallocate('prepare_fixedelec [Egrd_M2M]', 6, mm_atoms, eel%Egrd_M2M)
            end if
            
            if(do_gg .and. .not. allocated(eel%EHes_M2M)) then
                call mallocate('prepare_fixedelec [EHes_M2M]', 10, mm_atoms, eel%EHes_M2M)
            end if
            
            eel%V_M2M = 0.0_rp
            eel%E_M2M = 0.0_rp
            eel%Egrd_M2M = 0.0_rp
            if(do_gg) eel%EHes_M2M = 0.0_rp
            
            if(do_gg) then
                call elec_prop_M2M(eel, .true., .true., .true., .true.)
            else
                call elec_prop_M2M(eel, .true., .true., .true., .false.)
            end if
        else
            if(.not. allocated(eel%V_M2M)) then
                call mallocate('prepare_fixedelec [V_M2M]', mm_atoms, eel%V_M2M)
            end if

            if(do_gg .and. .not. allocated(eel%E_M2M)) then
                call mallocate('prepare_fixedelec [E_M2M]', 3, mm_atoms, eel%E_M2M)
            end if

            eel%V_M2M = 0.0_rp
            if(do_gg) eel%E_M2M = 0.0_rp

            if(do_gg) then
                call elec_prop_M2M(eel, .true., .true., .false., .false.)
            else
                call elec_prop_M2M(eel, .true., .false., .false., .false.)
            end if
        end if
        
        eel%M2M_done = .true.
        if(do_gg) eel%M2Mgg_done = .true.

    end subroutine prepare_fixedelec

    subroutine prepare_polelec(eel, arg_dogg)
        use mod_memory, only: mallocate
        implicit none

        type(ommp_electrostatics_type), intent(inout) :: eel
        logical, optional, intent(in) :: arg_dogg

        logical :: do_gg

        do_gg = .false.
        if(present(arg_dogg)) then
            if(arg_dogg) do_gg = .true.
        end if
        
        if(.not. do_gg .and. eel%M2D_done) return
        if(do_gg .and. eel%M2D_done .and. eel%M2Dgg_done) return
        
        if(.not. allocated(eel%E_M2D)) then
            call mallocate('prepare_polelec [E_M2D]', 3, eel%pol_atoms, &
                            eel%n_ipd, eel%E_M2D)
        end if
       
        if(do_gg) then
            if(.not. eel%ipd_done) call fatal_error("IPD should be computed &
                &before computing analytical geometrical gradients of &
                &polarization energy.")
            
            if(.not. allocated(eel%Egrd_M2D)) then
                call mallocate('prepare_polelec [Egrd_M2D]', 6, eel%pol_atoms, &
                                eel%n_ipd, eel%Egrd_M2D)
            end if
            
            if(.not. allocated(eel%Egrd_D2D)) then
                call mallocate('prepare_polelec [Egrd_D2D]', 6, eel%pol_atoms, &
                               eel%n_ipd, eel%Egrd_D2D)
            end if
            
            if(.not. allocated(eel%E_D2M)) then
                call mallocate('prepare_polelec [E_D2M]', 3, eel%top%mm_atoms, &
                               eel%E_D2M)
            end if
            if(.not. allocated(eel%Egrd_D2M) .and. eel%amoeba) then
                call mallocate('prepare_polelec [Egrd_D2M]', 6, eel%top%mm_atoms, &
                               eel%Egrd_D2M)
            end if
            if(.not. allocated(eel%EHes_D2M) .and. eel%amoeba) then
                call mallocate('prepare_polelec [E_D2M]', 10, eel%top%mm_atoms, &
                               eel%EHes_D2M)
            end if

        end if
        
        if(.not. do_gg) then
            eel%E_M2D = 0.0_rp
            call elec_prop_M2D(eel, .false., .true., .false., .false.)
        else
            eel%E_M2D = 0.0_rp
            eel%Egrd_M2D = 0.0_rp
            call elec_prop_M2D(eel, .false., .true., .true., .false.)

            eel%E_D2M = 0.0_rp
            eel%Egrd_D2D = 0.0_rp
            
            if(eel%amoeba) then
                eel%Egrd_D2M = 0.0_rp
                eel%EHes_D2M = 0.0_rp
                call elec_prop_D2M(eel, 'P', .false., .true., .true., .true.)
                call elec_prop_D2M(eel, 'D', .false., .true., .true., .true.)
        
                eel%E_D2M = eel%E_D2M * 0.5
                eel%Egrd_D2M = eel%Egrd_D2M * 0.5
                eel%EHes_D2M = eel%EHes_D2M * 0.5

                call elec_prop_D2D(eel, 'P', .false., .false., .true., .false.)
                call elec_prop_D2D(eel, 'D', .false., .false., .true., .false.)
            else
                call elec_prop_D2M(eel, '-', .false., .true., .false., .false.)
                call elec_prop_D2D(eel, '-', .false., .false., .true., .false.)
            end if
        end if
        
        if(do_gg) eel%M2Dgg_done = .true.
        eel%M2D_done = .true.
    end subroutine prepare_polelec

    subroutine elec_prop_M2M(eel, do_V, do_E, do_Egrd, do_EHes)
        !! Computes the electric potential, field and field gradients of 
        !! static multipoles at all sites (polarizable sites are a 
        !! subset of static ones)
        implicit none
        
        type(ommp_electrostatics_type), intent(inout) :: eel
        !! Electrostatics data structure
        logical, intent(in) :: do_V, do_E, do_Egrd, do_EHes
        !! Flags to enable/disable the calculation of different components


        real(rp) :: kernel(6), dr(3), tmpV, tmpE(3), tmpEgr(6), tmpHE(10), scalf
        integer(ip) :: i, j, idx, ikernel
        logical :: to_do, to_scale
        type(ommp_topology_type), pointer :: top

        top => eel%top

        if(do_EHes) then
            ikernel = 3 
        elseif(do_Egrd) then
            ikernel = 2
        elseif(do_E) then
            ikernel = 1
        elseif(do_V) then
            ikernel = 0
        else
            return
        end if
        if(eel%amoeba) ikernel = ikernel + 2

        if(eel%amoeba) then
            !$omp parallel do default(shared) schedule(dynamic) &
            !$omp private(i,j,idx,to_do,to_scale,scalf,dr,kernel,tmpV,tmpE,tmpEgr,tmpHE)
            do j=1, top%mm_atoms
                ! loop on sources
                do i=1, top%mm_atoms
                    if(j == i) cycle
                    !loop on target

                    to_do = .true.
                    to_scale = .false.
                    scalf = 1.0

                    ! Check if the element should be scaled
                    do idx=eel%list_S_S%ri(i), eel%list_S_S%ri(i+1)-1
                        if(eel%list_S_S%ci(idx) == j) then
                            to_scale = .true.
                            exit
                        end if
                    end do

                    !If it should set the correct variables
                    if(to_scale) then
                        to_do = eel%todo_S_S(idx)
                        scalf = eel%scalef_S_S(idx)
                    end if

                    if(to_do) then
                        dr = top%cmm(:,j) - top%cmm(:, i)
                        call coulomb_kernel(dr, ikernel, kernel)
                        
                        if(do_V) tmpV = 0.0_rp
                        if(do_E) tmpE = 0.0_rp
                        if(do_Egrd) tmpEgr = 0.0_rp
                        if(do_EHes) tmpHE = 0.0_rp

                        call q_elec_prop(eel%q(1,i), dr, kernel, &
                                         do_V, tmpV, & 
                                         do_E, tmpE, &
                                         do_Egrd, tmpEgr, &
                                         do_EHes, tmpHE)

                        call mu_elec_prop(eel%q(2:4,i), dr, kernel, &
                                          do_V, tmpV, & 
                                          do_E, tmpE, &
                                          do_Egrd, tmpEgr, &
                                          do_EHes, tmpHE)

                        call quad_elec_prop(eel%q(5:10,i), dr, kernel, &
                                            do_V, tmpV, & 
                                            do_E, tmpE, &
                                            do_Egrd, tmpEgr, &
                                            do_EHes, tmpHE)

                        if(to_scale) then
                            if(do_V) eel%V_M2M(j) = eel%V_M2M(j) + tmpV * scalf
                            if(do_E) eel%E_M2M(:,j) = eel%E_M2M(:,j) + tmpE * scalf
                            if(do_Egrd) eel%Egrd_M2M(:,j) = eel%Egrd_M2M(:,j) + tmpEgr * scalf
                            if(do_EHes) eel%EHes_M2M(:,j) = eel%EHes_M2M(:,j) + tmpHE * scalf
                        else
                            if(do_V) eel%V_M2M(j) = eel%V_M2M(j) + tmpV
                            if(do_E) eel%E_M2M(:,j) = eel%E_M2M(:,j) + tmpE
                            if(do_Egrd) eel%Egrd_M2M(:,j) = eel%Egrd_M2M(:,j) + tmpEgr
                            if(do_EHes) eel%EHes_M2M(:,j) = eel%EHes_M2M(:,j) + tmpHE
                        end if
                    end if
                end do
            end do
        else
            !$omp parallel do default(shared) schedule(dynamic) &
            !$omp private(i,j,idx,to_do,to_scale,scalf,dr,kernel,tmpV,tmpE,tmpEgr,tmpHE) 
            do j=1, top%mm_atoms
                ! loop on sources
                do i=1, top%mm_atoms
                    if(j == i) cycle
                    !loop on target
                    to_do = .true.
                    to_scale = .false.
                    scalf = 1.0

                    ! Check if the element should be scaled
                    do idx=eel%list_S_S%ri(i), eel%list_S_S%ri(i+1)-1
                        if(eel%list_S_S%ci(idx) == j) then
                            to_scale = .true.
                            exit
                        end if
                    end do

                    !If it should set the correct variables
                    if(to_scale) then
                        to_do = eel%todo_S_S(idx)
                        scalf = eel%scalef_S_S(idx)
                    end if
                    
                    if(to_do) then
                        dr = top%cmm(:,j) - top%cmm(:, i)
                        call coulomb_kernel(dr, ikernel, kernel)
                        
                        if(do_V) tmpV = 0.0_rp
                        if(do_E) tmpE = 0.0_rp
                        if(do_Egrd) tmpEgr = 0.0_rp
                        if(do_EHes) tmpHE = 0.0_rp

                        call q_elec_prop(eel%q(1,i), dr, kernel, &
                                         do_V, tmpV, & 
                                         do_E, tmpE, &
                                         do_Egrd, tmpEgr, &
                                         do_EHes, tmpHE)

                        if(to_scale) then
                            if(do_V) eel%V_M2M(j) = eel%V_M2M(j) + tmpV * scalf
                            if(do_E) eel%E_M2M(:,j) = eel%E_M2M(:,j) + tmpE * scalf
                            if(do_Egrd) eel%Egrd_M2M(:,j) = eel%Egrd_M2M(:,j) + tmpEgr * scalf
                            if(do_EHes) eel%EHes_M2M(:,j) = eel%EHes_M2M(:,j) + tmpHE * scalf
                        else
                            if(do_V) eel%V_M2M(j) = eel%V_M2M(j) + tmpV
                            if(do_E) eel%E_M2M(:,j) = eel%E_M2M(:,j) + tmpE
                            if(do_Egrd) eel%Egrd_M2M(:,j) = eel%Egrd_M2M(:,j) + tmpEgr
                            if(do_EHes) eel%EHes_M2M(:,j) = eel%EHes_M2M(:,j) + tmpHE
                        end if
                    end if
                end do
            end do
        end if
    end subroutine elec_prop_M2M
    
    subroutine field_extD2D(eel, ext_ipd, E)
        !! Computes the electric field of a trial set of induced point dipoles
        !! at polarizable sites. This is intended to be used as matrix-vector
        !! routine in the solution of the linear system.
        
        implicit none

        type(ommp_electrostatics_type), intent(in) :: eel
        !! Data structure for electrostatic part of the system
        real(rp), intent(in) :: ext_ipd(3, eel%pol_atoms)
        !! External induced point dipoles at polarizable sites
        real(rp), intent(inout) :: E(3, eel%pol_atoms)
        !! Electric field (results will be added)

        integer(ip) :: i, j, idx
        logical :: to_scale, to_do
        real(rp) :: kernel(5), dr(3), tmpV, tmpE(3), tmpEgr(6), tmpHE(10), scalf

        !$omp parallel do default(shared) schedule(dynamic) &
        !$omp private(i,j,to_do,to_scale,scalf,idx,tmpV,tmpE,tmpEgr,tmpHE,kernel,dr)
        do j=1, eel%pol_atoms
            do i=1, eel%pol_atoms
                if(j == i) cycle
                !loop on target
                to_do = .true.
                to_scale = .false.
                scalf = 1.0

                ! Check if the element should be scaled
                do idx=eel%list_P_P%ri(i), eel%list_P_P%ri(i+1)-1
                    if(eel%list_P_P%ci(idx) == j) then
                        to_scale = .true.
                        exit
                    end if
                end do

                !If it should set the correct variables
                if(to_scale) then
                    to_do = eel%todo_P_P(idx)
                    scalf = eel%scalef_P_P(idx)
                end if
                
                if(to_do) then
                    call damped_coulomb_kernel(eel, eel%polar_mm(i), &
                                               eel%polar_mm(j),& 
                                               2, kernel(1:3), dr)
                    
                    tmpE = 0.0_rp

                    call mu_elec_prop(ext_ipd(:,i), dr, kernel, .false., tmpV, &
                                      .true., tmpE, .false., tmpEgr, & 
                                      .false., tmpHE)
                    if(to_scale) then
                        E(:, j) = E(:, j) + tmpE * scalf
                    else
                        E(:, j) = E(:, j) + tmpE
                    end if
                end if
            end do
        end do
    end subroutine field_extD2D
    
    subroutine elec_prop_D2D(eel, in_kind, do_V, do_E, do_Egrd, do_EHes)
        !! Computes the electric field of a trial set of induced point dipoles
        !! at polarizable sites. This is intended to be used as matrix-vector
        !! routine in the solution of the linear system.
        implicit none

        type(ommp_electrostatics_type), intent(inout) :: eel
        !! Data structure for electrostatic part of the system
        logical, intent(in) :: do_V, do_E, do_Egrd, do_EHes
        !! Flag to control which properties have to be computed.
        character, intent(in) :: in_kind

        integer(ip) :: i, j, idx, ikernel, knd
        logical :: to_scale, to_do
        real(rp) :: kernel(5), dr(3), tmpV, tmpE(3), tmpEgr(6), tmpHE(10), scalf

        knd = 1 ! Default
        if(in_kind == 'P') then
            knd = _amoeba_P_
        elseif(in_kind == 'D') then 
            knd = _amoeba_D_
        elseif(eel%amoeba) then
            call fatal_error("Unrecognized interaction '"//in_kind//"' in elec&
                             &_prop_D2D.")
        else
            knd = 1
        end if

        if(do_EHes) then
            ikernel = 5 
        elseif(do_Egrd) then
            ikernel = 4
        elseif(do_E) then
            ikernel = 3
        elseif(do_V) then
            ikernel = 2
        else
            return
        end if

        !$omp parallel do default(shared) schedule(dynamic) &
        !$omp private(i,j,idx,dr,kernel,to_do,to_scale,scalf,tmpV,tmpE,tmpEgr,tmpHE) 
        do j=1, eel%pol_atoms
            do i=1, eel%pol_atoms
                if(j == i) cycle
                !loop on target
                to_do = .true.
                to_scale = .false.
                scalf = 1.0

                ! Check if the element should be scaled
                do idx=eel%list_P_P%ri(i), eel%list_P_P%ri(i+1)-1
                    if(eel%list_P_P%ci(idx) == j) then
                        to_scale = .true.
                        exit
                    end if
                end do

                !If it should set the correct variables
                if(to_scale) then
                    to_do = eel%todo_P_P(idx)
                    scalf = eel%scalef_P_P(idx)
                end if
                
                if(to_do) then
                    call damped_coulomb_kernel(eel, eel%polar_mm(i), &
                                               eel%polar_mm(j),& 
                                               ikernel, kernel, dr)
                    
                    if(do_V) tmpV = 0.0_rp
                    if(do_E) tmpE = 0.0_rp
                    if(do_Egrd) tmpEgr = 0.0_rp
                    if(do_EHes) tmpHE = 0.0_rp

                    call mu_elec_prop(eel%ipd(:,i,knd), dr, kernel, &
                                      do_V, tmpV, &
                                      do_E, tmpE, &
                                      do_Egrd, tmpEgr, & 
                                      do_EHes, tmpHE)
                    if(to_scale) then
                        if(do_V) eel%V_D2D(j,knd) = eel%V_D2D(j,knd) + tmpV * scalf
                        if(do_E) eel%E_D2D(:, j,knd) = eel%E_D2D(:, j,knd) + tmpE * scalf
                        if(do_Egrd) eel%Egrd_D2D(:, j,knd) = eel%Egrd_D2D(:, j,knd) + tmpEgr * scalf
                        if(do_EHes) eel%EHes_D2D(:, j,knd) = eel%EHes_D2D(:, j,knd) + tmpHE * scalf
                    else
                        if(do_V) eel%V_D2D(j,knd) = eel%V_D2D(j,knd) + tmpV
                        if(do_E) eel%E_D2D(:, j,knd) = eel%E_D2D(:, j,knd) + tmpE
                        if(do_Egrd) eel%Egrd_D2D(:, j,knd) = eel%Egrd_D2D(:, j,knd) + tmpEgr
                        if(do_EHes) eel%EHes_D2D(:, j,knd) = eel%EHes_D2D(:, j,knd) + tmpHE
                    end if
                end if
            end do
        end do
    end subroutine elec_prop_D2D
    
    subroutine elec_prop_M2D(eel, do_V, do_E, do_Egrd, do_EHes)
        !! Computes the electric field of static multipoles at induced dipoles
        !! sites. This is only intended to be used to build the RHS of the 
        !! linear system. This field is modified by the indroduction of the 
        !! damped kernels and by the connectivity-based screening rules.

        implicit none

        type(ommp_electrostatics_type), intent(inout) :: eel
        !! Electrostatics data structure
        logical, intent(in) :: do_V, do_E, do_Egrd, do_EHes
        !! Flag to control which properties have to be computed.

        integer(ip) :: i, j, idx, ikernel
        logical :: to_do_p, to_scale_p, to_do_d, to_scale_d, to_do, to_scale, &
                   amoeba
        real(rp) :: kernel(5), dr(3), tmpV, tmpE(3), tmpEgr(6), tmpHE(10), &
                    scalf_p, scalf_d, scalf
        type(ommp_topology_type), pointer :: top
        
        ! Shortcuts
        top => eel%top
        amoeba = eel%amoeba

        if(do_EHes) then
            ikernel = 3 
        elseif(do_Egrd) then
            ikernel = 2
        elseif(do_E) then
            ikernel = 1
        elseif(do_V) then
            ikernel = 0
        else
            return
        end if
        if(eel%amoeba) ikernel = ikernel + 2 
        
        if(amoeba) then
            !$omp parallel do default(shared) schedule(dynamic) &
            !$omp private(i,j,idx,dr,kernel,to_do_p,to_do_d,to_scale_p,to_scale_d,scalf_p,scalf_d,tmpV,tmpE,tmpEgr,tmpHE) 
            do j=1, eel%pol_atoms
                ! loop on sources
                do i=1, top%mm_atoms
                    if(eel%polar_mm(j) == i) cycle
                    !loop on target
                    to_do_p = .true.
                    to_scale_p = .false.
                    scalf_p = 1.0

                    ! Check if the element should be scaled
                    do idx=eel%list_S_P_P%ri(i), eel%list_S_P_P%ri(i+1)-1
                        if(eel%list_S_P_P%ci(idx) == j) then
                            to_scale_p = .true.
                            exit
                        end if
                    end do

                    !If it should set the correct variables
                    if(to_scale_p) then
                        to_do_p = eel%todo_S_P_P(idx)
                        scalf_p = eel%scalef_S_P_P(idx)
                    end if
                    
                    to_do_d = .true.
                    to_scale_d = .false.
                    scalf_d = 1.0

                    ! Check if the element should be scaled
                    do idx=eel%list_S_P_D%ri(i), eel%list_S_P_D%ri(i+1)-1
                        if(eel%list_S_P_D%ci(idx) == j) then
                            to_scale_d = .true.
                            exit
                        end if
                    end do

                    !If it should set the correct variables
                    if(to_scale_d) then
                        to_do_d = eel%todo_S_P_D(idx)
                        scalf_d = eel%scalef_S_P_D(idx)
                    end if
                    
                    if(to_do_p .or. to_do_d) then
                        call damped_coulomb_kernel(eel, i, eel%polar_mm(j), &
                                                   ikernel, kernel, dr)
                        
                        if(do_V) tmpV = 0.0_rp
                        if(do_E) tmpE = 0.0_rp
                        if(do_Egrd) tmpEgr = 0.0_rp
                        if(do_EHes) tmpHE = 0.0_rp

                        call q_elec_prop(eel%q(1,i), dr, kernel, &
                                         do_V, tmpV, &
                                         do_E, tmpE, &
                                         do_Egrd, tmpEgr, & 
                                         do_EHes, tmpHE)
                        call mu_elec_prop(eel%q(2:4,i), dr, kernel, &
                                          do_V, tmpV, &
                                          do_E, tmpE, &
                                          do_Egrd, tmpEgr, & 
                                          do_EHes, tmpHE)
                        call quad_elec_prop(eel%q(5:10,i), dr, kernel, &
                                            do_V, tmpV, &
                                            do_E, tmpE, &
                                            do_Egrd, tmpEgr, & 
                                            do_EHes, tmpHE)

                        if(to_do_p) then
                            if(to_scale_p) then
                                if(do_V) eel%V_M2D(j, _amoeba_P_) = eel%V_M2D(j, _amoeba_P_) + tmpV * scalf_p
                                if(do_E) eel%E_M2D(:, j, _amoeba_P_) = eel%E_M2D(:, j, _amoeba_P_) + tmpE * scalf_p
                                if(do_Egrd) eel%Egrd_M2D(:, j, _amoeba_P_) = eel%Egrd_M2D(:, j, _amoeba_P_) + tmpEgr * scalf_p
                                if(do_EHes) eel%EHes_M2D(:, j, _amoeba_P_) = eel%EHes_M2D(:, j, _amoeba_P_) + tmpHE * scalf_p
                            else
                                if(do_V) eel%V_M2D(j, _amoeba_P_) = eel%V_M2D(j, _amoeba_P_) + tmpV 
                                if(do_E) eel%E_M2D(:, j, _amoeba_P_) = eel%E_M2D(:, j, _amoeba_P_) + tmpE
                                if(do_Egrd) eel%Egrd_M2D(:, j, _amoeba_P_) = eel%Egrd_M2D(:, j, _amoeba_P_) + tmpEgr
                                if(do_EHes) eel%EHes_M2D(:, j, _amoeba_P_) = eel%EHes_M2D(:, j, _amoeba_P_) + tmpHE
                            end if
                        end if

                        if(to_do_d) then
                            if(to_scale_d) then
                                if(do_V) eel%V_M2D(j, _amoeba_D_) = eel%V_M2D(j, _amoeba_D_) + tmpV * scalf_d
                                if(do_E) eel%E_M2D(:, j, _amoeba_D_) = eel%E_M2D(:, j, _amoeba_D_) + tmpE * scalf_d
                                if(do_Egrd) eel%Egrd_M2D(:, j, _amoeba_D_) = eel%Egrd_M2D(:, j, _amoeba_D_) + tmpEgr * scalf_d
                                if(do_EHes) eel%EHes_M2D(:, j, _amoeba_D_) = eel%EHes_M2D(:, j, _amoeba_D_) + tmpHE * scalf_d
                            else
                                if(do_V) eel%V_M2D(j, _amoeba_D_) = eel%V_M2D(j, _amoeba_D_) + tmpV 
                                if(do_E) eel%E_M2D(:, j, _amoeba_D_) = eel%E_M2D(:, j, _amoeba_D_) + tmpE
                                if(do_Egrd) eel%Egrd_M2D(:, j, _amoeba_D_) = eel%Egrd_M2D(:, j, _amoeba_D_) + tmpEgr
                                if(do_EHes) eel%EHes_M2D(:, j, _amoeba_D_) = eel%EHes_M2D(:, j, _amoeba_D_) + tmpHE
                            end if
                        end if
                    end if
                end do
            end do
        else
            !$omp parallel do default(shared) schedule(dynamic) &
            !$omp private(i,j,idx,dr,kernel,to_do,to_scale,scalf,tmpV,tmpE,tmpEgr,tmpHE) 
            do j=1, eel%pol_atoms
                ! loop on sources
                do i=1, top%mm_atoms
                    if(eel%polar_mm(j) == i) cycle
                    !loop on target
                    to_do = .true.
                    to_scale = .false.
                    scalf = 1.0

                    ! Check if the element should be scaled
                    do idx=eel%list_S_P_P%ri(i), eel%list_S_P_P%ri(i+1)-1
                        if(eel%list_S_P_P%ci(idx) == j) then
                            to_scale = .true.
                            exit
                        end if
                    end do

                    !If it should set the correct variables
                    if(to_scale) then
                        to_do = eel%todo_S_P_P(idx)
                        scalf = eel%scalef_S_P_P(idx)
                    end if

                    if(to_do) then
                        call damped_coulomb_kernel(eel, i, eel%polar_mm(j), & 
                                                   ikernel, kernel, dr)
                       
                        if(do_V) tmpV = 0.0_rp
                        if(do_E) tmpE = 0.0_rp
                        if(do_Egrd) tmpEgr = 0.0_rp
                        if(do_EHes) tmpHE = 0.0_rp
                        
                        call q_elec_prop(eel%q(1,i), dr, kernel, & 
                                         do_V, tmpV, &
                                         do_E, tmpE, &
                                         do_Egrd, tmpEgr, &
                                         do_EHes, tmpHE)
                        if(to_scale) then
                            if(do_V) eel%V_M2D(j, 1) = eel%V_M2D(j, 1) + tmpV * scalf
                            if(do_E) eel%E_M2D(:, j, 1) = eel%E_M2D(:, j, 1) + tmpE * scalf
                            if(do_Egrd) eel%Egrd_M2D(:, j, 1) = eel%Egrd_M2D(:, j, 1) + tmpEgr * scalf
                            if(do_EHes) eel%EHes_M2D(:, j, 1) = eel%EHes_M2D(:, j, 1) + tmpHE * scalf
                        else
                            if(do_V) eel%V_M2D(j, 1) = eel%V_M2D(j, 1) + tmpV
                            if(do_E) eel%E_M2D(:, j, 1) = eel%E_M2D(:, j, 1) + tmpE
                            if(do_Egrd) eel%Egrd_M2D(:, j, 1) = eel%Egrd_M2D(:, j, 1) + tmpEgr
                            if(do_EHes) eel%EHes_M2D(:, j, 1) = eel%EHes_M2D(:, j, 1) + tmpHE
                        end if 
                    end if
                end do
            end do
        end if

    end subroutine
    
    subroutine elec_prop_D2M(eel, in_kind, do_V, do_E, do_Egrd, do_EHes)

        implicit none

        type(ommp_electrostatics_type), intent(inout) :: eel
        !! Electrostatics data structure
        character, intent(in) :: in_kind

        logical, intent(in) :: do_V, do_E, do_Egrd, do_EHes
        !! Flag to control which properties have to be computed.

        integer(ip) :: i, j, idx, ikernel, knd
        logical :: to_do, to_scale, amoeba
        real(rp) :: kernel(5), dr(3), tmpV, tmpE(3), tmpEgr(6), tmpHE(10), &
                    scalf
        type(ommp_topology_type), pointer :: top
        character :: screening_type
        
        ! Shortcuts
        top => eel%top
        amoeba = eel%amoeba
       
        knd =  1 ! Default
        if(eel%amoeba .and. in_kind /= 'P' .and. in_kind /= 'D') then
            call fatal_error("Unrecognized field '"//in_kind//"' for AMOEBA &
                             &force-field.")
        elseif(eel%amoeba .and. in_kind == 'P') then
            knd = _amoeba_P_
            screening_type = 'D'
        elseif(eel%amoeba .and. in_kind == 'D') then
            knd = _amoeba_D_
            screening_type = 'P'
        elseif(.not. eel%amoeba) then
            screening_type = '-'
        else
            call fatal_error("Unexpected error in elec_prop_D2M.")
        end if

        if(do_EHes) then
            ikernel = 4 
        elseif(do_Egrd) then
            ikernel = 3
        elseif(do_E) then
            ikernel = 2
        elseif(do_V) then 
            ikernel = 1
        else
            return
        end if
        
        if(amoeba) then
            !$omp parallel do default(shared) schedule(dynamic) &
            !$omp private(i,j,idx,dr,kernel,to_do,to_scale,scalf,tmpV,tmpE,tmpEgr,tmpHE) 
            do j=1, top%mm_atoms
                ! loop on sources
                do i=1, eel%pol_atoms
                    if(eel%polar_mm(i) == j) cycle
                    !loop on target
                    if(screening_type == 'P') then
                        to_do = .true.
                        to_scale = .false.
                        scalf = 1.0

                        ! Check if the element should be scaled
                        do idx=eel%list_S_P_P%ri(j), eel%list_S_P_P%ri(j+1)-1
                            if(eel%list_S_P_P%ci(idx) == i) then
                                to_scale = .true.
                                exit
                            end if
                        end do

                        !If it should set the correct variables
                        if(to_scale) then
                            to_do = eel%todo_S_P_P(idx)
                            scalf = eel%scalef_S_P_P(idx)
                        end if
                    else
                        to_do = .true.
                        to_scale = .false.
                        scalf = 1.0

                        ! Check if the element should be scaled
                        do idx=eel%list_S_P_D%ri(j), eel%list_S_P_D%ri(j+1)-1
                            if(eel%list_S_P_D%ci(idx) == i) then
                                to_scale = .true.
                                exit
                            end if
                        end do

                        !If it should set the correct variables
                        if(to_scale) then
                            to_do = eel%todo_S_P_D(idx)
                            scalf = eel%scalef_S_P_D(idx)
                        end if
                    end if

                    
                    if(to_do) then
                        call damped_coulomb_kernel(eel, eel%polar_mm(i), j, & 
                                                   ikernel, kernel, dr)
                       
                        if(do_V) tmpV = 0.0_rp
                        if(do_E) tmpE = 0.0_rp
                        if(do_Egrd) tmpEgr = 0.0_rp
                        if(do_EHes) tmpHE = 0.0_rp
                        
                        call mu_elec_prop(eel%ipd(:,i, knd), dr, kernel, & 
                                         do_V, tmpV, &
                                         do_E, tmpE, &
                                         do_Egrd, tmpEgr, &
                                         do_EHes, tmpHE)
                        if(to_scale) then
                            if(do_V) eel%V_D2M(j) = eel%V_D2M(j) + tmpV * scalf
                            if(do_E) eel%E_D2M(:, j) = eel%E_D2M(:, j) + tmpE * scalf
                            if(do_Egrd) eel%Egrd_D2M(:, j) = eel%Egrd_D2M(:, j) + tmpEgr * scalf
                            if(do_EHes) eel%EHes_D2M(:, j) = eel%EHes_D2M(:, j) + tmpHE * scalf
                        else
                            if(do_V) eel%V_D2M(j) = eel%V_D2M(j) + tmpV
                            if(do_E) eel%E_D2M(:, j) = eel%E_D2M(:, j) + tmpE
                            if(do_Egrd) eel%Egrd_D2M(:, j) = eel%Egrd_D2M(:, j) + tmpEgr
                            if(do_EHes) eel%EHes_D2M(:, j) = eel%EHes_D2M(:, j) + tmpHE
                        end if 
                    end if
                end do
            end do
        else
            !$omp parallel do default(shared) schedule(dynamic) &
            !$omp private(i,j,idx,dr,kernel,to_do,to_scale,scalf,tmpV,tmpE,tmpEgr,tmpHE) 
            do j=1, top%mm_atoms
                ! loop on sources
                do i=1, eel%pol_atoms
                    if(eel%polar_mm(i) == j) cycle
                    !loop on target
                    to_do = .true.
                    to_scale = .false.
                    scalf = 1.0

                    ! Check if the element should be scaled
                    do idx=eel%list_S_P_P%ri(j), eel%list_S_P_P%ri(j+1)-1
                        if(eel%list_S_P_P%ci(idx) == i) then
                            to_scale = .true.
                            exit
                        end if
                    end do

                    !If it should set the correct variables
                    if(to_scale) then
                        to_do = eel%todo_S_P_P(idx)
                        scalf = eel%scalef_S_P_P(idx)
                    end if
                    
                    if(to_do) then
                        call damped_coulomb_kernel(eel, eel%polar_mm(i), j, & 
                                                   ikernel, kernel, dr)
                       
                        if(do_V) tmpV = 0.0_rp
                        if(do_E) tmpE = 0.0_rp
                        if(do_Egrd) tmpEgr = 0.0_rp
                        if(do_EHes) tmpHE = 0.0_rp
                        
                        call mu_elec_prop(eel%ipd(:,i, knd), dr, kernel, & 
                                         do_V, tmpV, &
                                         do_E, tmpE, &
                                         do_Egrd, tmpEgr, &
                                         do_EHes, tmpHE)
                        if(to_scale) then
                            if(do_V) eel%V_D2M(j) = eel%V_D2M(j) + tmpV * scalf
                            if(do_E) eel%E_D2M(:, j) = eel%E_D2M(:, j) + tmpE * scalf
                            if(do_Egrd) eel%Egrd_D2M(:, j) = eel%Egrd_D2M(:, j) + tmpEgr * scalf
                            if(do_EHes) eel%EHes_D2M(:, j) = eel%EHes_D2M(:, j) + tmpHE * scalf
                        else
                            if(do_V) eel%V_D2M(j) = eel%V_D2M(j) + tmpV
                            if(do_E) eel%E_D2M(:, j) = eel%E_D2M(:, j) + tmpE
                            if(do_Egrd) eel%Egrd_D2M(:, j) = eel%Egrd_D2M(:, j) + tmpEgr
                            if(do_EHes) eel%EHes_D2M(:, j) = eel%EHes_D2M(:, j) + tmpHE
                        end if 
                    end if
                end do
            end do
        end if

    end subroutine

    subroutine potential_D2E(eel, cpt, V, amoeba_P_insted_of_D_)
        !! This subroutine computes the potential generated by the induced 
        !! point dipoles to a set of arbitrary coordinates, without applying
        !! any screening rules. Note: for AMOEBA D dipoles should be used. 
        
        implicit none

        type(ommp_electrostatics_type), intent(in) :: eel
        !! Electrostatics data structure
        real(rp), intent(inout) :: V(:)
        !! Electric field (results will be added)
        real(rp), intent(in) :: cpt(:,:)
        !! Coordinates at which the electric field is requested
        logical, optional, intent(in) :: amoeba_P_insted_of_D_
        !! For AMOEBA FF, if true the potential of P dipoles
        !! is computed, otherwise potential of D dipoles is computed

        integer(ip) :: i, j, n_cpt
        logical :: amoeba_P_insted_of_D
        real(rp) :: kernel(5), dr(3), tmpV, tmpE(3), tmpEgr(6), tmpHE(10)

        if(eel%pol_atoms < 1) return

        if(present(amoeba_P_insted_of_D_)) then
            amoeba_P_insted_of_D = amoeba_P_insted_of_D_
        else
            amoeba_P_insted_of_D = .false.
        end if

        if(.not. eel%ipd_done) call fatal_error("IPD should be computed before&
                                                & D2E potential.")
        n_cpt = size(cpt, 2)

        if(eel%amoeba) then
            if(.not. amoeba_P_insted_of_D) then
                !$omp parallel do default(shared) schedule(dynamic) collapse(2) &
                !$omp private(i,j,dr,kernel,tmpV,tmpE,tmpEgr,tmpHE) reduction(+:V)
                do i=1, eel%pol_atoms
                    do j=1, n_cpt
                        dr = cpt(:,j) - eel%cpol(:,i)
                        call coulomb_kernel(dr, 2, kernel(1:2))
                        tmpV = 0.0_rp
                        
                        call mu_elec_prop(eel%ipd(:,i,_amoeba_D_), &
                                        dr, kernel, .true., tmpV, &
                                        .false., tmpE, .false., tmpEgr, & 
                                        .false., tmpHE)

                        V(j) = V(j) + tmpV
                    end do
                end do
            else
                !$omp parallel do default(shared) schedule(dynamic) collapse(2) &
                !$omp private(i,j,dr,kernel,tmpV,tmpE,tmpEgr,tmpHE) reduction(+:V)
                do i=1, eel%pol_atoms
                    do j=1, n_cpt
                        dr = cpt(:,j) - eel%cpol(:,i)
                        call coulomb_kernel(dr, 2, kernel(1:2))
                        tmpV = 0.0_rp
                        
                        call mu_elec_prop(eel%ipd(:,i,_amoeba_P_), &
                                        dr, kernel, .true., tmpV, &
                                        .false., tmpE, .false., tmpEgr, & 
                                        .false., tmpHE)

                        V(j) = V(j) + tmpV
                    end do
                end do
            end if
        else
            !$omp parallel do default(shared) schedule(dynamic) collapse(2) &
            !$omp private(i,j,dr,kernel,tmpV,tmpE,tmpEgr,tmpHE) reduction(+:V)
            do i=1, eel%pol_atoms
                ! loop on sources
                do j=1, n_cpt
                    dr = cpt(:,j) - eel%cpol(:,i)
                    call coulomb_kernel(dr, 2, kernel(1:2))
                    tmpV = 0.0_rp
                    
                    call mu_elec_prop(eel%ipd(:,i,1), &
                                      dr, kernel, .true., tmpV, &
                                      .false., tmpE, .false., tmpEgr, & 
                                      .false., tmpHE)
                    
                    V(j) = V(j) + tmpV
                end do
            end do
        end if
    end subroutine potential_D2E

    subroutine potential_M2E(eel, cpt, V)
        !! This subroutine computes the potential generated by the static
        !! multipoles to a set of arbitrary coordinates, without applying
        !! any screening rules.
        
        implicit none

        type(ommp_electrostatics_type), intent(in) :: eel
        !! Electrostatics data structure
        real(rp), intent(inout) :: V(:)
        !! Electric field (results will be added)
        real(rp), intent(in) :: cpt(:,:)
        !! Coordinates at which the electric field is requested

        integer(ip) :: i, j, n_cpt
        real(rp) :: kernel(5), dr(3), tmpV, tmpE(3), tmpEgr(6), tmpHE(10)

        n_cpt = size(cpt, 2)

        if(eel%amoeba) then
            !$omp parallel do default(shared) schedule(dynamic) collapse(2) &
            !$omp private(i,j,dr,kernel,tmpV,tmpE,tmpEgr,tmpHE) reduction(+:V)
            do i=1, eel%top%mm_atoms
                do j=1, n_cpt
                    dr = cpt(:,j) - eel%top%cmm(:,i)
                    call coulomb_kernel(dr, 2, kernel(1:3))
                    tmpV = 0.0_rp
                    
                    call q_elec_prop(eel%q(1,i), dr, kernel, .true., tmpV, &
                                     .false., tmpE, .false., tmpEgr, & 
                                     .false., tmpHE)
                    call mu_elec_prop(eel%q(2:4,i), dr, kernel, .true., tmpV, &
                                      .false., tmpE, .false., tmpEgr, & 
                                      .false., tmpHE)
                    call quad_elec_prop(eel%q(5:10,i), dr, kernel, .true., tmpV, &
                                        .false., tmpE, .false., tmpEgr, & 
                                        .false., tmpHE)

                    V(j) = V(j) + tmpV
                end do
            end do
        else
            !$omp parallel do default(shared) schedule(dynamic) collapse(2) &
            !$omp private(i,j,dr,kernel,tmpV,tmpE,tmpEgr,tmpHE) reduction(+:V)
            do i=1, eel%top%mm_atoms
                ! loop on sources
                do j=1, n_cpt
                    dr = cpt(:,j) - eel%top%cmm(:,i)
                    call coulomb_kernel(dr, 1, kernel(1:2))
                    tmpV = 0.0_rp
                    
                    call q_elec_prop(eel%q(1,i), dr, kernel, .true., tmpV, &
                                     .false., tmpE, .false., tmpEgr, & 
                                     .false., tmpHE)
                    
                    V(j) = V(j) + tmpV
                end do
            end do
        end if
    end subroutine potential_M2E
    
    subroutine field_D2E(eel, cpt, E)
        !! This subroutine computes the potential generated by the static
        !! multipoles to a set of arbitrary coordinates, without applying
        !! any screening rules.
        
        implicit none

        type(ommp_electrostatics_type), intent(in) :: eel
        !! Electrostatics data structure
        real(rp), intent(inout) :: E(:,:)
        !! Electric field (results will be added)
        real(rp), intent(in) :: cpt(:,:)
        !! Coordinates at which the electric field is requested

        integer(ip) :: i, j, n_cpt
        real(rp) :: kernel(5), dr(3), tmpV, tmpE(3), tmpEgr(6), tmpHE(10)

        if(eel%pol_atoms < 1) return

        if(.not. eel%ipd_done) call fatal_error("IPD should be computed before&
                                                & computing D2E field.")
        n_cpt = size(cpt, 2)

        if(eel%amoeba) then
            !$omp parallel do default(shared) schedule(dynamic) &
            !$omp private(i,j,dr,kernel,tmpV,tmpE,tmpEgr,tmpHE) reduction(+:E)
            do i=1, eel%pol_atoms
                do j=1, n_cpt
                    dr = cpt(:,j) - eel%cpol(:,i)
                    call coulomb_kernel(dr, 3, kernel)
                    tmpE = 0.0_rp
                    !TODO
                    call mu_elec_prop(0.5*(eel%ipd(:,i, _amoeba_P_) + eel%ipd(:,i, _amoeba_D_)), dr, kernel, .false., tmpV, &
                                      .true., tmpE, .false., tmpEgr, & 
                                      .false., tmpHE)

                    E(:,j) = E(:,j) + tmpE 
                end do
            end do
        else
            !$omp parallel do default(shared) schedule(dynamic) &
            !$omp private(i,j,dr,kernel,tmpV,tmpE,tmpEgr,tmpHE) reduction(+:E)
            do i=1, eel%pol_atoms
                ! loop on sources
                do j=1, n_cpt
                    dr = cpt(:,j) - eel%cpol(:,i)
                    call coulomb_kernel(dr, 3, kernel)
                    tmpE = 0.0_rp
                    
                    call mu_elec_prop(eel%ipd(:,i,1), dr, kernel, .false., tmpV, &
                                     .true., tmpE, .false., tmpEgr, & 
                                     .false., tmpHE)
                    
                    E(:,j) = E(:,j) + tmpE
                end do
            end do
        end if
    end subroutine field_D2E

    subroutine field_M2E(eel, cpt, E)
        !! This subroutine computes the potential generated by the static
        !! multipoles to a set of arbitrary coordinates, without applying
        !! any screening rules.
        
        implicit none

        type(ommp_electrostatics_type), intent(in) :: eel
        !! Electrostatics data structure
        real(rp), intent(inout) :: E(:,:)
        !! Electric field (results will be added)
        real(rp), intent(in) :: cpt(:,:)
        !! Coordinates at which the electric field is requested

        integer(ip) :: i, j, n_cpt
        real(rp) :: kernel(5), dr(3), tmpV, tmpE(3), tmpEgr(6), tmpHE(10)

        n_cpt = size(cpt, 2)

        if(eel%amoeba) then
            !$omp parallel do default(shared) schedule(dynamic) &
            !$omp private(i,j,dr,kernel,tmpE,tmpV,tmpEgr,tmpHE) reduction(+:E)
            do i=1, eel%top%mm_atoms
                do j=1, n_cpt
                    dr = cpt(:,j) - eel%top%cmm(:,i)
                    call coulomb_kernel(dr, 4, kernel)
                    tmpE = 0.0_rp
                    
                    call q_elec_prop(eel%q(1,i), dr, kernel, .false., tmpV, &
                                     .true., tmpE, .false., tmpEgr, & 
                                     .false., tmpHE)
                    call mu_elec_prop(eel%q(2:4,i), dr, kernel, .false., tmpV, &
                                      .true., tmpE, .false., tmpEgr, & 
                                      .false., tmpHE)
                    call quad_elec_prop(eel%q(5:10,i), dr, kernel, .false., tmpV, &
                                        .true., tmpE, .false., tmpEgr, & 
                                        .false., tmpHE)

                    E(:,j) = E(:,j) + tmpE
                end do
            end do
        else
            !$omp parallel do default(shared) schedule(dynamic) &
            !$omp private(i,j,dr,kernel,tmpE,tmpV,tmpEgr,tmpHE) reduction(+:E)
            do i=1, eel%top%mm_atoms
                ! loop on sources
                do j=1, n_cpt
                    dr = cpt(:,j) - eel%top%cmm(:,i)
                    call coulomb_kernel(dr, 2, kernel)
                    tmpE = 0.0_rp
                    
                    call q_elec_prop(eel%q(1,i), dr, kernel, .false., tmpV, &
                                     .true., tmpE, .false., tmpEgr, & 
                                     .false., tmpHE)
                    
                    E(:,j) = E(:,j) + tmpE
                end do
            end do
        end if
    end subroutine field_M2E
    
    function screening_rules(eel, i, kind_i, j, kind_j, in_field) result(scalf)
        !! Utility function used to decide if an interaction between sites i and j
        !! should be computed and eventually scaled by a factor.
        !! This function is intended to be used in \(\mathcalO(n^2)\) code, for
        !! linear scaling code lists should be built.
        !! This is written to minimize code repetitions, all the screening rules
        !! are handled in two possible cases: 
        !! 1. rules based on adjacency matrix
        !! 2. rules based on AMOEBA polarization groups

        use mod_constants, only: eps_rp

        type(ommp_electrostatics_type), intent(in) :: eel
        !! Electrostatics object
        integer(ip), intent(in) :: i
        !! Index of source site (MM index is used for static sites, while
        !! Pol index is used for polarizable sites)
        integer(ip), intent(in) :: j
        !! Index of target site (MM index is used for static sites, while
        !! Pol index is used for polarizable sites)
        character, intent(in) :: in_field
        !! Which screening rules have to be applied? 'D' = screening rules
        !! for direct field; 'P' = screening rules for polarization field
        character, intent(in) :: kind_i, kind_j
        !! Type of sites i and j in the interaction for which the screening
        !! rules are required; possible choices are 'P' (polarizable site) or 
        !! 'S' (static site). Any other option will cause a fatal error.
        real(rp) :: scalf
        !! Scale factor for the interaction

        integer(ip) :: ineigh, grp, igrp, pg_i
        integer(ip) :: j_mm, i_mm
        character :: field, interaction
        logical :: amoeba
        type(ommp_topology_type), pointer :: top 
        real(rp) :: myscale(4), myscale_intra(4)

        ! Some shortcut
        top => eel%top
        amoeba = eel%amoeba

        ! Decide which kind of rule should be used
        if(kind_i == 'P') then
            i_mm = eel%polar_mm(i)
        else if(kind_i == 'S') then
            i_mm = i
        else
            call fatal_error('Unexpected value of kind_i in screening_rules')
            i_mm = 0
        end if

        if(kind_j == 'P') then
            j_mm = eel%polar_mm(j)
        else if(kind_j == 'S') then
            j_mm = j
        else
            call fatal_error('Unexpected value of kind_j in screening_rules')
            j_mm = 0
        end if

        myscale = 0.0_rp
        if(kind_j == 'P' .and. kind_i == 'P') then
            ! Use IPD-IPD screening rules
            myscale = eel%uscale
            interaction = 'P' !Pol-pol interaction is named P
        else if(kind_j == 'S' .and. kind_i == 'S') then
            ! Use static multipoles-static multipoles screening rules
            myscale = eel%mscale
            field = '-'
            interaction = 'S' !Static-static interaction is named S
        else
            ! Use static multipoles-IPD screening rules
            if(in_field == 'P' .and. amoeba) then
                myscale = eel%pscale
                myscale_intra = eel%pscale_intra
                field = 'P'
            else if(in_field == 'D' .and. amoeba) then
                myscale = eel%dscale
                field = 'D'
            else if(.not. amoeba) then
                myscale = eel%pscale
                field = 'P'
            else
                call fatal_error('Unexpected value of field in screening rules')
            end if
            interaction = 'M' !Mixed interaction is named M
        end if

        ! Default return value
        scalf = 1.0_rp
        
        if((.not. amoeba) .or. &
           (amoeba .and. interaction == 'M' .and. field == 'P') .or. &
           (amoeba .and. interaction == 'S')) then
            ! Screening rules based on connectivity matrix: 
            ! 1. all the ones of non-amoeba FF, 
            ! 2. the static to IPD polarization field
            ! 3. the static-static interaction of AMOEBA
            do ineigh=1,4
                ! Look if j is at distance ineigh from i
                if(any(top%conn(ineigh)%ci(top%conn(ineigh)%ri(i_mm): &
                                           top%conn(ineigh)%ri(i_mm+1)-1) == j_mm)) then
                   
                    if(interaction == 'M' .and. amoeba) then
                        ! Amoeba uses two different screening rules for atoms 
                        ! of the same group and for atoms of different group
                        if(eel%mmat_polgrp(i_mm) == eel%mmat_polgrp(j_mm)) then
                            scalf = myscale_intra(ineigh)
                        else
                            scalf = myscale(ineigh)
                        end if
                    else
                        ! In any other case just use the scale factor for this
                        ! distance
                        scalf = myscale(ineigh)
                    end if
                    
                    ! Exit the loop
                    exit 

                end if
            end do
        else if((amoeba .and. interaction == 'M' .and. field == 'D') .or. &
                (amoeba .and. interaction == 'P')) then
            ! Screening rules based on polarization groups: 
            ! 1. the static to IPD direct field
            ! 2. the IPD to IPD interaction of AMOEBA
            pg_i = eel%mmat_polgrp(i_mm)

            outer: do ineigh=1,4
                do igrp=eel%pg_conn(ineigh)%ri(pg_i), &
                        eel%pg_conn(ineigh)%ri(pg_i+1)-1
                    ! Indexes of groups at distance ineigh from group pg_i
                    grp = eel%pg_conn(ineigh)%ci(igrp)
                    
                    if(any(eel%polgrp_mmat%ci(eel%polgrp_mmat%ri(grp): &
                                              eel%polgrp_mmat%ri(grp+1)-1) == j_mm)) then
                        ! If atom j is in a group at distance ineigh from the 
                        ! one of atom i, the field is scaled according to dscale
                        scalf = myscale(ineigh)
                        exit outer
                    end if
                end do
            end do outer
        else 
            ! Unexpected error
            call fatal_error("Unexpected combination of parameter for screening_rules")
        end if

    end function screening_rules

end module mod_electrostatics