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