pw_poisson_types.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief functions related to the poisson solver on regular grids
!> \par History
!>      greens_fn: JGH (9-Mar-2001) : include influence_fn into
!>                         greens_fn_type add cell volume
!>                         as indicator for updates
!>      greens_fn: JGH (30-Mar-2001) : Added B-spline routines
!>      pws      : JGH (13-Mar-2001) : new pw_poisson_solver, delete
!>                         pw_greens_fn
!>      12.2004 condensed from pws, greens_fn and green_fns, by apsi and JGH,
!>              made thread safe, new input [fawzi]
!>      14-Mar-2006 : added short range screening function for SE codes
!> \author fawzi
! **************************************************************************************************
MODULE pw_poisson_types
   USE bessel_lib,                      ONLY: bessk0,&
                                              bessk1
   USE dielectric_types,                ONLY: dielectric_parameters
   USE dirichlet_bc_types,              ONLY: dirichlet_bc_parameters
   USE kinds,                           ONLY: dp
   USE mathconstants,                   ONLY: fourpi,&
                                              twopi,&
                                              z_zero
   USE mt_util,                         ONLY: mt0d,&
                                              mt1d,&
                                              mt2d,&
                                              mtin_create_screen_fn
   USE ps_implicit_types,               ONLY: mixed_bc,&
                                              neumann_bc,&
                                              ps_implicit_parameters,&
                                              ps_implicit_release,&
                                              ps_implicit_type
   USE ps_wavelet_types,                ONLY: wavelet0d,&
                                              ps_wavelet_release,&
                                              ps_wavelet_type
   USE pw_grid_types,                   ONLY: pw_grid_type
   USE pw_grids,                        ONLY: pw_grid_release
   USE pw_pool_types,                   ONLY: pw_pool_create,&
                                              pw_pool_p_type,&
                                              pw_pool_release,&
                                              pw_pool_type,&
                                              pw_pools_dealloc
   USE pw_types,                        ONLY: pw_c1d_gs_type,&
                                              pw_r1d_gs_type
   USE realspace_grid_types,            ONLY: realspace_grid_type,&
                                              rs_grid_release
#include "../base/base_uses.f90"
 
   IMPLICIT NONE
   PRIVATE
 
   LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .true.
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pw_poisson_types'
 
   PUBLIC :: pw_poisson_type
   PUBLIC :: greens_fn_type, pw_green_create, &
             pw_green_release
   PUBLIC :: pw_poisson_parameter_type
 
   INTEGER, PARAMETER, PUBLIC               :: pw_poisson_none = 0, &
                                               pw_poisson_periodic = 1, &
                                               pw_poisson_analytic = 2, &
                                               pw_poisson_mt = 3, &
                                               pw_poisson_hockney = 5, &
                                               pw_poisson_multipole = 4, &
                                               pw_poisson_wavelet = 6, &
                                               pw_poisson_implicit = 7
   ! EWALD methods
   INTEGER, PARAMETER, PUBLIC               :: do_ewald_none = 1, &
                                               do_ewald_ewald = 2, &
                                               do_ewald_pme = 3, &
                                               do_ewald_spme = 4
 
   INTEGER, PARAMETER, PUBLIC               :: periodic3d = 1000, &
                                               analytic2d = 1001, &
                                               analytic1d = 1002, &
                                               analytic0d = 1003, &
                                               hockney2d = 1201, &
                                               hockney1d = 1202, &
                                               hockney0d = 1203, &
                                               multipole2d = 1301, &
                                               multipole1d = 1302, &
                                               multipole0d = 1303, &
                                               ps_implicit = 1400
 
