cp_ddapc_forces.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief Density Derived atomic point charges from a QM calculation
!>      (see J. Chem. Phys. Vol. 103 pp. 7422-7428)
!> \par History
!>      08.2005 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
MODULE cp_ddapc_forces
 
   USE atomic_kind_types,               ONLY: atomic_kind_type,&
                                              get_atomic_kind_set
   USE cell_types,                      ONLY: cell_type
   USE cp_control_types,                ONLY: ddapc_restraint_type
   USE input_constants,                 ONLY: do_ddapc_constraint,&
                                              do_ddapc_restraint,&
                                              weight_type_mass,&
                                              weight_type_unit
   USE input_section_types,             ONLY: section_vals_type,&
                                              section_vals_val_get
   USE kinds,                           ONLY: dp
   USE mathconstants,                   ONLY: fourpi,&
                                              pi,&
                                              rootpi,&
                                              twopi
   USE message_passing,                 ONLY: mp_para_env_type
   USE particle_types,                  ONLY: particle_type
   USE pw_spline_utils,                 ONLY: eval_d_interp_spl3_pbc
   USE pw_types,                        ONLY: pw_r3d_rs_type
   USE qs_environment_types,            ONLY: get_qs_env,&
                                              qs_environment_type
   USE qs_force_types,                  ONLY: qs_force_type
   USE spherical_harmonics,             ONLY: dlegendre,&
                                              legendre
!NB for reducing results of calculations that use dq, which is now spread over nodes
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
   PRIVATE
 
   LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .true.
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_ddapc_forces'
   PUBLIC :: ewald_ddapc_force, &
             reset_ch_pulay, &
             evaluate_restraint_functional, &
             restraint_functional_force, &
             solvation_ddapc_force
 
CONTAINS
 
