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 : !> \par History
10 : !> 09.2005 created [fawzi]
11 : !> \author fawzi
12 : ! **************************************************************************************************
13 : MODULE pw_poisson_methods
14 :
15 : USE cp_log_handling, ONLY: cp_to_string
16 : USE dielectric_methods, ONLY: dielectric_compute
17 : USE kinds, ONLY: dp
18 : USE mathconstants, ONLY: fourpi
19 : USE mt_util, ONLY: MT0D, &
20 : MT1D, &
21 : MT2D
22 : USE ps_implicit_methods, ONLY: implicit_poisson_solver_mixed, &
23 : implicit_poisson_solver_mixed_periodic, &
24 : implicit_poisson_solver_neumann, &
25 : implicit_poisson_solver_periodic, &
26 : ps_implicit_create
27 : USE ps_implicit_types, ONLY: MIXED_BC, &
28 : MIXED_PERIODIC_BC, &
29 : NEUMANN_BC, &
30 : PERIODIC_BC
31 : USE ps_wavelet_methods, ONLY: cp2k_distribution_to_z_slices, &
32 : ps_wavelet_create, &
33 : ps_wavelet_solve, &
34 : z_slices_to_cp2k_distribution
35 : USE ps_wavelet_types, ONLY: WAVELET0D, &
36 : WAVELET1D, &
37 : WAVELET2D, &
38 : WAVELET3D, &
39 : ps_wavelet_type
40 : USE pw_grid_types, ONLY: pw_grid_type
41 : USE pw_grids, ONLY: pw_grid_compare, &
42 : pw_grid_release, &
43 : pw_grid_retain
44 : USE pw_methods, ONLY: pw_copy, &
45 : pw_derive, &
46 : pw_integral_ab, &
47 : pw_transfer, pw_multiply_with
48 : USE pw_poisson_types, ONLY: &
49 : ANALYTIC0D, ANALYTIC1D, ANALYTIC2D, MULTIPOLE0D, PERIODIC3D, PS_IMPLICIT, do_ewald_spme, &
50 : greens_fn_type, pw_green_create, pw_green_release, pw_poisson_analytic, &
51 : pw_poisson_implicit, pw_poisson_mt, pw_poisson_multipole, pw_poisson_none, &
52 : pw_poisson_parameter_type, pw_poisson_periodic, pw_poisson_type, pw_poisson_wavelet
53 : USE pw_pool_types, ONLY: pw_pool_p_type, &
54 : pw_pool_type, &
55 : pw_pools_copy, &
56 : pw_pools_dealloc
57 : USE pw_types, ONLY: &
58 : pw_r3d_rs_type, pw_c1d_gs_type, pw_r3d_rs_type
59 : #include "../base/base_uses.f90"
60 :
61 : IMPLICIT NONE
62 : PRIVATE
63 :
64 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
65 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pw_poisson_methods'
66 :
67 : PUBLIC :: pw_poisson_rebuild, &
68 : pw_poisson_solve, pw_poisson_set, pw_func_u_convolution
69 :
70 : INTEGER, PARAMETER :: use_rs_grid = 0, &
71 : use_gs_grid = 1
72 :
73 : INTERFACE pw_poisson_rebuild
74 : MODULE PROCEDURE pw_poisson_rebuild_nodens
75 : MODULE PROCEDURE pw_poisson_rebuild_c1d_gs, pw_poisson_rebuild_r3d_rs
76 : END INTERFACE
77 :
78 : INTERFACE pw_poisson_solve
79 : #:for kindd in ['r3d_rs', 'c1d_gs']
80 : MODULE PROCEDURE pw_poisson_solve_nov_nodv_${kindd}$
81 : #:for kindv in ['r3d_rs', 'c1d_gs']
82 : MODULE PROCEDURE pw_poisson_solve_v_nodv_${kindd}$_${kindv}$
83 : #:endfor
84 : #:for kindg in ['r3d_rs', 'c1d_gs']
85 : MODULE PROCEDURE pw_poisson_solve_nov_dv_${kindd}$_${kindg}$
86 : #:for kindv in ['r3d_rs', 'c1d_gs']
87 : MODULE PROCEDURE pw_poisson_solve_v_dv_${kindd}$_${kindv}$_${kindg}$
88 : #:endfor
89 : #:endfor
90 : #:endfor
91 : END INTERFACE
92 :
93 : CONTAINS
94 :
95 : ! **************************************************************************************************
96 : !> \brief removes all the object created from the parameters pw_pools and cell
97 : !> and used to solve the poisson equation like the green function and
98 : !> all the things allocated in pw_poisson_rebuild
99 : !> \param poisson_env ...
100 : !> \par History
101 : !> none
102 : ! **************************************************************************************************
103 71598 : SUBROUTINE pw_poisson_cleanup(poisson_env)
104 : TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
105 :
106 : TYPE(pw_pool_type), POINTER :: pw_pool
107 :
108 71598 : NULLIFY (pw_pool)
109 71598 : IF (ASSOCIATED(poisson_env%pw_pools)) THEN
110 59330 : pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
111 : END IF
112 71598 : IF (ASSOCIATED(poisson_env%green_fft)) THEN
113 10174 : CALL pw_green_release(poisson_env%green_fft, pw_pool=pw_pool)
114 10174 : DEALLOCATE (poisson_env%green_fft)
115 : END IF
116 71598 : poisson_env%rebuild = .TRUE.
117 :
118 71598 : END SUBROUTINE pw_poisson_cleanup
119 :
120 : ! **************************************************************************************************
121 : !> \brief checks if pw_poisson_rebuild has to be called and calls it if needed
122 : !> \param poisson_env the object to be checked
123 : !> \author fawzi
124 : ! **************************************************************************************************
125 26764 : SUBROUTINE pw_poisson_check(poisson_env)
126 : TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
127 :
128 : LOGICAL :: rebuild
129 : TYPE(greens_fn_type), POINTER :: green
130 : TYPE(ps_wavelet_type), POINTER :: wavelet
131 :
132 26764 : CPASSERT(ASSOCIATED(poisson_env%pw_pools))
133 53528 : CPASSERT(poisson_env%pw_level >= LBOUND(poisson_env%pw_pools, 1))
134 53528 : CPASSERT(poisson_env%pw_level <= UBOUND(poisson_env%pw_pools, 1))
135 26764 : green => poisson_env%green_fft
136 26764 : wavelet => poisson_env%wavelet
137 26764 : rebuild = poisson_env%rebuild
138 : rebuild = rebuild .OR. (poisson_env%method /= poisson_env%parameters%solver) &
139 26764 : .OR. .NOT. ASSOCIATED(green)
140 26764 : poisson_env%method = poisson_env%parameters%solver
141 :
142 26764 : IF (poisson_env%method == pw_poisson_wavelet) THEN
143 872 : poisson_env%used_grid = use_rs_grid
144 : ELSE
145 25892 : poisson_env%used_grid = use_gs_grid
146 : END IF
147 26764 : IF (.NOT. rebuild) THEN
148 0 : IF (poisson_env%parameters%ewald_type == do_ewald_spme) THEN
149 0 : rebuild = (poisson_env%parameters%ewald_alpha /= green%p3m_alpha) .OR. rebuild
150 0 : rebuild = (poisson_env%parameters%ewald_o_spline /= green%p3m_order) .OR. rebuild
151 : END IF
152 0 : SELECT CASE (poisson_env%method)
153 : CASE (pw_poisson_analytic)
154 0 : SELECT CASE (green%method)
155 : CASE (ANALYTIC0D, ANALYTIC1D, ANALYTIC2D, PERIODIC3D)
156 : CASE default
157 0 : rebuild = .TRUE.
158 : END SELECT
159 : CASE (pw_poisson_mt)
160 0 : SELECT CASE (green%method)
161 : CASE (MT0D, MT1D, MT2D)
162 : CASE default
163 0 : rebuild = .TRUE.
164 : END SELECT
165 0 : rebuild = (poisson_env%parameters%mt_alpha /= green%mt_alpha) .OR. rebuild
166 : CASE (pw_poisson_wavelet)
167 0 : rebuild = (poisson_env%parameters%wavelet_scf_type /= wavelet%itype_scf) .OR. rebuild
168 : CASE default
169 0 : CPABORT("Unknown Poisson solver")
170 : END SELECT
171 : END IF
172 0 : IF (rebuild) THEN
173 26764 : poisson_env%rebuild = .TRUE.
174 26764 : CALL pw_poisson_cleanup(poisson_env)
175 : END IF
176 26764 : END SUBROUTINE pw_poisson_check
177 :
178 : ! **************************************************************************************************
179 : !> \brief rebuilds all the internal values needed to use the poisson solver
180 : !> \param poisson_env the environment to rebuild
181 : !> \param density ...
182 : !> \author fawzi
183 : !> \note
184 : !> rebuilds if poisson_env%rebuild is true
185 : ! **************************************************************************************************
186 76714 : SUBROUTINE pw_poisson_rebuild_nodens(poisson_env)
187 : TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
188 :
189 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_rebuild'
190 :
191 : INTEGER :: handle
192 :
193 76714 : CALL timeset(routineN, handle)
194 :
195 76714 : CPASSERT(ASSOCIATED(poisson_env%pw_pools))
196 :
197 76714 : IF (poisson_env%rebuild) THEN
198 11832 : CALL pw_poisson_cleanup(poisson_env)
199 23664 : SELECT CASE (poisson_env%parameters%solver)
200 : CASE (pw_poisson_periodic, pw_poisson_analytic, pw_poisson_mt, pw_poisson_multipole)
201 11832 : ALLOCATE (poisson_env%green_fft)
202 : CALL pw_green_create(poisson_env%green_fft, cell_hmat=poisson_env%cell_hmat, &
203 : pw_pool=poisson_env%pw_pools(poisson_env%pw_level)%pool, &
204 : poisson_params=poisson_env%parameters, &
205 : mt_super_ref_pw_grid=poisson_env%mt_super_ref_pw_grid, &
206 11832 : dct_pw_grid=poisson_env%dct_pw_grid)
207 : CASE (pw_poisson_wavelet)
208 0 : CPABORT("Wavelet solver requires a density!")
209 : CASE (pw_poisson_implicit)
210 0 : ALLOCATE (poisson_env%green_fft)
211 : CALL pw_green_create(poisson_env%green_fft, cell_hmat=poisson_env%cell_hmat, &
212 : pw_pool=poisson_env%pw_pools(poisson_env%pw_level)%pool, &
213 : poisson_params=poisson_env%parameters, &
214 : mt_super_ref_pw_grid=poisson_env%mt_super_ref_pw_grid, &
215 0 : dct_pw_grid=poisson_env%dct_pw_grid)
216 : CALL ps_implicit_create(poisson_env%pw_pools(poisson_env%pw_level)%pool, &
217 : poisson_env%parameters, &
218 : poisson_env%dct_pw_grid, &
219 0 : poisson_env%green_fft, poisson_env%implicit_env)
220 : CASE (pw_poisson_none)
221 : CASE default
222 11832 : CPABORT("Unknown Poisson solver")
223 : END SELECT
224 11832 : poisson_env%rebuild = .FALSE.
225 : END IF
226 :
227 76714 : CALL timestop(handle)
228 :
229 76714 : END SUBROUTINE pw_poisson_rebuild_nodens
230 :
231 : #:for kindd in ["r3d_rs", "c1d_gs"]
232 : ! **************************************************************************************************
233 : !> \brief rebuilds all the internal values needed to use the poisson solver
234 : !> \param poisson_env the environment to rebuild
235 : !> \param density ...
236 : !> \author fawzi
237 : !> \note
238 : !> rebuilds if poisson_env%rebuild is true
239 : ! **************************************************************************************************
240 286708 : SUBROUTINE pw_poisson_rebuild_${kindd}$ (poisson_env, density)
241 : TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
242 : TYPE(pw_${kindd}$_type), INTENT(IN) :: density
243 :
244 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_rebuild'
245 :
246 : INTEGER :: handle
247 :
248 286708 : CALL timeset(routineN, handle)
249 :
250 286708 : CPASSERT(ASSOCIATED(poisson_env%pw_pools))
251 :
252 286708 : IF (poisson_env%rebuild) THEN
253 7908 : CALL pw_poisson_cleanup(poisson_env)
254 14906 : SELECT CASE (poisson_env%parameters%solver)
255 : CASE (pw_poisson_periodic, pw_poisson_analytic, pw_poisson_mt, pw_poisson_multipole)
256 6998 : ALLOCATE (poisson_env%green_fft)
257 : CALL pw_green_create(poisson_env%green_fft, cell_hmat=poisson_env%cell_hmat, &
258 : pw_pool=poisson_env%pw_pools(poisson_env%pw_level)%pool, &
259 : poisson_params=poisson_env%parameters, &
260 : mt_super_ref_pw_grid=poisson_env%mt_super_ref_pw_grid, &
261 6998 : dct_pw_grid=poisson_env%dct_pw_grid)
262 : CASE (pw_poisson_wavelet)
263 856 : CPASSERT(ASSOCIATED(density%pw_grid))
264 : CALL ps_wavelet_create(poisson_env%parameters, poisson_env%wavelet, &
265 856 : density%pw_grid)
266 : CASE (pw_poisson_implicit)
267 54 : ALLOCATE (poisson_env%green_fft)
268 : CALL pw_green_create(poisson_env%green_fft, cell_hmat=poisson_env%cell_hmat, &
269 : pw_pool=poisson_env%pw_pools(poisson_env%pw_level)%pool, &
270 : poisson_params=poisson_env%parameters, &
271 : mt_super_ref_pw_grid=poisson_env%mt_super_ref_pw_grid, &
272 54 : dct_pw_grid=poisson_env%dct_pw_grid)
273 : CALL ps_implicit_create(poisson_env%pw_pools(poisson_env%pw_level)%pool, &
274 : poisson_env%parameters, &
275 : poisson_env%dct_pw_grid, &
276 54 : poisson_env%green_fft, poisson_env%implicit_env)
277 : CASE (pw_poisson_none)
278 : CASE default
279 7908 : CPABORT("Unknown Poisson solver")
280 : END SELECT
281 7908 : poisson_env%rebuild = .FALSE.
282 : END IF
283 :
284 286708 : CALL timestop(handle)
285 :
286 286708 : END SUBROUTINE pw_poisson_rebuild_${kindd}$
287 : #:endfor
288 :
289 : #:for kindd in ['r3d_rs', 'c1d_gs']
290 : ! **************************************************************************************************
291 : !> \brief Solve Poisson equation in a plane wave basis set
292 : !> Obtains electrostatic potential and its derivatives with respect to r
293 : !> from the density
294 : !> \param poisson_env ...
295 : !> \param density ...
296 : !> \param ehartree ...
297 : !> \param h_stress ...
298 : !> \param rho_core ...
299 : !> \param greenfn ...
300 : !> \param aux_density Hartree energy and stress tensor between 2 different densities
301 : !> \par History
302 : !> JGH (13-Mar-2001) : completely revised
303 : !> \author apsi
304 : ! **************************************************************************************************
305 0 : SUBROUTINE pw_poisson_solve_nov_nodv_${kindd}$ (poisson_env, density, ehartree, &
306 : h_stress, rho_core, greenfn, aux_density)
307 :
308 : TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
309 : TYPE(pw_${kindd}$_type), INTENT(IN) :: density
310 : REAL(kind=dp), INTENT(out), OPTIONAL :: ehartree
311 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT), &
312 : OPTIONAL :: h_stress
313 : TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_core, greenfn
314 : TYPE(pw_${kindd}$_type), INTENT(IN), OPTIONAL :: aux_density
315 :
316 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_solve'
317 :
318 : INTEGER :: handle
319 : LOGICAL :: has_dielectric
320 : TYPE(pw_grid_type), POINTER :: pw_grid
321 : TYPE(pw_pool_type), POINTER :: pw_pool
322 : TYPE(pw_r3d_rs_type) :: rhor, vhartree_rs
323 : TYPE(pw_c1d_gs_type) :: influence_fn, rhog, rhog_aux, tmpg
324 :
325 0 : CALL timeset(routineN, handle)
326 :
327 0 : CALL pw_poisson_rebuild(poisson_env, density)
328 :
329 0 : has_dielectric = poisson_env%parameters%has_dielectric
330 :
331 : ! point pw
332 0 : pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
333 0 : pw_grid => pw_pool%pw_grid
334 : ! density in G space
335 0 : CALL pw_pool%create_pw(rhog)
336 0 : IF (PRESENT(aux_density)) THEN
337 0 : CALL pw_pool%create_pw(rhog_aux)
338 : END IF
339 :
340 0 : SELECT CASE (poisson_env%used_grid)
341 : CASE (use_gs_grid)
342 :
343 0 : SELECT CASE (poisson_env%green_fft%method)
344 : CASE (PERIODIC3D, ANALYTIC2D, ANALYTIC1D, ANALYTIC0D, MT2D, MT1D, MT0D, MULTIPOLE0D)
345 :
346 0 : CALL pw_transfer(density, rhog)
347 0 : IF (PRESENT(aux_density)) THEN
348 0 : CALL pw_transfer(aux_density, rhog_aux)
349 : END IF
350 0 : IF (PRESENT(ehartree)) THEN
351 0 : CALL pw_pool%create_pw(tmpg)
352 0 : CALL pw_copy(rhog, tmpg)
353 : END IF
354 0 : IF (PRESENT(greenfn)) THEN
355 0 : influence_fn = greenfn
356 : ELSE
357 0 : influence_fn = poisson_env%green_fft%influence_fn
358 : END IF
359 0 : CALL pw_multiply_with(rhog, influence_fn)
360 0 : IF (PRESENT(aux_density)) THEN
361 0 : CALL pw_multiply_with(rhog_aux, influence_fn)
362 : END IF
363 0 : IF (PRESENT(ehartree)) THEN
364 0 : IF (PRESENT(aux_density)) THEN
365 0 : ehartree = 0.5_dp*pw_integral_ab(rhog_aux, tmpg)
366 : ELSE
367 0 : ehartree = 0.5_dp*pw_integral_ab(rhog, tmpg)
368 : END IF
369 0 : CALL pw_pool%give_back_pw(tmpg)
370 : END IF
371 :
372 : CASE (PS_IMPLICIT)
373 :
374 0 : IF (PRESENT(h_stress)) THEN
375 0 : CPABORT("No stress tensor is implemented for the implicit Poisson solver.")
376 : END IF
377 :
378 0 : IF (has_dielectric .AND. PRESENT(rho_core)) THEN
379 0 : SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
380 : CASE (PERIODIC_BC, MIXED_PERIODIC_BC)
381 : CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
382 : poisson_env%diel_rs_grid, &
383 : poisson_env%pw_pools(poisson_env%pw_level)%pool, &
384 0 : density, rho_core=rho_core)
385 : CASE (NEUMANN_BC, MIXED_BC)
386 : CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
387 : poisson_env%diel_rs_grid, &
388 : poisson_env%pw_pools(poisson_env%pw_level)%pool, &
389 : poisson_env%dct_pw_grid, &
390 : poisson_env%parameters%ps_implicit_params%neumann_directions, &
391 : poisson_env%implicit_env%dct_env%recv_msgs_bnds, &
392 : poisson_env%implicit_env%dct_env%dests_expand, &
393 : poisson_env%implicit_env%dct_env%srcs_expand, &
394 : poisson_env%implicit_env%dct_env%flipg_stat, &
395 : poisson_env%implicit_env%dct_env%bounds_shftd, &
396 0 : density, rho_core=rho_core)
397 : END SELECT
398 : END IF
399 :
400 0 : CALL pw_pool%create_pw(rhor)
401 0 : CALL pw_pool%create_pw(vhartree_rs)
402 0 : CALL pw_transfer(density, rhor)
403 :
404 0 : SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
405 : CASE (PERIODIC_BC)
406 : CALL implicit_poisson_solver_periodic(poisson_env, rhor, vhartree_rs, &
407 0 : ehartree=ehartree)
408 : CASE (NEUMANN_BC)
409 : CALL implicit_poisson_solver_neumann(poisson_env, rhor, vhartree_rs, &
410 0 : ehartree=ehartree)
411 : CASE (MIXED_PERIODIC_BC)
412 : CALL implicit_poisson_solver_mixed_periodic(poisson_env, rhor, vhartree_rs, &
413 0 : electric_enthalpy=ehartree)
414 : CASE (MIXED_BC)
415 : CALL implicit_poisson_solver_mixed(poisson_env, rhor, vhartree_rs, &
416 0 : electric_enthalpy=ehartree)
417 : END SELECT
418 :
419 0 : IF (PRESENT(aux_density)) THEN
420 0 : CALL pw_transfer(aux_density, rhor)
421 0 : ehartree = 0.5_dp*pw_integral_ab(rhor, vhartree_rs)
422 : END IF
423 :
424 0 : CALL pw_pool%give_back_pw(rhor)
425 0 : CALL pw_pool%give_back_pw(vhartree_rs)
426 :
427 : CASE DEFAULT
428 : CALL cp_abort(__LOCATION__, &
429 : "unknown poisson method "// &
430 0 : cp_to_string(poisson_env%green_fft%method))
431 : END SELECT
432 :
433 : CASE (use_rs_grid)
434 :
435 0 : CALL pw_pool%create_pw(rhor)
436 0 : CALL pw_transfer(density, rhor)
437 0 : CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
438 0 : CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
439 0 : CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
440 0 : IF (PRESENT(ehartree)) THEN
441 0 : IF (PRESENT(aux_density)) THEN
442 : #:if kindd=="r3d_rs"
443 0 : ehartree = 0.5_dp*pw_integral_ab(aux_density, rhor)
444 : #:else
445 0 : IF (.NOT. PRESENT(h_stress)) CALL pw_pool%create_pw(rhog)
446 0 : CALL pw_transfer(rhor, rhog)
447 0 : ehartree = 0.5_dp*pw_integral_ab(aux_density, rhog)
448 0 : IF (.NOT. PRESENT(h_stress)) CALL pw_pool%give_back_pw(rhog)
449 : #:endif
450 : ELSE
451 : #:if kindd=="r3d_rs"
452 0 : ehartree = 0.5_dp*pw_integral_ab(density, rhor)
453 : #:else
454 0 : IF (.NOT. PRESENT(h_stress)) CALL pw_pool%create_pw(rhog)
455 0 : CALL pw_transfer(rhor, rhog)
456 0 : ehartree = 0.5_dp*pw_integral_ab(density, rhog)
457 0 : IF (.NOT. PRESENT(h_stress)) CALL pw_pool%give_back_pw(rhog)
458 : #:endif
459 : END IF
460 : END IF
461 0 : IF (PRESENT(h_stress)) THEN
462 0 : CALL pw_transfer(rhor, rhog)
463 0 : IF (PRESENT(aux_density)) THEN
464 0 : CALL pw_transfer(aux_density, rhor)
465 0 : CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
466 0 : CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
467 0 : CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
468 0 : CALL pw_transfer(rhor, rhog_aux)
469 : END IF
470 : END IF
471 0 : CALL pw_pool%give_back_pw(rhor)
472 :
473 : END SELECT
474 :
475 0 : IF (PRESENT(aux_density)) THEN
476 0 : CALL calc_stress_and_gradient_${kindd}$ (poisson_env, rhog, ehartree, rhog_aux, h_stress)
477 : ELSE
478 0 : CALL calc_stress_and_gradient_${kindd}$ (poisson_env, rhog, ehartree, h_stress=h_stress)
479 : END IF
480 :
481 0 : CALL pw_pool%give_back_pw(rhog)
482 0 : IF (PRESENT(aux_density)) THEN
483 0 : CALL pw_pool%give_back_pw(rhog_aux)
484 : END IF
485 :
486 0 : CALL timestop(handle)
487 :
488 0 : END SUBROUTINE pw_poisson_solve_nov_nodv_${kindd}$
489 : #:endfor
490 :
491 : #:for kindd in ['r3d_rs', 'c1d_gs']
492 : #:for kindv in ['r3d_rs', 'c1d_gs']
493 : ! **************************************************************************************************
494 : !> \brief Solve Poisson equation in a plane wave basis set
495 : !> Obtains electrostatic potential and its derivatives with respect to r
496 : !> from the density
497 : !> \param poisson_env ...
498 : !> \param density ...
499 : !> \param ehartree ...
500 : !> \param vhartree ...
501 : !> \param h_stress ...
502 : !> \param rho_core ...
503 : !> \param greenfn ...
504 : !> \param aux_density Hartree energy and stress tensor between 2 different densities
505 : !> \par History
506 : !> JGH (13-Mar-2001) : completely revised
507 : !> \author apsi
508 : ! **************************************************************************************************
509 221667 : SUBROUTINE pw_poisson_solve_v_nodv_${kindd}$_${kindv}$ (poisson_env, density, ehartree, vhartree, &
510 : h_stress, rho_core, greenfn, aux_density)
511 :
512 : TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
513 : TYPE(pw_${kindd}$_type), INTENT(IN) :: density
514 : REAL(kind=dp), INTENT(out), OPTIONAL :: ehartree
515 : TYPE(pw_${kindv}$_type), INTENT(INOUT) :: vhartree
516 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT), &
517 : OPTIONAL :: h_stress
518 : TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_core, greenfn
519 : TYPE(pw_${kindd}$_type), INTENT(IN), OPTIONAL :: aux_density
520 :
521 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_solve'
522 :
523 : INTEGER :: handle
524 : LOGICAL :: has_dielectric
525 : TYPE(pw_grid_type), POINTER :: pw_grid
526 : TYPE(pw_pool_type), POINTER :: pw_pool
527 : TYPE(pw_r3d_rs_type) :: &
528 : rhor, vhartree_rs
529 : TYPE(pw_c1d_gs_type) :: influence_fn, rhog, rhog_aux
530 :
531 221667 : CALL timeset(routineN, handle)
532 :
533 221667 : CALL pw_poisson_rebuild(poisson_env, density)
534 :
535 221667 : has_dielectric = poisson_env%parameters%has_dielectric
536 :
537 : ! point pw
538 221667 : pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
539 221667 : pw_grid => pw_pool%pw_grid
540 221667 : IF (.NOT. pw_grid_compare(pw_pool%pw_grid, vhartree%pw_grid)) THEN
541 0 : CPABORT("vhartree has a different grid than the poisson solver")
542 : END IF
543 : ! density in G space
544 221667 : CALL pw_pool%create_pw(rhog)
545 221667 : IF (PRESENT(aux_density)) THEN
546 404 : CALL pw_pool%create_pw(rhog_aux)
547 : END IF
548 :
549 410507 : SELECT CASE (poisson_env%used_grid)
550 : CASE (use_gs_grid)
551 :
552 410055 : SELECT CASE (poisson_env%green_fft%method)
553 : CASE (PERIODIC3D, ANALYTIC2D, ANALYTIC1D, ANALYTIC0D, MT2D, MT1D, MT0D, MULTIPOLE0D)
554 :
555 188388 : CALL pw_transfer(density, rhog)
556 188388 : IF (PRESENT(aux_density)) THEN
557 404 : CALL pw_transfer(aux_density, rhog_aux)
558 : END IF
559 188388 : IF (PRESENT(greenfn)) THEN
560 302 : influence_fn = greenfn
561 : ELSE
562 188086 : influence_fn = poisson_env%green_fft%influence_fn
563 : END IF
564 188388 : CALL pw_multiply_with(rhog, influence_fn)
565 188388 : IF (PRESENT(aux_density)) THEN
566 404 : CALL pw_multiply_with(rhog_aux, influence_fn)
567 : END IF
568 188388 : CALL pw_transfer(rhog, vhartree)
569 188388 : IF (PRESENT(ehartree)) THEN
570 139312 : IF (PRESENT(aux_density)) THEN
571 : #:if kindd==kindv
572 404 : ehartree = 0.5_dp*pw_integral_ab(aux_density, vhartree)
573 : #:elif kindd=="c1d_gs"
574 0 : ehartree = 0.5_dp*pw_integral_ab(aux_density, rhog)
575 : #:else
576 0 : CALL pw_transfer(aux_density, rhog)
577 0 : ehartree = 0.5_dp*pw_integral_ab(rhog, vhartree)
578 : #:endif
579 : ELSE
580 : #:if kindd==kindv
581 138908 : ehartree = 0.5_dp*pw_integral_ab(density, vhartree)
582 : #:elif kindd=="c1d_gs"
583 0 : ehartree = 0.5_dp*pw_integral_ab(density, rhog)
584 : #:else
585 0 : CALL pw_transfer(density, rhog)
586 0 : ehartree = 0.5_dp*pw_integral_ab(rhog, vhartree)
587 : #:endif
588 : END IF
589 : END IF
590 :
591 : CASE (PS_IMPLICIT)
592 452 : IF (PRESENT(h_stress)) THEN
593 0 : CPABORT("No stress tensor is implemented for the implicit Poisson solver.")
594 : END IF
595 :
596 452 : IF (has_dielectric .AND. PRESENT(rho_core)) THEN
597 748 : SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
598 : CASE (PERIODIC_BC, MIXED_PERIODIC_BC)
599 : CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
600 : poisson_env%diel_rs_grid, &
601 : poisson_env%pw_pools(poisson_env%pw_level)%pool, &
602 296 : density, rho_core=rho_core)
603 : CASE (NEUMANN_BC, MIXED_BC)
604 : CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
605 : poisson_env%diel_rs_grid, &
606 : poisson_env%pw_pools(poisson_env%pw_level)%pool, &
607 : poisson_env%dct_pw_grid, &
608 : poisson_env%parameters%ps_implicit_params%neumann_directions, &
609 : poisson_env%implicit_env%dct_env%recv_msgs_bnds, &
610 : poisson_env%implicit_env%dct_env%dests_expand, &
611 : poisson_env%implicit_env%dct_env%srcs_expand, &
612 : poisson_env%implicit_env%dct_env%flipg_stat, &
613 : poisson_env%implicit_env%dct_env%bounds_shftd, &
614 452 : density, rho_core=rho_core)
615 : END SELECT
616 : END IF
617 :
618 452 : CALL pw_pool%create_pw(rhor)
619 452 : CALL pw_pool%create_pw(vhartree_rs)
620 452 : CALL pw_transfer(density, rhor)
621 :
622 556 : SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
623 : CASE (PERIODIC_BC)
624 : CALL implicit_poisson_solver_periodic(poisson_env, rhor, vhartree_rs, &
625 128 : ehartree=ehartree)
626 : CASE (NEUMANN_BC)
627 : CALL implicit_poisson_solver_neumann(poisson_env, rhor, vhartree_rs, &
628 216 : ehartree=ehartree)
629 : CASE (MIXED_PERIODIC_BC)
630 : CALL implicit_poisson_solver_mixed_periodic(poisson_env, rhor, vhartree_rs, &
631 324 : electric_enthalpy=ehartree)
632 : CASE (MIXED_BC)
633 : CALL implicit_poisson_solver_mixed(poisson_env, rhor, vhartree_rs, &
634 452 : electric_enthalpy=ehartree)
635 : END SELECT
636 :
637 452 : IF (PRESENT(aux_density)) THEN
638 0 : CALL pw_transfer(aux_density, rhor)
639 0 : ehartree = 0.5_dp*pw_integral_ab(rhor, vhartree_rs)
640 : END IF
641 :
642 452 : CALL pw_transfer(vhartree_rs, vhartree)
643 :
644 452 : CALL pw_pool%give_back_pw(rhor)
645 452 : CALL pw_pool%give_back_pw(vhartree_rs)
646 :
647 : CASE DEFAULT
648 : CALL cp_abort(__LOCATION__, &
649 : "unknown poisson method "// &
650 188840 : cp_to_string(poisson_env%green_fft%method))
651 : END SELECT
652 :
653 : CASE (use_rs_grid)
654 :
655 32827 : CALL pw_pool%create_pw(rhor)
656 32827 : CALL pw_transfer(density, rhor)
657 32827 : CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
658 32827 : CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
659 32827 : CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
660 32827 : CALL pw_transfer(rhor, vhartree)
661 32827 : IF (PRESENT(ehartree)) THEN
662 8948 : IF (PRESENT(aux_density)) THEN
663 : #:if kindd==kindv
664 0 : ehartree = 0.5_dp*pw_integral_ab(aux_density, vhartree)
665 : #:elif kindd=="r3d_rs"
666 0 : ehartree = 0.5_dp*pw_integral_ab(aux_density, rhor)
667 : #:else
668 0 : CALL pw_transfer(vhartree, rhog)
669 0 : ehartree = 0.5_dp*pw_integral_ab(aux_density, rhog)
670 : #:endif
671 : ELSE
672 : #:if kindd==kindv
673 8948 : ehartree = 0.5_dp*pw_integral_ab(density, vhartree)
674 : #:elif kindd=="r3d_rs"
675 0 : ehartree = 0.5_dp*pw_integral_ab(density, rhor)
676 : #:else
677 0 : CALL pw_transfer(vhartree, rhog)
678 0 : ehartree = 0.5_dp*pw_integral_ab(density, rhog)
679 : #:endif
680 : END IF
681 : END IF
682 32827 : IF (PRESENT(h_stress)) THEN
683 12 : CALL pw_transfer(rhor, rhog)
684 12 : IF (PRESENT(aux_density)) THEN
685 0 : CALL pw_transfer(aux_density, rhor)
686 0 : CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
687 0 : CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
688 0 : CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
689 0 : CALL pw_transfer(rhor, rhog_aux)
690 : END IF
691 : END IF
692 254494 : CALL pw_pool%give_back_pw(rhor)
693 :
694 : END SELECT
695 :
696 221667 : IF (PRESENT(aux_density)) THEN
697 404 : CALL calc_stress_and_gradient_${kindd}$ (poisson_env, rhog, ehartree, rhog_aux, h_stress)
698 : ELSE
699 221263 : CALL calc_stress_and_gradient_${kindd}$ (poisson_env, rhog, ehartree, h_stress=h_stress)
700 : END IF
701 :
702 221667 : CALL pw_pool%give_back_pw(rhog)
703 221667 : IF (PRESENT(aux_density)) THEN
704 404 : CALL pw_pool%give_back_pw(rhog_aux)
705 : END IF
706 :
707 221667 : CALL timestop(handle)
708 :
709 221667 : END SUBROUTINE pw_poisson_solve_v_nodv_${kindd}$_${kindv}$
710 : #:endfor
711 : #:endfor
712 :
713 : #:for kindd in ['r3d_rs', 'c1d_gs']
714 : #:for kindg in ['r3d_rs', 'c1d_gs']
715 : ! **************************************************************************************************
716 : !> \brief Solve Poisson equation in a plane wave basis set
717 : !> Obtains electrostatic potential and its derivatives with respect to r
718 : !> from the density
719 : !> \param poisson_env ...
720 : !> \param density ...
721 : !> \param ehartree ...
722 : !> \param dvhartree ...
723 : !> \param h_stress ...
724 : !> \param rho_core ...
725 : !> \param greenfn ...
726 : !> \param aux_density Hartree energy and stress tensor between 2 different densities
727 : !> \par History
728 : !> JGH (13-Mar-2001) : completely revised
729 : !> \author apsi
730 : ! **************************************************************************************************
731 0 : SUBROUTINE pw_poisson_solve_nov_dv_${kindd}$_${kindg}$ (poisson_env, density, ehartree, &
732 : dvhartree, h_stress, rho_core, greenfn, aux_density)
733 :
734 : TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
735 : TYPE(pw_${kindd}$_type), INTENT(IN) :: density
736 : REAL(kind=dp), INTENT(out), OPTIONAL :: ehartree
737 : TYPE(pw_${kindg}$_type), DIMENSION(3) :: dvhartree
738 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT), &
739 : OPTIONAL :: h_stress
740 : TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_core, greenfn
741 : TYPE(pw_${kindd}$_type), INTENT(IN), OPTIONAL :: aux_density
742 :
743 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_solve'
744 :
745 : INTEGER :: handle
746 : LOGICAL :: has_dielectric
747 : TYPE(pw_grid_type), POINTER :: pw_grid
748 : TYPE(pw_pool_type), POINTER :: pw_pool
749 : TYPE(pw_r3d_rs_type) :: &
750 : rhor, vhartree_rs
751 : TYPE(pw_c1d_gs_type) :: influence_fn, rhog, rhog_aux, tmpg
752 :
753 0 : CALL timeset(routineN, handle)
754 :
755 0 : CALL pw_poisson_rebuild(poisson_env, density)
756 :
757 0 : has_dielectric = poisson_env%parameters%has_dielectric
758 :
759 : ! point pw
760 0 : pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
761 0 : pw_grid => pw_pool%pw_grid
762 : ! density in G space
763 0 : CALL pw_pool%create_pw(rhog)
764 0 : IF (PRESENT(aux_density)) THEN
765 0 : CALL pw_pool%create_pw(rhog_aux)
766 : END IF
767 :
768 0 : SELECT CASE (poisson_env%used_grid)
769 : CASE (use_gs_grid)
770 :
771 0 : SELECT CASE (poisson_env%green_fft%method)
772 : CASE (PERIODIC3D, ANALYTIC2D, ANALYTIC1D, ANALYTIC0D, MT2D, MT1D, MT0D, MULTIPOLE0D)
773 :
774 0 : CALL pw_transfer(density, rhog)
775 0 : IF (PRESENT(aux_density)) THEN
776 0 : CALL pw_transfer(aux_density, rhog_aux)
777 : END IF
778 0 : IF (PRESENT(ehartree)) THEN
779 0 : CALL pw_pool%create_pw(tmpg)
780 0 : CALL pw_copy(rhog, tmpg)
781 : END IF
782 0 : IF (PRESENT(greenfn)) THEN
783 0 : influence_fn = greenfn
784 : ELSE
785 0 : influence_fn = poisson_env%green_fft%influence_fn
786 : END IF
787 0 : CALL pw_multiply_with(rhog, influence_fn)
788 0 : IF (PRESENT(aux_density)) THEN
789 0 : CALL pw_multiply_with(rhog_aux, influence_fn)
790 0 : rhog_aux%array(:) = rhog_aux%array(:)*influence_fn%array(:)
791 : END IF
792 0 : IF (PRESENT(ehartree)) THEN
793 0 : IF (PRESENT(aux_density)) THEN
794 0 : ehartree = 0.5_dp*pw_integral_ab(rhog_aux, tmpg)
795 : ELSE
796 0 : ehartree = 0.5_dp*pw_integral_ab(rhog, tmpg)
797 : END IF
798 0 : CALL pw_pool%give_back_pw(tmpg)
799 : END IF
800 :
801 : CASE (PS_IMPLICIT)
802 0 : IF (PRESENT(h_stress)) THEN
803 0 : CPABORT("No stress tensor is implemented for the implicit Poisson solver.")
804 : END IF
805 :
806 0 : IF (has_dielectric .AND. PRESENT(rho_core)) THEN
807 0 : SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
808 : CASE (PERIODIC_BC, MIXED_PERIODIC_BC)
809 : CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
810 : poisson_env%diel_rs_grid, &
811 : poisson_env%pw_pools(poisson_env%pw_level)%pool, &
812 0 : density, rho_core=rho_core)
813 : CASE (NEUMANN_BC, MIXED_BC)
814 : CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
815 : poisson_env%diel_rs_grid, &
816 : poisson_env%pw_pools(poisson_env%pw_level)%pool, &
817 : poisson_env%dct_pw_grid, &
818 : poisson_env%parameters%ps_implicit_params%neumann_directions, &
819 : poisson_env%implicit_env%dct_env%recv_msgs_bnds, &
820 : poisson_env%implicit_env%dct_env%dests_expand, &
821 : poisson_env%implicit_env%dct_env%srcs_expand, &
822 : poisson_env%implicit_env%dct_env%flipg_stat, &
823 : poisson_env%implicit_env%dct_env%bounds_shftd, &
824 0 : density, rho_core=rho_core)
825 : END SELECT
826 : END IF
827 :
828 0 : CALL pw_pool%create_pw(rhor)
829 0 : CALL pw_pool%create_pw(vhartree_rs)
830 0 : CALL pw_transfer(density, rhor)
831 :
832 0 : SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
833 : CASE (PERIODIC_BC)
834 : CALL implicit_poisson_solver_periodic(poisson_env, rhor, vhartree_rs, &
835 0 : ehartree=ehartree)
836 : CASE (NEUMANN_BC)
837 : CALL implicit_poisson_solver_neumann(poisson_env, rhor, vhartree_rs, &
838 0 : ehartree=ehartree)
839 : CASE (MIXED_PERIODIC_BC)
840 : CALL implicit_poisson_solver_mixed_periodic(poisson_env, rhor, vhartree_rs, &
841 0 : electric_enthalpy=ehartree)
842 : CASE (MIXED_BC)
843 : CALL implicit_poisson_solver_mixed(poisson_env, rhor, vhartree_rs, &
844 0 : electric_enthalpy=ehartree)
845 : END SELECT
846 :
847 0 : CALL pw_transfer(rhor, rhog)
848 :
849 0 : IF (PRESENT(aux_density)) THEN
850 0 : CALL pw_transfer(aux_density, rhor)
851 0 : ehartree = 0.5_dp*pw_integral_ab(rhor, vhartree_rs)
852 : END IF
853 :
854 0 : CALL pw_pool%give_back_pw(rhor)
855 0 : CALL pw_pool%give_back_pw(vhartree_rs)
856 :
857 : CASE DEFAULT
858 : CALL cp_abort(__LOCATION__, &
859 : "unknown poisson method "// &
860 0 : cp_to_string(poisson_env%green_fft%method))
861 : END SELECT
862 :
863 : CASE (use_rs_grid)
864 :
865 0 : CALL pw_pool%create_pw(rhor)
866 0 : CALL pw_transfer(density, rhor)
867 0 : CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
868 0 : CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
869 0 : CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
870 0 : CALL pw_transfer(rhor, rhog)
871 0 : IF (PRESENT(ehartree)) THEN
872 : #! This is actually not consequent but to keep it working, I leave it that way
873 : #! Correctly, one checks the spaces but in CP2K, there is a separation in r-space/3D and g-space/1D in most cases
874 : #:if kindd=="r3d_rs"
875 0 : ehartree = 0.5_dp*pw_integral_ab(density, rhor)
876 : #:else
877 0 : ehartree = 0.5_dp*pw_integral_ab(density, rhog)
878 : #:endif
879 : END IF
880 0 : CALL pw_pool%give_back_pw(rhor)
881 :
882 : END SELECT
883 :
884 0 : IF (PRESENT(aux_density)) THEN
885 0 : CALL calc_stress_and_gradient_${kindg}$ (poisson_env, rhog, ehartree, rhog_aux, h_stress, dvhartree=dvhartree)
886 : ELSE
887 0 : CALL calc_stress_and_gradient_${kindg}$ (poisson_env, rhog, ehartree, h_stress=h_stress, dvhartree=dvhartree)
888 : END IF
889 :
890 0 : CALL pw_pool%give_back_pw(rhog)
891 0 : IF (PRESENT(aux_density)) THEN
892 0 : CALL pw_pool%give_back_pw(rhog_aux)
893 : END IF
894 :
895 0 : CALL timestop(handle)
896 :
897 0 : END SUBROUTINE pw_poisson_solve_nov_dv_${kindd}$_${kindg}$
898 : #:endfor
899 : #:endfor
900 :
901 : #:for kindd in ['r3d_rs', 'c1d_gs']
902 : #:for kindg in ['r3d_rs', 'c1d_gs']
903 : #:for kindv in ['r3d_rs', 'c1d_gs']
904 : ! **************************************************************************************************
905 : !> \brief Solve Poisson equation in a plane wave basis set
906 : !> Obtains electrostatic potential and its derivatives with respect to r
907 : !> from the density
908 : !> \param poisson_env ...
909 : !> \param density ...
910 : !> \param ehartree ...
911 : !> \param vhartree ...
912 : !> \param dvhartree ...
913 : !> \param h_stress ...
914 : !> \param rho_core ...
915 : !> \param greenfn ...
916 : !> \param aux_density Hartree energy and stress tensor between 2 different densities
917 : !> \par History
918 : !> JGH (13-Mar-2001) : completely revised
919 : !> \author apsi
920 : ! **************************************************************************************************
921 65041 : SUBROUTINE pw_poisson_solve_v_dv_${kindd}$_${kindv}$_${kindg}$ (poisson_env, density, ehartree, vhartree, &
922 : dvhartree, h_stress, rho_core, greenfn, aux_density)
923 :
924 : TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
925 : TYPE(pw_${kindd}$_type), INTENT(IN) :: density
926 : REAL(kind=dp), INTENT(out), OPTIONAL :: ehartree
927 : TYPE(pw_${kindv}$_type), INTENT(INOUT) :: vhartree
928 : TYPE(pw_${kindg}$_type), DIMENSION(3) :: dvhartree
929 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT), &
930 : OPTIONAL :: h_stress
931 : TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_core, greenfn
932 : TYPE(pw_${kindd}$_type), INTENT(IN), OPTIONAL :: aux_density
933 :
934 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_solve'
935 :
936 : INTEGER :: handle
937 : LOGICAL :: has_dielectric
938 : TYPE(pw_grid_type), POINTER :: pw_grid
939 : TYPE(pw_pool_type), POINTER :: pw_pool
940 : TYPE(pw_r3d_rs_type) :: rhor, vhartree_rs
941 : TYPE(pw_c1d_gs_type) :: influence_fn, rhog, rhog_aux
942 :
943 65041 : CALL timeset(routineN, handle)
944 :
945 65041 : CALL pw_poisson_rebuild(poisson_env, density)
946 :
947 65041 : has_dielectric = poisson_env%parameters%has_dielectric
948 :
949 : ! point pw
950 65041 : pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
951 65041 : pw_grid => pw_pool%pw_grid
952 65041 : IF (.NOT. pw_grid_compare(pw_pool%pw_grid, vhartree%pw_grid)) THEN
953 0 : CPABORT("vhartree has a different grid than the poisson solver")
954 : END IF
955 : ! density in G space
956 65041 : CALL pw_pool%create_pw(rhog)
957 65041 : IF (PRESENT(aux_density)) THEN
958 0 : CALL pw_pool%create_pw(rhog_aux)
959 : END IF
960 :
961 129614 : SELECT CASE (poisson_env%used_grid)
962 : CASE (use_gs_grid)
963 :
964 129614 : SELECT CASE (poisson_env%green_fft%method)
965 : CASE (PERIODIC3D, ANALYTIC2D, ANALYTIC1D, ANALYTIC0D, MT2D, MT1D, MT0D, MULTIPOLE0D)
966 :
967 64573 : CALL pw_transfer(density, rhog)
968 64573 : IF (PRESENT(aux_density)) THEN
969 0 : CALL pw_transfer(aux_density, rhog_aux)
970 : END IF
971 64573 : IF (PRESENT(greenfn)) THEN
972 0 : influence_fn = greenfn
973 : ELSE
974 64573 : influence_fn = poisson_env%green_fft%influence_fn
975 : END IF
976 64573 : CALL pw_multiply_with(rhog, influence_fn)
977 64573 : IF (PRESENT(aux_density)) THEN
978 0 : CALL pw_multiply_with(rhog_aux, influence_fn)
979 : END IF
980 64573 : CALL pw_transfer(rhog, vhartree)
981 64573 : IF (PRESENT(ehartree)) THEN
982 61632 : IF (PRESENT(aux_density)) THEN
983 : #:if kindd==kindv
984 0 : ehartree = 0.5_dp*pw_integral_ab(aux_density, vhartree)
985 : #:elif kindd=="c1d_gs"
986 0 : ehartree = 0.5_dp*pw_integral_ab(aux_density, rhog)
987 : #:else
988 0 : CALL pw_transfer(aux_density, rhog)
989 0 : ehartree = 0.5_dp*pw_integral_ab(rhog, vhartree)
990 : #:endif
991 : ELSE
992 : #:if kindd==kindv
993 61632 : ehartree = 0.5_dp*pw_integral_ab(density, vhartree)
994 : #:elif kindd=="c1d_gs"
995 0 : ehartree = 0.5_dp*pw_integral_ab(density, rhog)
996 : #:else
997 0 : CALL pw_transfer(density, rhog)
998 0 : ehartree = 0.5_dp*pw_integral_ab(rhog, vhartree)
999 : #:endif
1000 : END IF
1001 : END IF
1002 :
1003 : CASE (PS_IMPLICIT)
1004 0 : IF (PRESENT(h_stress)) THEN
1005 : CALL cp_abort(__LOCATION__, &
1006 0 : "No stress tensor is implemented for the implicit Poisson solver.")
1007 : END IF
1008 :
1009 0 : IF (has_dielectric .AND. PRESENT(rho_core)) THEN
1010 0 : SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
1011 : CASE (PERIODIC_BC, MIXED_PERIODIC_BC)
1012 : CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
1013 : poisson_env%diel_rs_grid, &
1014 : poisson_env%pw_pools(poisson_env%pw_level)%pool, &
1015 0 : density, rho_core=rho_core)
1016 : CASE (NEUMANN_BC, MIXED_BC)
1017 : CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
1018 : poisson_env%diel_rs_grid, &
1019 : poisson_env%pw_pools(poisson_env%pw_level)%pool, &
1020 : poisson_env%dct_pw_grid, &
1021 : poisson_env%parameters%ps_implicit_params%neumann_directions, &
1022 : poisson_env%implicit_env%dct_env%recv_msgs_bnds, &
1023 : poisson_env%implicit_env%dct_env%dests_expand, &
1024 : poisson_env%implicit_env%dct_env%srcs_expand, &
1025 : poisson_env%implicit_env%dct_env%flipg_stat, &
1026 : poisson_env%implicit_env%dct_env%bounds_shftd, &
1027 0 : density, rho_core=rho_core)
1028 : END SELECT
1029 : END IF
1030 :
1031 0 : CALL pw_pool%create_pw(rhor)
1032 0 : CALL pw_pool%create_pw(vhartree_rs)
1033 0 : CALL pw_transfer(density, rhor)
1034 :
1035 0 : SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
1036 : CASE (PERIODIC_BC)
1037 : CALL implicit_poisson_solver_periodic(poisson_env, rhor, vhartree_rs, &
1038 0 : ehartree=ehartree)
1039 : CASE (NEUMANN_BC)
1040 : CALL implicit_poisson_solver_neumann(poisson_env, rhor, vhartree_rs, &
1041 0 : ehartree=ehartree)
1042 : CASE (MIXED_PERIODIC_BC)
1043 : CALL implicit_poisson_solver_mixed_periodic(poisson_env, rhor, vhartree_rs, &
1044 0 : electric_enthalpy=ehartree)
1045 : CASE (MIXED_BC)
1046 : CALL implicit_poisson_solver_mixed(poisson_env, rhor, vhartree_rs, &
1047 0 : electric_enthalpy=ehartree)
1048 : END SELECT
1049 :
1050 0 : CALL pw_transfer(vhartree_rs, vhartree)
1051 0 : CALL pw_transfer(rhor, rhog)
1052 :
1053 0 : IF (PRESENT(aux_density)) THEN
1054 0 : CALL pw_transfer(aux_density, rhor)
1055 0 : ehartree = 0.5_dp*pw_integral_ab(rhor, vhartree_rs)
1056 : END IF
1057 :
1058 0 : CALL pw_pool%give_back_pw(rhor)
1059 0 : CALL pw_pool%give_back_pw(vhartree_rs)
1060 :
1061 : CASE DEFAULT
1062 : CALL cp_abort(__LOCATION__, &
1063 : "unknown poisson method "// &
1064 64573 : cp_to_string(poisson_env%green_fft%method))
1065 : END SELECT
1066 :
1067 : CASE (use_rs_grid)
1068 :
1069 468 : CALL pw_pool%create_pw(rhor)
1070 468 : CALL pw_transfer(density, rhor)
1071 468 : CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
1072 468 : CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
1073 468 : CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
1074 468 : CALL pw_transfer(rhor, vhartree)
1075 468 : CALL pw_transfer(rhor, rhog)
1076 468 : IF (PRESENT(ehartree)) THEN
1077 0 : IF (PRESENT(aux_density)) THEN
1078 : #:if kindd==kindv
1079 0 : ehartree = 0.5_dp*pw_integral_ab(aux_density, vhartree)
1080 : #:elif kindd=="r3d_rs"
1081 0 : ehartree = 0.5_dp*pw_integral_ab(aux_density, rhor)
1082 : #:else
1083 0 : ehartree = 0.5_dp*pw_integral_ab(aux_density, rhog)
1084 : #:endif
1085 : ELSE
1086 : #:if kindd==kindv
1087 0 : ehartree = 0.5_dp*pw_integral_ab(density, vhartree)
1088 : #:elif kindd=="r3d_rs"
1089 0 : ehartree = 0.5_dp*pw_integral_ab(density, rhor)
1090 : #:else
1091 0 : ehartree = 0.5_dp*pw_integral_ab(density, rhog)
1092 : #:endif
1093 : END IF
1094 : END IF
1095 468 : CALL pw_transfer(rhor, rhog)
1096 468 : IF (PRESENT(aux_density)) THEN
1097 0 : CALL pw_transfer(aux_density, rhor)
1098 0 : CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
1099 0 : CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
1100 0 : CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
1101 0 : CALL pw_transfer(rhor, rhog_aux)
1102 : END IF
1103 65509 : CALL pw_pool%give_back_pw(rhor)
1104 :
1105 : END SELECT
1106 :
1107 65041 : IF (PRESENT(aux_density)) THEN
1108 0 : CALL calc_stress_and_gradient_${kindg}$ (poisson_env, rhog, ehartree, rhog_aux, h_stress, dvhartree=dvhartree)
1109 : ELSE
1110 65041 : CALL calc_stress_and_gradient_${kindg}$ (poisson_env, rhog, ehartree, h_stress=h_stress, dvhartree=dvhartree)
1111 : END IF
1112 :
1113 65041 : CALL pw_pool%give_back_pw(rhog)
1114 65041 : IF (PRESENT(aux_density)) THEN
1115 0 : CALL pw_pool%give_back_pw(rhog_aux)
1116 : END IF
1117 :
1118 65041 : CALL timestop(handle)
1119 :
1120 65041 : END SUBROUTINE pw_poisson_solve_v_dv_${kindd}$_${kindv}$_${kindg}$
1121 : #:endfor
1122 : #:endfor
1123 : #:endfor
1124 :
1125 : #:for kind in ["c1d_gs", "r3d_rs"]
1126 286708 : SUBROUTINE calc_stress_and_gradient_${kind}$ (poisson_env, rhog, ehartree, rhog_aux, h_stress, dvhartree)
1127 : TYPE(pw_poisson_type), INTENT(IN) :: poisson_env
1128 : TYPE(pw_c1d_gs_type), INTENT(IN) :: rhog
1129 : REAL(KIND=dp) :: ehartree
1130 : TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rhog_aux
1131 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT), OPTIONAL :: h_stress
1132 : TYPE(pw_${kind}$_type), DIMENSION(3), INTENT(INOUT), OPTIONAL :: dvhartree
1133 :
1134 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_set'
1135 :
1136 : REAL(KIND=dp) :: ffa
1137 : INTEGER :: alpha, beta, n(3), handle, i
1138 2293664 : TYPE(pw_c1d_gs_type) :: dvg(3), dvg_aux(3)
1139 : TYPE(pw_pool_type), POINTER :: pw_pool
1140 :
1141 286708 : CALL timeset(routineN, handle)
1142 :
1143 286708 : pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
1144 :
1145 1146832 : DO i = 1, 3
1146 860124 : CALL pw_pool%create_pw(dvg(i))
1147 860124 : n = 0
1148 860124 : n(i) = 1
1149 860124 : CALL pw_copy(rhog, dvg(i))
1150 860124 : CALL pw_derive(dvg(i), n)
1151 1146832 : IF (PRESENT(rhog_aux)) THEN
1152 1212 : CALL pw_pool%create_pw(dvg_aux(i))
1153 1212 : CALL pw_copy(rhog_aux, dvg_aux(i))
1154 1212 : CALL pw_derive(dvg_aux(i), n)
1155 : END IF
1156 : END DO
1157 : ! save the derivatives
1158 286708 : IF (PRESENT(dvhartree)) THEN
1159 260164 : DO i = 1, 3
1160 260164 : CALL pw_transfer(dvg(i), dvhartree(i))
1161 : END DO
1162 : END IF
1163 : ! Calculate the contribution to the stress tensor this is only the contribution from
1164 : ! the Greens FUNCTION and the volume factor of the plane waves
1165 286708 : IF (PRESENT(h_stress)) THEN
1166 31724 : ffa = -1.0_dp/fourpi
1167 31724 : h_stress = 0.0_dp
1168 126896 : DO alpha = 1, 3
1169 95172 : h_stress(alpha, alpha) = ehartree
1170 126896 : IF (PRESENT(rhog_aux)) THEN
1171 3636 : DO beta = alpha, 3
1172 : h_stress(alpha, beta) = h_stress(alpha, beta) &
1173 2424 : + ffa*pw_integral_ab(dvg_aux(alpha), dvg(beta))
1174 3636 : h_stress(beta, alpha) = h_stress(alpha, beta)
1175 : END DO
1176 : ELSE
1177 281880 : DO beta = alpha, 3
1178 : h_stress(alpha, beta) = h_stress(alpha, beta) &
1179 187920 : + ffa*pw_integral_ab(dvg(alpha), dvg(beta))
1180 281880 : h_stress(beta, alpha) = h_stress(alpha, beta)
1181 : END DO
1182 : END IF
1183 : END DO
1184 :
1185 : ! Handle the periodicity cases for the Stress Tensor
1186 63436 : SELECT CASE (poisson_env%used_grid)
1187 : CASE (use_gs_grid)
1188 :
1189 : ! FFT based Poisson-Solver
1190 31740 : SELECT CASE (poisson_env%green_fft%method)
1191 : CASE (PERIODIC3D, PS_IMPLICIT)
1192 : ! Do Nothing
1193 : CASE (ANALYTIC2D, MT2D)
1194 : ! Zero the 1 non-periodic component. Only the periodic-subspace
1195 : ! stress components are defined for 2D systems.
1196 16 : alpha = poisson_env%green_fft%special_dimension
1197 64 : h_stress(:, alpha) = 0.0_dp
1198 64 : h_stress(alpha, :) = 0.0_dp
1199 : CASE (ANALYTIC1D, MT1D)
1200 : ! Zero the 2 non-periodic components
1201 0 : DO alpha = 1, 3
1202 0 : DO beta = alpha, 3
1203 0 : IF ((alpha /= poisson_env%green_fft%special_dimension) .OR. &
1204 0 : (beta /= poisson_env%green_fft%special_dimension)) THEN
1205 0 : h_stress(alpha, beta) = 0.0_dp
1206 0 : h_stress(beta, alpha) = 0.0_dp
1207 : END IF
1208 : END DO
1209 : END DO
1210 0 : CPABORT("Stress Tensor not tested for 1D systems.")
1211 : CASE (ANALYTIC0D, MT0D, MULTIPOLE0D)
1212 : ! Zero the full stress tensor
1213 302 : h_stress = 0.0_dp
1214 : CASE DEFAULT
1215 : CALL cp_abort(__LOCATION__, &
1216 : "unknown poisson method"// &
1217 31712 : cp_to_string(poisson_env%green_fft%method))
1218 : END SELECT
1219 :
1220 : CASE (use_rs_grid)
1221 :
1222 : ! Wavelet based Poisson-Solver
1223 31724 : SELECT CASE (poisson_env%wavelet%method)
1224 : CASE (WAVELET3D)
1225 : ! Do Nothing
1226 : CASE (WAVELET2D)
1227 : ! Zero the 1 non-periodic component
1228 0 : alpha = poisson_env%wavelet%special_dimension
1229 0 : h_stress(:, alpha) = 0.0_dp
1230 0 : h_stress(alpha, :) = 0.0_dp
1231 0 : CPABORT("Stress Tensor not tested for WAVELET 2D.")
1232 : CASE (WAVELET1D)
1233 : ! Zero the 2 non-periodic components
1234 0 : CPABORT("WAVELET 1D not implemented!")
1235 : CASE (WAVELET0D)
1236 : ! Zero the full stress tensor
1237 12 : h_stress = 0.0_dp
1238 : END SELECT
1239 :
1240 : END SELECT
1241 : END IF
1242 :
1243 1146832 : DO i = 1, 3
1244 860124 : CALL pw_pool%give_back_pw(dvg(i))
1245 1146832 : IF (PRESENT(rhog_aux)) THEN
1246 1212 : CALL pw_pool%give_back_pw(dvg_aux(i))
1247 : END IF
1248 : END DO
1249 :
1250 286708 : CALL timestop(handle)
1251 :
1252 286708 : END SUBROUTINE calc_stress_and_gradient_${kind}$
1253 : #:endfor
1254 :
1255 : ! **************************************************************************************************
1256 : !> \brief sets cell, grids and parameters used by the poisson solver
1257 : !> You should call this at least once (and set everything)
1258 : !> before using the poisson solver.
1259 : !> Smart, doesn't set the thing twice to the same value
1260 : !> Keeps track of the need to rebuild the poisson_env
1261 : !> \param poisson_env ...
1262 : !> \param cell_hmat ...
1263 : !> \param parameters ...
1264 : !> \param pw_pools ...
1265 : !> \param use_level ...
1266 : !> \param mt_super_ref_pw_grid ...
1267 : !> \param dct_pw_grid ...
1268 : !> \param force_rebuild ...
1269 : !> \author fawzi
1270 : !> \note
1271 : !> Checks everything at the end. This means that after *each* call to
1272 : !> this method the poisson env must be fully ready, so the first time
1273 : !> you have to set everything at once. Change this behaviour?
1274 : ! **************************************************************************************************
1275 26764 : SUBROUTINE pw_poisson_set(poisson_env, cell_hmat, parameters, pw_pools, use_level, &
1276 : mt_super_ref_pw_grid, dct_pw_grid, force_rebuild)
1277 :
1278 : TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
1279 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN), &
1280 : OPTIONAL :: cell_hmat
1281 : TYPE(pw_poisson_parameter_type), INTENT(IN), &
1282 : OPTIONAL :: parameters
1283 : TYPE(pw_pool_p_type), DIMENSION(:), OPTIONAL, &
1284 : POINTER :: pw_pools
1285 : INTEGER, INTENT(in), OPTIONAL :: use_level
1286 : TYPE(pw_grid_type), OPTIONAL, POINTER :: mt_super_ref_pw_grid, dct_pw_grid
1287 : LOGICAL, INTENT(in), OPTIONAL :: force_rebuild
1288 :
1289 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_set'
1290 :
1291 : INTEGER :: handle, i
1292 : LOGICAL :: same
1293 26764 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: tmp_pools
1294 :
1295 26764 : CALL timeset(routineN, handle)
1296 :
1297 26764 : IF (PRESENT(parameters)) poisson_env%parameters = parameters
1298 :
1299 26764 : IF (PRESENT(cell_hmat)) THEN
1300 63280 : IF (ANY(poisson_env%cell_hmat /= cell_hmat)) THEN
1301 25094 : CALL pw_poisson_cleanup(poisson_env)
1302 : END IF
1303 347932 : poisson_env%cell_hmat(:, :) = cell_hmat(:, :)
1304 26764 : poisson_env%rebuild = .TRUE.
1305 : END IF
1306 :
1307 26764 : IF (PRESENT(pw_pools)) THEN
1308 26764 : CPASSERT(ASSOCIATED(pw_pools))
1309 26764 : same = .FALSE.
1310 26764 : IF (ASSOCIATED(poisson_env%pw_pools)) THEN
1311 14496 : same = SIZE(poisson_env%pw_pools) == SIZE(pw_pools)
1312 14496 : IF (same) THEN
1313 31170 : DO i = 1, SIZE(pw_pools)
1314 16674 : IF (.NOT. ASSOCIATED(poisson_env%pw_pools(i)%pool, &
1315 19714 : pw_pools(i)%pool)) same = .FALSE.
1316 : END DO
1317 : END IF
1318 : END IF
1319 14496 : IF (.NOT. same) THEN
1320 15464 : poisson_env%rebuild = .TRUE.
1321 15464 : CALL pw_pools_copy(pw_pools, tmp_pools)
1322 15464 : CALL pw_pools_dealloc(poisson_env%pw_pools)
1323 15464 : poisson_env%pw_pools => tmp_pools
1324 : END IF
1325 : END IF
1326 :
1327 26764 : IF (PRESENT(use_level)) poisson_env%pw_level = use_level
1328 :
1329 26764 : IF (PRESENT(dct_pw_grid)) THEN
1330 13040 : IF (ASSOCIATED(dct_pw_grid)) THEN
1331 0 : CALL pw_grid_retain(dct_pw_grid)
1332 : END IF
1333 13040 : CALL pw_grid_release(poisson_env%dct_pw_grid)
1334 13040 : poisson_env%dct_pw_grid => dct_pw_grid
1335 : END IF
1336 :
1337 26764 : IF (PRESENT(mt_super_ref_pw_grid)) THEN
1338 13040 : IF (ASSOCIATED(mt_super_ref_pw_grid)) THEN
1339 1542 : CALL pw_grid_retain(mt_super_ref_pw_grid)
1340 : END IF
1341 13040 : CALL pw_grid_release(poisson_env%mt_super_ref_pw_grid)
1342 13040 : poisson_env%mt_super_ref_pw_grid => mt_super_ref_pw_grid
1343 : END IF
1344 :
1345 26764 : IF (PRESENT(force_rebuild)) THEN
1346 0 : IF (force_rebuild) poisson_env%rebuild = .TRUE.
1347 : END IF
1348 :
1349 26764 : CALL pw_poisson_check(poisson_env)
1350 :
1351 26764 : CALL timestop(handle)
1352 :
1353 26764 : END SUBROUTINE pw_poisson_set
1354 :
1355 : ! **************************************************************************************************
1356 : !> \brief computes the convolution f * u using reciprocal-space multiplication
1357 : !> \par
1358 : !> Given two real-space functions f and u, this routine evaluates their
1359 : !> convolution by transforming both to G space, multiplying their Fourier
1360 : !> components, and transforming the product back to real space.
1361 : !> \param poisson_env poisson environment
1362 : !> \param func input function f in real space
1363 : !> \param convolution convolution result (f * u) in real space
1364 : !> \param u input function/kernel u in real space
1365 : !> \author Ziwei Chai
1366 : !> \date 12.03.2026
1367 : ! **************************************************************************************************
1368 0 : SUBROUTINE pw_func_u_convolution(poisson_env, func, convolution, u)
1369 :
1370 : TYPE(pw_poisson_type), POINTER :: poisson_env
1371 : TYPE(pw_r3d_rs_type), INTENT(IN) :: func, u
1372 : TYPE(pw_r3d_rs_type), INTENT(INOUT), OPTIONAL :: convolution
1373 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_func_u_convolution'
1374 : INTEGER :: handle
1375 : TYPE(pw_grid_type), POINTER :: pw_grid
1376 : TYPE(pw_pool_type), POINTER :: pw_pool
1377 : TYPE(pw_c1d_gs_type) :: funcg, ug
1378 :
1379 0 : CALL timeset(routineN, handle)
1380 :
1381 0 : CPASSERT(ASSOCIATED(poisson_env))
1382 0 : CALL pw_poisson_rebuild(poisson_env, func)
1383 :
1384 : ! point pw
1385 0 : pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
1386 0 : pw_grid => pw_pool%pw_grid
1387 : ! func and u in G space
1388 0 : CALL pw_pool%create_pw(funcg)
1389 0 : CALL pw_pool%create_pw(ug)
1390 :
1391 0 : CALL pw_transfer(func, funcg)
1392 0 : CALL pw_transfer(u, ug)
1393 :
1394 0 : funcg%array(:) = funcg%array(:)*ug%array(:)
1395 :
1396 0 : CALL pw_transfer(funcg, convolution)
1397 :
1398 0 : CALL pw_pool%give_back_pw(funcg)
1399 0 : CALL pw_pool%give_back_pw(ug)
1400 :
1401 0 : CALL timestop(handle)
1402 :
1403 0 : END SUBROUTINE pw_func_u_convolution
1404 :
1405 : END MODULE pw_poisson_methods
|
