cp_ddapc_methods.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 contains information regarding the decoupling/recoupling method of Bloechl
!> \author Teodoro Laino
! **************************************************************************************************
MODULE cp_ddapc_methods
   USE cell_types,                      ONLY: cell_type
   USE cp_log_handling,                 ONLY: cp_logger_get_default_io_unit
   USE input_constants,                 ONLY: weight_type_mass,&
                                              weight_type_unit
   USE input_section_types,             ONLY: section_vals_type,&
                                              section_vals_val_get,&
                                              section_vals_val_set
   USE kahan_sum,                       ONLY: accurate_sum
   USE kinds,                           ONLY: dp
   USE mathconstants,                   ONLY: fourpi,&
                                              oorootpi,&
                                              pi,&
                                              twopi
   USE mathlib,                         ONLY: diamat_all,&
                                              invert_matrix
   USE message_passing,                 ONLY: mp_para_env_type
   USE particle_types,                  ONLY: particle_type
   USE pw_spline_utils,                 ONLY: eval_interp_spl3_pbc
   USE pw_types,                        ONLY: pw_c1d_gs_type,&
                                              pw_r3d_rs_type
   USE spherical_harmonics,             ONLY: legendre
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
   PRIVATE
   LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .false.
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_ddapc_methods'
   PUBLIC :: ddapc_eval_gfunc, &
             build_b_vector, &
             build_der_b_vector, &
             build_a_matrix, &
             build_der_a_matrix_rows, &
             prep_g_dot_rvec_sin_cos, &
             cleanup_g_dot_rvec_sin_cos, &
             ddapc_eval_ami, &
             ewald_ddapc_pot, &
             solvation_ddapc_pot
 
CONTAINS
 
! **************************************************************************************************
!> \brief ...
!> \param gfunc ...
!> \param w ...
!> \param gcut ...
!> \param rho_tot_g ...
!> \param radii ...
! **************************************************************************************************
   SUBROUTINE ddapc_eval_gfunc(gfunc, w, gcut, rho_tot_g, radii)
      REAL(kind=dp), DIMENSION(:, :), POINTER            :: gfunc
      REAL(kind=dp), DIMENSION(:), POINTER               :: w
      REAL(kind=dp), INTENT(IN)                          :: gcut
      TYPE(pw_c1d_gs_type), INTENT(IN)                   :: rho_tot_g
      REAL(kind=dp), DIMENSION(:), POINTER               :: radii
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'ddapc_eval_gfunc'
 
      INTEGER                                            :: e_dim, handle, ig, igauss, s_dim
      REAL(kind=dp)                                      :: g2, gcut2, rc, rc2
 
      CALL timeset(routinen, handle)
      gcut2 = gcut*gcut
      !
      s_dim = rho_tot_g%pw_grid%first_gne0
      e_dim = rho_tot_g%pw_grid%ngpts_cut_local
      ALLOCATE (gfunc(s_dim:e_dim, SIZE(radii)))
      ALLOCATE (w(s_dim:e_dim))
      gfunc = 0.0_dp
      w = 0.0_dp
      DO igauss = 1, SIZE(radii)
         rc = radii(igauss)
         rc2 = rc*rc
         DO ig = s_dim, e_dim
            g2 = rho_tot_g%pw_grid%gsq(ig)
            IF (g2 > gcut2) EXIT
            gfunc(ig, igauss) = exp(-g2*rc2/4.0_dp)
         END DO
      END DO
      DO ig = s_dim, e_dim
         g2 = rho_tot_g%pw_grid%gsq(ig)
         IF (g2 > gcut2) EXIT
         w(ig) = fourpi*(g2 - gcut2)**2/(g2*gcut2)
      END DO
      CALL timestop(handle)
   END SUBROUTINE ddapc_eval_gfunc
 