! **************************************************************************************************
!> \brief Evaluates the Ewald term E2 and E3 energy term for the decoupling/coupling
!>      of periodic images
!> \param qs_env ...
!> \param coeff ...
!> \param apply_qmmm_periodic ...
!> \param factor ...
!> \param multipole_section ...
!> \param cell ...
!> \param particle_set ...
!> \param radii ...
!> \param dq ...
!> \param charges ...
!> \par History
!>      08.2005 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   RECURSIVE SUBROUTINE ewald_ddapc_force(qs_env, coeff, apply_qmmm_periodic, &
                                          factor, multipole_section, cell, particle_set, radii, dq, charges)
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(pw_r3d_rs_type), INTENT(IN), POINTER          :: coeff
      LOGICAL, INTENT(IN)                                :: apply_qmmm_periodic
      REAL(kind=dp), INTENT(IN)                          :: factor
      TYPE(section_vals_type), POINTER                   :: multipole_section
      TYPE(cell_type), POINTER                           :: cell
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      REAL(kind=dp), DIMENSION(:), POINTER               :: radii
      REAL(kind=dp), DIMENSION(:, :, :), OPTIONAL, &
         POINTER                                         :: dq
      REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER     :: charges
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'ewald_ddapc_force'
 
      INTEGER                                            :: handle, idim, ip1, ip2, iparticle1, &
                                                            iparticle2, k1, k2, k3, n_rep, r1, r2, &
                                                            r3
      INTEGER, DIMENSION(3)                              :: gmax, image_cell, rmax, rmin
      LOGICAL                                            :: analyt
      REAL(kind=dp)                                      :: alpha, eps, fac, frac_radius, fs, &
                                                            galpha, gsq, gsqi, ij_fac, q1t, q2t, &
                                                            r, r2tmp, rcut, rcut2, t1, t2, tol, &
                                                            tol1
      REAL(kind=dp), DIMENSION(3)                        :: drvec, g_index, gvec, ra, rvec, svec
      REAL(kind=dp), DIMENSION(:, :), POINTER            :: d_el, m
      TYPE(mp_para_env_type), POINTER                    :: para_env
 
      NULLIFY (d_el, m, para_env)
      CALL timeset(routinen, handle)
      CALL get_qs_env(qs_env, para_env=para_env)
      cpassert(PRESENT(charges))
      cpassert(ASSOCIATED(radii))
      rcut = min(norm2(cell%hmat(:, 1)), norm2(cell%hmat(:, 2)), norm2(cell%hmat(:, 3)))/2.0_dp
      CALL section_vals_val_get(multipole_section, "RCUT", n_rep_val=n_rep)
      IF (n_rep == 1) CALL section_vals_val_get(multipole_section, "RCUT", r_val=rcut)
      CALL section_vals_val_get(multipole_section, "EWALD_PRECISION", r_val=eps)
      CALL section_vals_val_get(multipole_section, "ANALYTICAL_GTERM", l_val=analyt)
      ! The spline interpolation path is only valid for orthorhombic grids.
      analyt = analyt .OR. .NOT. cell%orthorhombic .OR. .NOT. ASSOCIATED(coeff)
      rcut2 = rcut**2
      !
      ! Setting-up parameters for Ewald summation
      !
      eps = min(abs(eps), 0.5_dp)
      tol = sqrt(abs(log(eps*rcut)))
      alpha = sqrt(abs(log(eps*rcut*tol)))/rcut
      galpha = 1.0_dp/(4.0_dp*alpha*alpha)
      tol1 = sqrt(-log(eps*rcut*(2.0_dp*tol*alpha)**2))
      IF (cell%orthorhombic) THEN
         DO idim = 1, 3
            gmax(idim) = nint(0.25_dp + cell%hmat(idim, idim)*alpha*tol1/pi)
         END DO
      ELSE
         DO idim = 1, 3
            gmax(idim) = ceiling(alpha*tol1*norm2(cell%hmat(:, idim))/pi)
         END DO
      END IF
 
      ALLOCATE (d_el(3, SIZE(particle_set)))
      d_el = 0.0_dp
      fac = 1.e0_dp/cell%deth
      !
      DO iparticle1 = 1, SIZE(particle_set)
         !NB parallelization
         IF (mod(iparticle1, para_env%num_pe) /= para_env%mepos) cycle
         ip1 = (iparticle1 - 1)*SIZE(radii)
         q1t = sum(charges(ip1 + 1:ip1 + SIZE(radii)))
         DO iparticle2 = 1, iparticle1
            ij_fac = 1.0_dp
            IF (iparticle1 == iparticle2) ij_fac = 0.5_dp
 
            ip2 = (iparticle2 - 1)*SIZE(radii)
            q2t = sum(charges(ip2 + 1:ip2 + SIZE(radii)))
            !
            ! Real-Space Contribution
            !
            rvec = particle_set(iparticle1)%r - particle_set(iparticle2)%r
            IF (iparticle1 /= iparticle2) THEN
               ra = rvec
               r2tmp = dot_product(ra, ra)
               IF (r2tmp <= rcut2) THEN
                  r = sqrt(r2tmp)
                  t1 = erfc(alpha*r)/r
                  drvec = ra/r*q1t*q2t*factor
                  t2 = -2.0_dp*alpha*exp(-alpha*alpha*r*r)/(r*rootpi) - t1/r
                  d_el(1:3, iparticle1) = d_el(1:3, iparticle1) - t2*drvec
                  d_el(1:3, iparticle2) = d_el(1:3, iparticle2) + t2*drvec
               END IF
            END IF
            svec = matmul(cell%h_inv, rvec)
            DO idim = 1, 3
               frac_radius = rcut*norm2(cell%h_inv(idim, :))
               rmin(idim) = floor(-svec(idim) - frac_radius)
               rmax(idim) = ceiling(-svec(idim) + frac_radius)
            END DO
            DO r1 = rmin(1), rmax(1)
               DO r2 = rmin(2), rmax(2)
                  DO r3 = rmin(3), rmax(3)
                     IF ((r1 == 0) .AND. (r2 == 0) .AND. (r3 == 0)) cycle
                     image_cell = [r1, r2, r3]
                     ra = rvec + matmul(cell%hmat, real(image_cell, kind=dp))
                     r2tmp = dot_product(ra, ra)
                     IF (r2tmp <= rcut2) THEN
                        r = sqrt(r2tmp)
                        t1 = erfc(alpha*r)/r
                        drvec = ra/r*q1t*q2t*factor*ij_fac
                        t2 = -2.0_dp*alpha*exp(-alpha*alpha*r*r)/(r*rootpi) - t1/r
                        d_el(1, iparticle1) = d_el(1, iparticle1) - t2*drvec(1)
                        d_el(2, iparticle1) = d_el(2, iparticle1) - t2*drvec(2)
                        d_el(3, iparticle1) = d_el(3, iparticle1) - t2*drvec(3)
                        d_el(1, iparticle2) = d_el(1, iparticle2) + t2*drvec(1)
                        d_el(2, iparticle2) = d_el(2, iparticle2) + t2*drvec(2)
                        d_el(3, iparticle2) = d_el(3, iparticle2) + t2*drvec(3)
                     END IF
                  END DO
               END DO
            END DO
            !
            ! G-space Contribution
            !
            IF (analyt) THEN
               DO k1 = 0, gmax(1)
                  DO k2 = -gmax(2), gmax(2)
                     DO k3 = -gmax(3), gmax(3)
                        IF (k1 == 0 .AND. k2 == 0 .AND. k3 == 0) cycle
                        fs = 2.0_dp; IF (k1 == 0) fs = 1.0_dp
                        g_index = [real(k1, kind=dp), real(k2, kind=dp), real(k3, kind=dp)]
                        gvec = twopi*matmul(transpose(cell%h_inv), g_index)
                        gsq = dot_product(gvec, gvec)
                        gsqi = fs/gsq
                        t1 = fac*gsqi*exp(-galpha*gsq)
                        t2 = -sin(dot_product(gvec, rvec))*t1*q1t*q2t*factor*fourpi
                        d_el(1:3, iparticle1) = d_el(1:3, iparticle1) - t2*gvec
                        d_el(1:3, iparticle2) = d_el(1:3, iparticle2) + t2*gvec
                     END DO
                  END DO
               END DO
            ELSE
               gvec = eval_d_interp_spl3_pbc(rvec, coeff)*q1t*q2t*factor*fourpi
               d_el(1:3, iparticle1) = d_el(1:3, iparticle1) - gvec
               d_el(1:3, iparticle2) = d_el(1:3, iparticle2) + gvec
            END IF
            IF (iparticle1 /= iparticle2) THEN
               ra = rvec
               r = sqrt(dot_product(ra, ra))
               t2 = -1.0_dp/(r*r)*factor
               drvec = ra/r*q1t*q2t
               d_el(1:3, iparticle1) = d_el(1:3, iparticle1) + t2*drvec
               d_el(1:3, iparticle2) = d_el(1:3, iparticle2) - t2*drvec
            END IF
         END DO ! iparticle2
      END DO ! iparticle1
      !NB parallelization
      CALL para_env%sum(d_el)
      m => qs_env%cp_ddapc_env%Md
      IF (apply_qmmm_periodic) m => qs_env%cp_ddapc_env%Mr
      CALL cp_decpl_ddapc_forces(qs_env, m, charges, dq, d_el, particle_set)
      DEALLOCATE (d_el)
      CALL timestop(handle)
   END SUBROUTINE ewald_ddapc_force
 
