Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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
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 53870 : 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 53870 : NULLIFY (pw_pool)
109 53870 : IF (ASSOCIATED(poisson_env%pw_pools)) THEN
110 43958 : pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
111 : END IF
112 53870 : IF (ASSOCIATED(poisson_env%green_fft)) THEN
113 7688 : CALL pw_green_release(poisson_env%green_fft, pw_pool=pw_pool)
114 7688 : DEALLOCATE (poisson_env%green_fft)
115 : END IF
116 53870 : poisson_env%rebuild = .TRUE.
117 :
118 53870 : 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 19732 : 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 19732 : CPASSERT(ASSOCIATED(poisson_env%pw_pools))
133 39464 : CPASSERT(poisson_env%pw_level >= LBOUND(poisson_env%pw_pools, 1))
134 39464 : CPASSERT(poisson_env%pw_level <= UBOUND(poisson_env%pw_pools, 1))
135 19732 : green => poisson_env%green_fft
136 19732 : wavelet => poisson_env%wavelet
137 19732 : rebuild = poisson_env%rebuild
138 : rebuild = rebuild .OR. (poisson_env%method /= poisson_env%parameters%solver) &
139 19732 : .OR. .NOT. ASSOCIATED(green)
140 19732 : poisson_env%method = poisson_env%parameters%solver
141 :
142 19732 : IF (poisson_env%method == pw_poisson_wavelet) THEN
143 768 : poisson_env%used_grid = use_rs_grid
144 : ELSE
145 18964 : poisson_env%used_grid = use_gs_grid
146 : END IF
147 19732 : 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("")
170 : END SELECT
171 : END IF
172 19732 : IF (rebuild) THEN
173 19732 : poisson_env%rebuild = .TRUE.
174 19732 : CALL pw_poisson_cleanup(poisson_env)
175 : END IF
176 19732 : 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 49438 : 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 49438 : CALL timeset(routineN, handle)
194 :
195 49438 : CPASSERT(ASSOCIATED(poisson_env%pw_pools))
196 :
197 49438 : IF (poisson_env%rebuild) THEN
198 9536 : CALL pw_poisson_cleanup(poisson_env)
199 19072 : SELECT CASE (poisson_env%parameters%solver)
200 : CASE (pw_poisson_periodic, pw_poisson_analytic, pw_poisson_mt, pw_poisson_multipole)
201 9536 : 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 9536 : 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 9536 : CPABORT("Unknown Poisson solver")
223 : END SELECT
224 9536 : poisson_env%rebuild = .FALSE.
225 : END IF
226 :
227 49438 : CALL timestop(handle)
228 :
229 49438 : 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 226125 : 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 226125 : CALL timeset(routineN, handle)
249 :
250 226125 : CPASSERT(ASSOCIATED(poisson_env%pw_pools))
251 :
252 226125 : IF (poisson_env%rebuild) THEN
253 6394 : CALL pw_poisson_cleanup(poisson_env)
254 11982 : SELECT CASE (poisson_env%parameters%solver)
255 : CASE (pw_poisson_periodic, pw_poisson_analytic, pw_poisson_mt, pw_poisson_multipole)
256 5588 : 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 5588 : dct_pw_grid=poisson_env%dct_pw_grid)
262 : CASE (pw_poisson_wavelet)
263 752 : CPASSERT(ASSOCIATED(density%pw_grid))
264 : CALL ps_wavelet_create(poisson_env%parameters, poisson_env%wavelet, &
265 752 : 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 6394 : CPABORT("Unknown Poisson solver")
280 : END SELECT
281 6394 : poisson_env%rebuild = .FALSE.
282 : END IF
283 :
284 226125 : CALL timestop(handle)
285 :
286 226125 : 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 173600 : 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 173600 : CALL timeset(routineN, handle)
532 :
533 173600 : CALL pw_poisson_rebuild(poisson_env, density)
534 :
535 173600 : has_dielectric = poisson_env%parameters%has_dielectric
536 :
537 : ! point pw
538 173600 : pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
539 173600 : pw_grid => pw_pool%pw_grid
540 173600 : IF (.NOT. pw_grid_compare(pw_pool%pw_grid, vhartree%pw_grid)) &
541 0 : CPABORT("vhartree has a different grid than the poisson solver")
542 : ! density in G space
543 173600 : CALL pw_pool%create_pw(rhog)
544 173600 : IF (PRESENT(aux_density)) THEN
545 392 : CALL pw_pool%create_pw(rhog_aux)
546 : END IF
547 :
548 317579 : SELECT CASE (poisson_env%used_grid)
549 : CASE (use_gs_grid)
550 :
551 317131 : SELECT CASE (poisson_env%green_fft%method)
552 : CASE (PERIODIC3D, ANALYTIC2D, ANALYTIC1D, ANALYTIC0D, MT2D, MT1D, MT0D, MULTIPOLE0D)
553 :
554 143531 : CALL pw_transfer(density, rhog)
555 143531 : IF (PRESENT(aux_density)) THEN
556 392 : CALL pw_transfer(aux_density, rhog_aux)
557 : END IF
558 143531 : IF (PRESENT(greenfn)) THEN
559 182 : influence_fn = greenfn
560 : ELSE
561 143349 : influence_fn = poisson_env%green_fft%influence_fn
562 : END IF
563 143531 : CALL pw_multiply_with(rhog, influence_fn)
564 143531 : IF (PRESENT(aux_density)) THEN
565 392 : CALL pw_multiply_with(rhog_aux, influence_fn)
566 : END IF
567 143531 : CALL pw_transfer(rhog, vhartree)
568 143531 : IF (PRESENT(ehartree)) THEN
569 111511 : IF (PRESENT(aux_density)) THEN
570 : #:if kindd==kindv
571 392 : ehartree = 0.5_dp*pw_integral_ab(aux_density, vhartree)
572 : #:elif kindd=="c1d_gs"
573 0 : ehartree = 0.5_dp*pw_integral_ab(aux_density, rhog)
574 : #:else
575 0 : CALL pw_transfer(aux_density, rhog)
576 0 : ehartree = 0.5_dp*pw_integral_ab(rhog, vhartree)
577 : #:endif
578 : ELSE
579 : #:if kindd==kindv
580 111119 : ehartree = 0.5_dp*pw_integral_ab(density, vhartree)
581 : #:elif kindd=="c1d_gs"
582 0 : ehartree = 0.5_dp*pw_integral_ab(density, rhog)
583 : #:else
584 0 : CALL pw_transfer(density, rhog)
585 0 : ehartree = 0.5_dp*pw_integral_ab(rhog, vhartree)
586 : #:endif
587 : END IF
588 : END IF
589 :
590 : CASE (PS_IMPLICIT)
591 448 : IF (PRESENT(h_stress)) THEN
592 0 : CPABORT("No stress tensor is implemented for the implicit Poisson solver.")
593 : END IF
594 :
595 448 : IF (has_dielectric .AND. PRESENT(rho_core)) THEN
596 448 : SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
597 : CASE (PERIODIC_BC, MIXED_PERIODIC_BC)
598 : CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
599 : poisson_env%diel_rs_grid, &
600 : poisson_env%pw_pools(poisson_env%pw_level)%pool, &
601 288 : density, rho_core=rho_core)
602 : CASE (NEUMANN_BC, MIXED_BC)
603 : CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
604 : poisson_env%diel_rs_grid, &
605 : poisson_env%pw_pools(poisson_env%pw_level)%pool, &
606 : poisson_env%dct_pw_grid, &
607 : poisson_env%parameters%ps_implicit_params%neumann_directions, &
608 : poisson_env%implicit_env%dct_env%recv_msgs_bnds, &
609 : poisson_env%implicit_env%dct_env%dests_expand, &
610 : poisson_env%implicit_env%dct_env%srcs_expand, &
611 : poisson_env%implicit_env%dct_env%flipg_stat, &
612 : poisson_env%implicit_env%dct_env%bounds_shftd, &
613 160 : density, rho_core=rho_core)
614 : END SELECT
615 : END IF
616 :
617 448 : CALL pw_pool%create_pw(rhor)
618 448 : CALL pw_pool%create_pw(vhartree_rs)
619 448 : CALL pw_transfer(density, rhor)
620 :
621 448 : SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
622 : CASE (PERIODIC_BC)
623 : CALL implicit_poisson_solver_periodic(poisson_env, rhor, vhartree_rs, &
624 104 : ehartree=ehartree)
625 : CASE (NEUMANN_BC)
626 : CALL implicit_poisson_solver_neumann(poisson_env, rhor, vhartree_rs, &
627 24 : ehartree=ehartree)
628 : CASE (MIXED_PERIODIC_BC)
629 : CALL implicit_poisson_solver_mixed_periodic(poisson_env, rhor, vhartree_rs, &
630 184 : electric_enthalpy=ehartree)
631 : CASE (MIXED_BC)
632 : CALL implicit_poisson_solver_mixed(poisson_env, rhor, vhartree_rs, &
633 136 : electric_enthalpy=ehartree)
634 : END SELECT
635 :
636 448 : IF (PRESENT(aux_density)) THEN
637 0 : CALL pw_transfer(aux_density, rhor)
638 0 : ehartree = 0.5_dp*pw_integral_ab(rhor, vhartree_rs)
639 : END IF
640 :
641 448 : CALL pw_transfer(vhartree_rs, vhartree)
642 :
643 448 : CALL pw_pool%give_back_pw(rhor)
644 448 : CALL pw_pool%give_back_pw(vhartree_rs)
645 :
646 : CASE DEFAULT
647 : CALL cp_abort(__LOCATION__, &
648 : "unknown poisson method "// &
649 143979 : cp_to_string(poisson_env%green_fft%method))
650 : END SELECT
651 :
652 : CASE (use_rs_grid)
653 :
654 29621 : CALL pw_pool%create_pw(rhor)
655 29621 : CALL pw_transfer(density, rhor)
656 29621 : CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
657 29621 : CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
658 29621 : CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
659 29621 : CALL pw_transfer(rhor, vhartree)
660 29621 : IF (PRESENT(ehartree)) THEN
661 7490 : IF (PRESENT(aux_density)) THEN
662 : #:if kindd==kindv
663 0 : ehartree = 0.5_dp*pw_integral_ab(aux_density, vhartree)
664 : #:elif kindd=="r3d_rs"
665 0 : ehartree = 0.5_dp*pw_integral_ab(aux_density, rhor)
666 : #:else
667 0 : CALL pw_transfer(vhartree, rhog)
668 0 : ehartree = 0.5_dp*pw_integral_ab(aux_density, rhog)
669 : #:endif
670 : ELSE
671 : #:if kindd==kindv
672 7490 : ehartree = 0.5_dp*pw_integral_ab(density, vhartree)
673 : #:elif kindd=="r3d_rs"
674 0 : ehartree = 0.5_dp*pw_integral_ab(density, rhor)
675 : #:else
676 0 : CALL pw_transfer(vhartree, rhog)
677 0 : ehartree = 0.5_dp*pw_integral_ab(density, rhog)
678 : #:endif
679 : END IF
680 : END IF
681 29621 : IF (PRESENT(h_stress)) THEN
682 12 : CALL pw_transfer(rhor, rhog)
683 12 : IF (PRESENT(aux_density)) THEN
684 0 : CALL pw_transfer(aux_density, rhor)
685 0 : CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
686 0 : CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
687 0 : CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
688 0 : CALL pw_transfer(rhor, rhog_aux)
689 : END IF
690 : END IF
691 203221 : CALL pw_pool%give_back_pw(rhor)
692 :
693 : END SELECT
694 :
695 173600 : IF (PRESENT(aux_density)) THEN
696 392 : CALL calc_stress_and_gradient_${kindd}$ (poisson_env, rhog, ehartree, rhog_aux, h_stress)
697 : ELSE
698 173208 : CALL calc_stress_and_gradient_${kindd}$ (poisson_env, rhog, ehartree, h_stress=h_stress)
699 : END IF
700 :
701 173600 : CALL pw_pool%give_back_pw(rhog)
702 173600 : IF (PRESENT(aux_density)) THEN
703 392 : CALL pw_pool%give_back_pw(rhog_aux)
704 : END IF
705 :
706 173600 : CALL timestop(handle)
707 :
708 173600 : END SUBROUTINE pw_poisson_solve_v_nodv_${kindd}$_${kindv}$
709 : #:endfor
710 : #:endfor
711 :
712 : #:for kindd in ['r3d_rs', 'c1d_gs']
713 : #:for kindg in ['r3d_rs', 'c1d_gs']
714 : ! **************************************************************************************************
715 : !> \brief Solve Poisson equation in a plane wave basis set
716 : !> Obtains electrostatic potential and its derivatives with respect to r
717 : !> from the density
718 : !> \param poisson_env ...
719 : !> \param density ...
720 : !> \param ehartree ...
721 : !> \param dvhartree ...
722 : !> \param h_stress ...
723 : !> \param rho_core ...
724 : !> \param greenfn ...
725 : !> \param aux_density Hartree energy and stress tensor between 2 different densities
726 : !> \par History
727 : !> JGH (13-Mar-2001) : completely revised
728 : !> \author apsi
729 : ! **************************************************************************************************
730 0 : SUBROUTINE pw_poisson_solve_nov_dv_${kindd}$_${kindg}$ (poisson_env, density, ehartree, &
731 : dvhartree, h_stress, rho_core, greenfn, aux_density)
732 :
733 : TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
734 : TYPE(pw_${kindd}$_type), INTENT(IN) :: density
735 : REAL(kind=dp), INTENT(out), OPTIONAL :: ehartree
736 : TYPE(pw_${kindg}$_type), DIMENSION(3) :: dvhartree
737 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT), &
738 : OPTIONAL :: h_stress
739 : TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_core, greenfn
740 : TYPE(pw_${kindd}$_type), INTENT(IN), OPTIONAL :: aux_density
741 :
742 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_solve'
743 :
744 : INTEGER :: handle
745 : LOGICAL :: has_dielectric
746 : TYPE(pw_grid_type), POINTER :: pw_grid
747 : TYPE(pw_pool_type), POINTER :: pw_pool
748 : TYPE(pw_r3d_rs_type) :: &
749 : rhor, vhartree_rs
750 : TYPE(pw_c1d_gs_type) :: influence_fn, rhog, rhog_aux, tmpg
751 :
752 0 : CALL timeset(routineN, handle)
753 :
754 0 : CALL pw_poisson_rebuild(poisson_env, density)
755 :
756 0 : has_dielectric = poisson_env%parameters%has_dielectric
757 :
758 : ! point pw
759 0 : pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
760 0 : pw_grid => pw_pool%pw_grid
761 : ! density in G space
762 0 : CALL pw_pool%create_pw(rhog)
763 0 : IF (PRESENT(aux_density)) THEN
764 0 : CALL pw_pool%create_pw(rhog_aux)
765 : END IF
766 :
767 0 : SELECT CASE (poisson_env%used_grid)
768 : CASE (use_gs_grid)
769 :
770 0 : SELECT CASE (poisson_env%green_fft%method)
771 : CASE (PERIODIC3D, ANALYTIC2D, ANALYTIC1D, ANALYTIC0D, MT2D, MT1D, MT0D, MULTIPOLE0D)
772 :
773 0 : CALL pw_transfer(density, rhog)
774 0 : IF (PRESENT(aux_density)) THEN
775 0 : CALL pw_transfer(aux_density, rhog_aux)
776 : END IF
777 0 : IF (PRESENT(ehartree)) THEN
778 0 : CALL pw_pool%create_pw(tmpg)
779 0 : CALL pw_copy(rhog, tmpg)
780 : END IF
781 0 : IF (PRESENT(greenfn)) THEN
782 0 : influence_fn = greenfn
783 : ELSE
784 0 : influence_fn = poisson_env%green_fft%influence_fn
785 : END IF
786 0 : CALL pw_multiply_with(rhog, influence_fn)
787 0 : IF (PRESENT(aux_density)) THEN
788 0 : CALL pw_multiply_with(rhog_aux, influence_fn)
789 0 : rhog_aux%array(:) = rhog_aux%array(:)*influence_fn%array(:)
790 : END IF
791 0 : IF (PRESENT(ehartree)) THEN
792 0 : IF (PRESENT(aux_density)) THEN
793 0 : ehartree = 0.5_dp*pw_integral_ab(rhog_aux, tmpg)
794 : ELSE
795 0 : ehartree = 0.5_dp*pw_integral_ab(rhog, tmpg)
796 : END IF
797 0 : CALL pw_pool%give_back_pw(tmpg)
798 : END IF
799 :
800 : CASE (PS_IMPLICIT)
801 0 : IF (PRESENT(h_stress)) THEN
802 0 : CPABORT("No stress tensor is implemented for the implicit Poisson solver.")
803 : END IF
804 :
805 0 : IF (has_dielectric .AND. PRESENT(rho_core)) THEN
806 0 : SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
807 : CASE (PERIODIC_BC, MIXED_PERIODIC_BC)
808 : CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
809 : poisson_env%diel_rs_grid, &
810 : poisson_env%pw_pools(poisson_env%pw_level)%pool, &
811 0 : density, rho_core=rho_core)
812 : CASE (NEUMANN_BC, MIXED_BC)
813 : CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
814 : poisson_env%diel_rs_grid, &
815 : poisson_env%pw_pools(poisson_env%pw_level)%pool, &
816 : poisson_env%dct_pw_grid, &
817 : poisson_env%parameters%ps_implicit_params%neumann_directions, &
818 : poisson_env%implicit_env%dct_env%recv_msgs_bnds, &
819 : poisson_env%implicit_env%dct_env%dests_expand, &
820 : poisson_env%implicit_env%dct_env%srcs_expand, &
821 : poisson_env%implicit_env%dct_env%flipg_stat, &
822 : poisson_env%implicit_env%dct_env%bounds_shftd, &
823 0 : density, rho_core=rho_core)
824 : END SELECT
825 : END IF
826 :
827 0 : CALL pw_pool%create_pw(rhor)
828 0 : CALL pw_pool%create_pw(vhartree_rs)
829 0 : CALL pw_transfer(density, rhor)
830 :
831 0 : SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
832 : CASE (PERIODIC_BC)
833 : CALL implicit_poisson_solver_periodic(poisson_env, rhor, vhartree_rs, &
834 0 : ehartree=ehartree)
835 : CASE (NEUMANN_BC)
836 : CALL implicit_poisson_solver_neumann(poisson_env, rhor, vhartree_rs, &
837 0 : ehartree=ehartree)
838 : CASE (MIXED_PERIODIC_BC)
839 : CALL implicit_poisson_solver_mixed_periodic(poisson_env, rhor, vhartree_rs, &
840 0 : electric_enthalpy=ehartree)
841 : CASE (MIXED_BC)
842 : CALL implicit_poisson_solver_mixed(poisson_env, rhor, vhartree_rs, &
843 0 : electric_enthalpy=ehartree)
844 : END SELECT
845 :
846 0 : CALL pw_transfer(rhor, rhog)
847 :
848 0 : IF (PRESENT(aux_density)) THEN
849 0 : CALL pw_transfer(aux_density, rhor)
850 0 : ehartree = 0.5_dp*pw_integral_ab(rhor, vhartree_rs)
851 : END IF
852 :
853 0 : CALL pw_pool%give_back_pw(rhor)
854 0 : CALL pw_pool%give_back_pw(vhartree_rs)
855 :
856 : CASE DEFAULT
857 : CALL cp_abort(__LOCATION__, &
858 : "unknown poisson method "// &
859 0 : cp_to_string(poisson_env%green_fft%method))
860 : END SELECT
861 :
862 : CASE (use_rs_grid)
863 :
864 0 : CALL pw_pool%create_pw(rhor)
865 0 : CALL pw_transfer(density, rhor)
866 0 : CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
867 0 : CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
868 0 : CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
869 0 : CALL pw_transfer(rhor, rhog)
870 0 : IF (PRESENT(ehartree)) THEN
871 : #! This is actually not consequent but to keep it working, I leave it that way
872 : #! Correctly, one checks the spaces but in CP2K, there is a separation in r-space/3D and g-space/1D in most cases
873 : #:if kindd=="r3d_rs"
874 0 : ehartree = 0.5_dp*pw_integral_ab(density, rhor)
875 : #:else
876 0 : ehartree = 0.5_dp*pw_integral_ab(density, rhog)
877 : #:endif
878 : END IF
879 0 : CALL pw_pool%give_back_pw(rhor)
880 :
881 : END SELECT
882 :
883 0 : IF (PRESENT(aux_density)) THEN
884 0 : CALL calc_stress_and_gradient_${kindg}$ (poisson_env, rhog, ehartree, rhog_aux, h_stress, dvhartree=dvhartree)
885 : ELSE
886 0 : CALL calc_stress_and_gradient_${kindg}$ (poisson_env, rhog, ehartree, h_stress=h_stress, dvhartree=dvhartree)
887 : END IF
888 :
889 0 : CALL pw_pool%give_back_pw(rhog)
890 0 : IF (PRESENT(aux_density)) THEN
891 0 : CALL pw_pool%give_back_pw(rhog_aux)
892 : END IF
893 :
894 0 : CALL timestop(handle)
895 :
896 0 : END SUBROUTINE pw_poisson_solve_nov_dv_${kindd}$_${kindg}$
897 : #:endfor
898 : #:endfor
899 :
900 : #:for kindd in ['r3d_rs', 'c1d_gs']
901 : #:for kindg in ['r3d_rs', 'c1d_gs']
902 : #:for kindv in ['r3d_rs', 'c1d_gs']
903 : ! **************************************************************************************************
904 : !> \brief Solve Poisson equation in a plane wave basis set
905 : !> Obtains electrostatic potential and its derivatives with respect to r
906 : !> from the density
907 : !> \param poisson_env ...
908 : !> \param density ...
909 : !> \param ehartree ...
910 : !> \param vhartree ...
911 : !> \param dvhartree ...
912 : !> \param h_stress ...
913 : !> \param rho_core ...
914 : !> \param greenfn ...
915 : !> \param aux_density Hartree energy and stress tensor between 2 different densities
916 : !> \par History
917 : !> JGH (13-Mar-2001) : completely revised
918 : !> \author apsi
919 : ! **************************************************************************************************
920 52525 : SUBROUTINE pw_poisson_solve_v_dv_${kindd}$_${kindv}$_${kindg}$ (poisson_env, density, ehartree, vhartree, &
921 : dvhartree, h_stress, rho_core, greenfn, aux_density)
922 :
923 : TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
924 : TYPE(pw_${kindd}$_type), INTENT(IN) :: density
925 : REAL(kind=dp), INTENT(out), OPTIONAL :: ehartree
926 : TYPE(pw_${kindv}$_type), INTENT(INOUT) :: vhartree
927 : TYPE(pw_${kindg}$_type), DIMENSION(3) :: dvhartree
928 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT), &
929 : OPTIONAL :: h_stress
930 : TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_core, greenfn
931 : TYPE(pw_${kindd}$_type), INTENT(IN), OPTIONAL :: aux_density
932 :
933 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_solve'
934 :
935 : INTEGER :: handle
936 : LOGICAL :: has_dielectric
937 : TYPE(pw_grid_type), POINTER :: pw_grid
938 : TYPE(pw_pool_type), POINTER :: pw_pool
939 : TYPE(pw_r3d_rs_type) :: rhor, vhartree_rs
940 : TYPE(pw_c1d_gs_type) :: influence_fn, rhog, rhog_aux
941 :
942 52525 : CALL timeset(routineN, handle)
943 :
944 52525 : CALL pw_poisson_rebuild(poisson_env, density)
945 :
946 52525 : has_dielectric = poisson_env%parameters%has_dielectric
947 :
948 : ! point pw
949 52525 : pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
950 52525 : pw_grid => pw_pool%pw_grid
951 52525 : IF (.NOT. pw_grid_compare(pw_pool%pw_grid, vhartree%pw_grid)) &
952 0 : CPABORT("vhartree has a different grid than the poisson solver")
953 : ! density in G space
954 52525 : CALL pw_pool%create_pw(rhog)
955 52525 : IF (PRESENT(aux_density)) THEN
956 0 : CALL pw_pool%create_pw(rhog_aux)
957 : END IF
958 :
959 104582 : SELECT CASE (poisson_env%used_grid)
960 : CASE (use_gs_grid)
961 :
962 104582 : SELECT CASE (poisson_env%green_fft%method)
963 : CASE (PERIODIC3D, ANALYTIC2D, ANALYTIC1D, ANALYTIC0D, MT2D, MT1D, MT0D, MULTIPOLE0D)
964 :
965 52057 : CALL pw_transfer(density, rhog)
966 52057 : IF (PRESENT(aux_density)) THEN
967 0 : CALL pw_transfer(aux_density, rhog_aux)
968 : END IF
969 52057 : IF (PRESENT(greenfn)) THEN
970 0 : influence_fn = greenfn
971 : ELSE
972 52057 : influence_fn = poisson_env%green_fft%influence_fn
973 : END IF
974 52057 : CALL pw_multiply_with(rhog, influence_fn)
975 52057 : IF (PRESENT(aux_density)) THEN
976 0 : CALL pw_multiply_with(rhog_aux, influence_fn)
977 : END IF
978 52057 : CALL pw_transfer(rhog, vhartree)
979 52057 : IF (PRESENT(ehartree)) THEN
980 49802 : IF (PRESENT(aux_density)) THEN
981 : #:if kindd==kindv
982 0 : ehartree = 0.5_dp*pw_integral_ab(aux_density, vhartree)
983 : #:elif kindd=="c1d_gs"
984 0 : ehartree = 0.5_dp*pw_integral_ab(aux_density, rhog)
985 : #:else
986 0 : CALL pw_transfer(aux_density, rhog)
987 0 : ehartree = 0.5_dp*pw_integral_ab(rhog, vhartree)
988 : #:endif
989 : ELSE
990 : #:if kindd==kindv
991 49802 : ehartree = 0.5_dp*pw_integral_ab(density, vhartree)
992 : #:elif kindd=="c1d_gs"
993 0 : ehartree = 0.5_dp*pw_integral_ab(density, rhog)
994 : #:else
995 0 : CALL pw_transfer(density, rhog)
996 0 : ehartree = 0.5_dp*pw_integral_ab(rhog, vhartree)
997 : #:endif
998 : END IF
999 : END IF
1000 :
1001 : CASE (PS_IMPLICIT)
1002 0 : IF (PRESENT(h_stress)) THEN
1003 : CALL cp_abort(__LOCATION__, &
1004 0 : "No stress tensor is implemented for the implicit Poisson solver.")
1005 : END IF
1006 :
1007 0 : IF (has_dielectric .AND. PRESENT(rho_core)) THEN
1008 0 : SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
1009 : CASE (PERIODIC_BC, MIXED_PERIODIC_BC)
1010 : CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
1011 : poisson_env%diel_rs_grid, &
1012 : poisson_env%pw_pools(poisson_env%pw_level)%pool, &
1013 0 : density, rho_core=rho_core)
1014 : CASE (NEUMANN_BC, MIXED_BC)
1015 : CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
1016 : poisson_env%diel_rs_grid, &
1017 : poisson_env%pw_pools(poisson_env%pw_level)%pool, &
1018 : poisson_env%dct_pw_grid, &
1019 : poisson_env%parameters%ps_implicit_params%neumann_directions, &
1020 : poisson_env%implicit_env%dct_env%recv_msgs_bnds, &
1021 : poisson_env%implicit_env%dct_env%dests_expand, &
1022 : poisson_env%implicit_env%dct_env%srcs_expand, &
1023 : poisson_env%implicit_env%dct_env%flipg_stat, &
1024 : poisson_env%implicit_env%dct_env%bounds_shftd, &
1025 0 : density, rho_core=rho_core)
1026 : END SELECT
1027 : END IF
1028 :
1029 0 : CALL pw_pool%create_pw(rhor)
1030 0 : CALL pw_pool%create_pw(vhartree_rs)
1031 0 : CALL pw_transfer(density, rhor)
1032 :
1033 0 : SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
1034 : CASE (PERIODIC_BC)
1035 : CALL implicit_poisson_solver_periodic(poisson_env, rhor, vhartree_rs, &
1036 0 : ehartree=ehartree)
1037 : CASE (NEUMANN_BC)
1038 : CALL implicit_poisson_solver_neumann(poisson_env, rhor, vhartree_rs, &
1039 0 : ehartree=ehartree)
1040 : CASE (MIXED_PERIODIC_BC)
1041 : CALL implicit_poisson_solver_mixed_periodic(poisson_env, rhor, vhartree_rs, &
1042 0 : electric_enthalpy=ehartree)
1043 : CASE (MIXED_BC)
1044 : CALL implicit_poisson_solver_mixed(poisson_env, rhor, vhartree_rs, &
1045 0 : electric_enthalpy=ehartree)
1046 : END SELECT
1047 :
1048 0 : CALL pw_transfer(vhartree_rs, vhartree)
1049 0 : CALL pw_transfer(rhor, rhog)
1050 :
1051 0 : IF (PRESENT(aux_density)) THEN
1052 0 : CALL pw_transfer(aux_density, rhor)
1053 0 : ehartree = 0.5_dp*pw_integral_ab(rhor, vhartree_rs)
1054 : END IF
1055 :
1056 0 : CALL pw_pool%give_back_pw(rhor)
1057 0 : CALL pw_pool%give_back_pw(vhartree_rs)
1058 :
1059 : CASE DEFAULT
1060 : CALL cp_abort(__LOCATION__, &
1061 : "unknown poisson method "// &
1062 52057 : cp_to_string(poisson_env%green_fft%method))
1063 : END SELECT
1064 :
1065 : CASE (use_rs_grid)
1066 :
1067 468 : CALL pw_pool%create_pw(rhor)
1068 468 : CALL pw_transfer(density, rhor)
1069 468 : CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
1070 468 : CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
1071 468 : CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
1072 468 : CALL pw_transfer(rhor, vhartree)
1073 468 : CALL pw_transfer(rhor, rhog)
1074 468 : IF (PRESENT(ehartree)) THEN
1075 0 : IF (PRESENT(aux_density)) THEN
1076 : #:if kindd==kindv
1077 0 : ehartree = 0.5_dp*pw_integral_ab(aux_density, vhartree)
1078 : #:elif kindd=="r3d_rs"
1079 0 : ehartree = 0.5_dp*pw_integral_ab(aux_density, rhor)
1080 : #:else
1081 0 : ehartree = 0.5_dp*pw_integral_ab(aux_density, rhog)
1082 : #:endif
1083 : ELSE
1084 : #:if kindd==kindv
1085 0 : ehartree = 0.5_dp*pw_integral_ab(density, vhartree)
1086 : #:elif kindd=="r3d_rs"
1087 0 : ehartree = 0.5_dp*pw_integral_ab(density, rhor)
1088 : #:else
1089 0 : ehartree = 0.5_dp*pw_integral_ab(density, rhog)
1090 : #:endif
1091 : END IF
1092 : END IF
1093 468 : CALL pw_transfer(rhor, rhog)
1094 468 : IF (PRESENT(aux_density)) THEN
1095 0 : CALL pw_transfer(aux_density, rhor)
1096 0 : CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
1097 0 : CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
1098 0 : CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
1099 0 : CALL pw_transfer(rhor, rhog_aux)
1100 : END IF
1101 52993 : CALL pw_pool%give_back_pw(rhor)
1102 :
1103 : END SELECT
1104 :
1105 52525 : IF (PRESENT(aux_density)) THEN
1106 0 : CALL calc_stress_and_gradient_${kindg}$ (poisson_env, rhog, ehartree, rhog_aux, h_stress, dvhartree=dvhartree)
1107 : ELSE
1108 52525 : CALL calc_stress_and_gradient_${kindg}$ (poisson_env, rhog, ehartree, h_stress=h_stress, dvhartree=dvhartree)
1109 : END IF
1110 :
1111 52525 : CALL pw_pool%give_back_pw(rhog)
1112 52525 : IF (PRESENT(aux_density)) THEN
1113 0 : CALL pw_pool%give_back_pw(rhog_aux)
1114 : END IF
1115 :
1116 52525 : CALL timestop(handle)
1117 :
1118 52525 : END SUBROUTINE pw_poisson_solve_v_dv_${kindd}$_${kindv}$_${kindg}$
1119 : #:endfor
1120 : #:endfor
1121 : #:endfor
1122 :
1123 : #:for kind in ["c1d_gs", "r3d_rs"]
1124 226125 : SUBROUTINE calc_stress_and_gradient_${kind}$ (poisson_env, rhog, ehartree, rhog_aux, h_stress, dvhartree)
1125 : TYPE(pw_poisson_type), INTENT(IN) :: poisson_env
1126 : TYPE(pw_c1d_gs_type), INTENT(IN) :: rhog
1127 : REAL(KIND=dp) :: ehartree
1128 : TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rhog_aux
1129 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT), OPTIONAL :: h_stress
1130 : TYPE(pw_${kind}$_type), DIMENSION(3), INTENT(INOUT), OPTIONAL :: dvhartree
1131 :
1132 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_set'
1133 :
1134 : REAL(KIND=dp) :: ffa
1135 : INTEGER :: alpha, beta, n(3), handle, i
1136 1809000 : TYPE(pw_c1d_gs_type) :: dvg(3), dvg_aux(3)
1137 : TYPE(pw_pool_type), POINTER :: pw_pool
1138 :
1139 226125 : CALL timeset(routineN, handle)
1140 :
1141 226125 : pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
1142 :
1143 904500 : DO i = 1, 3
1144 678375 : CALL pw_pool%create_pw(dvg(i))
1145 678375 : n = 0
1146 678375 : n(i) = 1
1147 678375 : CALL pw_copy(rhog, dvg(i))
1148 678375 : CALL pw_derive(dvg(i), n)
1149 904500 : IF (PRESENT(rhog_aux)) THEN
1150 1176 : CALL pw_pool%create_pw(dvg_aux(i))
1151 1176 : CALL pw_copy(rhog_aux, dvg_aux(i))
1152 1176 : CALL pw_derive(dvg_aux(i), n)
1153 : END IF
1154 : END DO
1155 : ! save the derivatives
1156 226125 : IF (PRESENT(dvhartree)) THEN
1157 210100 : DO i = 1, 3
1158 210100 : CALL pw_transfer(dvg(i), dvhartree(i))
1159 : END DO
1160 : END IF
1161 : ! Calculate the contribution to the stress tensor this is only the contribution from
1162 : ! the Greens FUNCTION and the volume factor of the plane waves
1163 226125 : IF (PRESENT(h_stress)) THEN
1164 32332 : ffa = -1.0_dp/fourpi
1165 32332 : h_stress = 0.0_dp
1166 129328 : DO alpha = 1, 3
1167 96996 : h_stress(alpha, alpha) = ehartree
1168 129328 : IF (PRESENT(rhog_aux)) THEN
1169 3528 : DO beta = alpha, 3
1170 : h_stress(alpha, beta) = h_stress(alpha, beta) &
1171 2352 : + ffa*pw_integral_ab(dvg_aux(alpha), dvg(beta))
1172 3528 : h_stress(beta, alpha) = h_stress(alpha, beta)
1173 : END DO
1174 : ELSE
1175 287460 : DO beta = alpha, 3
1176 : h_stress(alpha, beta) = h_stress(alpha, beta) &
1177 191640 : + ffa*pw_integral_ab(dvg(alpha), dvg(beta))
1178 287460 : h_stress(beta, alpha) = h_stress(alpha, beta)
1179 : END DO
1180 : END IF
1181 : END DO
1182 :
1183 : ! Handle the periodicity cases for the Stress Tensor
1184 64652 : SELECT CASE (poisson_env%used_grid)
1185 : CASE (use_gs_grid)
1186 :
1187 : ! FFT based Poisson-Solver
1188 32332 : SELECT CASE (poisson_env%green_fft%method)
1189 : CASE (PERIODIC3D, PS_IMPLICIT)
1190 : ! Do Nothing
1191 : CASE (ANALYTIC2D, MT2D)
1192 : ! Zero the 1 non-periodic component
1193 0 : alpha = poisson_env%green_fft%special_dimension
1194 0 : h_stress(:, alpha) = 0.0_dp
1195 0 : h_stress(alpha, :) = 0.0_dp
1196 0 : CPABORT("Stress Tensor not tested for 2D systems.")
1197 : CASE (ANALYTIC1D, MT1D)
1198 : ! Zero the 2 non-periodic components
1199 0 : DO alpha = 1, 3
1200 0 : DO beta = alpha, 3
1201 0 : IF ((alpha /= poisson_env%green_fft%special_dimension) .OR. &
1202 0 : (beta /= poisson_env%green_fft%special_dimension)) THEN
1203 0 : h_stress(alpha, beta) = 0.0_dp
1204 0 : h_stress(beta, alpha) = 0.0_dp
1205 : END IF
1206 : END DO
1207 : END DO
1208 0 : CPABORT("Stress Tensor not tested for 1D systems.")
1209 : CASE (ANALYTIC0D, MT0D, MULTIPOLE0D)
1210 : ! Zero the full stress tensor
1211 302 : h_stress = 0.0_dp
1212 : CASE DEFAULT
1213 : CALL cp_abort(__LOCATION__, &
1214 : "unknown poisson method"// &
1215 32320 : cp_to_string(poisson_env%green_fft%method))
1216 : END SELECT
1217 :
1218 : CASE (use_rs_grid)
1219 :
1220 : ! Wavelet based Poisson-Solver
1221 32332 : SELECT CASE (poisson_env%wavelet%method)
1222 : CASE (WAVELET3D)
1223 : ! Do Nothing
1224 : CASE (WAVELET2D)
1225 : ! Zero the 1 non-periodic component
1226 0 : alpha = poisson_env%wavelet%special_dimension
1227 0 : h_stress(:, alpha) = 0.0_dp
1228 0 : h_stress(alpha, :) = 0.0_dp
1229 0 : CPABORT("Stress Tensor not tested for 2D systems.")
1230 : CASE (WAVELET1D)
1231 : ! Zero the 2 non-periodic components
1232 0 : CPABORT("WAVELET 1D not implemented!")
1233 : CASE (WAVELET0D)
1234 : ! Zero the full stress tensor
1235 0 : h_stress = 0.0_dp
1236 : END SELECT
1237 :
1238 : END SELECT
1239 : END IF
1240 :
1241 904500 : DO i = 1, 3
1242 678375 : CALL pw_pool%give_back_pw(dvg(i))
1243 904500 : IF (PRESENT(rhog_aux)) THEN
1244 1176 : CALL pw_pool%give_back_pw(dvg_aux(i))
1245 : END IF
1246 : END DO
1247 :
1248 226125 : CALL timestop(handle)
1249 :
1250 226125 : END SUBROUTINE calc_stress_and_gradient_${kind}$
1251 : #:endfor
1252 :
1253 : ! **************************************************************************************************
1254 : !> \brief sets cell, grids and parameters used by the poisson solver
1255 : !> You should call this at least once (and set everything)
1256 : !> before using the poisson solver.
1257 : !> Smart, doesn't set the thing twice to the same value
1258 : !> Keeps track of the need to rebuild the poisson_env
1259 : !> \param poisson_env ...
1260 : !> \param cell_hmat ...
1261 : !> \param parameters ...
1262 : !> \param pw_pools ...
1263 : !> \param use_level ...
1264 : !> \param mt_super_ref_pw_grid ...
1265 : !> \param dct_pw_grid ...
1266 : !> \param force_rebuild ...
1267 : !> \author fawzi
1268 : !> \note
1269 : !> Checks everything at the end. This means that after *each* call to
1270 : !> this method the poisson env must be fully ready, so the first time
1271 : !> you have to set everything at once. Change this behaviour?
1272 : ! **************************************************************************************************
1273 19732 : SUBROUTINE pw_poisson_set(poisson_env, cell_hmat, parameters, pw_pools, use_level, &
1274 : mt_super_ref_pw_grid, dct_pw_grid, force_rebuild)
1275 :
1276 : TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
1277 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN), &
1278 : OPTIONAL :: cell_hmat
1279 : TYPE(pw_poisson_parameter_type), INTENT(IN), &
1280 : OPTIONAL :: parameters
1281 : TYPE(pw_pool_p_type), DIMENSION(:), OPTIONAL, &
1282 : POINTER :: pw_pools
1283 : INTEGER, INTENT(in), OPTIONAL :: use_level
1284 : TYPE(pw_grid_type), OPTIONAL, POINTER :: mt_super_ref_pw_grid, dct_pw_grid
1285 : LOGICAL, INTENT(in), OPTIONAL :: force_rebuild
1286 :
1287 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_set'
1288 :
1289 : INTEGER :: handle, i
1290 : LOGICAL :: same
1291 19732 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: tmp_pools
1292 :
1293 19732 : CALL timeset(routineN, handle)
1294 :
1295 19732 : IF (PRESENT(parameters)) &
1296 19732 : poisson_env%parameters = parameters
1297 :
1298 19732 : IF (PRESENT(cell_hmat)) THEN
1299 43494 : IF (ANY(poisson_env%cell_hmat /= cell_hmat)) &
1300 18208 : CALL pw_poisson_cleanup(poisson_env)
1301 256516 : poisson_env%cell_hmat(:, :) = cell_hmat(:, :)
1302 19732 : poisson_env%rebuild = .TRUE.
1303 : END IF
1304 :
1305 19732 : IF (PRESENT(pw_pools)) THEN
1306 19732 : CPASSERT(ASSOCIATED(pw_pools))
1307 19732 : same = .FALSE.
1308 19732 : IF (ASSOCIATED(poisson_env%pw_pools)) THEN
1309 9820 : same = SIZE(poisson_env%pw_pools) == SIZE(pw_pools)
1310 9820 : IF (same) THEN
1311 20788 : DO i = 1, SIZE(pw_pools)
1312 10968 : IF (.NOT. ASSOCIATED(poisson_env%pw_pools(i)%pool, &
1313 11712 : pw_pools(i)%pool)) same = .FALSE.
1314 : END DO
1315 : END IF
1316 : END IF
1317 9820 : IF (.NOT. same) THEN
1318 10796 : poisson_env%rebuild = .TRUE.
1319 10796 : CALL pw_pools_copy(pw_pools, tmp_pools)
1320 10796 : CALL pw_pools_dealloc(poisson_env%pw_pools)
1321 10796 : poisson_env%pw_pools => tmp_pools
1322 : END IF
1323 : END IF
1324 :
1325 19732 : IF (PRESENT(use_level)) poisson_env%pw_level = use_level
1326 :
1327 19732 : IF (PRESENT(dct_pw_grid)) THEN
1328 8554 : IF (ASSOCIATED(dct_pw_grid)) THEN
1329 0 : CALL pw_grid_retain(dct_pw_grid)
1330 : END IF
1331 8554 : CALL pw_grid_release(poisson_env%dct_pw_grid)
1332 8554 : poisson_env%dct_pw_grid => dct_pw_grid
1333 : END IF
1334 :
1335 19732 : IF (PRESENT(mt_super_ref_pw_grid)) THEN
1336 8554 : IF (ASSOCIATED(mt_super_ref_pw_grid)) THEN
1337 1104 : CALL pw_grid_retain(mt_super_ref_pw_grid)
1338 : END IF
1339 8554 : CALL pw_grid_release(poisson_env%mt_super_ref_pw_grid)
1340 8554 : poisson_env%mt_super_ref_pw_grid => mt_super_ref_pw_grid
1341 : END IF
1342 :
1343 19732 : IF (PRESENT(force_rebuild)) THEN
1344 0 : IF (force_rebuild) poisson_env%rebuild = .TRUE.
1345 : END IF
1346 :
1347 19732 : CALL pw_poisson_check(poisson_env)
1348 :
1349 19732 : CALL timestop(handle)
1350 :
1351 19732 : END SUBROUTINE pw_poisson_set
1352 :
1353 : END MODULE pw_poisson_methods
|