! **************************************************************************************************
!> \brief parameters for the poisson solver independet of input_section
!> \author Ole Schuett
! **************************************************************************************************
   TYPE pw_poisson_parameter_type
      INTEGER                        :: solver = pw_poisson_none
 
      INTEGER, DIMENSION(3)          :: periodic = 0
      INTEGER                        :: ewald_type = do_ewald_none
      INTEGER                        :: ewald_o_spline = 0
      REAL(kind=dp)                 :: ewald_alpha = 0.0_dp
 
      REAL(kind=dp)                 :: mt_rel_cutoff = 0.0_dp
      REAL(kind=dp)                 :: mt_alpha = 0.0_dp
 
      INTEGER                        :: wavelet_scf_type = 0
      INTEGER                        :: wavelet_method = wavelet0d
      INTEGER                        :: wavelet_special_dimension = 0
      CHARACTER(LEN=1)               :: wavelet_geocode = "S"
 
      LOGICAL                        :: has_dielectric = .false.
      TYPE(dielectric_parameters)    :: dielectric_params = dielectric_parameters()
      TYPE(ps_implicit_parameters)   :: ps_implicit_params = ps_implicit_parameters()
      TYPE(dirichlet_bc_parameters)  :: dbc_params = dirichlet_bc_parameters()
   END TYPE pw_poisson_parameter_type
 
! **************************************************************************************************
!> \brief environment for the poisson solver
!> \author fawzi
! **************************************************************************************************
   TYPE pw_poisson_type
      INTEGER :: pw_level = 0
      INTEGER :: method = pw_poisson_none
      INTEGER :: used_grid = 0
      LOGICAL :: rebuild = .true.
      TYPE(greens_fn_type), POINTER               :: green_fft => null()
      TYPE(ps_wavelet_type), POINTER               :: wavelet => null()
      TYPE(pw_poisson_parameter_type)             :: parameters = pw_poisson_parameter_type()
      REAL(kind=dp), DIMENSION(3, 3)             :: cell_hmat = 0.0_dp
      TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools => null()
      TYPE(pw_grid_type), POINTER                 :: mt_super_ref_pw_grid => null()
      TYPE(ps_implicit_type), POINTER             :: implicit_env => null()
      TYPE(pw_grid_type), POINTER                 :: dct_pw_grid => null()
      TYPE(realspace_grid_type), POINTER          :: diel_rs_grid => null()
   CONTAINS
      PROCEDURE, PUBLIC, non_overridable :: create => pw_poisson_create
      PROCEDURE, PUBLIC, non_overridable :: release => pw_poisson_release
   END TYPE pw_poisson_type
 
! **************************************************************************************************
!> \brief contains all the informations needed by the fft based poisson solvers
!> \author JGH,Teo,fawzi
! **************************************************************************************************
   TYPE greens_fn_type
      INTEGER :: method = periodic3d
      INTEGER :: special_dimension = 0
      REAL(kind=dp) :: radius = 0.0_dp
      REAL(kind=dp) :: mt_alpha = 1.0_dp
      REAL(kind=dp) :: mt_rel_cutoff = 1.0_dp
      REAL(kind=dp) :: slab_size = 0.0_dp
      REAL(kind=dp) :: alpha = 0.0_dp
      LOGICAL :: p3m = .false.
      INTEGER :: p3m_order = 0
      REAL(kind=dp) :: p3m_alpha = 0.0_dp
      REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: p3m_coeff
      REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: p3m_bm2
      LOGICAL :: sr_screening = .false.
      REAL(kind=dp) :: sr_alpha = 1.0_dp
      REAL(kind=dp) :: sr_rc = 0.0_dp
      TYPE(pw_c1d_gs_type) :: influence_fn = pw_c1d_gs_type()
      TYPE(pw_c1d_gs_type), POINTER :: dct_influence_fn => null()
      TYPE(pw_c1d_gs_type), POINTER :: screen_fn => null()
      TYPE(pw_r1d_gs_type), POINTER :: p3m_charge => null()
   END TYPE greens_fn_type
 
CONTAINS
 