! **************************************************************************************************
!> \brief Evaluation of the pulay forces due to the fitted charge density
!> \param qs_env ...
!> \param M ...
!> \param charges ...
!> \param dq ...
!> \param d_el ...
!> \param particle_set ...
!> \par History
!>      08.2005 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE cp_decpl_ddapc_forces(qs_env, M, charges, dq, d_el, particle_set)
      TYPE(qs_environment_type), POINTER                 :: qs_env
      REAL(kind=dp), DIMENSION(:, :), POINTER            :: m
      REAL(kind=dp), DIMENSION(:), POINTER               :: charges
      REAL(kind=dp), DIMENSION(:, :, :), POINTER         :: dq
      REAL(kind=dp), DIMENSION(:, :), POINTER            :: d_el
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
 
      CHARACTER(len=*), PARAMETER :: routinen = 'cp_decpl_ddapc_forces'
 
      INTEGER                                            :: handle, i, iatom, ikind, j, k, natom
      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind, kind_of
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: uv
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: chf
      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
 
      NULLIFY (para_env)
      CALL timeset(routinen, handle)
      natom = SIZE(particle_set)
      CALL get_qs_env(qs_env=qs_env, &
                      atomic_kind_set=atomic_kind_set, &
                      para_env=para_env, &
                      force=force)
      ALLOCATE (chf(3, natom))
      CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
                               atom_of_kind=atom_of_kind, &
                               kind_of=kind_of)
 
      ALLOCATE (uv(SIZE(m, 1)))
      uv(:) = matmul(m, charges)
      DO k = 1, natom
         DO j = 1, 3
            chf(j, k) = dot_product(uv, dq(:, k, j))
         END DO
      END DO
      !NB now that get_ddapc returns dq that's spread over nodes, must reduce chf here
      CALL para_env%sum(chf)
      DO iatom = 1, natom
         ikind = kind_of(iatom)
         i = atom_of_kind(iatom)
         force(ikind)%ch_pulay(1:3, i) = force(ikind)%ch_pulay(1:3, i) + chf(1:3, iatom) + d_el(1:3, iatom)
      END DO
      DEALLOCATE (chf)
      DEALLOCATE (uv)
      CALL timestop(handle)
   END SUBROUTINE cp_decpl_ddapc_forces
 