! **************************************************************************************************
!> \brief Computes the B vector for the solution of the linear system
!> \param bv ...
!> \param gfunc ...
!> \param w ...
!> \param particle_set ...
!> \param radii ...
!> \param rho_tot_g ...
!> \param gcut ...
!> \par History
!>      08.2005 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE build_b_vector(bv, gfunc, w, particle_set, radii, rho_tot_g, gcut)
      REAL(kind=dp), DIMENSION(:), INTENT(INOUT)         :: bv
      REAL(kind=dp), DIMENSION(:, :), POINTER            :: gfunc
      REAL(kind=dp), DIMENSION(:), POINTER               :: w
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      REAL(kind=dp), DIMENSION(:), POINTER               :: radii
      TYPE(pw_c1d_gs_type), INTENT(IN)                   :: rho_tot_g
      REAL(kind=dp), INTENT(IN)                          :: gcut
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'build_b_vector'
 
      COMPLEX(KIND=dp)                                   :: phase
      INTEGER                                            :: e_dim, handle, idim, ig, igauss, igmax, &
                                                            iparticle, s_dim
      REAL(kind=dp)                                      :: arg, g2, gcut2, gvec(3), rvec(3)
      REAL(kind=dp), DIMENSION(:), POINTER               :: my_bv, my_bvw
 
      CALL timeset(routinen, handle)
      NULLIFY (my_bv, my_bvw)
      gcut2 = gcut*gcut
      s_dim = rho_tot_g%pw_grid%first_gne0
      e_dim = rho_tot_g%pw_grid%ngpts_cut_local
      igmax = 0
      DO ig = s_dim, e_dim
         g2 = rho_tot_g%pw_grid%gsq(ig)
         IF (g2 > gcut2) EXIT
         igmax = ig
      END DO
      IF (igmax >= s_dim) THEN
         ALLOCATE (my_bv(s_dim:igmax))
         ALLOCATE (my_bvw(s_dim:igmax))
         !
         DO iparticle = 1, SIZE(particle_set)
            rvec = particle_set(iparticle)%r
            my_bv = 0.0_dp
            DO ig = s_dim, igmax
               gvec = rho_tot_g%pw_grid%g(:, ig)
               arg = dot_product(gvec, rvec)
               phase = cmplx(cos(arg), -sin(arg), kind=dp)
               my_bv(ig) = w(ig)*real(conjg(rho_tot_g%array(ig))*phase, kind=dp)
            END DO
            DO igauss = 1, SIZE(radii)
               idim = (iparticle - 1)*SIZE(radii) + igauss
               DO ig = s_dim, igmax
                  my_bvw(ig) = my_bv(ig)*gfunc(ig, igauss)
               END DO
               bv(idim) = accurate_sum(my_bvw)
            END DO
         END DO
         DEALLOCATE (my_bvw)
         DEALLOCATE (my_bv)
      ELSE
         DO iparticle = 1, SIZE(particle_set)
            DO igauss = 1, SIZE(radii)
               idim = (iparticle - 1)*SIZE(radii) + igauss
               bv(idim) = 0.0_dp
            END DO
         END DO
      END IF
      CALL timestop(handle)
   END SUBROUTINE build_b_vector
 
! **************************************************************************************************
!> \brief Computes the A matrix for the solution of the linear system
!> \param Am ...
!> \param gfunc ...
!> \param w ...
!> \param particle_set ...
!> \param radii ...
!> \param rho_tot_g ...
!> \param gcut ...
!> \param g_dot_rvec_sin ...
!> \param g_dot_rvec_cos ...
!> \par History
!>      08.2005 created [tlaino]
!> \author Teodoro Laino
!> \note NB accept g_dot_rvec_* arrays
! **************************************************************************************************
   SUBROUTINE build_a_matrix(Am, gfunc, w, particle_set, radii, rho_tot_g, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
      REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT)      :: am
      REAL(kind=dp), DIMENSION(:, :), POINTER            :: gfunc
      REAL(kind=dp), DIMENSION(:), POINTER               :: w
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      REAL(kind=dp), DIMENSION(:), POINTER               :: radii
      TYPE(pw_c1d_gs_type), INTENT(IN)                   :: rho_tot_g
      REAL(kind=dp), INTENT(IN)                          :: gcut
      REAL(kind=dp), DIMENSION(:, :), INTENT(IN)         :: g_dot_rvec_sin, g_dot_rvec_cos
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'build_A_matrix'
 
      INTEGER                                            :: e_dim, handle, idim1, idim2, ig, &
                                                            igauss1, igauss2, igmax, iparticle1, &
                                                            iparticle2, istart_g, s_dim
      REAL(kind=dp)                                      :: g2, gcut2, tmp
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: my_am, my_amw
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: gfunc_sq(:, :, :)
 
!NB precalculate as many things outside of the innermost loop as possible, in particular w(ig)*gfunc(ig,igauus1)*gfunc(ig,igauss2)
 
      CALL timeset(routinen, handle)
      gcut2 = gcut*gcut
      s_dim = rho_tot_g%pw_grid%first_gne0
      e_dim = rho_tot_g%pw_grid%ngpts_cut_local
      igmax = 0
      DO ig = s_dim, e_dim
         g2 = rho_tot_g%pw_grid%gsq(ig)
         IF (g2 > gcut2) EXIT
         igmax = ig
      END DO
      IF (igmax >= s_dim) THEN
         ALLOCATE (my_am(s_dim:igmax))
         ALLOCATE (my_amw(s_dim:igmax))
         ALLOCATE (gfunc_sq(s_dim:igmax, SIZE(radii), SIZE(radii)))
 
         DO igauss1 = 1, SIZE(radii)
            DO igauss2 = 1, SIZE(radii)
               gfunc_sq(s_dim:igmax, igauss1, igauss2) = w(s_dim:igmax)*gfunc(s_dim:igmax, igauss1)*gfunc(s_dim:igmax, igauss2)
            END DO
         END DO
 
         DO iparticle1 = 1, SIZE(particle_set)
            DO iparticle2 = iparticle1, SIZE(particle_set)
               DO ig = s_dim, igmax
                  !NB replace explicit dot product and cosine with cos(A+B) formula - much faster
                  my_am(ig) = (g_dot_rvec_cos(ig - s_dim + 1, iparticle1)*g_dot_rvec_cos(ig - s_dim + 1, iparticle2) + &
                               g_dot_rvec_sin(ig - s_dim + 1, iparticle1)*g_dot_rvec_sin(ig - s_dim + 1, iparticle2))
               END DO
               DO igauss1 = 1, SIZE(radii)
                  idim1 = (iparticle1 - 1)*SIZE(radii) + igauss1
                  istart_g = 1
                  IF (iparticle2 == iparticle1) istart_g = igauss1
                  DO igauss2 = istart_g, SIZE(radii)
                     idim2 = (iparticle2 - 1)*SIZE(radii) + igauss2
                     my_amw(s_dim:igmax) = my_am(s_dim:igmax)*gfunc_sq(s_dim:igmax, igauss1, igauss2)
                     !NB no loss of accuracy in my test cases
                     !tmp = accurate_sum(my_Amw)
                     tmp = sum(my_amw)
                     am(idim2, idim1) = tmp
                     am(idim1, idim2) = tmp
                  END DO
               END DO
            END DO
         END DO
         DEALLOCATE (gfunc_sq)
         DEALLOCATE (my_amw)
         DEALLOCATE (my_am)
      END IF
      CALL timestop(handle)
   END SUBROUTINE build_a_matrix
 
