Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief functions related to the poisson solver on regular grids
10 : !> \par History
11 : !> greens_fn: JGH (9-Mar-2001) : include influence_fn into
12 : !> greens_fn_type add cell volume
13 : !> as indicator for updates
14 : !> greens_fn: JGH (30-Mar-2001) : Added B-spline routines
15 : !> pws : JGH (13-Mar-2001) : new pw_poisson_solver, delete
16 : !> pw_greens_fn
17 : !> 12.2004 condensed from pws, greens_fn and green_fns, by apsi and JGH,
18 : !> made thread safe, new input [fawzi]
19 : !> 14-Mar-2006 : added short range screening function for SE codes
20 : !> \author fawzi
21 : ! **************************************************************************************************
22 : MODULE pw_poisson_types
23 : USE bessel_lib, ONLY: bessk0,&
24 : bessk1
25 : USE dielectric_types, ONLY: dielectric_parameters
26 : USE dirichlet_bc_types, ONLY: dirichlet_bc_parameters
27 : USE kinds, ONLY: dp
28 : USE mathconstants, ONLY: fourpi,&
29 : twopi,&
30 : z_zero
31 : USE mt_util, ONLY: MT0D,&
32 : MT1D,&
33 : MT2D,&
34 : MTin_create_screen_fn
35 : USE ps_implicit_types, ONLY: MIXED_BC,&
36 : NEUMANN_BC,&
37 : ps_implicit_parameters,&
38 : ps_implicit_release,&
39 : ps_implicit_type
40 : USE ps_wavelet_types, ONLY: WAVELET0D,&
41 : ps_wavelet_release,&
42 : ps_wavelet_type
43 : USE pw_grid_types, ONLY: pw_grid_type
44 : USE pw_grids, ONLY: pw_grid_release
45 : USE pw_pool_types, ONLY: pw_pool_create,&
46 : pw_pool_p_type,&
47 : pw_pool_release,&
48 : pw_pool_type,&
49 : pw_pools_dealloc
50 : USE pw_types, ONLY: pw_c1d_gs_type,&
51 : pw_r1d_gs_type
52 : USE realspace_grid_types, ONLY: realspace_grid_type,&
53 : rs_grid_release
54 : #include "../base/base_uses.f90"
55 :
56 : IMPLICIT NONE
57 : PRIVATE
58 :
59 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
60 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pw_poisson_types'
61 :
62 : PUBLIC :: pw_poisson_type
63 : PUBLIC :: greens_fn_type, pw_green_create, &
64 : pw_green_release
65 : PUBLIC :: pw_poisson_parameter_type
66 :
67 : INTEGER, PARAMETER, PUBLIC :: pw_poisson_none = 0, &
68 : pw_poisson_periodic = 1, &
69 : pw_poisson_analytic = 2, &
70 : pw_poisson_mt = 3, &
71 : pw_poisson_hockney = 5, &
72 : pw_poisson_multipole = 4, &
73 : pw_poisson_wavelet = 6, &
74 : pw_poisson_implicit = 7
75 : ! EWALD methods
76 : INTEGER, PARAMETER, PUBLIC :: do_ewald_none = 1, &
77 : do_ewald_ewald = 2, &
78 : do_ewald_pme = 3, &
79 : do_ewald_spme = 4
80 :
81 : INTEGER, PARAMETER, PUBLIC :: PERIODIC3D = 1000, &
82 : ANALYTIC2D = 1001, &
83 : ANALYTIC1D = 1002, &
84 : ANALYTIC0D = 1003, &
85 : HOCKNEY2D = 1201, &
86 : HOCKNEY1D = 1202, &
87 : HOCKNEY0D = 1203, &
88 : MULTIPOLE2D = 1301, &
89 : MULTIPOLE1D = 1302, &
90 : MULTIPOLE0D = 1303, &
91 : PS_IMPLICIT = 1400
92 :
93 : ! **************************************************************************************************
94 : !> \brief parameters for the poisson solver independet of input_section
95 : !> \author Ole Schuett
96 : ! **************************************************************************************************
97 : TYPE pw_poisson_parameter_type
98 : INTEGER :: solver = pw_poisson_none
99 :
100 : INTEGER, DIMENSION(3) :: periodic = 0
101 : INTEGER :: ewald_type = do_ewald_none
102 : INTEGER :: ewald_o_spline = 0
103 : REAL(KIND=dp) :: ewald_alpha = 0.0_dp
104 :
105 : REAL(KIND=dp) :: mt_rel_cutoff = 0.0_dp
106 : REAL(KIND=dp) :: mt_alpha = 0.0_dp
107 :
108 : INTEGER :: wavelet_scf_type = 0
109 : INTEGER :: wavelet_method = WAVELET0D
110 : INTEGER :: wavelet_special_dimension = 0
111 : CHARACTER(LEN=1) :: wavelet_geocode = "S"
112 :
113 : LOGICAL :: has_dielectric = .FALSE.
114 : TYPE(dielectric_parameters) :: dielectric_params = dielectric_parameters()
115 : TYPE(ps_implicit_parameters) :: ps_implicit_params = ps_implicit_parameters()
116 : TYPE(dirichlet_bc_parameters) :: dbc_params = dirichlet_bc_parameters()
117 : END TYPE pw_poisson_parameter_type
118 :
119 : ! **************************************************************************************************
120 : !> \brief environment for the poisson solver
121 : !> \author fawzi
122 : ! **************************************************************************************************
123 : TYPE pw_poisson_type
124 : INTEGER :: pw_level = 0
125 : INTEGER :: method = pw_poisson_none
126 : INTEGER :: used_grid = 0
127 : LOGICAL :: rebuild = .TRUE.
128 : TYPE(greens_fn_type), POINTER :: green_fft => NULL()
129 : TYPE(ps_wavelet_type), POINTER :: wavelet => NULL()
130 : TYPE(pw_poisson_parameter_type) :: parameters = pw_poisson_parameter_type()
131 : REAL(KIND=dp), DIMENSION(3, 3) :: cell_hmat = 0.0_dp
132 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools => NULL()
133 : TYPE(pw_grid_type), POINTER :: mt_super_ref_pw_grid => NULL()
134 : TYPE(ps_implicit_type), POINTER :: implicit_env => NULL()
135 : TYPE(pw_grid_type), POINTER :: dct_pw_grid => NULL()
136 : TYPE(realspace_grid_type), POINTER :: diel_rs_grid => NULL()
137 : CONTAINS
138 : PROCEDURE, PUBLIC, NON_OVERRIDABLE :: create => pw_poisson_create
139 : PROCEDURE, PUBLIC, NON_OVERRIDABLE :: release => pw_poisson_release
140 : END TYPE pw_poisson_type
141 :
142 : ! **************************************************************************************************
143 : !> \brief contains all the informations needed by the fft based poisson solvers
144 : !> \author JGH,Teo,fawzi
145 : ! **************************************************************************************************
146 : TYPE greens_fn_type
147 : INTEGER :: method = PERIODIC3D
148 : INTEGER :: special_dimension = 0
149 : REAL(KIND=dp) :: radius = 0.0_dp
150 : REAL(KIND=dp) :: MT_alpha = 1.0_dp
151 : REAL(KIND=dp) :: MT_rel_cutoff = 1.0_dp
152 : REAL(KIND=dp) :: slab_size = 0.0_dp
153 : REAL(KIND=dp) :: alpha = 0.0_dp
154 : LOGICAL :: p3m = .FALSE.
155 : INTEGER :: p3m_order = 0
156 : REAL(KIND=dp) :: p3m_alpha = 0.0_dp
157 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: p3m_coeff
158 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: p3m_bm2
159 : LOGICAL :: sr_screening = .FALSE.
160 : REAL(KIND=dp) :: sr_alpha = 1.0_dp
161 : REAL(KIND=dp) :: sr_rc = 0.0_dp
162 : TYPE(pw_c1d_gs_type) :: influence_fn = pw_c1d_gs_type()
163 : TYPE(pw_c1d_gs_type), POINTER :: dct_influence_fn => NULL()
164 : TYPE(pw_c1d_gs_type), POINTER :: screen_fn => NULL()
165 : TYPE(pw_r1d_gs_type), POINTER :: p3m_charge => NULL()
166 : END TYPE greens_fn_type
167 :
168 : CONTAINS
169 :
170 : ! **************************************************************************************************
171 : !> \brief Allocates and sets up the green functions for the fft based poisson
172 : !> solvers
173 : !> \param green ...
174 : !> \param poisson_params ...
175 : !> \param cell_hmat ...
176 : !> \param pw_pool ...
177 : !> \param mt_super_ref_pw_grid ...
178 : !> \param dct_pw_grid ...
179 : !> \author Fawzi, based on previous functions by JGH and Teo
180 : ! **************************************************************************************************
181 15794 : SUBROUTINE pw_green_create(green, poisson_params, cell_hmat, pw_pool, &
182 : mt_super_ref_pw_grid, dct_pw_grid)
183 : TYPE(greens_fn_type), INTENT(OUT) :: green
184 : TYPE(pw_poisson_parameter_type), INTENT(IN) :: poisson_params
185 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: cell_hmat
186 : TYPE(pw_pool_type), POINTER :: pw_pool
187 : TYPE(pw_grid_type), POINTER :: mt_super_ref_pw_grid, dct_pw_grid
188 :
189 : INTEGER :: dim, i, ig, iz, n, nz
190 : REAL(KIND=dp) :: g2, g3d, gg, gxy, gz, j0g, j1g, k0g, &
191 : k1g, rlength, zlength
192 : REAL(KIND=dp), DIMENSION(3) :: abc
193 : TYPE(pw_c1d_gs_type), POINTER :: dct_gf
194 : TYPE(pw_grid_type), POINTER :: dct_grid
195 : TYPE(pw_pool_type), POINTER :: pw_pool_xpndd
196 :
197 : !CPASSERT(cell%orthorhombic)
198 63176 : DO i = 1, 3
199 63176 : abc(i) = cell_hmat(i, i)
200 : END DO
201 63176 : dim = COUNT(poisson_params%periodic == 1)
202 :
203 30116 : SELECT CASE (poisson_params%solver)
204 : CASE (pw_poisson_periodic)
205 : green%method = PERIODIC3D
206 14322 : IF (dim /= 3) THEN
207 0 : CPABORT("Illegal combination of periodicity and Poisson solver periodic3d")
208 : END IF
209 : CASE (pw_poisson_multipole)
210 20 : green%method = MULTIPOLE0D
211 20 : IF (dim /= 0) THEN
212 0 : CPABORT("Illegal combination of periodicity and Poisson solver mulipole0d")
213 : END IF
214 : CASE (pw_poisson_analytic)
215 1382 : SELECT CASE (dim)
216 : CASE (0)
217 170 : green%method = ANALYTIC0D
218 850 : green%radius = 0.5_dp*MINVAL(abc)
219 : CASE (1)
220 6 : green%method = ANALYTIC1D
221 24 : green%special_dimension = MAXLOC(poisson_params%periodic(1:3), 1)
222 30 : green%radius = MAXVAL(abc(1:3))
223 24 : DO i = 1, 3
224 18 : IF (i == green%special_dimension) CYCLE
225 24 : green%radius = MIN(green%radius, 0.5_dp*abc(i))
226 : END DO
227 : CASE (2)
228 8 : green%method = ANALYTIC2D
229 32 : i = MINLOC(poisson_params%periodic, 1)
230 8 : green%special_dimension = i
231 8 : green%slab_size = abc(i)
232 : CASE (3)
233 0 : green%method = PERIODIC3D
234 : CASE DEFAULT
235 186 : CPABORT("")
236 : END SELECT
237 : CASE (pw_poisson_mt)
238 1212 : green%MT_rel_cutoff = poisson_params%mt_rel_cutoff
239 1212 : green%MT_alpha = poisson_params%mt_alpha
240 6060 : green%MT_alpha = green%MT_alpha/MINVAL(abc)
241 1264 : SELECT CASE (dim)
242 : CASE (0)
243 1210 : green%method = MT0D
244 6050 : green%radius = 0.5_dp*MINVAL(abc)
245 : CASE (1)
246 0 : green%method = MT1D
247 0 : green%special_dimension = MAXLOC(poisson_params%periodic(1:3), 1)
248 0 : green%radius = MAXVAL(abc(1:3))
249 0 : DO i = 1, 3
250 0 : IF (i == green%special_dimension) CYCLE
251 0 : green%radius = MIN(green%radius, 0.5_dp*abc(i))
252 : END DO
253 : CASE (2)
254 2 : green%method = MT2D
255 8 : i = MINLOC(poisson_params%periodic, 1)
256 2 : green%special_dimension = i
257 2 : green%slab_size = abc(i)
258 : CASE (3)
259 0 : CPABORT("Illegal combination of periodicity and Poisson solver (MT)")
260 : CASE DEFAULT
261 1212 : CPABORT("")
262 : END SELECT
263 : CASE (pw_poisson_implicit)
264 54 : green%method = PS_IMPLICIT
265 : CASE DEFAULT
266 15794 : CPABORT("An unknown Poisson solver was specified")
267 : END SELECT
268 :
269 : ! allocate influence function,...
270 31588 : SELECT CASE (green%method)
271 : CASE (PERIODIC3D, ANALYTIC2D, ANALYTIC1D, ANALYTIC0D, MT2D, MT1D, MT0D, MULTIPOLE0D, PS_IMPLICIT)
272 15794 : CALL pw_pool%create_pw(green%influence_fn)
273 :
274 15794 : IF (poisson_params%ewald_type == do_ewald_spme) THEN
275 9870 : green%p3m = .TRUE.
276 9870 : green%p3m_order = poisson_params%ewald_o_spline
277 9870 : green%p3m_alpha = poisson_params%ewald_alpha
278 9870 : n = green%p3m_order
279 39480 : ALLOCATE (green%p3m_coeff(-(n - 1):n - 1, 0:n - 1))
280 9870 : CALL spme_coeff_calculate(n, green%p3m_coeff)
281 : NULLIFY (green%p3m_charge)
282 9870 : ALLOCATE (green%p3m_charge)
283 9870 : CALL pw_pool%create_pw(green%p3m_charge)
284 9870 : CALL influence_factor(green)
285 9870 : CALL calc_p3m_charge(green)
286 : ELSE
287 5924 : green%p3m = .FALSE.
288 : END IF
289 : !
290 17006 : SELECT CASE (green%method)
291 : CASE (MT0D, MT1D, MT2D)
292 : CALL MTin_create_screen_fn(green%screen_fn, pw_pool=pw_pool, method=green%method, &
293 : alpha=green%MT_alpha, &
294 : special_dimension=green%special_dimension, slab_size=green%slab_size, &
295 1212 : super_ref_pw_grid=mt_super_ref_pw_grid)
296 : CASE (PS_IMPLICIT)
297 54 : IF ((poisson_params%ps_implicit_params%boundary_condition == MIXED_BC) .OR. &
298 15794 : (poisson_params%ps_implicit_params%boundary_condition == NEUMANN_BC)) THEN
299 22 : CALL pw_pool_create(pw_pool_xpndd, pw_grid=dct_pw_grid)
300 : NULLIFY (green%dct_influence_fn)
301 22 : ALLOCATE (green%dct_influence_fn)
302 22 : CALL pw_pool_xpndd%create_pw(green%dct_influence_fn)
303 22 : CALL pw_pool_release(pw_pool_xpndd)
304 : END IF
305 : END SELECT
306 :
307 : CASE DEFAULT
308 15794 : CPABORT("")
309 : END SELECT
310 :
311 : ! initialize influence function
312 : ASSOCIATE (gf => green%influence_fn, grid => green%influence_fn%pw_grid)
313 30138 : SELECT CASE (green%method)
314 : CASE (PERIODIC3D, MULTIPOLE0D)
315 :
316 345671671 : DO ig = grid%first_gne0, grid%ngpts_cut_local
317 345657327 : g2 = grid%gsq(ig)
318 345671671 : gf%array(ig) = fourpi/g2
319 : END DO
320 14344 : IF (grid%have_g0) gf%array(1) = 0.0_dp
321 :
322 : CASE (ANALYTIC2D)
323 :
324 8 : iz = green%special_dimension ! iz is the direction with NO PBC
325 8 : zlength = green%slab_size ! zlength is the thickness of the cell
326 634372 : DO ig = grid%first_gne0, grid%ngpts_cut_local
327 634364 : nz = grid%g_hat(iz, ig)
328 634364 : g2 = grid%gsq(ig)
329 634364 : g3d = fourpi/g2
330 634364 : gxy = MAX(0.0_dp, g2 - grid%g(iz, ig)*grid%g(iz, ig))
331 634364 : gg = 0.5_dp*SQRT(gxy)
332 634372 : gf%array(ig) = g3d*(1.0_dp - (-1.0_dp)**nz*EXP(-gg*zlength))
333 : END DO
334 8 : IF (grid%have_g0) gf%array(1) = 0.0_dp
335 :
336 : CASE (ANALYTIC1D)
337 : ! see 'ab initio molecular dynamics' table 3.1
338 : ! iz is the direction of the PBC ( can be 1,2,3 -> x,y,z )
339 6 : iz = green%special_dimension
340 : ! rlength is the radius of the tube
341 6 : rlength = green%radius
342 192003 : DO ig = grid%first_gne0, grid%ngpts_cut_local
343 191997 : g2 = grid%gsq(ig)
344 191997 : g3d = fourpi/g2
345 191997 : gxy = SQRT(MAX(0.0_dp, g2 - grid%g(iz, ig)*grid%g(iz, ig)))
346 191997 : gz = ABS(grid%g(iz, ig))
347 191997 : j0g = BESSEL_J0(rlength*gxy)
348 191997 : j1g = BESSEL_J1(rlength*gxy)
349 191997 : IF (gz > 0) THEN
350 187200 : k0g = bessk0(rlength*gz)
351 187200 : k1g = bessk1(rlength*gz)
352 : ELSE
353 : k0g = 0
354 : k1g = 0
355 : END IF
356 : gf%array(ig) = g3d*(1.0_dp + rlength* &
357 192003 : (gxy*j1g*k0g - gz*j0g*k1g))
358 : END DO
359 6 : IF (grid%have_g0) gf%array(1) = 0.0_dp
360 :
361 : CASE (ANALYTIC0D)
362 :
363 170 : rlength = green%radius ! rlength is the radius of the sphere
364 39004540 : DO ig = grid%first_gne0, grid%ngpts_cut_local
365 39004370 : g2 = grid%gsq(ig)
366 39004370 : gg = SQRT(g2)
367 39004370 : g3d = fourpi/g2
368 39004540 : gf%array(ig) = g3d*(1.0_dp - COS(rlength*gg))
369 : END DO
370 170 : IF (grid%have_g0) &
371 85 : gf%array(1) = 0.5_dp*fourpi*rlength*rlength
372 :
373 : CASE (MT2D, MT1D, MT0D)
374 :
375 81587146 : DO ig = grid%first_gne0, grid%ngpts_cut_local
376 81585934 : g2 = grid%gsq(ig)
377 81585934 : g3d = fourpi/g2
378 81587146 : gf%array(ig) = g3d + green%screen_fn%array(ig)
379 : END DO
380 1212 : IF (grid%have_g0) &
381 612 : gf%array(1) = green%screen_fn%array(1)
382 :
383 : CASE (PS_IMPLICIT)
384 :
385 6554816 : DO ig = grid%first_gne0, grid%ngpts_cut_local
386 6554762 : g2 = grid%gsq(ig)
387 6554816 : gf%array(ig) = fourpi/g2
388 : END DO
389 54 : IF (grid%have_g0) gf%array(1) = 0.0_dp
390 :
391 54 : IF (ASSOCIATED(green%dct_influence_fn)) THEN
392 22 : dct_gf => green%dct_influence_fn
393 22 : dct_grid => green%dct_influence_fn%pw_grid
394 :
395 9851771 : DO ig = dct_grid%first_gne0, dct_grid%ngpts_cut_local
396 9851749 : g2 = dct_grid%gsq(ig)
397 9851771 : dct_gf%array(ig) = fourpi/g2
398 : END DO
399 22 : IF (dct_grid%have_g0) dct_gf%array(1) = 0.0_dp
400 : END IF
401 :
402 : CASE DEFAULT
403 15794 : CPABORT("")
404 : END SELECT
405 : END ASSOCIATE
406 :
407 15794 : END SUBROUTINE pw_green_create
408 :
409 : ! **************************************************************************************************
410 : !> \brief destroys the type (deallocates data)
411 : !> \param gftype ...
412 : !> \param pw_pool ...
413 : !> \par History
414 : !> none
415 : !> \author Joost VandeVondele
416 : !> Teodoro Laino
417 : ! **************************************************************************************************
418 15794 : SUBROUTINE pw_green_release(gftype, pw_pool)
419 : TYPE(greens_fn_type), INTENT(INOUT) :: gftype
420 : TYPE(pw_pool_type), OPTIONAL, POINTER :: pw_pool
421 :
422 : LOGICAL :: can_give_back
423 :
424 15794 : can_give_back = PRESENT(pw_pool)
425 15794 : IF (can_give_back) can_give_back = ASSOCIATED(pw_pool)
426 15794 : IF (can_give_back) THEN
427 8054 : CALL pw_pool%give_back_pw(gftype%influence_fn)
428 8054 : IF (ASSOCIATED(gftype%dct_influence_fn)) THEN
429 0 : CALL pw_pool%give_back_pw(gftype%dct_influence_fn)
430 0 : DEALLOCATE (gftype%dct_influence_fn)
431 : END IF
432 8054 : IF (ASSOCIATED(gftype%screen_fn)) THEN
433 56 : CALL pw_pool%give_back_pw(gftype%screen_fn)
434 56 : DEALLOCATE (gftype%screen_fn)
435 : END IF
436 8054 : IF (ASSOCIATED(gftype%p3m_charge)) THEN
437 7682 : CALL pw_pool%give_back_pw(gftype%p3m_charge)
438 7682 : DEALLOCATE (gftype%p3m_charge)
439 : END IF
440 : ELSE
441 7740 : CALL gftype%influence_fn%release()
442 7740 : IF (ASSOCIATED(gftype%dct_influence_fn)) THEN
443 22 : CALL gftype%dct_influence_fn%release()
444 22 : DEALLOCATE (gftype%dct_influence_fn)
445 : END IF
446 7740 : IF (ASSOCIATED(gftype%screen_fn)) THEN
447 1156 : CALL gftype%screen_fn%release()
448 1156 : DEALLOCATE (gftype%screen_fn)
449 : END IF
450 7740 : IF (ASSOCIATED(gftype%p3m_charge)) THEN
451 2188 : CALL gftype%p3m_charge%release()
452 2188 : DEALLOCATE (gftype%p3m_charge)
453 : END IF
454 : END IF
455 15794 : IF (ALLOCATED(gftype%p3m_bm2)) &
456 9870 : DEALLOCATE (gftype%p3m_bm2)
457 15794 : IF (ALLOCATED(gftype%p3m_coeff)) &
458 9870 : DEALLOCATE (gftype%p3m_coeff)
459 15794 : END SUBROUTINE pw_green_release
460 :
461 : ! **************************************************************************************************
462 : !> \brief Calculates the influence_factor for the
463 : !> SPME Green's function in reciprocal space'''
464 : !> \param gftype ...
465 : !> \par History
466 : !> none
467 : !> \author DH (29-Mar-2001)
468 : ! **************************************************************************************************
469 9870 : SUBROUTINE influence_factor(gftype)
470 : TYPE(greens_fn_type), INTENT(INOUT) :: gftype
471 :
472 : COMPLEX(KIND=dp) :: b_m, exp_m, sum_m
473 : INTEGER :: dim, j, k, l, n, pt
474 : INTEGER, DIMENSION(3) :: npts
475 9870 : INTEGER, DIMENSION(:), POINTER :: lb, ub
476 : REAL(KIND=dp) :: l_arg, prod_arg, val
477 9870 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: m_assign
478 :
479 9870 : n = gftype%p3m_order
480 :
481 : ! calculate the assignment function values
482 :
483 9870 : lb => gftype%influence_fn%pw_grid%bounds(1, :)
484 9870 : ub => gftype%influence_fn%pw_grid%bounds(2, :)
485 9870 : IF (ALLOCATED(gftype%p3m_bm2)) THEN
486 0 : IF (LBOUND(gftype%p3m_bm2, 2) /= MINVAL(lb(:)) .OR. &
487 : UBOUND(gftype%p3m_bm2, 2) /= MAXVAL(ub(:))) THEN
488 0 : DEALLOCATE (gftype%p3m_bm2)
489 : END IF
490 : END IF
491 9870 : IF (.NOT. ALLOCATED(gftype%p3m_bm2)) THEN
492 88830 : ALLOCATE (gftype%p3m_bm2(3, MINVAL(lb(:)):MAXVAL(ub(:))))
493 : END IF
494 :
495 29610 : ALLOCATE (m_assign(0:n - 2))
496 58400 : m_assign = 0.0_dp
497 58400 : DO k = 0, n - 2
498 48530 : j = -(n - 1) + 2*k
499 347200 : DO l = 0, n - 1
500 288800 : l_arg = 0.5_dp**l
501 288800 : prod_arg = gftype%p3m_coeff(j, l)*l_arg
502 337330 : m_assign(k) = m_assign(k) + prod_arg
503 : END DO
504 : END DO
505 :
506 : ! calculate the absolute b values
507 :
508 39480 : npts(:) = ub(:) - lb(:) + 1
509 39480 : DO dim = 1, 3
510 943190 : DO pt = lb(dim), ub(dim)
511 903710 : val = twopi*(REAL(pt, KIND=dp)/REAL(npts(dim), KIND=dp))
512 903710 : exp_m = CMPLX(COS(val), -SIN(val), KIND=dp)
513 903710 : sum_m = z_zero
514 5349240 : DO k = 0, n - 2
515 5349240 : sum_m = sum_m + m_assign(k)*exp_m**k
516 : END DO
517 903710 : b_m = exp_m**(n - 1)/sum_m
518 933320 : gftype%p3m_bm2(dim, pt) = SQRT(REAL(b_m*CONJG(b_m), KIND=dp))
519 : END DO
520 : END DO
521 :
522 9870 : DEALLOCATE (m_assign)
523 9870 : END SUBROUTINE influence_factor
524 :
525 : ! **************************************************************************************************
526 : !> \brief ...
527 : !> \param gf ...
528 : ! **************************************************************************************************
529 9870 : PURE SUBROUTINE calc_p3m_charge(gf)
530 :
531 : TYPE(greens_fn_type), INTENT(INOUT) :: gf
532 :
533 : INTEGER :: ig, l, m, n
534 : REAL(KIND=dp) :: arg, novol
535 :
536 : ! check if charge function is consistent with current box volume
537 :
538 : ASSOCIATE (grid => gf%influence_fn%pw_grid, bm2 => gf%p3m_bm2)
539 9870 : arg = 1.0_dp/(8.0_dp*gf%p3m_alpha**2)
540 9870 : novol = REAL(grid%ngpts, KIND=dp)/grid%vol
541 57442488 : DO ig = 1, grid%ngpts_cut_local
542 57432618 : l = grid%g_hat(1, ig)
543 57432618 : m = grid%g_hat(2, ig)
544 57432618 : n = grid%g_hat(3, ig)
545 : gf%p3m_charge%array(ig) = novol*EXP(-arg*grid%gsq(ig))* &
546 57442488 : bm2(1, l)*bm2(2, m)*bm2(3, n)
547 : END DO
548 : END ASSOCIATE
549 :
550 9870 : END SUBROUTINE calc_p3m_charge
551 :
552 : ! **************************************************************************************************
553 : !> \brief Initialize the poisson solver
554 : !> You should call this just before calling the work routine
555 : !> pw_poisson_solver
556 : !> Call pw_poisson_release when you have finished
557 : !> \param poisson_env ...
558 : !> \par History
559 : !> none
560 : !> \author JGH (12-Mar-2001)
561 : ! **************************************************************************************************
562 10944 : SUBROUTINE pw_poisson_create(poisson_env)
563 :
564 : CLASS(pw_poisson_type), INTENT(INOUT) :: poisson_env
565 :
566 : ! Cleanup a potential previous poisson_env
567 10944 : CALL poisson_env%release()
568 :
569 10944 : END SUBROUTINE pw_poisson_create
570 :
571 : ! **************************************************************************************************
572 : !> \brief releases the poisson solver
573 : !> \param poisson_env ...
574 : !> \par History
575 : !> none
576 : !> \author fawzi (11.2002)
577 : ! **************************************************************************************************
578 21888 : SUBROUTINE pw_poisson_release(poisson_env)
579 :
580 : CLASS(pw_poisson_type), INTENT(INOUT) :: poisson_env
581 :
582 21888 : IF (ASSOCIATED(poisson_env%pw_pools)) THEN
583 10944 : CALL pw_pools_dealloc(poisson_env%pw_pools)
584 : END IF
585 :
586 21888 : IF (ASSOCIATED(poisson_env%green_fft)) THEN
587 7740 : CALL pw_green_release(poisson_env%green_fft)
588 7740 : DEALLOCATE (poisson_env%green_fft)
589 : END IF
590 21888 : CALL pw_grid_release(poisson_env%mt_super_ref_pw_grid)
591 21888 : CALL ps_wavelet_release(poisson_env%wavelet)
592 : CALL ps_implicit_release(poisson_env%implicit_env, &
593 21888 : poisson_env%parameters%ps_implicit_params)
594 21888 : CALL pw_grid_release(poisson_env%dct_pw_grid)
595 21888 : IF (ASSOCIATED(poisson_env%diel_rs_grid)) THEN
596 50 : CALL rs_grid_release(poisson_env%diel_rs_grid)
597 50 : DEALLOCATE (poisson_env%diel_rs_grid)
598 : END IF
599 :
600 21888 : END SUBROUTINE pw_poisson_release
601 :
602 : ! **************************************************************************************************
603 : !> \brief Calculates the coefficients for the charge assignment function
604 : !> \param n ...
605 : !> \param coeff ...
606 : !> \par History
607 : !> none
608 : !> \author DG (29-Mar-2001)
609 : ! **************************************************************************************************
610 9870 : SUBROUTINE spme_coeff_calculate(n, coeff)
611 :
612 : INTEGER, INTENT(IN) :: n
613 : REAL(KIND=dp), DIMENSION(-(n-1):n-1, 0:n-1), &
614 : INTENT(OUT) :: coeff
615 :
616 : INTEGER :: i, j, l, m
617 : REAL(KIND=dp) :: b
618 9870 : REAL(KIND=dp), DIMENSION(n, -n:n, 0:n-1) :: a
619 :
620 5311598 : a = 0.0_dp
621 9870 : a(1, 0, 0) = 1.0_dp
622 :
623 58400 : DO i = 2, n
624 48530 : m = i - 1
625 251330 : DO j = -m, m, 2
626 864018 : DO l = 0, m - 1
627 : b = (a(m, j - 1, l) + &
628 : REAL((-1)**l, KIND=dp)*a(m, j + 1, l))/ &
629 671088 : REAL((l + 1)*2**(l + 1), KIND=dp)
630 864018 : a(i, j, 0) = a(i, j, 0) + b
631 : END DO
632 912548 : DO l = 0, m - 1
633 : a(i, j, l + 1) = (a(m, j + 1, l) - &
634 864018 : a(m, j - 1, l))/REAL(l + 1, KIND=dp)
635 : END DO
636 : END DO
637 : END DO
638 :
639 704270 : coeff = 0.0_dp
640 68270 : DO i = 0, n - 1
641 68270 : DO j = -(n - 1), n - 1, 2
642 347200 : coeff(j, i) = a(n, j, i)
643 : END DO
644 : END DO
645 :
646 9870 : END SUBROUTINE spme_coeff_calculate
647 :
648 0 : END MODULE pw_poisson_types
|