! **************************************************************************************************
!> \brief Evaluation of the pulay forces due to the fitted charge density
!> \param qs_env ...
!> \par History
!>      08.2005 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE reset_ch_pulay(qs_env)
      TYPE(qs_environment_type), POINTER                 :: qs_env
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'reset_ch_pulay'
 
      INTEGER                                            :: handle, ind
      TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
 
      CALL timeset(routinen, handle)
      CALL get_qs_env(qs_env=qs_env, &
                      force=force)
      DO ind = 1, SIZE(force)
         force(ind)%ch_pulay = 0.0_dp
      END DO
      CALL timestop(handle)
   END SUBROUTINE reset_ch_pulay
 
! **************************************************************************************************
!> \brief computes energy and derivatives given a set of charges
!> \param ddapc_restraint_control ...
!> \param n_gauss ...
!> \param uv derivate of energy wrt the corresponding charge entry
!> \param charges current value of the charges (one number for each gaussian used)
!>
!> \param energy_res energy due to the restraint
!> \par History
!>      02.2006 [Joost VandeVondele]
!>               modified [Teo]
!> \note
!>       should be easy to adapt for other specialized cases
! **************************************************************************************************
   SUBROUTINE evaluate_restraint_functional(ddapc_restraint_control, n_gauss, uv, &
                                            charges, energy_res)
      TYPE(ddapc_restraint_type), INTENT(INOUT)          :: ddapc_restraint_control
      INTEGER, INTENT(in)                                :: n_gauss
      REAL(kind=dp), DIMENSION(:)                        :: uv
      REAL(kind=dp), DIMENSION(:), POINTER               :: charges
      REAL(kind=dp), INTENT(INOUT)                       :: energy_res
 
      INTEGER                                            :: i, ind
      REAL(kind=dp)                                      :: de, order_p
 
! order parameter (i.e. the sum of the charges of the selected atoms)
 
      order_p = 0.0_dp
      DO i = 1, ddapc_restraint_control%natoms
         ind = (ddapc_restraint_control%atoms(i) - 1)*n_gauss
         order_p = order_p + ddapc_restraint_control%coeff(i)*sum(charges(ind + 1:ind + n_gauss))
      END DO
      ddapc_restraint_control%ddapc_order_p = order_p
 
      SELECT CASE (ddapc_restraint_control%functional_form)
      CASE (do_ddapc_restraint)
         ! the restraint energy
         energy_res = ddapc_restraint_control%strength*(order_p - ddapc_restraint_control%target)**2.0_dp
 
         ! derivative of the energy
         de = 2.0_dp*ddapc_restraint_control%strength*(order_p - ddapc_restraint_control%target)
         DO i = 1, ddapc_restraint_control%natoms
            ind = (ddapc_restraint_control%atoms(i) - 1)*n_gauss
            uv(ind + 1:ind + n_gauss) = de*ddapc_restraint_control%coeff(i)
         END DO
      CASE (do_ddapc_constraint)
         energy_res = ddapc_restraint_control%strength*(order_p - ddapc_restraint_control%target)
 
         ! derivative of the energy
         DO i = 1, ddapc_restraint_control%natoms
            ind = (ddapc_restraint_control%atoms(i) - 1)*n_gauss
            uv(ind + 1:ind + n_gauss) = ddapc_restraint_control%strength*ddapc_restraint_control%coeff(i)
         END DO
 
      CASE DEFAULT
         cpabort("")
      END SELECT
 
   END SUBROUTINE evaluate_restraint_functional
 