! **************************************************************************************************
!> \brief Computes the derivative of B vector for the evaluation of the Pulay forces
!> \param dbv ...
!> \param gfunc ...
!> \param w ...
!> \param particle_set ...
!> \param radii ...
!> \param rho_tot_g ...
!> \param gcut ...
!> \param iparticle0 ...
!> \par History
!>      08.2005 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE build_der_b_vector(dbv, gfunc, w, particle_set, radii, rho_tot_g, gcut, iparticle0)
      REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT)      :: dbv
      REAL(kind=dp), DIMENSION(:, :), POINTER            :: gfunc
      REAL(kind=dp), DIMENSION(:), POINTER               :: w
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      REAL(kind=dp), DIMENSION(:), INTENT(IN)            :: radii
      TYPE(pw_c1d_gs_type), INTENT(IN)                   :: rho_tot_g
      REAL(kind=dp), INTENT(IN)                          :: gcut
      INTEGER, INTENT(IN)                                :: iparticle0
 
      CHARACTER(len=*), PARAMETER :: routinen = 'build_der_b_vector'
 
      COMPLEX(KIND=dp)                                   :: dphase
      INTEGER                                            :: e_dim, handle, idim, ig, igauss, igmax, &
                                                            iparticle, s_dim
      REAL(kind=dp)                                      :: arg, g2, gcut2
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: my_dbvw
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: my_dbv
      REAL(kind=dp), DIMENSION(3)                        :: gvec, rvec
 
      CALL timeset(routinen, handle)
      gcut2 = gcut*gcut
      s_dim = rho_tot_g%pw_grid%first_gne0
      e_dim = rho_tot_g%pw_grid%ngpts_cut_local
      igmax = 0
      DO ig = s_dim, e_dim
         g2 = rho_tot_g%pw_grid%gsq(ig)
         IF (g2 > gcut2) EXIT
         igmax = ig
      END DO
      IF (igmax >= s_dim) THEN
         ALLOCATE (my_dbv(3, s_dim:igmax))
         ALLOCATE (my_dbvw(s_dim:igmax))
         DO iparticle = 1, SIZE(particle_set)
            IF (iparticle /= iparticle0) cycle
            rvec = particle_set(iparticle)%r
            DO ig = s_dim, igmax
               gvec = rho_tot_g%pw_grid%g(:, ig)
               arg = dot_product(gvec, rvec)
               dphase = -cmplx(sin(arg), cos(arg), kind=dp)
               my_dbv(:, ig) = w(ig)*real(conjg(rho_tot_g%array(ig))*dphase, kind=dp)*gvec(:)
            END DO
            DO igauss = 1, SIZE(radii)
               idim = (iparticle - 1)*SIZE(radii) + igauss
               DO ig = s_dim, igmax
                  my_dbvw(ig) = my_dbv(1, ig)*gfunc(ig, igauss)
               END DO
               dbv(idim, 1) = accurate_sum(my_dbvw)
               DO ig = s_dim, igmax
                  my_dbvw(ig) = my_dbv(2, ig)*gfunc(ig, igauss)
               END DO
               dbv(idim, 2) = accurate_sum(my_dbvw)
               DO ig = s_dim, igmax
                  my_dbvw(ig) = my_dbv(3, ig)*gfunc(ig, igauss)
               END DO
               dbv(idim, 3) = accurate_sum(my_dbvw)
            END DO
         END DO
         DEALLOCATE (my_dbvw)
         DEALLOCATE (my_dbv)
      ELSE
         DO iparticle = 1, SIZE(particle_set)
            IF (iparticle /= iparticle0) cycle
            DO igauss = 1, SIZE(radii)
               idim = (iparticle - 1)*SIZE(radii) + igauss
               dbv(idim, 1:3) = 0.0_dp
            END DO
         END DO
      END IF
      CALL timestop(handle)
   END SUBROUTINE build_der_b_vector
 