! **************************************************************************************************
!> \brief Allocates and sets up the green functions for the fft based poisson
!>      solvers
!> \param green ...
!> \param poisson_params ...
!> \param cell_hmat ...
!> \param pw_pool ...
!> \param mt_super_ref_pw_grid ...
!> \param dct_pw_grid ...
!> \author Fawzi, based on previous functions by JGH and Teo
! **************************************************************************************************
   SUBROUTINE pw_green_create(green, poisson_params, cell_hmat, pw_pool, &
                              mt_super_ref_pw_grid, dct_pw_grid)
      TYPE(greens_fn_type), INTENT(OUT)                  :: green
      TYPE(pw_poisson_parameter_type), INTENT(IN)        :: poisson_params
      REAL(kind=dp), DIMENSION(3, 3), INTENT(IN)         :: cell_hmat
      TYPE(pw_pool_type), POINTER                        :: pw_pool
      TYPE(pw_grid_type), POINTER                        :: mt_super_ref_pw_grid, dct_pw_grid
 
      INTEGER                                            :: dim, i, ig, iz, n, nz
      REAL(kind=dp)                                      :: g2, g3d, gg, gxy, gz, j0g, j1g, k0g, &
                                                            k1g, rlength, zlength
      REAL(kind=dp), DIMENSION(3)                        :: abc
      TYPE(pw_c1d_gs_type), POINTER                      :: dct_gf
      TYPE(pw_grid_type), POINTER                        :: dct_grid
      TYPE(pw_pool_type), POINTER                        :: pw_pool_xpndd
 
      !CPASSERT(cell%orthorhombic)
      DO i = 1, 3
         abc(i) = cell_hmat(i, i)
      END DO
      dim = count(poisson_params%periodic == 1)
 
      SELECT CASE (poisson_params%solver)
      CASE (pw_poisson_periodic)
         green%method = periodic3d
         IF (dim /= 3) THEN
            cpabort("Illegal combination of periodicity and Poisson solver periodic3d")
         END IF
      CASE (pw_poisson_multipole)
         green%method = multipole0d
         IF (dim /= 0) THEN
            cpabort("Illegal combination of periodicity and Poisson solver mulipole0d")
         END IF
      CASE (pw_poisson_analytic)
         SELECT CASE (dim)
         CASE (0)
            green%method = analytic0d
            green%radius = 0.5_dp*minval(abc)
         CASE (1)
            green%method = analytic1d
            green%special_dimension = maxloc(poisson_params%periodic(1:3), 1)
            green%radius = maxval(abc(1:3))
            DO i = 1, 3
               IF (i == green%special_dimension) cycle
               green%radius = min(green%radius, 0.5_dp*abc(i))
            END DO
         CASE (2)
            green%method = analytic2d
            i = minloc(poisson_params%periodic, 1)
            green%special_dimension = i
            green%slab_size = abc(i)
         CASE (3)
            green%method = periodic3d
         CASE DEFAULT
            cpabort("")
         END SELECT
      CASE (pw_poisson_mt)
         green%MT_rel_cutoff = poisson_params%mt_rel_cutoff
         green%MT_alpha = poisson_params%mt_alpha
         green%MT_alpha = green%MT_alpha/minval(abc)
         SELECT CASE (dim)
         CASE (0)
            green%method = mt0d
            green%radius = 0.5_dp*minval(abc)
         CASE (1)
            green%method = mt1d
            green%special_dimension = maxloc(poisson_params%periodic(1:3), 1)
            green%radius = maxval(abc(1:3))
            DO i = 1, 3
               IF (i == green%special_dimension) cycle
               green%radius = min(green%radius, 0.5_dp*abc(i))
            END DO
         CASE (2)
            green%method = mt2d
            i = minloc(poisson_params%periodic, 1)
            green%special_dimension = i
            green%slab_size = abc(i)
         CASE (3)
            cpabort("Illegal combination of periodicity and Poisson solver (MT)")
         CASE DEFAULT
            cpabort("")
         END SELECT
      CASE (pw_poisson_implicit)
         green%method = ps_implicit
      CASE DEFAULT
         cpabort("An unknown Poisson solver was specified")
      END SELECT
 
      ! allocate influence function,...
      SELECT CASE (green%method)
      CASE (periodic3d, analytic2d, analytic1d, analytic0d, mt2d, mt1d, mt0d, multipole0d, ps_implicit)
         CALL pw_pool%create_pw(green%influence_fn)
 
         IF (poisson_params%ewald_type == do_ewald_spme) THEN
            green%p3m = .true.
            green%p3m_order = poisson_params%ewald_o_spline
            green%p3m_alpha = poisson_params%ewald_alpha
            n = green%p3m_order
            ALLOCATE (green%p3m_coeff(-(n - 1):n - 1, 0:n - 1))
            CALL spme_coeff_calculate(n, green%p3m_coeff)
            NULLIFY (green%p3m_charge)
            ALLOCATE (green%p3m_charge)
            CALL pw_pool%create_pw(green%p3m_charge)
            CALL influence_factor(green)
            CALL calc_p3m_charge(green)
         ELSE
            green%p3m = .false.
         END IF
         !
         SELECT CASE (green%method)
         CASE (mt0d, mt1d, mt2d)
            CALL mtin_create_screen_fn(green%screen_fn, pw_pool=pw_pool, method=green%method, &
                                       alpha=green%MT_alpha, &
                                       special_dimension=green%special_dimension, slab_size=green%slab_size, &
                                       super_ref_pw_grid=mt_super_ref_pw_grid)
         CASE (ps_implicit)
            IF ((poisson_params%ps_implicit_params%boundary_condition == mixed_bc) .OR. &
                (poisson_params%ps_implicit_params%boundary_condition == neumann_bc)) THEN
               CALL pw_pool_create(pw_pool_xpndd, pw_grid=dct_pw_grid)
               NULLIFY (green%dct_influence_fn)
               ALLOCATE (green%dct_influence_fn)
               CALL pw_pool_xpndd%create_pw(green%dct_influence_fn)
               CALL pw_pool_release(pw_pool_xpndd)
            END IF
         END SELECT
 
      CASE DEFAULT
         cpabort("")
      END SELECT
 
      ! initialize influence function
      associate(gf => green%influence_fn, grid => green%influence_fn%pw_grid)
         SELECT CASE (green%method)
         CASE (periodic3d, multipole0d)
 
            DO ig = grid%first_gne0, grid%ngpts_cut_local
               g2 = grid%gsq(ig)
               gf%array(ig) = fourpi/g2
            END DO
            IF (grid%have_g0) gf%array(1) = 0.0_dp
 
         CASE (analytic2d)
 
            iz = green%special_dimension ! iz is the direction with NO PBC
            zlength = green%slab_size ! zlength is the thickness of the cell
            DO ig = grid%first_gne0, grid%ngpts_cut_local
               nz = grid%g_hat(iz, ig)
               g2 = grid%gsq(ig)
               g3d = fourpi/g2
               gxy = max(0.0_dp, g2 - grid%g(iz, ig)*grid%g(iz, ig))
               gg = 0.5_dp*sqrt(gxy)
               gf%array(ig) = g3d*(1.0_dp - (-1.0_dp)**nz*exp(-gg*zlength))
            END DO
            IF (grid%have_g0) gf%array(1) = 0.0_dp
 
         CASE (analytic1d)
            ! see 'ab initio molecular dynamics' table 3.1
            ! iz is the direction of the PBC ( can be 1,2,3 -> x,y,z )
            iz = green%special_dimension
            ! rlength is the radius of the tube
            rlength = green%radius
            DO ig = grid%first_gne0, grid%ngpts_cut_local
               g2 = grid%gsq(ig)
               g3d = fourpi/g2
               gxy = sqrt(max(0.0_dp, g2 - grid%g(iz, ig)*grid%g(iz, ig)))
               gz = abs(grid%g(iz, ig))
               j0g = bessel_j0(rlength*gxy)
               j1g = bessel_j1(rlength*gxy)
               IF (gz > 0) THEN
                  k0g = bessk0(rlength*gz)
                  k1g = bessk1(rlength*gz)
               ELSE
                  k0g = 0
                  k1g = 0
               END IF
               gf%array(ig) = g3d*(1.0_dp + rlength* &
                                   (gxy*j1g*k0g - gz*j0g*k1g))
            END DO
            IF (grid%have_g0) gf%array(1) = 0.0_dp
 
         CASE (analytic0d)
 
            rlength = green%radius ! rlength is the radius of the sphere
            DO ig = grid%first_gne0, grid%ngpts_cut_local
               g2 = grid%gsq(ig)
               gg = sqrt(g2)
               g3d = fourpi/g2
               gf%array(ig) = g3d*(1.0_dp - cos(rlength*gg))
            END DO
            IF (grid%have_g0) &
               gf%array(1) = 0.5_dp*fourpi*rlength*rlength
 
         CASE (mt2d, mt1d, mt0d)
 
            DO ig = grid%first_gne0, grid%ngpts_cut_local
               g2 = grid%gsq(ig)
               g3d = fourpi/g2
               gf%array(ig) = g3d + green%screen_fn%array(ig)
            END DO
            IF (grid%have_g0) &
               gf%array(1) = green%screen_fn%array(1)
 
         CASE (ps_implicit)
 
            DO ig = grid%first_gne0, grid%ngpts_cut_local
               g2 = grid%gsq(ig)
               gf%array(ig) = fourpi/g2
            END DO
            IF (grid%have_g0) gf%array(1) = 0.0_dp
 
            IF (ASSOCIATED(green%dct_influence_fn)) THEN
               dct_gf => green%dct_influence_fn
               dct_grid => green%dct_influence_fn%pw_grid
 
               DO ig = dct_grid%first_gne0, dct_grid%ngpts_cut_local
                  g2 = dct_grid%gsq(ig)
                  dct_gf%array(ig) = fourpi/g2
               END DO
               IF (dct_grid%have_g0) dct_gf%array(1) = 0.0_dp
            END IF
 
         CASE DEFAULT
            cpabort("")
         END SELECT
      END associate
 
   END SUBROUTINE pw_green_create
 
