36 ps_implicit_parameters,&
53 #include "../base/base_uses.f90"
58 LOGICAL,
PRIVATE,
PARAMETER :: debug_this_module = .true.
59 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'pw_poisson_types'
61 PUBLIC :: pw_poisson_type
64 PUBLIC :: pw_poisson_parameter_type
96 TYPE pw_poisson_parameter_type
99 INTEGER,
DIMENSION(3) :: periodic = 0
101 INTEGER :: ewald_o_spline = 0
102 REAL(kind=
dp) :: ewald_alpha = 0.0_dp
104 REAL(kind=
dp) :: mt_rel_cutoff = 0.0_dp
105 REAL(kind=
dp) :: mt_alpha = 0.0_dp
107 INTEGER :: wavelet_scf_type = 0
109 INTEGER :: wavelet_special_dimension = 0
110 CHARACTER(LEN=1) :: wavelet_geocode =
"S"
112 LOGICAL :: has_dielectric = .false.
113 TYPE(dielectric_parameters) :: dielectric_params = dielectric_parameters()
114 TYPE(ps_implicit_parameters) :: ps_implicit_params = ps_implicit_parameters()
115 TYPE(dirichlet_bc_parameters) :: dbc_params = dirichlet_bc_parameters()
116 END TYPE pw_poisson_parameter_type
123 INTEGER :: pw_level = 0
125 INTEGER :: used_grid = 0
126 LOGICAL :: rebuild = .true.
127 TYPE(greens_fn_type),
POINTER :: green_fft => null()
128 TYPE(ps_wavelet_type),
POINTER :: wavelet => null()
129 TYPE(pw_poisson_parameter_type) :: parameters = pw_poisson_parameter_type()
130 REAL(kind=
dp),
DIMENSION(3, 3) :: cell_hmat = 0.0_dp
131 TYPE(pw_pool_p_type),
DIMENSION(:),
POINTER :: pw_pools => null()
132 TYPE(pw_grid_type),
POINTER :: mt_super_ref_pw_grid => null()
133 TYPE(ps_implicit_type),
POINTER :: implicit_env => null()
134 TYPE(pw_grid_type),
POINTER :: dct_pw_grid => null()
135 TYPE(realspace_grid_type),
POINTER :: diel_rs_grid => null()
137 PROCEDURE,
PUBLIC, non_overridable :: create => pw_poisson_create
138 PROCEDURE,
PUBLIC, non_overridable :: release => pw_poisson_release
139 END TYPE pw_poisson_type
147 INTEGER :: special_dimension = 0
148 REAL(kind=
dp) :: radius = 0.0_dp
149 REAL(kind=
dp) :: mt_alpha = 1.0_dp
150 REAL(kind=
dp) :: mt_rel_cutoff = 1.0_dp
151 REAL(kind=
dp) :: slab_size = 0.0_dp
152 REAL(kind=
dp) :: alpha = 0.0_dp
153 LOGICAL :: p3m = .false.
154 INTEGER :: p3m_order = 0
155 REAL(kind=
dp) :: p3m_alpha = 0.0_dp
156 REAL(kind=
dp),
DIMENSION(:, :),
ALLOCATABLE :: p3m_coeff
157 REAL(kind=
dp),
DIMENSION(:, :),
ALLOCATABLE :: p3m_bm2
158 LOGICAL :: sr_screening = .false.
159 REAL(kind=
dp) :: sr_alpha = 1.0_dp
160 REAL(kind=
dp) :: sr_rc = 0.0_dp
161 TYPE(pw_c1d_gs_type) :: influence_fn = pw_c1d_gs_type()
162 TYPE(pw_c1d_gs_type),
POINTER :: dct_influence_fn => null()
163 TYPE(pw_c1d_gs_type),
POINTER :: screen_fn => null()
164 TYPE(pw_r1d_gs_type),
POINTER :: p3m_charge => null()
165 END TYPE greens_fn_type
181 mt_super_ref_pw_grid, dct_pw_grid)
182 TYPE(greens_fn_type),
INTENT(OUT) :: green
183 TYPE(pw_poisson_parameter_type),
INTENT(IN) :: poisson_params
184 REAL(kind=
dp),
DIMENSION(3, 3),
INTENT(IN) :: cell_hmat
185 TYPE(pw_pool_type),
POINTER :: pw_pool
186 TYPE(pw_grid_type),
POINTER :: mt_super_ref_pw_grid, dct_pw_grid
188 INTEGER :: dim, i, ig, iz, n, nz
189 REAL(kind=
dp) :: g2, g3d, gg, gxy, gz, j0g, j1g, k0g, &
190 k1g, rlength, zlength
191 REAL(kind=
dp),
DIMENSION(3) :: abc
192 TYPE(pw_c1d_gs_type),
POINTER :: dct_gf
193 TYPE(pw_grid_type),
POINTER :: dct_grid
194 TYPE(pw_pool_type),
POINTER :: pw_pool_xpndd
198 abc(i) = cell_hmat(i, i)
200 dim = count(poisson_params%periodic == 1)
202 SELECT CASE (poisson_params%solver)
206 cpabort(
"Illegal combination of periodicity and Poisson solver periodic3d")
211 cpabort(
"Illegal combination of periodicity and Poisson solver mulipole0d")
217 green%radius = 0.5_dp*minval(abc)
220 green%special_dimension = maxloc(poisson_params%periodic(1:3), 1)
221 green%radius = maxval(abc(1:3))
223 IF (i == green%special_dimension) cycle
224 green%radius = min(green%radius, 0.5_dp*abc(i))
228 i = minloc(poisson_params%periodic, 1)
229 green%special_dimension = i
230 green%slab_size = abc(i)
237 green%MT_rel_cutoff = poisson_params%mt_rel_cutoff
238 green%MT_alpha = poisson_params%mt_alpha
239 green%MT_alpha = green%MT_alpha/minval(abc)
243 green%radius = 0.5_dp*minval(abc)
246 green%special_dimension = maxloc(poisson_params%periodic(1:3), 1)
247 green%radius = maxval(abc(1:3))
249 IF (i == green%special_dimension) cycle
250 green%radius = min(green%radius, 0.5_dp*abc(i))
254 i = minloc(poisson_params%periodic, 1)
255 green%special_dimension = i
256 green%slab_size = abc(i)
258 cpabort(
"Illegal combination of periodicity and Poisson solver (MT)")
265 cpabort(
"An unknown Poisson solver was specified")
269 SELECT CASE (green%method)
271 CALL pw_pool%create_pw(green%influence_fn)
275 green%p3m_order = poisson_params%ewald_o_spline
276 green%p3m_alpha = poisson_params%ewald_alpha
278 ALLOCATE (green%p3m_coeff(-(n - 1):n - 1, 0:n - 1))
279 CALL spme_coeff_calculate(n, green%p3m_coeff)
280 NULLIFY (green%p3m_charge)
281 ALLOCATE (green%p3m_charge)
282 CALL pw_pool%create_pw(green%p3m_charge)
283 CALL influence_factor(green)
284 CALL calc_p3m_charge(green)
289 SELECT CASE (green%method)
292 alpha=green%MT_alpha, &
293 special_dimension=green%special_dimension, slab_size=green%slab_size, &
294 super_ref_pw_grid=mt_super_ref_pw_grid)
296 IF ((poisson_params%ps_implicit_params%boundary_condition .EQ.
mixed_bc) .OR. &
297 (poisson_params%ps_implicit_params%boundary_condition .EQ.
neumann_bc))
THEN
299 NULLIFY (green%dct_influence_fn)
300 ALLOCATE (green%dct_influence_fn)
301 CALL pw_pool_xpndd%create_pw(green%dct_influence_fn)
311 associate(gf => green%influence_fn, grid => green%influence_fn%pw_grid)
312 SELECT CASE (green%method)
315 DO ig = grid%first_gne0, grid%ngpts_cut_local
319 IF (grid%have_g0) gf%array(1) = 0.0_dp
323 iz = green%special_dimension
324 zlength = green%slab_size
325 DO ig = grid%first_gne0, grid%ngpts_cut_local
326 nz = grid%g_hat(iz, ig)
329 gxy = max(0.0_dp, g2 - grid%g(iz, ig)*grid%g(iz, ig))
330 gg = 0.5_dp*sqrt(gxy)
331 gf%array(ig) = g3d*(1.0_dp - (-1.0_dp)**nz*exp(-gg*zlength))
333 IF (grid%have_g0) gf%array(1) = 0.0_dp
338 iz = green%special_dimension
340 rlength = green%radius
341 DO ig = grid%first_gne0, grid%ngpts_cut_local
344 gxy = sqrt(max(0.0_dp, g2 - grid%g(iz, ig)*grid%g(iz, ig)))
345 gz = abs(grid%g(iz, ig))
346 j0g = bessel_j0(rlength*gxy)
347 j1g = bessel_j1(rlength*gxy)
355 gf%array(ig) = g3d*(1.0_dp + rlength* &
356 (gxy*j1g*k0g - gz*j0g*k1g))
358 IF (grid%have_g0) gf%array(1) = 0.0_dp
362 rlength = green%radius
363 DO ig = grid%first_gne0, grid%ngpts_cut_local
367 gf%array(ig) = g3d*(1.0_dp - cos(rlength*gg))
370 gf%array(1) = 0.5_dp*
fourpi*rlength*rlength
374 DO ig = grid%first_gne0, grid%ngpts_cut_local
377 gf%array(ig) = g3d + green%screen_fn%array(ig)
380 gf%array(1) = green%screen_fn%array(1)
384 DO ig = grid%first_gne0, grid%ngpts_cut_local
388 IF (grid%have_g0) gf%array(1) = 0.0_dp
390 IF (
ASSOCIATED(green%dct_influence_fn))
THEN
391 dct_gf => green%dct_influence_fn
392 dct_grid => green%dct_influence_fn%pw_grid
394 DO ig = dct_grid%first_gne0, dct_grid%ngpts_cut_local
395 g2 = dct_grid%gsq(ig)
396 dct_gf%array(ig) =
fourpi/g2
398 IF (dct_grid%have_g0) dct_gf%array(1) = 0.0_dp
418 TYPE(greens_fn_type),
INTENT(INOUT) :: gftype
419 TYPE(pw_pool_type),
OPTIONAL,
POINTER :: pw_pool
421 LOGICAL :: can_give_back
423 can_give_back =
PRESENT(pw_pool)
424 IF (can_give_back) can_give_back =
ASSOCIATED(pw_pool)
425 IF (can_give_back)
THEN
426 CALL pw_pool%give_back_pw(gftype%influence_fn)
427 IF (
ASSOCIATED(gftype%dct_influence_fn))
THEN
428 CALL pw_pool%give_back_pw(gftype%dct_influence_fn)
429 DEALLOCATE (gftype%dct_influence_fn)
431 IF (
ASSOCIATED(gftype%screen_fn))
THEN
432 CALL pw_pool%give_back_pw(gftype%screen_fn)
433 DEALLOCATE (gftype%screen_fn)
435 IF (
ASSOCIATED(gftype%p3m_charge))
THEN
436 CALL pw_pool%give_back_pw(gftype%p3m_charge)
437 DEALLOCATE (gftype%p3m_charge)
440 CALL gftype%influence_fn%release()
441 IF (
ASSOCIATED(gftype%dct_influence_fn))
THEN
442 CALL gftype%dct_influence_fn%release()
443 DEALLOCATE (gftype%dct_influence_fn)
445 IF (
ASSOCIATED(gftype%screen_fn))
THEN
446 CALL gftype%screen_fn%release()
447 DEALLOCATE (gftype%screen_fn)
449 IF (
ASSOCIATED(gftype%p3m_charge))
THEN
450 CALL gftype%p3m_charge%release()
451 DEALLOCATE (gftype%p3m_charge)
454 IF (
ALLOCATED(gftype%p3m_bm2)) &
455 DEALLOCATE (gftype%p3m_bm2)
456 IF (
ALLOCATED(gftype%p3m_coeff)) &
457 DEALLOCATE (gftype%p3m_coeff)
468 SUBROUTINE influence_factor(gftype)
469 TYPE(greens_fn_type),
INTENT(INOUT) :: gftype
471 COMPLEX(KIND=dp) :: b_m, exp_m, sum_m
472 INTEGER :: dim, j, k, l, n, pt
473 INTEGER,
DIMENSION(3) :: npts
474 INTEGER,
DIMENSION(:),
POINTER :: lb, ub
475 REAL(kind=dp) :: l_arg, prod_arg, val
476 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:) :: m_assign
482 lb => gftype%influence_fn%pw_grid%bounds(1, :)
483 ub => gftype%influence_fn%pw_grid%bounds(2, :)
484 IF (
ALLOCATED(gftype%p3m_bm2))
THEN
485 IF (lbound(gftype%p3m_bm2, 2) /= minval(lb(:)) .OR. &
486 ubound(gftype%p3m_bm2, 2) /= maxval(ub(:)))
THEN
487 DEALLOCATE (gftype%p3m_bm2)
490 IF (.NOT.
ALLOCATED(gftype%p3m_bm2))
THEN
491 ALLOCATE (gftype%p3m_bm2(3, minval(lb(:)):maxval(ub(:))))
494 ALLOCATE (m_assign(0:n - 2))
500 prod_arg = gftype%p3m_coeff(j, l)*l_arg
501 m_assign(k) = m_assign(k) + prod_arg
507 npts(:) = ub(:) - lb(:) + 1
509 DO pt = lb(dim), ub(dim)
510 val = twopi*(real(pt, kind=dp)/real(npts(dim), kind=dp))
511 exp_m = cmplx(cos(val), -sin(val), kind=dp)
512 sum_m = cmplx(0.0_dp, 0.0_dp, kind=dp)
514 sum_m = sum_m + m_assign(k)*exp_m**k
516 b_m = exp_m**(n - 1)/sum_m
517 gftype%p3m_bm2(dim, pt) = sqrt(real(b_m*conjg(b_m), kind=dp))
521 DEALLOCATE (m_assign)
522 END SUBROUTINE influence_factor
528 PURE SUBROUTINE calc_p3m_charge(gf)
530 TYPE(greens_fn_type),
INTENT(INOUT) :: gf
532 INTEGER :: ig, l, m, n
533 REAL(kind=dp) :: arg, novol
537 associate(grid => gf%influence_fn%pw_grid, bm2 => gf%p3m_bm2)
538 arg = 1.0_dp/(8.0_dp*gf%p3m_alpha**2)
539 novol = real(grid%ngpts, kind=dp)/grid%vol
540 DO ig = 1, grid%ngpts_cut_local
541 l = grid%g_hat(1, ig)
542 m = grid%g_hat(2, ig)
543 n = grid%g_hat(3, ig)
544 gf%p3m_charge%array(ig) = novol*exp(-arg*grid%gsq(ig))* &
545 bm2(1, l)*bm2(2, m)*bm2(3, n)
549 END SUBROUTINE calc_p3m_charge
561 SUBROUTINE pw_poisson_create(poisson_env)
563 CLASS(pw_poisson_type),
INTENT(INOUT) :: poisson_env
566 CALL poisson_env%release()
568 END SUBROUTINE pw_poisson_create
577 SUBROUTINE pw_poisson_release(poisson_env)
579 CLASS(pw_poisson_type),
INTENT(INOUT) :: poisson_env
581 IF (
ASSOCIATED(poisson_env%pw_pools))
THEN
582 CALL pw_pools_dealloc(poisson_env%pw_pools)
585 IF (
ASSOCIATED(poisson_env%green_fft))
THEN
587 DEALLOCATE (poisson_env%green_fft)
589 CALL pw_grid_release(poisson_env%mt_super_ref_pw_grid)
590 CALL ps_wavelet_release(poisson_env%wavelet)
591 CALL ps_implicit_release(poisson_env%implicit_env, &
592 poisson_env%parameters%ps_implicit_params)
593 CALL pw_grid_release(poisson_env%dct_pw_grid)
594 IF (
ASSOCIATED(poisson_env%diel_rs_grid))
THEN
595 CALL rs_grid_release(poisson_env%diel_rs_grid)
596 DEALLOCATE (poisson_env%diel_rs_grid)
599 END SUBROUTINE pw_poisson_release
609 SUBROUTINE spme_coeff_calculate(n, coeff)
611 INTEGER,
INTENT(IN) :: n
612 REAL(kind=dp),
DIMENSION(-(n-1):n-1, 0:n-1), &
615 INTEGER :: i, j, l, m
617 REAL(kind=dp),
DIMENSION(n, -n:n, 0:n-1) :: a
626 b = (a(m, j - 1, l) + &
627 REAL((-1)**l, kind=dp)*a(m, j + 1, l))/ &
628 REAL((l + 1)*2**(l + 1), kind=dp)
629 a(i, j, 0) = a(i, j, 0) + b
632 a(i, j, l + 1) = (a(m, j + 1, l) - &
633 a(m, j - 1, l))/real(l + 1, kind=dp)
640 DO j = -(n - 1), n - 1, 2
641 coeff(j, i) = a(n, j, i)
645 END SUBROUTINE spme_coeff_calculate
Calculates Bessel functions.
elemental real(kind=dp) function, public bessk1(x)
...
elemental real(kind=dp) function, public bessk0(x)
...
dielectric constant data type
Dirichlet boundary condition data types.
Defines the basic variable types.
integer, parameter, public dp
Definition of mathematical constants and functions.
real(kind=dp), parameter, public fourpi
real(kind=dp), parameter, public twopi
integer, parameter, public mt0d
integer, parameter, public mt2d
subroutine, public mtin_create_screen_fn(screen_function, pw_pool, method, alpha, special_dimension, slab_size, super_ref_pw_grid)
Initialize the Martyna && Tuckerman Poisson Solver.
integer, parameter, public mt1d
Types containing essential information for running implicit (iterative) Poisson solver.
integer, parameter, public neumann_bc
subroutine, public ps_implicit_release(ps_implicit_env, ps_implicit_params, pw_pool)
Deallocates ps_implicit.
integer, parameter, public mixed_bc
Definition and initialisation of the ps_wavelet data type.
integer, parameter, public wavelet0d
subroutine, public ps_wavelet_release(wavelet)
...
This module defines the grid data type and some basic operations on it.
subroutine, public pw_grid_release(pw_grid)
releases the given pw grid
functions related to the poisson solver on regular grids
integer, parameter, public pw_poisson_wavelet
subroutine, public pw_green_create(green, poisson_params, cell_hmat, pw_pool, mt_super_ref_pw_grid, dct_pw_grid)
Allocates and sets up the green functions for the fft based poisson solvers.
integer, parameter, public pw_poisson_periodic
integer, parameter, public multipole1d
subroutine, public pw_green_release(gftype, pw_pool)
destroys the type (deallocates data)
integer, parameter, public pw_poisson_none
integer, parameter, public pw_poisson_mt
integer, parameter, public hockney0d
integer, parameter, public multipole2d
integer, parameter, public pw_poisson_hockney
integer, parameter, public hockney2d
integer, parameter, public do_ewald_pme
integer, parameter, public do_ewald_ewald
integer, parameter, public analytic2d
integer, parameter, public analytic1d
integer, parameter, public periodic3d
integer, parameter, public hockney1d
integer, parameter, public ps_implicit
integer, parameter, public analytic0d
integer, parameter, public multipole0d
integer, parameter, public pw_poisson_implicit
integer, parameter, public do_ewald_none
integer, parameter, public pw_poisson_analytic
integer, parameter, public do_ewald_spme
integer, parameter, public pw_poisson_multipole
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
subroutine, public pw_pools_dealloc(pools)
deallocates the given pools (releasing each of the underlying pools)
subroutine, public pw_pool_release(pool)
releases the given pool (see cp2k/doc/ReferenceCounting.html)
subroutine, public pw_pool_create(pool, pw_grid, max_cache)
creates a pool for pw
subroutine, public rs_grid_release(rs_grid)
releases the given rs grid (see doc/ReferenceCounting.html)