! **************************************************************************************************
!> \brief computes derivatives for DDAPC restraint
!> \param qs_env ...
!> \param ddapc_restraint_control ...
!> \param dq ...
!> \param charges ...
!> \param n_gauss ...
!> \param particle_set ...
!> \par History
!>      02.2006 [Joost VandeVondele]
!>              modified [Teo]
!> \note
!>       should be easy to adapt for other specialized cases
! **************************************************************************************************
   SUBROUTINE restraint_functional_force(qs_env, ddapc_restraint_control, dq, charges, &
                                         n_gauss, particle_set)
 
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(ddapc_restraint_type), INTENT(INOUT)          :: ddapc_restraint_control
      REAL(kind=dp), DIMENSION(:, :, :), POINTER         :: dq
      REAL(kind=dp), DIMENSION(:), POINTER               :: charges
      INTEGER, INTENT(in)                                :: n_gauss
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
 
      CHARACTER(len=*), PARAMETER :: routinen = 'restraint_functional_force'
 
      INTEGER                                            :: handle, i, iatom, ikind, j, k, natom
      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind, kind_of
      REAL(kind=dp)                                      :: dum
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: uv
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: chf
      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
 
      NULLIFY (para_env)
      CALL timeset(routinen, handle)
      natom = SIZE(particle_set)
      CALL get_qs_env(qs_env=qs_env, &
                      atomic_kind_set=atomic_kind_set, &
                      para_env=para_env, &
                      force=force)
      ALLOCATE (chf(3, natom))
      CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
                               atom_of_kind=atom_of_kind, &
                               kind_of=kind_of)
 
      ALLOCATE (uv(SIZE(dq, 1)))
      uv = 0.0_dp
      CALL evaluate_restraint_functional(ddapc_restraint_control, n_gauss, uv, &
                                         charges, dum)
      DO k = 1, natom
         DO j = 1, 3
            chf(j, k) = dot_product(uv, dq(:, k, j))
         END DO
      END DO
      !NB now that get_ddapc returns dq that's spread over nodes, must reduce chf here
      CALL para_env%sum(chf)
      DO iatom = 1, natom
         ikind = kind_of(iatom)
         i = atom_of_kind(iatom)
         force(ikind)%ch_pulay(1:3, i) = force(ikind)%ch_pulay(1:3, i) + chf(1:3, iatom)
      END DO
      DEALLOCATE (chf)
      DEALLOCATE (uv)
      CALL timestop(handle)
 
   END SUBROUTINE restraint_functional_force
 