! **************************************************************************************************
!> \brief destroys the type (deallocates data)
!> \param gftype ...
!> \param pw_pool ...
!> \par History
!>      none
!> \author Joost VandeVondele
!>      Teodoro Laino
! **************************************************************************************************
   SUBROUTINE pw_green_release(gftype, pw_pool)
      TYPE(greens_fn_type), INTENT(INOUT)                :: gftype
      TYPE(pw_pool_type), OPTIONAL, POINTER              :: pw_pool
 
      LOGICAL                                            :: can_give_back
 
      can_give_back = PRESENT(pw_pool)
      IF (can_give_back) can_give_back = ASSOCIATED(pw_pool)
      IF (can_give_back) THEN
         CALL pw_pool%give_back_pw(gftype%influence_fn)
         IF (ASSOCIATED(gftype%dct_influence_fn)) THEN
            CALL pw_pool%give_back_pw(gftype%dct_influence_fn)
            DEALLOCATE (gftype%dct_influence_fn)
         END IF
         IF (ASSOCIATED(gftype%screen_fn)) THEN
            CALL pw_pool%give_back_pw(gftype%screen_fn)
            DEALLOCATE (gftype%screen_fn)
         END IF
         IF (ASSOCIATED(gftype%p3m_charge)) THEN
            CALL pw_pool%give_back_pw(gftype%p3m_charge)
            DEALLOCATE (gftype%p3m_charge)
         END IF
      ELSE
         CALL gftype%influence_fn%release()
         IF (ASSOCIATED(gftype%dct_influence_fn)) THEN
            CALL gftype%dct_influence_fn%release()
            DEALLOCATE (gftype%dct_influence_fn)
         END IF
         IF (ASSOCIATED(gftype%screen_fn)) THEN
            CALL gftype%screen_fn%release()
            DEALLOCATE (gftype%screen_fn)
         END IF
         IF (ASSOCIATED(gftype%p3m_charge)) THEN
            CALL gftype%p3m_charge%release()
            DEALLOCATE (gftype%p3m_charge)
         END IF
      END IF
      IF (ALLOCATED(gftype%p3m_bm2)) &
         DEALLOCATE (gftype%p3m_bm2)
      IF (ALLOCATED(gftype%p3m_coeff)) &
         DEALLOCATE (gftype%p3m_coeff)
   END SUBROUTINE pw_green_release
 