! **************************************************************************************************
!> \brief Computes the derivative of the A matrix for the evaluation of the
!>      Pulay forces
!> \param dAm ...
!> \param gfunc ...
!> \param w ...
!> \param particle_set ...
!> \param radii ...
!> \param rho_tot_g ...
!> \param gcut ...
!> \param iparticle0 ...
!> \param nparticles ...
!> \param g_dot_rvec_sin ...
!> \param g_dot_rvec_cos ...
!> \par History
!>      08.2005 created [tlaino]
!> \author Teodoro Laino
!> \note NB accept g_dot_rvec_* arrays
! **************************************************************************************************
   SUBROUTINE build_der_a_matrix_rows(dAm, gfunc, w, particle_set, radii, &
                                      rho_tot_g, gcut, iparticle0, nparticles, g_dot_rvec_sin, g_dot_rvec_cos)
      REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: dam
      REAL(kind=dp), DIMENSION(:, :), POINTER            :: gfunc
      REAL(kind=dp), DIMENSION(:), POINTER               :: w
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      REAL(kind=dp), DIMENSION(:), POINTER               :: radii
      TYPE(pw_c1d_gs_type), INTENT(IN)                   :: rho_tot_g
      REAL(kind=dp), INTENT(IN)                          :: gcut
      INTEGER, INTENT(IN)                                :: iparticle0, nparticles
      REAL(kind=dp), DIMENSION(:, :), INTENT(IN)         :: g_dot_rvec_sin, g_dot_rvec_cos
 
      CHARACTER(len=*), PARAMETER :: routinen = 'build_der_A_matrix_rows'
 
      INTEGER                                            :: e_dim, handle, ig, igauss2, igmax, &
                                                            iparticle1, iparticle2, s_dim
      REAL(kind=dp)                                      :: g2, gcut2
 
!NB calculate derivatives for a block of particles, just the row parts (since derivative matrix is symmetric)
!NB Use DGEMM to speed up calculation, can't do accurate_sum() anymore because dgemm does the sum over g
 
      EXTERNAL dgemm
      REAL(kind=dp), ALLOCATABLE :: lhs(:, :), rhs(:, :)
      INTEGER :: nr, np, ng, icomp, ipp
 
      CALL timeset(routinen, handle)
      gcut2 = gcut*gcut
      s_dim = rho_tot_g%pw_grid%first_gne0
      e_dim = rho_tot_g%pw_grid%ngpts_cut_local
      igmax = 0
      DO ig = s_dim, e_dim
         g2 = rho_tot_g%pw_grid%gsq(ig)
         IF (g2 > gcut2) EXIT
         igmax = ig
      END DO
 
      nr = SIZE(radii)
      np = SIZE(particle_set)
      ng = igmax - s_dim + 1
      IF (igmax >= s_dim) THEN
         ALLOCATE (lhs(nparticles*nr, ng))
         ALLOCATE (rhs(ng, np*nr))
 
         ! rhs with first term of sin(g.(rvec1-rvec2))
         ! rhs has all parts that depend on iparticle2
         DO iparticle2 = 1, np
            DO igauss2 = 1, nr
               rhs(1:ng, (iparticle2 - 1)*nr + igauss2) = g_dot_rvec_sin(1:ng, iparticle2)*gfunc(s_dim:igmax, igauss2)
            END DO
         END DO
         DO icomp = 1, 3
            ! create lhs, which has all parts that depend on iparticle1
            DO ipp = 1, nparticles
               iparticle1 = iparticle0 + ipp - 1
               DO ig = s_dim, igmax
                  lhs((ipp - 1)*nr + 1:(ipp - 1)*nr + nr, ig - s_dim + 1) = w(ig)*rho_tot_g%pw_grid%g(icomp, ig)* &
                                                                          gfunc(ig, 1:nr)*g_dot_rvec_cos(ig - s_dim + 1, iparticle1)
               END DO
            END DO ! ipp
            ! do main multiply
            CALL dgemm('N', 'N', nparticles*nr, np*nr, ng, 1.0d0, lhs(1, 1), nparticles*nr, rhs(1, 1), &
                       ng, 0.0d0, dam((iparticle0 - 1)*nr + 1, 1, icomp), np*nr)
            ! do extra multiplies to compensate for missing factor of 2
            DO ipp = 1, nparticles
               iparticle1 = iparticle0 + ipp - 1
               CALL dgemm('N', 'N', nr, nr, ng, 1.0d0, lhs((ipp - 1)*nr + 1, 1), nparticles*nr, rhs(1, (iparticle1 - 1)*nr + 1), &
                          ng, 1.0d0, dam((iparticle1 - 1)*nr + 1, (iparticle1 - 1)*nr + 1, icomp), np*nr)
            END DO
            ! now extra columns to account for factor of 2 in some rhs columns
         END DO ! icomp
 
         ! rhs with second term of sin(g.(rvec1-rvec2))
         ! rhs has all parts that depend on iparticle2
         DO iparticle2 = 1, np
            DO igauss2 = 1, nr
               rhs(1:ng, (iparticle2 - 1)*nr + igauss2) = -g_dot_rvec_cos(1:ng, iparticle2)*gfunc(s_dim:igmax, igauss2)
            END DO
         END DO
         DO icomp = 1, 3
            ! create lhs, which has all parts that depend on iparticle1
            DO ipp = 1, nparticles
               iparticle1 = iparticle0 + ipp - 1
               DO ig = s_dim, igmax
                  lhs((ipp - 1)*nr + 1:(ipp - 1)*nr + nr, ig - s_dim + 1) = w(ig)*rho_tot_g%pw_grid%g(icomp, ig)*gfunc(ig, 1:nr)* &
                                                                            g_dot_rvec_sin(ig - s_dim + 1, iparticle1)
               END DO
            END DO
            ! do main multiply
            CALL dgemm('N', 'N', nparticles*nr, np*nr, ng, 1.0d0, lhs(1, 1), nparticles*nr, rhs(1, 1), &
                       ng, 1.0d0, dam((iparticle0 - 1)*nr + 1, 1, icomp), np*nr)
            ! do extra multiples to compensate for missing factor of 2
            DO ipp = 1, nparticles
               iparticle1 = iparticle0 + ipp - 1
               CALL dgemm('N', 'N', nr, nr, ng, 1.0d0, lhs((ipp - 1)*nr + 1, 1), nparticles*nr, rhs(1, (iparticle1 - 1)*nr + 1), &
                          ng, 1.0d0, dam((iparticle1 - 1)*nr + 1, (iparticle1 - 1)*nr + 1, icomp), np*nr)
            END DO
         END DO
 
         DEALLOCATE (rhs)
         DEALLOCATE (lhs)
      END IF
      CALL timestop(handle)
   END SUBROUTINE build_der_a_matrix_rows
 