! **************************************************************************************************
!> \brief Evaluates the electrostatic potential due to a simple solvation model
!>      Spherical cavity in a dieletric medium
!> \param qs_env ...
!> \param solvation_section ...
!> \param particle_set ...
!> \param radii ...
!> \param dq ...
!> \param charges ...
!> \par History
!>      08.2006 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE solvation_ddapc_force(qs_env, solvation_section, particle_set, &
                                    radii, dq, charges)
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(section_vals_type), POINTER                   :: solvation_section
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      REAL(kind=dp), DIMENSION(:), POINTER               :: radii
      REAL(kind=dp), DIMENSION(:, :, :), OPTIONAL, &
         POINTER                                         :: dq
      REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER     :: charges
 
      INTEGER                                            :: i, ip1, ip2, iparticle1, iparticle2, l, &
                                                            lmax, n_rep1, n_rep2, weight
      INTEGER, DIMENSION(:), POINTER                     :: list
      LOGICAL                                            :: fixed_center
      REAL(kind=dp) :: center(3), dcos1(3), dcos2(3), dpos1(3), dpos2(3), eps_in, eps_out, &
         factor1(3), factor2(3), lr, mass, mycos, pos1, pos1i, pos2, pos2i, ptcos, q1t, q2t, &
         r1(3), r1s, r2(3), r2s, rs, rvec(3)
      REAL(kind=dp), DIMENSION(:), POINTER               :: pos, r0
      REAL(kind=dp), DIMENSION(:, :), POINTER            :: d_el, locp, m
 
      fixed_center = .false.
      NULLIFY (d_el, m)
      eps_in = 1.0_dp
      CALL section_vals_val_get(solvation_section, "EPS_OUT", r_val=eps_out)
      CALL section_vals_val_get(solvation_section, "LMAX", i_val=lmax)
      CALL section_vals_val_get(solvation_section, "SPHERE%RADIUS", r_val=rs)
      CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%XYZ", n_rep_val=n_rep1)
      IF (n_rep1 /= 0) THEN
         CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%XYZ", r_vals=r0)
         center = r0
      ELSE
         CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%ATOM_LIST", &
                                   n_rep_val=n_rep2)
         IF (n_rep2 /= 0) THEN
            CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%ATOM_LIST", i_vals=list)
            CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%WEIGHT_TYPE", i_val=weight)
            ALLOCATE (r0(SIZE(list)))
            SELECT CASE (weight)
            CASE (weight_type_unit)
               r0 = 0.0_dp
               DO i = 1, SIZE(list)
                  r0 = r0 + particle_set(list(i))%r
               END DO
               r0 = r0/real(SIZE(list), kind=dp)
            CASE (weight_type_mass)
               r0 = 0.0_dp
               mass = 0.0_dp
               DO i = 1, SIZE(list)
                  r0 = r0 + particle_set(list(i))%r*particle_set(list(i))%atomic_kind%mass
                  mass = mass + particle_set(list(i))%atomic_kind%mass
               END DO
               r0 = r0/mass
            END SELECT
            center = r0
            DEALLOCATE (r0)
         END IF
      END IF
      cpassert(n_rep1 /= 0 .OR. n_rep2 /= 0)
      ! Potential calculation
      ALLOCATE (locp(0:lmax, SIZE(particle_set)))
      ALLOCATE (pos(SIZE(particle_set)))
      ALLOCATE (d_el(3, SIZE(particle_set)))
      d_el = 0.0_dp
      ! Determining the single atomic contribution to the dielectric dipole
      DO i = 1, SIZE(particle_set)
         rvec = particle_set(i)%r - center
         r2s = dot_product(rvec, rvec)
         r1s = sqrt(r2s)
         locp(:, i) = 0.0_dp
         IF (r1s /= 0.0_dp) THEN
            DO l = 0, lmax
               locp(l, i) = (r1s**l*real(l + 1, kind=dp)*(eps_in - eps_out))/ &
                            (rs**(2*l + 1)*eps_in*(real(l, kind=dp)*eps_in + real(l + 1, kind=dp)*eps_out))
            END DO
         END IF
         pos(i) = r1s
      END DO
      ! Computes the full derivatives of the interaction energy
      DO iparticle1 = 1, SIZE(particle_set)
         ip1 = (iparticle1 - 1)*SIZE(radii)
         q1t = sum(charges(ip1 + 1:ip1 + SIZE(radii)))
         DO iparticle2 = 1, iparticle1
            ip2 = (iparticle2 - 1)*SIZE(radii)
            q2t = sum(charges(ip2 + 1:ip2 + SIZE(radii)))
            !
            r1 = particle_set(iparticle1)%r - center
            r2 = particle_set(iparticle2)%r - center
            pos1 = pos(iparticle1)
            pos2 = pos(iparticle2)
            factor1 = 0.0_dp
            factor2 = 0.0_dp
            IF (pos1*pos2 /= 0.0_dp) THEN
               pos1i = 1.0_dp/pos1
               pos2i = 1.0_dp/pos2
               dpos1 = pos1i*r1
               dpos2 = pos2i*r2
               ptcos = dot_product(r1, r2)
               mycos = ptcos/(pos1*pos2)
               IF (abs(mycos) > 1.0_dp) mycos = sign(1.0_dp, mycos)
               dcos1 = (r2*(pos1*pos2) - pos2*dpos1*ptcos)/(pos1*pos2)**2
               dcos2 = (r1*(pos1*pos2) - pos1*dpos2*ptcos)/(pos1*pos2)**2
 
               DO l = 1, lmax
                  lr = real(l, kind=dp)
                  factor1 = factor1 + lr*locp(l, iparticle2)*pos1**(l - 1)*legendre(mycos, l, 0)*dpos1 &
                            + locp(l, iparticle2)*pos1**l*dlegendre(mycos, l, 0)*dcos1
                  factor2 = factor2 + lr*locp(l, iparticle1)*pos2**(l - 1)*legendre(mycos, l, 0)*dpos2 &
                            + locp(l, iparticle1)*pos2**l*dlegendre(mycos, l, 0)*dcos2
               END DO
            END IF
            factor1 = factor1*q1t*q2t
            factor2 = factor2*q1t*q2t
            d_el(1:3, iparticle1) = d_el(1:3, iparticle1) + 0.5_dp*factor1
            d_el(1:3, iparticle2) = d_el(1:3, iparticle2) + 0.5_dp*factor2
         END DO
      END DO
      DEALLOCATE (pos)
      DEALLOCATE (locp)
      m => qs_env%cp_ddapc_env%Ms
      CALL cp_decpl_ddapc_forces(qs_env, m, charges, dq, d_el, particle_set)
      DEALLOCATE (d_el)
   END SUBROUTINE solvation_ddapc_force
 
END MODULE cp_ddapc_forces