! **************************************************************************************************
!> \brief Calculates the influence_factor for the
!>      SPME Green's function in reciprocal space'''
!> \param gftype ...
!> \par History
!>      none
!> \author DH (29-Mar-2001)
! **************************************************************************************************
   SUBROUTINE influence_factor(gftype)
      TYPE(greens_fn_type), INTENT(INOUT)                :: gftype
 
      COMPLEX(KIND=dp)                                   :: b_m, exp_m, sum_m
      INTEGER                                            :: dim, j, k, l, n, pt
      INTEGER, DIMENSION(3)                              :: npts
      INTEGER, DIMENSION(:), POINTER                     :: lb, ub
      REAL(kind=dp)                                      :: l_arg, prod_arg, val
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: m_assign
 
      n = gftype%p3m_order
 
      ! calculate the assignment function values
 
      lb => gftype%influence_fn%pw_grid%bounds(1, :)
      ub => gftype%influence_fn%pw_grid%bounds(2, :)
      IF (ALLOCATED(gftype%p3m_bm2)) THEN
         IF (lbound(gftype%p3m_bm2, 2) /= minval(lb(:)) .OR. &
             ubound(gftype%p3m_bm2, 2) /= maxval(ub(:))) THEN
            DEALLOCATE (gftype%p3m_bm2)
         END IF
      END IF
      IF (.NOT. ALLOCATED(gftype%p3m_bm2)) THEN
         ALLOCATE (gftype%p3m_bm2(3, minval(lb(:)):maxval(ub(:))))
      END IF
 
      ALLOCATE (m_assign(0:n - 2))
      m_assign = 0.0_dp
      DO k = 0, n - 2
         j = -(n - 1) + 2*k
         DO l = 0, n - 1
            l_arg = 0.5_dp**l
            prod_arg = gftype%p3m_coeff(j, l)*l_arg
            m_assign(k) = m_assign(k) + prod_arg
         END DO
      END DO
 
      ! calculate the absolute b values
 
      npts(:) = ub(:) - lb(:) + 1
      DO dim = 1, 3
         DO pt = lb(dim), ub(dim)
            val = twopi*(real(pt, kind=dp)/real(npts(dim), kind=dp))
            exp_m = cmplx(cos(val), -sin(val), kind=dp)
            sum_m = z_zero
            DO k = 0, n - 2
               sum_m = sum_m + m_assign(k)*exp_m**k
            END DO
            b_m = exp_m**(n - 1)/sum_m
            gftype%p3m_bm2(dim, pt) = sqrt(real(b_m*conjg(b_m), kind=dp))
         END DO
      END DO
 
      DEALLOCATE (m_assign)
   END SUBROUTINE influence_factor
 
