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