Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 19222 : 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 76888 : DO i = 1, 3
199 76888 : abc(i) = cell_hmat(i, i)
200 : END DO
201 76888 : dim = COUNT(poisson_params%periodic == 1)
202 :
203 36444 : SELECT CASE (poisson_params%solver)
204 : CASE (pw_poisson_periodic)
205 : green%method = PERIODIC3D
206 17222 : IF (dim /= 3) THEN
207 0 : CPABORT("Illegal combination of periodicity and Poisson solver periodic3d")
208 : END IF
209 : CASE (pw_poisson_multipole)
210 46 : green%method = MULTIPOLE0D
211 46 : IF (dim /= 0) THEN
212 0 : CPABORT("Illegal combination of periodicity and Poisson solver mulipole0d")
213 : END IF
214 : CASE (pw_poisson_analytic)
215 1796 : SELECT CASE (dim)
216 : CASE (0)
217 242 : green%method = ANALYTIC0D
218 968 : green%radius = 0.5_dp*MINVAL(abc)
219 : CASE (1)
220 8 : green%method = ANALYTIC1D
221 32 : green%special_dimension = MAXLOC(poisson_params%periodic(1:3), 1)
222 32 : green%radius = MAXVAL(abc(1:3))
223 32 : DO i = 1, 3
224 24 : IF (i == green%special_dimension) CYCLE
225 32 : green%radius = MIN(green%radius, 0.5_dp*abc(i))
226 : END DO
227 : CASE (2)
228 82 : green%method = ANALYTIC2D
229 328 : i = MINLOC(poisson_params%periodic, 1)
230 82 : green%special_dimension = i
231 82 : green%slab_size = abc(i)
232 : CASE (3)
233 0 : green%method = PERIODIC3D
234 : CASE DEFAULT
235 346 : CPABORT("Periodicity that is not 0, 1, 2, 3 is invalid")
236 : END SELECT
237 : CASE (pw_poisson_mt)
238 1554 : green%MT_rel_cutoff = poisson_params%mt_rel_cutoff
239 1554 : green%MT_alpha = poisson_params%mt_alpha
240 6216 : green%MT_alpha = green%MT_alpha/MINVAL(abc)
241 1542 : SELECT CASE (dim)
242 : CASE (0)
243 1488 : green%method = MT0D
244 5952 : 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 66 : green%method = MT2D
255 264 : i = MINLOC(poisson_params%periodic, 1)
256 66 : green%special_dimension = i
257 66 : green%slab_size = abc(i)
258 : CASE (3)
259 0 : CPABORT("Illegal combination of periodicity and Poisson solver (MT)")
260 : CASE DEFAULT
261 1554 : CPABORT("Periodicity that is not 0, 1, 2, 3 is invalid")
262 : END SELECT
263 : CASE (pw_poisson_implicit)
264 54 : green%method = PS_IMPLICIT
265 : CASE DEFAULT
266 19222 : CPABORT("An unknown Poisson solver was specified")
267 : END SELECT
268 :
269 : ! allocate influence function,...
270 38444 : SELECT CASE (green%method)
271 : CASE (PERIODIC3D, ANALYTIC2D, ANALYTIC1D, ANALYTIC0D, MT2D, MT1D, MT0D, MULTIPOLE0D, PS_IMPLICIT)
272 19222 : CALL pw_pool%create_pw(green%influence_fn)
273 :
274 19222 : IF (poisson_params%ewald_type == do_ewald_spme) THEN
275 11970 : green%p3m = .TRUE.
276 11970 : green%p3m_order = poisson_params%ewald_o_spline
277 11970 : green%p3m_alpha = poisson_params%ewald_alpha
278 11970 : n = green%p3m_order
279 47880 : ALLOCATE (green%p3m_coeff(-(n - 1):n - 1, 0:n - 1))
280 11970 : CALL spme_coeff_calculate(n, green%p3m_coeff)
281 : NULLIFY (green%p3m_charge)
282 11970 : ALLOCATE (green%p3m_charge)
283 11970 : CALL pw_pool%create_pw(green%p3m_charge)
284 11970 : CALL influence_factor(green)
285 11970 : CALL calc_p3m_charge(green)
286 : ELSE
287 7252 : green%p3m = .FALSE.
288 : END IF
289 : !
290 20776 : 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 1554 : 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 19222 : (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 19222 : CPABORT("Unknown method for green function")
309 : END SELECT
310 :
311 : ! initialize influence function
312 : ASSOCIATE (gf => green%influence_fn, grid => green%influence_fn%pw_grid)
313 36504 : SELECT CASE (green%method)
314 : CASE (PERIODIC3D, MULTIPOLE0D)
315 :
316 398918412 : DO ig = grid%first_gne0, grid%ngpts_cut_local
317 398901130 : g2 = grid%gsq(ig)
318 398918412 : gf%array(ig) = fourpi/g2
319 : END DO
320 17282 : IF (grid%have_g0) gf%array(1) = 0.0_dp
321 :
322 : CASE (ANALYTIC2D)
323 :
324 82 : iz = green%special_dimension ! iz is the direction with NO PBC
325 82 : zlength = green%slab_size ! zlength is the thickness of the cell
326 9304649 : DO ig = grid%first_gne0, grid%ngpts_cut_local
327 9304567 : nz = grid%g_hat(iz, ig)
328 9304567 : g2 = grid%gsq(ig)
329 9304567 : g3d = fourpi/g2
330 9304567 : gxy = MAX(0.0_dp, g2 - grid%g(iz, ig)*grid%g(iz, ig))
331 9304567 : gg = 0.5_dp*SQRT(gxy)
332 9304649 : gf%array(ig) = g3d*(1.0_dp - (-1.0_dp)**nz*EXP(-gg*zlength))
333 : END DO
334 82 : 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 8 : iz = green%special_dimension
340 : ! rlength is the radius of the tube
341 8 : rlength = green%radius
342 483604 : DO ig = grid%first_gne0, grid%ngpts_cut_local
343 483596 : g2 = grid%gsq(ig)
344 483596 : g3d = fourpi/g2
345 483596 : gxy = SQRT(MAX(0.0_dp, g2 - grid%g(iz, ig)*grid%g(iz, ig)))
346 483596 : gz = ABS(grid%g(iz, ig))
347 483596 : j0g = BESSEL_J0(rlength*gxy)
348 483596 : j1g = BESSEL_J1(rlength*gxy)
349 483596 : IF (gz > 0) THEN
350 470700 : k0g = bessk0(rlength*gz)
351 470700 : 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 483604 : (gxy*j1g*k0g - gz*j0g*k1g))
358 : END DO
359 8 : IF (grid%have_g0) gf%array(1) = 0.0_dp
360 :
361 : CASE (ANALYTIC0D)
362 :
363 242 : rlength = green%radius ! rlength is the radius of the sphere
364 47735002 : DO ig = grid%first_gne0, grid%ngpts_cut_local
365 47734760 : g2 = grid%gsq(ig)
366 47734760 : gg = SQRT(g2)
367 47734760 : g3d = fourpi/g2
368 47735002 : gf%array(ig) = g3d*(1.0_dp - COS(rlength*gg))
369 : END DO
370 242 : IF (grid%have_g0) THEN
371 124 : gf%array(1) = 0.5_dp*fourpi*rlength*rlength
372 : END IF
373 :
374 : CASE (MT2D, MT1D, MT0D)
375 :
376 111182694 : DO ig = grid%first_gne0, grid%ngpts_cut_local
377 111181140 : g2 = grid%gsq(ig)
378 111181140 : g3d = fourpi/g2
379 111182694 : gf%array(ig) = g3d + green%screen_fn%array(ig)
380 : END DO
381 1554 : IF (grid%have_g0) THEN
382 783 : gf%array(1) = green%screen_fn%array(1)
383 : END IF
384 :
385 : CASE (PS_IMPLICIT)
386 :
387 6554816 : DO ig = grid%first_gne0, grid%ngpts_cut_local
388 6554762 : g2 = grid%gsq(ig)
389 6554816 : gf%array(ig) = fourpi/g2
390 : END DO
391 54 : IF (grid%have_g0) gf%array(1) = 0.0_dp
392 :
393 54 : IF (ASSOCIATED(green%dct_influence_fn)) THEN
394 22 : dct_gf => green%dct_influence_fn
395 22 : dct_grid => green%dct_influence_fn%pw_grid
396 :
397 9851771 : DO ig = dct_grid%first_gne0, dct_grid%ngpts_cut_local
398 9851749 : g2 = dct_grid%gsq(ig)
399 9851771 : dct_gf%array(ig) = fourpi/g2
400 : END DO
401 22 : IF (dct_grid%have_g0) dct_gf%array(1) = 0.0_dp
402 : END IF
403 :
404 : CASE DEFAULT
405 19222 : CPABORT("Unknown method for green function")
406 : END SELECT
407 : END ASSOCIATE
408 :
409 19222 : END SUBROUTINE pw_green_create
410 :
411 : ! **************************************************************************************************
412 : !> \brief destroys the type (deallocates data)
413 : !> \param gftype ...
414 : !> \param pw_pool ...
415 : !> \par History
416 : !> none
417 : !> \author Joost VandeVondele
418 : !> Teodoro Laino
419 : ! **************************************************************************************************
420 19222 : SUBROUTINE pw_green_release(gftype, pw_pool)
421 : TYPE(greens_fn_type), INTENT(INOUT) :: gftype
422 : TYPE(pw_pool_type), OPTIONAL, POINTER :: pw_pool
423 :
424 : LOGICAL :: can_give_back
425 :
426 19222 : can_give_back = PRESENT(pw_pool)
427 19222 : IF (can_give_back) can_give_back = ASSOCIATED(pw_pool)
428 19222 : IF (can_give_back) THEN
429 10324 : CALL pw_pool%give_back_pw(gftype%influence_fn)
430 10324 : IF (ASSOCIATED(gftype%dct_influence_fn)) THEN
431 0 : CALL pw_pool%give_back_pw(gftype%dct_influence_fn)
432 0 : DEALLOCATE (gftype%dct_influence_fn)
433 : END IF
434 10324 : IF (ASSOCIATED(gftype%screen_fn)) THEN
435 112 : CALL pw_pool%give_back_pw(gftype%screen_fn)
436 112 : DEALLOCATE (gftype%screen_fn)
437 : END IF
438 10324 : IF (ASSOCIATED(gftype%p3m_charge)) THEN
439 9590 : CALL pw_pool%give_back_pw(gftype%p3m_charge)
440 9590 : DEALLOCATE (gftype%p3m_charge)
441 : END IF
442 : ELSE
443 8898 : CALL gftype%influence_fn%release()
444 8898 : IF (ASSOCIATED(gftype%dct_influence_fn)) THEN
445 22 : CALL gftype%dct_influence_fn%release()
446 22 : DEALLOCATE (gftype%dct_influence_fn)
447 : END IF
448 8898 : IF (ASSOCIATED(gftype%screen_fn)) THEN
449 1442 : CALL gftype%screen_fn%release()
450 1442 : DEALLOCATE (gftype%screen_fn)
451 : END IF
452 8898 : IF (ASSOCIATED(gftype%p3m_charge)) THEN
453 2380 : CALL gftype%p3m_charge%release()
454 2380 : DEALLOCATE (gftype%p3m_charge)
455 : END IF
456 : END IF
457 19222 : IF (ALLOCATED(gftype%p3m_bm2)) THEN
458 11970 : DEALLOCATE (gftype%p3m_bm2)
459 : END IF
460 19222 : IF (ALLOCATED(gftype%p3m_coeff)) THEN
461 11970 : DEALLOCATE (gftype%p3m_coeff)
462 : END IF
463 19222 : END SUBROUTINE pw_green_release
464 :
465 : ! **************************************************************************************************
466 : !> \brief Calculates the influence_factor for the
467 : !> SPME Green's function in reciprocal space'''
468 : !> \param gftype ...
469 : !> \par History
470 : !> none
471 : !> \author DH (29-Mar-2001)
472 : ! **************************************************************************************************
473 11970 : SUBROUTINE influence_factor(gftype)
474 : TYPE(greens_fn_type), INTENT(INOUT) :: gftype
475 :
476 : COMPLEX(KIND=dp) :: b_m, exp_m, sum_m
477 : INTEGER :: dim, j, k, l, n, pt
478 : INTEGER, DIMENSION(3) :: npts
479 11970 : INTEGER, DIMENSION(:), POINTER :: lb, ub
480 : REAL(KIND=dp) :: l_arg, prod_arg, val
481 11970 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: m_assign
482 :
483 11970 : n = gftype%p3m_order
484 :
485 : ! calculate the assignment function values
486 :
487 11970 : lb => gftype%influence_fn%pw_grid%bounds(1, :)
488 11970 : ub => gftype%influence_fn%pw_grid%bounds(2, :)
489 11970 : IF (ALLOCATED(gftype%p3m_bm2)) THEN
490 0 : IF (LBOUND(gftype%p3m_bm2, 2) /= MINVAL(lb(:)) .OR. &
491 : UBOUND(gftype%p3m_bm2, 2) /= MAXVAL(ub(:))) THEN
492 0 : DEALLOCATE (gftype%p3m_bm2)
493 : END IF
494 : END IF
495 11970 : IF (.NOT. ALLOCATED(gftype%p3m_bm2)) THEN
496 107730 : ALLOCATE (gftype%p3m_bm2(3, MINVAL(lb(:)):MAXVAL(ub(:))))
497 : END IF
498 :
499 35910 : ALLOCATE (m_assign(0:n - 2))
500 11970 : m_assign = 0.0_dp
501 70828 : DO k = 0, n - 2
502 58858 : j = -(n - 1) + 2*k
503 421080 : DO l = 0, n - 1
504 350252 : l_arg = 0.5_dp**l
505 350252 : prod_arg = gftype%p3m_coeff(j, l)*l_arg
506 409110 : m_assign(k) = m_assign(k) + prod_arg
507 : END DO
508 : END DO
509 :
510 : ! calculate the absolute b values
511 :
512 47880 : npts(:) = ub(:) - lb(:) + 1
513 47880 : DO dim = 1, 3
514 949196 : DO pt = lb(dim), ub(dim)
515 901316 : val = twopi*(REAL(pt, KIND=dp)/REAL(npts(dim), KIND=dp))
516 901316 : exp_m = CMPLX(COS(val), -SIN(val), KIND=dp)
517 901316 : sum_m = z_zero
518 5315136 : DO k = 0, n - 2
519 5315136 : sum_m = sum_m + m_assign(k)*exp_m**k
520 : END DO
521 901316 : b_m = exp_m**(n - 1)/sum_m
522 937226 : gftype%p3m_bm2(dim, pt) = SQRT(REAL(b_m*CONJG(b_m), KIND=dp))
523 : END DO
524 : END DO
525 :
526 11970 : DEALLOCATE (m_assign)
527 11970 : END SUBROUTINE influence_factor
528 :
529 : ! **************************************************************************************************
530 : !> \brief ...
531 : !> \param gf ...
532 : ! **************************************************************************************************
533 11970 : PURE SUBROUTINE calc_p3m_charge(gf)
534 :
535 : TYPE(greens_fn_type), INTENT(INOUT) :: gf
536 :
537 : INTEGER :: ig, l, m, n
538 : REAL(KIND=dp) :: arg, novol
539 :
540 : ! check if charge function is consistent with current box volume
541 :
542 : ASSOCIATE (grid => gf%influence_fn%pw_grid, bm2 => gf%p3m_bm2)
543 11970 : arg = 1.0_dp/(8.0_dp*gf%p3m_alpha**2)
544 11970 : novol = REAL(grid%ngpts, KIND=dp)/grid%vol
545 54965590 : DO ig = 1, grid%ngpts_cut_local
546 54953620 : l = grid%g_hat(1, ig)
547 54953620 : m = grid%g_hat(2, ig)
548 54953620 : n = grid%g_hat(3, ig)
549 : gf%p3m_charge%array(ig) = novol*EXP(-arg*grid%gsq(ig))* &
550 54965590 : bm2(1, l)*bm2(2, m)*bm2(3, n)
551 : END DO
552 : END ASSOCIATE
553 :
554 11970 : END SUBROUTINE calc_p3m_charge
555 :
556 : ! **************************************************************************************************
557 : !> \brief Initialize the poisson solver
558 : !> You should call this just before calling the work routine
559 : !> pw_poisson_solver
560 : !> Call pw_poisson_release when you have finished
561 : !> \param poisson_env ...
562 : !> \par History
563 : !> none
564 : !> \author JGH (12-Mar-2001)
565 : ! **************************************************************************************************
566 12524 : SUBROUTINE pw_poisson_create(poisson_env)
567 :
568 : CLASS(pw_poisson_type), INTENT(INOUT) :: poisson_env
569 :
570 : ! Cleanup a potential previous poisson_env
571 12524 : CALL poisson_env%release()
572 :
573 12524 : END SUBROUTINE pw_poisson_create
574 :
575 : ! **************************************************************************************************
576 : !> \brief releases the poisson solver
577 : !> \param poisson_env ...
578 : !> \par History
579 : !> none
580 : !> \author fawzi (11.2002)
581 : ! **************************************************************************************************
582 25048 : SUBROUTINE pw_poisson_release(poisson_env)
583 :
584 : CLASS(pw_poisson_type), INTENT(INOUT) :: poisson_env
585 :
586 25048 : IF (ASSOCIATED(poisson_env%pw_pools)) THEN
587 12524 : CALL pw_pools_dealloc(poisson_env%pw_pools)
588 : END IF
589 :
590 25048 : IF (ASSOCIATED(poisson_env%green_fft)) THEN
591 8898 : CALL pw_green_release(poisson_env%green_fft)
592 8898 : DEALLOCATE (poisson_env%green_fft)
593 : END IF
594 25048 : CALL pw_grid_release(poisson_env%mt_super_ref_pw_grid)
595 25048 : CALL ps_wavelet_release(poisson_env%wavelet)
596 : CALL ps_implicit_release(poisson_env%implicit_env, &
597 25048 : poisson_env%parameters%ps_implicit_params)
598 25048 : CALL pw_grid_release(poisson_env%dct_pw_grid)
599 25048 : IF (ASSOCIATED(poisson_env%diel_rs_grid)) THEN
600 50 : CALL rs_grid_release(poisson_env%diel_rs_grid)
601 50 : DEALLOCATE (poisson_env%diel_rs_grid)
602 : END IF
603 :
604 25048 : END SUBROUTINE pw_poisson_release
605 :
606 : ! **************************************************************************************************
607 : !> \brief Calculates the coefficients for the charge assignment function
608 : !> \param n ...
609 : !> \param coeff ...
610 : !> \par History
611 : !> none
612 : !> \author DG (29-Mar-2001)
613 : ! **************************************************************************************************
614 11970 : SUBROUTINE spme_coeff_calculate(n, coeff)
615 :
616 : INTEGER, INTENT(IN) :: n
617 : REAL(KIND=dp), DIMENSION(-(n-1):n-1, 0:n-1), &
618 : INTENT(OUT) :: coeff
619 :
620 : INTEGER :: i, j, l, m
621 : REAL(KIND=dp) :: b
622 11970 : REAL(KIND=dp), DIMENSION(n, -n:n, 0:n-1) :: a
623 :
624 6441250 : a = 0.0_dp
625 11970 : a(1, 0, 0) = 1.0_dp
626 :
627 70828 : DO i = 2, n
628 58858 : m = i - 1
629 304812 : DO j = -m, m, 2
630 1047772 : DO l = 0, m - 1
631 : b = (a(m, j - 1, l) + &
632 : REAL((-1)**l, KIND=dp)*a(m, j + 1, l))/ &
633 813788 : REAL((l + 1)*2**(l + 1), KIND=dp)
634 1047772 : a(i, j, 0) = a(i, j, 0) + b
635 : END DO
636 1106630 : DO l = 0, m - 1
637 : a(i, j, l + 1) = (a(m, j + 1, l) - &
638 1047772 : a(m, j - 1, l))/REAL(l + 1, KIND=dp)
639 : END DO
640 : END DO
641 : END DO
642 :
643 854130 : coeff = 0.0_dp
644 82798 : DO i = 0, n - 1
645 82798 : DO j = -(n - 1), n - 1, 2
646 421080 : coeff(j, i) = a(n, j, i)
647 : END DO
648 : END DO
649 :
650 11970 : END SUBROUTINE spme_coeff_calculate
651 :
652 0 : END MODULE pw_poisson_types
|