! **************************************************************************************************
!> \brief ...
!> \param gf ...
! **************************************************************************************************
   PURE SUBROUTINE calc_p3m_charge(gf)
 
      TYPE(greens_fn_type), INTENT(INOUT)                :: gf
 
      INTEGER                                            :: ig, l, m, n
      REAL(kind=dp)                                      :: arg, novol
 
      ! check if charge function is consistent with current box volume
 
      associate(grid => gf%influence_fn%pw_grid, bm2 => gf%p3m_bm2)
         arg = 1.0_dp/(8.0_dp*gf%p3m_alpha**2)
         novol = real(grid%ngpts, kind=dp)/grid%vol
         DO ig = 1, grid%ngpts_cut_local
            l = grid%g_hat(1, ig)
            m = grid%g_hat(2, ig)
            n = grid%g_hat(3, ig)
            gf%p3m_charge%array(ig) = novol*exp(-arg*grid%gsq(ig))* &
                                      bm2(1, l)*bm2(2, m)*bm2(3, n)
         END DO
      END associate
 
   END SUBROUTINE calc_p3m_charge
 
! **************************************************************************************************
!> \brief Initialize the poisson solver
!>      You should call this just before calling the work routine
!>      pw_poisson_solver
!>      Call pw_poisson_release when you have finished
!> \param poisson_env ...
!> \par History
!>      none
!> \author JGH (12-Mar-2001)
! **************************************************************************************************
   SUBROUTINE pw_poisson_create(poisson_env)
 
      CLASS(pw_poisson_type), INTENT(INOUT)              :: poisson_env
 
      ! Cleanup a potential previous poisson_env
      CALL poisson_env%release()
 
   END SUBROUTINE pw_poisson_create
 