! **************************************************************************************************
!> \brief deallocate g_dot_rvec_* arrays
!> \param g_dot_rvec_sin ...
!> \param g_dot_rvec_cos ...
! **************************************************************************************************
   SUBROUTINE cleanup_g_dot_rvec_sin_cos(g_dot_rvec_sin, g_dot_rvec_cos)
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: g_dot_rvec_sin, g_dot_rvec_cos
 
      IF (ALLOCATED(g_dot_rvec_sin)) DEALLOCATE (g_dot_rvec_sin)
      IF (ALLOCATED(g_dot_rvec_cos)) DEALLOCATE (g_dot_rvec_cos)
   END SUBROUTINE cleanup_g_dot_rvec_sin_cos
 
! **************************************************************************************************
!> \brief precompute sin(g.r) and cos(g.r) for quicker evaluations of sin(g.(r1-r2)) and cos(g.(r1-r2))
!> \param rho_tot_g ...
!> \param particle_set ...
!> \param gcut ...
!> \param g_dot_rvec_sin ...
!> \param g_dot_rvec_cos ...
! **************************************************************************************************
   SUBROUTINE prep_g_dot_rvec_sin_cos(rho_tot_g, particle_set, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
      TYPE(pw_c1d_gs_type), INTENT(IN)                   :: rho_tot_g
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      REAL(kind=dp), INTENT(IN)                          :: gcut
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: g_dot_rvec_sin, g_dot_rvec_cos
 
      INTEGER                                            :: e_dim, ig, igmax, iparticle, s_dim
      REAL(kind=dp)                                      :: g2, g_dot_rvec, gcut2, rvec(3)
 
      gcut2 = gcut*gcut
      s_dim = rho_tot_g%pw_grid%first_gne0
      e_dim = rho_tot_g%pw_grid%ngpts_cut_local
      igmax = 0
      DO ig = s_dim, e_dim
         g2 = rho_tot_g%pw_grid%gsq(ig)
         IF (g2 > gcut2) EXIT
         igmax = ig
      END DO
 
      IF (igmax >= s_dim) THEN
         ALLOCATE (g_dot_rvec_sin(1:igmax - s_dim + 1, SIZE(particle_set)))
         ALLOCATE (g_dot_rvec_cos(1:igmax - s_dim + 1, SIZE(particle_set)))
 
         DO iparticle = 1, SIZE(particle_set)
            rvec = particle_set(iparticle)%r
            DO ig = s_dim, igmax
               g_dot_rvec = dot_product(rho_tot_g%pw_grid%g(:, ig), rvec)
               g_dot_rvec_sin(ig - s_dim + 1, iparticle) = sin(g_dot_rvec)
               g_dot_rvec_cos(ig - s_dim + 1, iparticle) = cos(g_dot_rvec)
            END DO
         END DO
      END IF
 
   END SUBROUTINE prep_g_dot_rvec_sin_cos
 
! **************************************************************************************************
!> \brief Computes the inverse AmI of the Am matrix
!> \param GAmI ...
!> \param c0 ...
!> \param gfunc ...
!> \param w ...
!> \param particle_set ...
!> \param gcut ...
!> \param rho_tot_g ...
!> \param radii ...
!> \param iw ...
!> \param Vol ...
!> \par History
!>      12.2005 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE ddapc_eval_ami(GAmI, c0, gfunc, w, particle_set, gcut, &
                             rho_tot_g, radii, iw, Vol)
      REAL(kind=dp), DIMENSION(:, :), POINTER            :: gami
      REAL(kind=dp), INTENT(OUT)                         :: c0
      REAL(kind=dp), DIMENSION(:, :), POINTER            :: gfunc
      REAL(kind=dp), DIMENSION(:), POINTER               :: w
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      REAL(kind=dp), INTENT(IN)                          :: gcut
      TYPE(pw_c1d_gs_type), INTENT(IN)                   :: rho_tot_g
      REAL(kind=dp), DIMENSION(:), POINTER               :: radii
      INTEGER, INTENT(IN)                                :: iw
      REAL(kind=dp), INTENT(IN)                          :: vol
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'ddapc_eval_AmI'
 
      INTEGER                                            :: handle, ndim
      REAL(kind=dp)                                      :: condition_number, inv_error
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: ame, cv
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: am, ami, amw, g_dot_rvec_cos, &
                                                            g_dot_rvec_sin
 
!NB for precomputation of sin(g.r) and cos(g.r)
 
      CALL timeset(routinen, handle)
      ndim = SIZE(particle_set)*SIZE(radii)
      ALLOCATE (am(ndim, ndim))
      ALLOCATE (ami(ndim, ndim))
      ALLOCATE (gami(ndim, ndim))
      ALLOCATE (cv(ndim))
      am = 0.0_dp
      ami = 0.0_dp
      cv = 1.0_dp/vol
      !NB precompute sin(g.r) and cos(g.r) for faster evaluation of cos(g.(r1-r2)) in build_A_matrix()
      CALL prep_g_dot_rvec_sin_cos(rho_tot_g, particle_set, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
      CALL build_a_matrix(am, gfunc, w, particle_set, radii, rho_tot_g, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
      CALL cleanup_g_dot_rvec_sin_cos(g_dot_rvec_sin, g_dot_rvec_cos)
      am(:, :) = am(:, :)/(vol*vol)
      CALL rho_tot_g%pw_grid%para%group%sum(am)
      IF (iw > 0) THEN
         ! Checking conditions numbers and eigenvalues
         ALLOCATE (amw(ndim, ndim))
         ALLOCATE (ame(ndim))
         amw(:, :) = am
         CALL diamat_all(amw, ame)
         condition_number = maxval(abs(ame))/minval(abs(ame))
         WRITE (iw, '(T3,A)') " Eigenvalues of Matrix A:"
         WRITE (iw, '(T3,4E15.8)') ame
         WRITE (iw, '(T3,A,1E15.9)') " Condition number:", condition_number
         IF (condition_number > 1.0e12_dp) THEN
            WRITE (iw, fmt="(/,T2,A)") &
               "WARNING: high condition number => possibly ill-conditioned matrix"
         END IF
         DEALLOCATE (amw)
         DEALLOCATE (ame)
      END IF
      CALL invert_matrix(am, ami, inv_error, "N", improve=.false.)
      IF (iw > 0) THEN
         WRITE (iw, '(T3,A,F15.9)') " Error inverting the A matrix: ", inv_error
      END IF
      c0 = dot_product(cv, matmul(ami, cv))
      DEALLOCATE (am)
      DEALLOCATE (cv)
      gami = ami
      DEALLOCATE (ami)
      CALL timestop(handle)
   END SUBROUTINE ddapc_eval_ami
 
! **************************************************************************************************
!> \brief Evaluates the Ewald term E2 and E3 energy term for the decoupling/coupling
!>      of periodic images
!> \param cp_para_env ...
!> \param coeff ...
!> \param factor ...
!> \param cell ...
!> \param multipole_section ...
!> \param particle_set ...
!> \param M ...
!> \param radii ...
!> \par History
!>      08.2005 created [tlaino]
!> \author Teodoro Laino
!> \note NB receive cp_para_env for parallelization
! **************************************************************************************************
   RECURSIVE SUBROUTINE ewald_ddapc_pot(cp_para_env, coeff, factor, cell, multipole_section, &
                                        particle_set, M, radii)
      TYPE(mp_para_env_type), INTENT(IN)                 :: cp_para_env
      TYPE(pw_r3d_rs_type), INTENT(IN), POINTER          :: coeff
      REAL(kind=dp), INTENT(IN)                          :: factor
      TYPE(cell_type), POINTER                           :: cell
      TYPE(section_vals_type), POINTER                   :: multipole_section
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      REAL(kind=dp), DIMENSION(:, :), POINTER            :: m
      REAL(kind=dp), DIMENSION(:), POINTER               :: radii
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'ewald_ddapc_pot'
 
      INTEGER :: ewmdim, handle, iaxis, idim, idim1, idim2, idimo, igauss1, igauss2, ip1, ip2, &
         iparticle1, iparticle2, istart_g, k1, k2, k3, n_rep, ndim, r1, r2, r3
      INTEGER, DIMENSION(3)                              :: gmax, image_cell, rmax, rmin
      LOGICAL                                            :: analyt
      REAL(kind=dp) :: alpha, eps, ew_neut, fac, fac3, frac_radius, fs, g_ewald, galpha, gsq, &
         gsqi, ij_fac, my_val, r, r2tmp, r_ewald, rc1, rc12, rc2, rc22, rcut, rcut2, t1, tol, tol1
      REAL(kind=dp), DIMENSION(3)                        :: g_index, gvec, ra, rvec, svec
      REAL(kind=dp), DIMENSION(:), POINTER               :: ewm
 
      NULLIFY (ewm)
      CALL timeset(routinen, handle)
      cpassert(.NOT. ASSOCIATED(m))
      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 iaxis = 1, 3
            gmax(iaxis) = nint(0.25_dp + cell%hmat(iaxis, iaxis)*alpha*tol1/pi)
         END DO
      ELSE
         DO iaxis = 1, 3
            gmax(iaxis) = ceiling(alpha*tol1*norm2(cell%hmat(:, iaxis))/pi)
         END DO
      END IF
      fac = 1.e0_dp/cell%deth
      fac3 = fac*pi
      ew_neut = -fac*pi/alpha**2
      !
      ewmdim = SIZE(particle_set)*(SIZE(particle_set) + 1)/2
      ndim = SIZE(particle_set)*SIZE(radii)
      ALLOCATE (ewm(ewmdim))
      ALLOCATE (m(ndim, ndim))
      m = 0.0_dp
      !
      idim = 0
      ewm = 0.0_dp
      DO iparticle1 = 1, SIZE(particle_set)
         ip1 = (iparticle1 - 1)*SIZE(radii)
         DO iparticle2 = 1, iparticle1
            ij_fac = 1.0_dp
            IF (iparticle1 == iparticle2) ij_fac = 0.5_dp
 
            ip2 = (iparticle2 - 1)*SIZE(radii)
            idim = idim + 1
            !NB parallelization, done here so indexing is right
            IF (mod(iparticle1, cp_para_env%num_pe) /= cp_para_env%mepos) cycle
            !
            ! Real-Space Contribution
            !
            my_val = 0.0_dp
            rvec = particle_set(iparticle1)%r - particle_set(iparticle2)%r
            r_ewald = 0.0_dp
            IF (iparticle1 /= iparticle2) THEN
               ra = rvec
               r2tmp = dot_product(ra, ra)
               IF (r2tmp <= rcut2) THEN
                  r = sqrt(r2tmp)
                  t1 = erfc(alpha*r)/r
                  r_ewald = t1
               END IF
            END IF
            svec = matmul(cell%h_inv, rvec)
            DO iaxis = 1, 3
               frac_radius = rcut*norm2(cell%h_inv(iaxis, :))
               rmin(iaxis) = floor(-svec(iaxis) - frac_radius)
               rmax(iaxis) = ceiling(-svec(iaxis) + 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
                        r_ewald = r_ewald + t1*ij_fac
                     END IF
                  END DO
               END DO
            END DO
            !
            ! G-space Contribution
            !
            IF (analyt) THEN
               g_ewald = 0.0_dp
               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)
                        g_ewald = g_ewald + t1*cos(dot_product(gvec, rvec))
                     END DO
                  END DO
               END DO
            ELSE
               g_ewald = eval_interp_spl3_pbc(rvec, coeff)
            END IF
            !
            ! G-EWALD, R-EWALD
            !
            g_ewald = r_ewald + fourpi*g_ewald
            !
            ! Self Contribution
            !
            IF (iparticle1 == iparticle2) THEN
               g_ewald = g_ewald - 2.0_dp*alpha*oorootpi
            END IF
            !
            IF (iparticle1 /= iparticle2) THEN
               ra = rvec
               r = sqrt(dot_product(ra, ra))
               my_val = factor/r
            END IF
            ewm(idim) = my_val - factor*g_ewald
         END DO ! iparticle2
      END DO ! iparticle1
      !NB sum over parallelized contributions of different nodes
      CALL cp_para_env%sum(ewm)
      idim = 0
      DO iparticle2 = 1, SIZE(particle_set)
         ip2 = (iparticle2 - 1)*SIZE(radii)
         idimo = (iparticle2 - 1)
         idimo = idimo*(idimo + 1)/2
         DO igauss2 = 1, SIZE(radii)
            idim2 = ip2 + igauss2
            rc2 = radii(igauss2)
            rc22 = rc2*rc2
            DO iparticle1 = 1, iparticle2
               ip1 = (iparticle1 - 1)*SIZE(radii)
               idim = idimo + iparticle1
               istart_g = 1
               IF (iparticle1 == iparticle2) istart_g = igauss2
               DO igauss1 = istart_g, SIZE(radii)
                  idim1 = ip1 + igauss1
                  rc1 = radii(igauss1)
                  rc12 = rc1*rc1
                  m(idim1, idim2) = ewm(idim) - factor*ew_neut - factor*fac3*(rc12 + rc22)
                  m(idim2, idim1) = m(idim1, idim2)
               END DO
            END DO
         END DO ! iparticle2
      END DO ! iparticle1
      DEALLOCATE (ewm)
      CALL timestop(handle)
   END SUBROUTINE ewald_ddapc_pot
 
! **************************************************************************************************
!> \brief Evaluates the electrostatic potential due to a simple solvation model
!>      Spherical cavity in a dieletric medium
!> \param solvation_section ...
!> \param particle_set ...
!> \param M ...
!> \param radii ...
!> \par History
!>      08.2006 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE solvation_ddapc_pot(solvation_section, particle_set, M, radii)
      TYPE(section_vals_type), POINTER                   :: solvation_section
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      REAL(kind=dp), DIMENSION(:, :), POINTER            :: m
      REAL(kind=dp), DIMENSION(:), POINTER               :: radii
 
      INTEGER :: i, idim, idim1, idim2, igauss1, igauss2, ip1, ip2, iparticle1, iparticle2, &
         istart_g, j, l, lmax, n_rep1, n_rep2, ndim, output_unit, weight
      INTEGER, DIMENSION(:), POINTER                     :: list
      LOGICAL                                            :: fixed_center
      REAL(kind=dp)                                      :: center(3), eps_in, eps_out, factor, &
                                                            mass, mycos, r1, r2, rs, rvec(3)
      REAL(kind=dp), DIMENSION(:), POINTER               :: pos, r0
      REAL(kind=dp), DIMENSION(:, :), POINTER            :: cost, locp
 
      fixed_center = .false.
      output_unit = cp_logger_get_default_io_unit()
      ndim = SIZE(particle_set)*SIZE(radii)
      ALLOCATE (m(ndim, ndim))
      m = 0.0_dp
      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(3))
            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
            CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%FIXED", l_val=fixed_center)
            IF (fixed_center) THEN
               CALL section_vals_val_set(solvation_section, "SPHERE%CENTER%XYZ", &
                                         r_vals_ptr=r0)
            ELSE
               DEALLOCATE (r0)
            END IF
         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 (cost(SIZE(particle_set), SIZE(particle_set)))
      ! Determining the single atomic contribution to the dielectric dipole
      DO i = 1, SIZE(particle_set)
         rvec = particle_set(i)%r - center
         r2 = dot_product(rvec, rvec)
         r1 = sqrt(r2)
         IF (r1 >= rs) THEN
            IF (output_unit > 0) THEN
               WRITE (output_unit, '(A,I6,A)') "Atom number :: ", i, " is out of the solvation sphere"
               WRITE (output_unit, '(2(A,F12.6))') "Distance from the center::", r1, " Radius of the sphere::", rs
            END IF
            cpabort("")
         END IF
         locp(:, i) = 0.0_dp
         IF (r1 /= 0.0_dp) THEN
            DO l = 0, lmax
               locp(l, i) = (r1**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
         ELSE
            ! limit for r->0
            locp(0, i) = (eps_in - eps_out)/(rs*eps_in*eps_out)
         END IF
         pos(i) = r1
      END DO
      ! Particle-Particle potential energy matrix
      cost = 0.0_dp
      DO i = 1, SIZE(particle_set)
         DO j = 1, i
            factor = 0.0_dp
            IF (pos(i)*pos(j) /= 0.0_dp) THEN
               mycos = dot_product(particle_set(i)%r - center, particle_set(j)%r - center)/(pos(i)*pos(j))
               IF (abs(mycos) > 1.0_dp) mycos = sign(1.0_dp, mycos)
               DO l = 0, lmax
                  factor = factor + locp(l, i)*pos(j)**l*legendre(mycos, l, 0)
               END DO
            ELSE
               factor = locp(0, i)
            END IF
            cost(i, j) = factor
            cost(j, i) = factor
         END DO
      END DO
      ! Computes the full potential energy matrix
      idim = 0
      DO iparticle2 = 1, SIZE(particle_set)
         ip2 = (iparticle2 - 1)*SIZE(radii)
         DO igauss2 = 1, SIZE(radii)
            idim2 = ip2 + igauss2
            DO iparticle1 = 1, iparticle2
               ip1 = (iparticle1 - 1)*SIZE(radii)
               istart_g = 1
               IF (iparticle1 == iparticle2) istart_g = igauss2
               DO igauss1 = istart_g, SIZE(radii)
                  idim1 = ip1 + igauss1
                  m(idim1, idim2) = cost(iparticle1, iparticle2)
                  m(idim2, idim1) = m(idim1, idim2)
               END DO
            END DO
         END DO
      END DO
      DEALLOCATE (cost)
      DEALLOCATE (pos)
      DEALLOCATE (locp)
   END SUBROUTINE solvation_ddapc_pot
 
END MODULE cp_ddapc_methods