! **************************************************************************************************
!> \brief releases the poisson solver
!> \param poisson_env ...
!> \par History
!>      none
!> \author fawzi (11.2002)
! **************************************************************************************************
   SUBROUTINE pw_poisson_release(poisson_env)
 
      CLASS(pw_poisson_type), INTENT(INOUT)               :: poisson_env
 
      IF (ASSOCIATED(poisson_env%pw_pools)) THEN
         CALL pw_pools_dealloc(poisson_env%pw_pools)
      END IF
 
      IF (ASSOCIATED(poisson_env%green_fft)) THEN
         CALL pw_green_release(poisson_env%green_fft)
         DEALLOCATE (poisson_env%green_fft)
      END IF
      CALL pw_grid_release(poisson_env%mt_super_ref_pw_grid)
      CALL ps_wavelet_release(poisson_env%wavelet)
      CALL ps_implicit_release(poisson_env%implicit_env, &
                               poisson_env%parameters%ps_implicit_params)
      CALL pw_grid_release(poisson_env%dct_pw_grid)
      IF (ASSOCIATED(poisson_env%diel_rs_grid)) THEN
         CALL rs_grid_release(poisson_env%diel_rs_grid)
         DEALLOCATE (poisson_env%diel_rs_grid)
      END IF
 
   END SUBROUTINE pw_poisson_release
 
! **************************************************************************************************
!> \brief Calculates the coefficients for the charge assignment function
!> \param n ...
!> \param coeff ...
!> \par History
!>      none
!> \author DG (29-Mar-2001)
! **************************************************************************************************
   SUBROUTINE spme_coeff_calculate(n, coeff)
 
      INTEGER, INTENT(IN)                                :: n
      REAL(kind=dp), DIMENSION(-(n-1):n-1, 0:n-1), &
         INTENT(OUT)                                     :: coeff
 
      INTEGER                                            :: i, j, l, m
      REAL(kind=dp)                                      :: b
      REAL(kind=dp), DIMENSION(n, -n:n, 0:n-1)           :: a
 
      a = 0.0_dp
      a(1, 0, 0) = 1.0_dp
 
      DO i = 2, n
         m = i - 1
         DO j = -m, m, 2
            DO l = 0, m - 1
               b = (a(m, j - 1, l) + &
                    REAL((-1)**l, kind=dp)*a(m, j + 1, l))/ &
                   REAL((l + 1)*2**(l + 1), kind=dp)
               a(i, j, 0) = a(i, j, 0) + b
            END DO
            DO l = 0, m - 1
               a(i, j, l + 1) = (a(m, j + 1, l) - &
                                 a(m, j - 1, l))/real(l + 1, kind=dp)
            END DO
         END DO
      END DO
 
      coeff = 0.0_dp
      DO i = 0, n - 1
         DO j = -(n - 1), n - 1, 2
            coeff(j, i) = a(n, j, i)
         END DO
      END DO
 
   END SUBROUTINE spme_coeff_calculate
 
END MODULE pw_poisson_types