(git:b6ef100)
Loading...
Searching...
No Matches
pw_poisson_methods.F
Go to the documentation of this file.
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! **************************************************************************************************
13MODULE pw_poisson_methods
14
15 USE cp_log_handling, ONLY: cp_to_string
16 USE dielectric_methods, ONLY: dielectric_compute
17 USE kinds, ONLY: dp
18 USE mathconstants, ONLY: fourpi
19 USE mt_util, ONLY: mt0d, &
20 mt1d, &
21 mt2d
22 USE ps_implicit_methods, ONLY: implicit_poisson_solver_mixed, &
23 implicit_poisson_solver_mixed_periodic, &
24 implicit_poisson_solver_neumann, &
25 implicit_poisson_solver_periodic, &
26 ps_implicit_create
27 USE ps_implicit_types, ONLY: mixed_bc, &
28 mixed_periodic_bc, &
29 neumann_bc, &
30 periodic_bc
31 USE ps_wavelet_methods, ONLY: cp2k_distribution_to_z_slices, &
32 ps_wavelet_create, &
33 ps_wavelet_solve, &
34 z_slices_to_cp2k_distribution
35 USE ps_wavelet_types, ONLY: wavelet0d, &
36 wavelet1d, &
37 wavelet2d, &
38 wavelet3d, &
39 ps_wavelet_type
40 USE pw_grid_types, ONLY: pw_grid_type
41 USE pw_grids, ONLY: pw_grid_compare, &
42 pw_grid_release, &
43 pw_grid_retain
44 USE pw_methods, ONLY: pw_copy, &
45 pw_derive, &
46 pw_integral_ab, &
47 pw_transfer, pw_multiply_with
48 USE pw_poisson_types, ONLY: &
49 analytic0d, analytic1d, analytic2d, multipole0d, periodic3d, ps_implicit, do_ewald_spme, &
50 greens_fn_type, pw_green_create, pw_green_release, pw_poisson_analytic, &
51 pw_poisson_implicit, pw_poisson_mt, pw_poisson_multipole, pw_poisson_none, &
52 pw_poisson_parameter_type, pw_poisson_periodic, pw_poisson_type, pw_poisson_wavelet
53 USE pw_pool_types, ONLY: pw_pool_p_type, &
54 pw_pool_type, &
55 pw_pools_copy, &
56 pw_pools_dealloc
57 USE pw_types, ONLY: &
58 pw_r3d_rs_type, pw_c1d_gs_type, pw_r3d_rs_type
59#include "../base/base_uses.f90"
60
61 IMPLICIT NONE
62 PRIVATE
63
64 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .true.
65 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pw_poisson_methods'
66
67 PUBLIC :: pw_poisson_rebuild, &
68 pw_poisson_solve, pw_poisson_set, pw_func_u_convolution
69
70 INTEGER, PARAMETER :: use_rs_grid = 0, &
71 use_gs_grid = 1
72
73 INTERFACE pw_poisson_rebuild
74 MODULE PROCEDURE pw_poisson_rebuild_nodens
75 MODULE PROCEDURE pw_poisson_rebuild_c1d_gs, pw_poisson_rebuild_r3d_rs
76 END INTERFACE
77
78 INTERFACE pw_poisson_solve
79 MODULE PROCEDURE pw_poisson_solve_nov_nodv_r3d_rs
80 MODULE PROCEDURE pw_poisson_solve_v_nodv_r3d_rs_r3d_rs
81 MODULE PROCEDURE pw_poisson_solve_v_nodv_r3d_rs_c1d_gs
82 MODULE PROCEDURE pw_poisson_solve_nov_dv_r3d_rs_r3d_rs
83 MODULE PROCEDURE pw_poisson_solve_v_dv_r3d_rs_r3d_rs_r3d_rs
84 MODULE PROCEDURE pw_poisson_solve_v_dv_r3d_rs_c1d_gs_r3d_rs
85 MODULE PROCEDURE pw_poisson_solve_nov_dv_r3d_rs_c1d_gs
86 MODULE PROCEDURE pw_poisson_solve_v_dv_r3d_rs_r3d_rs_c1d_gs
87 MODULE PROCEDURE pw_poisson_solve_v_dv_r3d_rs_c1d_gs_c1d_gs
88 MODULE PROCEDURE pw_poisson_solve_nov_nodv_c1d_gs
89 MODULE PROCEDURE pw_poisson_solve_v_nodv_c1d_gs_r3d_rs
90 MODULE PROCEDURE pw_poisson_solve_v_nodv_c1d_gs_c1d_gs
91 MODULE PROCEDURE pw_poisson_solve_nov_dv_c1d_gs_r3d_rs
92 MODULE PROCEDURE pw_poisson_solve_v_dv_c1d_gs_r3d_rs_r3d_rs
93 MODULE PROCEDURE pw_poisson_solve_v_dv_c1d_gs_c1d_gs_r3d_rs
94 MODULE PROCEDURE pw_poisson_solve_nov_dv_c1d_gs_c1d_gs
95 MODULE PROCEDURE pw_poisson_solve_v_dv_c1d_gs_r3d_rs_c1d_gs
96 MODULE PROCEDURE pw_poisson_solve_v_dv_c1d_gs_c1d_gs_c1d_gs
97 END INTERFACE
98
99CONTAINS
100
101! **************************************************************************************************
102!> \brief removes all the object created from the parameters pw_pools and cell
103!> and used to solve the poisson equation like the green function and
104!> all the things allocated in pw_poisson_rebuild
105!> \param poisson_env ...
106!> \par History
107!> none
108! **************************************************************************************************
109 SUBROUTINE pw_poisson_cleanup(poisson_env)
110 TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
111
112 TYPE(pw_pool_type), POINTER :: pw_pool
113
114 NULLIFY (pw_pool)
115 IF (ASSOCIATED(poisson_env%pw_pools)) THEN
116 pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
117 END IF
118 IF (ASSOCIATED(poisson_env%green_fft)) THEN
119 CALL pw_green_release(poisson_env%green_fft, pw_pool=pw_pool)
120 DEALLOCATE (poisson_env%green_fft)
121 END IF
122 poisson_env%rebuild = .true.
123
124 END SUBROUTINE pw_poisson_cleanup
125
126! **************************************************************************************************
127!> \brief checks if pw_poisson_rebuild has to be called and calls it if needed
128!> \param poisson_env the object to be checked
129!> \author fawzi
130! **************************************************************************************************
131 SUBROUTINE pw_poisson_check(poisson_env)
132 TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
133
134 LOGICAL :: rebuild
135 TYPE(greens_fn_type), POINTER :: green
136 TYPE(ps_wavelet_type), POINTER :: wavelet
137
138 cpassert(ASSOCIATED(poisson_env%pw_pools))
139 cpassert(poisson_env%pw_level >= lbound(poisson_env%pw_pools, 1))
140 cpassert(poisson_env%pw_level <= ubound(poisson_env%pw_pools, 1))
141 green => poisson_env%green_fft
142 wavelet => poisson_env%wavelet
143 rebuild = poisson_env%rebuild
144 rebuild = rebuild .OR. (poisson_env%method /= poisson_env%parameters%solver) &
145 .OR. .NOT. ASSOCIATED(green)
146 poisson_env%method = poisson_env%parameters%solver
147
148 IF (poisson_env%method == pw_poisson_wavelet) THEN
149 poisson_env%used_grid = use_rs_grid
150 ELSE
151 poisson_env%used_grid = use_gs_grid
152 END IF
153 IF (.NOT. rebuild) THEN
154 IF (poisson_env%parameters%ewald_type == do_ewald_spme) THEN
155 rebuild = (poisson_env%parameters%ewald_alpha /= green%p3m_alpha) .OR. rebuild
156 rebuild = (poisson_env%parameters%ewald_o_spline /= green%p3m_order) .OR. rebuild
157 END IF
158 SELECT CASE (poisson_env%method)
159 CASE (pw_poisson_analytic)
160 SELECT CASE (green%method)
161 CASE (analytic0d, analytic1d, analytic2d, periodic3d)
162 CASE default
163 rebuild = .true.
164 END SELECT
165 CASE (pw_poisson_mt)
166 SELECT CASE (green%method)
167 CASE (mt0d, mt1d, mt2d)
168 CASE default
169 rebuild = .true.
170 END SELECT
171 rebuild = (poisson_env%parameters%mt_alpha /= green%mt_alpha) .OR. rebuild
172 CASE (pw_poisson_wavelet)
173 rebuild = (poisson_env%parameters%wavelet_scf_type /= wavelet%itype_scf) .OR. rebuild
174 CASE default
175 cpabort("Unknown Poisson solver")
176 END SELECT
177 END IF
178 IF (rebuild) THEN
179 poisson_env%rebuild = .true.
180 CALL pw_poisson_cleanup(poisson_env)
181 END IF
182 END SUBROUTINE pw_poisson_check
183
184! **************************************************************************************************
185!> \brief rebuilds all the internal values needed to use the poisson solver
186!> \param poisson_env the environment to rebuild
187!> \param density ...
188!> \author fawzi
189!> \note
190!> rebuilds if poisson_env%rebuild is true
191! **************************************************************************************************
192 SUBROUTINE pw_poisson_rebuild_nodens(poisson_env)
193 TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
194
195 CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_rebuild'
196
197 INTEGER :: handle
198
199 CALL timeset(routinen, handle)
200
201 cpassert(ASSOCIATED(poisson_env%pw_pools))
202
203 IF (poisson_env%rebuild) THEN
204 CALL pw_poisson_cleanup(poisson_env)
205 SELECT CASE (poisson_env%parameters%solver)
206 CASE (pw_poisson_periodic, pw_poisson_analytic, pw_poisson_mt, pw_poisson_multipole)
207 ALLOCATE (poisson_env%green_fft)
208 CALL pw_green_create(poisson_env%green_fft, cell_hmat=poisson_env%cell_hmat, &
209 pw_pool=poisson_env%pw_pools(poisson_env%pw_level)%pool, &
210 poisson_params=poisson_env%parameters, &
211 mt_super_ref_pw_grid=poisson_env%mt_super_ref_pw_grid, &
212 dct_pw_grid=poisson_env%dct_pw_grid)
213 CASE (pw_poisson_wavelet)
214 cpabort("Wavelet solver requires a density!")
215 CASE (pw_poisson_implicit)
216 ALLOCATE (poisson_env%green_fft)
217 CALL pw_green_create(poisson_env%green_fft, cell_hmat=poisson_env%cell_hmat, &
218 pw_pool=poisson_env%pw_pools(poisson_env%pw_level)%pool, &
219 poisson_params=poisson_env%parameters, &
220 mt_super_ref_pw_grid=poisson_env%mt_super_ref_pw_grid, &
221 dct_pw_grid=poisson_env%dct_pw_grid)
222 CALL ps_implicit_create(poisson_env%pw_pools(poisson_env%pw_level)%pool, &
223 poisson_env%parameters, &
224 poisson_env%dct_pw_grid, &
225 poisson_env%green_fft, poisson_env%implicit_env)
226 CASE (pw_poisson_none)
227 CASE default
228 cpabort("Unknown Poisson solver")
229 END SELECT
230 poisson_env%rebuild = .false.
231 END IF
232
233 CALL timestop(handle)
234
235 END SUBROUTINE pw_poisson_rebuild_nodens
236
237! **************************************************************************************************
238!> \brief rebuilds all the internal values needed to use the poisson solver
239!> \param poisson_env the environment to rebuild
240!> \param density ...
241!> \author fawzi
242!> \note
243!> rebuilds if poisson_env%rebuild is true
244! **************************************************************************************************
245 SUBROUTINE pw_poisson_rebuild_r3d_rs (poisson_env, density)
246 TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
247 TYPE(pw_r3d_rs_type), INTENT(IN) :: density
248
249 CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_rebuild'
250
251 INTEGER :: handle
252
253 CALL timeset(routinen, handle)
254
255 cpassert(ASSOCIATED(poisson_env%pw_pools))
256
257 IF (poisson_env%rebuild) THEN
258 CALL pw_poisson_cleanup(poisson_env)
259 SELECT CASE (poisson_env%parameters%solver)
260 CASE (pw_poisson_periodic, pw_poisson_analytic, pw_poisson_mt, pw_poisson_multipole)
261 ALLOCATE (poisson_env%green_fft)
262 CALL pw_green_create(poisson_env%green_fft, cell_hmat=poisson_env%cell_hmat, &
263 pw_pool=poisson_env%pw_pools(poisson_env%pw_level)%pool, &
264 poisson_params=poisson_env%parameters, &
265 mt_super_ref_pw_grid=poisson_env%mt_super_ref_pw_grid, &
266 dct_pw_grid=poisson_env%dct_pw_grid)
267 CASE (pw_poisson_wavelet)
268 cpassert(ASSOCIATED(density%pw_grid))
269 CALL ps_wavelet_create(poisson_env%parameters, poisson_env%wavelet, &
270 density%pw_grid)
271 CASE (pw_poisson_implicit)
272 ALLOCATE (poisson_env%green_fft)
273 CALL pw_green_create(poisson_env%green_fft, cell_hmat=poisson_env%cell_hmat, &
274 pw_pool=poisson_env%pw_pools(poisson_env%pw_level)%pool, &
275 poisson_params=poisson_env%parameters, &
276 mt_super_ref_pw_grid=poisson_env%mt_super_ref_pw_grid, &
277 dct_pw_grid=poisson_env%dct_pw_grid)
278 CALL ps_implicit_create(poisson_env%pw_pools(poisson_env%pw_level)%pool, &
279 poisson_env%parameters, &
280 poisson_env%dct_pw_grid, &
281 poisson_env%green_fft, poisson_env%implicit_env)
282 CASE (pw_poisson_none)
283 CASE default
284 cpabort("Unknown Poisson solver")
285 END SELECT
286 poisson_env%rebuild = .false.
287 END IF
288
289 CALL timestop(handle)
290
291 END SUBROUTINE pw_poisson_rebuild_r3d_rs
292! **************************************************************************************************
293!> \brief rebuilds all the internal values needed to use the poisson solver
294!> \param poisson_env the environment to rebuild
295!> \param density ...
296!> \author fawzi
297!> \note
298!> rebuilds if poisson_env%rebuild is true
299! **************************************************************************************************
300 SUBROUTINE pw_poisson_rebuild_c1d_gs (poisson_env, density)
301 TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
302 TYPE(pw_c1d_gs_type), INTENT(IN) :: density
303
304 CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_rebuild'
305
306 INTEGER :: handle
307
308 CALL timeset(routinen, handle)
309
310 cpassert(ASSOCIATED(poisson_env%pw_pools))
311
312 IF (poisson_env%rebuild) THEN
313 CALL pw_poisson_cleanup(poisson_env)
314 SELECT CASE (poisson_env%parameters%solver)
315 CASE (pw_poisson_periodic, pw_poisson_analytic, pw_poisson_mt, pw_poisson_multipole)
316 ALLOCATE (poisson_env%green_fft)
317 CALL pw_green_create(poisson_env%green_fft, cell_hmat=poisson_env%cell_hmat, &
318 pw_pool=poisson_env%pw_pools(poisson_env%pw_level)%pool, &
319 poisson_params=poisson_env%parameters, &
320 mt_super_ref_pw_grid=poisson_env%mt_super_ref_pw_grid, &
321 dct_pw_grid=poisson_env%dct_pw_grid)
322 CASE (pw_poisson_wavelet)
323 cpassert(ASSOCIATED(density%pw_grid))
324 CALL ps_wavelet_create(poisson_env%parameters, poisson_env%wavelet, &
325 density%pw_grid)
326 CASE (pw_poisson_implicit)
327 ALLOCATE (poisson_env%green_fft)
328 CALL pw_green_create(poisson_env%green_fft, cell_hmat=poisson_env%cell_hmat, &
329 pw_pool=poisson_env%pw_pools(poisson_env%pw_level)%pool, &
330 poisson_params=poisson_env%parameters, &
331 mt_super_ref_pw_grid=poisson_env%mt_super_ref_pw_grid, &
332 dct_pw_grid=poisson_env%dct_pw_grid)
333 CALL ps_implicit_create(poisson_env%pw_pools(poisson_env%pw_level)%pool, &
334 poisson_env%parameters, &
335 poisson_env%dct_pw_grid, &
336 poisson_env%green_fft, poisson_env%implicit_env)
337 CASE (pw_poisson_none)
338 CASE default
339 cpabort("Unknown Poisson solver")
340 END SELECT
341 poisson_env%rebuild = .false.
342 END IF
343
344 CALL timestop(handle)
345
346 END SUBROUTINE pw_poisson_rebuild_c1d_gs
347
348! **************************************************************************************************
349!> \brief Solve Poisson equation in a plane wave basis set
350!> Obtains electrostatic potential and its derivatives with respect to r
351!> from the density
352!> \param poisson_env ...
353!> \param density ...
354!> \param ehartree ...
355!> \param h_stress ...
356!> \param rho_core ...
357!> \param greenfn ...
358!> \param aux_density Hartree energy and stress tensor between 2 different densities
359!> \par History
360!> JGH (13-Mar-2001) : completely revised
361!> \author apsi
362! **************************************************************************************************
363 SUBROUTINE pw_poisson_solve_nov_nodv_r3d_rs (poisson_env, density, ehartree, &
364 h_stress, rho_core, greenfn, aux_density)
365
366 TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
367 TYPE(pw_r3d_rs_type), INTENT(IN) :: density
368 REAL(kind=dp), INTENT(out), OPTIONAL :: ehartree
369 REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT), &
370 OPTIONAL :: h_stress
371 TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_core, greenfn
372 TYPE(pw_r3d_rs_type), INTENT(IN), OPTIONAL :: aux_density
373
374 CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_solve'
375
376 INTEGER :: handle
377 LOGICAL :: has_dielectric
378 TYPE(pw_grid_type), POINTER :: pw_grid
379 TYPE(pw_pool_type), POINTER :: pw_pool
380 TYPE(pw_r3d_rs_type) :: rhor, vhartree_rs
381 TYPE(pw_c1d_gs_type) :: influence_fn, rhog, rhog_aux, tmpg
382
383 CALL timeset(routinen, handle)
384
385 CALL pw_poisson_rebuild(poisson_env, density)
386
387 has_dielectric = poisson_env%parameters%has_dielectric
388
389 ! point pw
390 pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
391 pw_grid => pw_pool%pw_grid
392 ! density in G space
393 CALL pw_pool%create_pw(rhog)
394 IF (PRESENT(aux_density)) THEN
395 CALL pw_pool%create_pw(rhog_aux)
396 END IF
397
398 SELECT CASE (poisson_env%used_grid)
399 CASE (use_gs_grid)
400
401 SELECT CASE (poisson_env%green_fft%method)
402 CASE (periodic3d, analytic2d, analytic1d, analytic0d, mt2d, mt1d, mt0d, multipole0d)
403
404 CALL pw_transfer(density, rhog)
405 IF (PRESENT(aux_density)) THEN
406 CALL pw_transfer(aux_density, rhog_aux)
407 END IF
408 IF (PRESENT(ehartree)) THEN
409 CALL pw_pool%create_pw(tmpg)
410 CALL pw_copy(rhog, tmpg)
411 END IF
412 IF (PRESENT(greenfn)) THEN
413 influence_fn = greenfn
414 ELSE
415 influence_fn = poisson_env%green_fft%influence_fn
416 END IF
417 CALL pw_multiply_with(rhog, influence_fn)
418 IF (PRESENT(aux_density)) THEN
419 CALL pw_multiply_with(rhog_aux, influence_fn)
420 END IF
421 IF (PRESENT(ehartree)) THEN
422 IF (PRESENT(aux_density)) THEN
423 ehartree = 0.5_dp*pw_integral_ab(rhog_aux, tmpg)
424 ELSE
425 ehartree = 0.5_dp*pw_integral_ab(rhog, tmpg)
426 END IF
427 CALL pw_pool%give_back_pw(tmpg)
428 END IF
429
430 CASE (ps_implicit)
431
432 IF (PRESENT(h_stress)) THEN
433 cpabort("No stress tensor is implemented for the implicit Poisson solver.")
434 END IF
435
436 IF (has_dielectric .AND. PRESENT(rho_core)) THEN
437 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
438 CASE (periodic_bc, mixed_periodic_bc)
439 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
440 poisson_env%diel_rs_grid, &
441 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
442 density, rho_core=rho_core)
443 CASE (neumann_bc, mixed_bc)
444 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
445 poisson_env%diel_rs_grid, &
446 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
447 poisson_env%dct_pw_grid, &
448 poisson_env%parameters%ps_implicit_params%neumann_directions, &
449 poisson_env%implicit_env%dct_env%recv_msgs_bnds, &
450 poisson_env%implicit_env%dct_env%dests_expand, &
451 poisson_env%implicit_env%dct_env%srcs_expand, &
452 poisson_env%implicit_env%dct_env%flipg_stat, &
453 poisson_env%implicit_env%dct_env%bounds_shftd, &
454 density, rho_core=rho_core)
455 END SELECT
456 END IF
457
458 CALL pw_pool%create_pw(rhor)
459 CALL pw_pool%create_pw(vhartree_rs)
460 CALL pw_transfer(density, rhor)
461
462 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
463 CASE (periodic_bc)
464 CALL implicit_poisson_solver_periodic(poisson_env, rhor, vhartree_rs, &
465 ehartree=ehartree)
466 CASE (neumann_bc)
467 CALL implicit_poisson_solver_neumann(poisson_env, rhor, vhartree_rs, &
468 ehartree=ehartree)
469 CASE (mixed_periodic_bc)
470 CALL implicit_poisson_solver_mixed_periodic(poisson_env, rhor, vhartree_rs, &
471 electric_enthalpy=ehartree)
472 CASE (mixed_bc)
473 CALL implicit_poisson_solver_mixed(poisson_env, rhor, vhartree_rs, &
474 electric_enthalpy=ehartree)
475 END SELECT
476
477 IF (PRESENT(aux_density)) THEN
478 CALL pw_transfer(aux_density, rhor)
479 ehartree = 0.5_dp*pw_integral_ab(rhor, vhartree_rs)
480 END IF
481
482 CALL pw_pool%give_back_pw(rhor)
483 CALL pw_pool%give_back_pw(vhartree_rs)
484
485 CASE DEFAULT
486 CALL cp_abort(__location__, &
487 "unknown poisson method "// &
488 cp_to_string(poisson_env%green_fft%method))
489 END SELECT
490
491 CASE (use_rs_grid)
492
493 CALL pw_pool%create_pw(rhor)
494 CALL pw_transfer(density, rhor)
495 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
496 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
497 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
498 IF (PRESENT(ehartree)) THEN
499 IF (PRESENT(aux_density)) THEN
500 ehartree = 0.5_dp*pw_integral_ab(aux_density, rhor)
501 ELSE
502 ehartree = 0.5_dp*pw_integral_ab(density, rhor)
503 END IF
504 END IF
505 IF (PRESENT(h_stress)) THEN
506 CALL pw_transfer(rhor, rhog)
507 IF (PRESENT(aux_density)) THEN
508 CALL pw_transfer(aux_density, rhor)
509 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
510 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
511 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
512 CALL pw_transfer(rhor, rhog_aux)
513 END IF
514 END IF
515 CALL pw_pool%give_back_pw(rhor)
516
517 END SELECT
518
519 IF (PRESENT(aux_density)) THEN
520 CALL calc_stress_and_gradient_r3d_rs (poisson_env, rhog, ehartree, rhog_aux, h_stress)
521 ELSE
522 CALL calc_stress_and_gradient_r3d_rs (poisson_env, rhog, ehartree, h_stress=h_stress)
523 END IF
524
525 CALL pw_pool%give_back_pw(rhog)
526 IF (PRESENT(aux_density)) THEN
527 CALL pw_pool%give_back_pw(rhog_aux)
528 END IF
529
530 CALL timestop(handle)
531
532 END SUBROUTINE pw_poisson_solve_nov_nodv_r3d_rs
533! **************************************************************************************************
534!> \brief Solve Poisson equation in a plane wave basis set
535!> Obtains electrostatic potential and its derivatives with respect to r
536!> from the density
537!> \param poisson_env ...
538!> \param density ...
539!> \param ehartree ...
540!> \param h_stress ...
541!> \param rho_core ...
542!> \param greenfn ...
543!> \param aux_density Hartree energy and stress tensor between 2 different densities
544!> \par History
545!> JGH (13-Mar-2001) : completely revised
546!> \author apsi
547! **************************************************************************************************
548 SUBROUTINE pw_poisson_solve_nov_nodv_c1d_gs (poisson_env, density, ehartree, &
549 h_stress, rho_core, greenfn, aux_density)
550
551 TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
552 TYPE(pw_c1d_gs_type), INTENT(IN) :: density
553 REAL(kind=dp), INTENT(out), OPTIONAL :: ehartree
554 REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT), &
555 OPTIONAL :: h_stress
556 TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_core, greenfn
557 TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: aux_density
558
559 CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_solve'
560
561 INTEGER :: handle
562 LOGICAL :: has_dielectric
563 TYPE(pw_grid_type), POINTER :: pw_grid
564 TYPE(pw_pool_type), POINTER :: pw_pool
565 TYPE(pw_r3d_rs_type) :: rhor, vhartree_rs
566 TYPE(pw_c1d_gs_type) :: influence_fn, rhog, rhog_aux, tmpg
567
568 CALL timeset(routinen, handle)
569
570 CALL pw_poisson_rebuild(poisson_env, density)
571
572 has_dielectric = poisson_env%parameters%has_dielectric
573
574 ! point pw
575 pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
576 pw_grid => pw_pool%pw_grid
577 ! density in G space
578 CALL pw_pool%create_pw(rhog)
579 IF (PRESENT(aux_density)) THEN
580 CALL pw_pool%create_pw(rhog_aux)
581 END IF
582
583 SELECT CASE (poisson_env%used_grid)
584 CASE (use_gs_grid)
585
586 SELECT CASE (poisson_env%green_fft%method)
587 CASE (periodic3d, analytic2d, analytic1d, analytic0d, mt2d, mt1d, mt0d, multipole0d)
588
589 CALL pw_transfer(density, rhog)
590 IF (PRESENT(aux_density)) THEN
591 CALL pw_transfer(aux_density, rhog_aux)
592 END IF
593 IF (PRESENT(ehartree)) THEN
594 CALL pw_pool%create_pw(tmpg)
595 CALL pw_copy(rhog, tmpg)
596 END IF
597 IF (PRESENT(greenfn)) THEN
598 influence_fn = greenfn
599 ELSE
600 influence_fn = poisson_env%green_fft%influence_fn
601 END IF
602 CALL pw_multiply_with(rhog, influence_fn)
603 IF (PRESENT(aux_density)) THEN
604 CALL pw_multiply_with(rhog_aux, influence_fn)
605 END IF
606 IF (PRESENT(ehartree)) THEN
607 IF (PRESENT(aux_density)) THEN
608 ehartree = 0.5_dp*pw_integral_ab(rhog_aux, tmpg)
609 ELSE
610 ehartree = 0.5_dp*pw_integral_ab(rhog, tmpg)
611 END IF
612 CALL pw_pool%give_back_pw(tmpg)
613 END IF
614
615 CASE (ps_implicit)
616
617 IF (PRESENT(h_stress)) THEN
618 cpabort("No stress tensor is implemented for the implicit Poisson solver.")
619 END IF
620
621 IF (has_dielectric .AND. PRESENT(rho_core)) THEN
622 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
623 CASE (periodic_bc, mixed_periodic_bc)
624 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
625 poisson_env%diel_rs_grid, &
626 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
627 density, rho_core=rho_core)
628 CASE (neumann_bc, mixed_bc)
629 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
630 poisson_env%diel_rs_grid, &
631 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
632 poisson_env%dct_pw_grid, &
633 poisson_env%parameters%ps_implicit_params%neumann_directions, &
634 poisson_env%implicit_env%dct_env%recv_msgs_bnds, &
635 poisson_env%implicit_env%dct_env%dests_expand, &
636 poisson_env%implicit_env%dct_env%srcs_expand, &
637 poisson_env%implicit_env%dct_env%flipg_stat, &
638 poisson_env%implicit_env%dct_env%bounds_shftd, &
639 density, rho_core=rho_core)
640 END SELECT
641 END IF
642
643 CALL pw_pool%create_pw(rhor)
644 CALL pw_pool%create_pw(vhartree_rs)
645 CALL pw_transfer(density, rhor)
646
647 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
648 CASE (periodic_bc)
649 CALL implicit_poisson_solver_periodic(poisson_env, rhor, vhartree_rs, &
650 ehartree=ehartree)
651 CASE (neumann_bc)
652 CALL implicit_poisson_solver_neumann(poisson_env, rhor, vhartree_rs, &
653 ehartree=ehartree)
654 CASE (mixed_periodic_bc)
655 CALL implicit_poisson_solver_mixed_periodic(poisson_env, rhor, vhartree_rs, &
656 electric_enthalpy=ehartree)
657 CASE (mixed_bc)
658 CALL implicit_poisson_solver_mixed(poisson_env, rhor, vhartree_rs, &
659 electric_enthalpy=ehartree)
660 END SELECT
661
662 IF (PRESENT(aux_density)) THEN
663 CALL pw_transfer(aux_density, rhor)
664 ehartree = 0.5_dp*pw_integral_ab(rhor, vhartree_rs)
665 END IF
666
667 CALL pw_pool%give_back_pw(rhor)
668 CALL pw_pool%give_back_pw(vhartree_rs)
669
670 CASE DEFAULT
671 CALL cp_abort(__location__, &
672 "unknown poisson method "// &
673 cp_to_string(poisson_env%green_fft%method))
674 END SELECT
675
676 CASE (use_rs_grid)
677
678 CALL pw_pool%create_pw(rhor)
679 CALL pw_transfer(density, rhor)
680 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
681 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
682 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
683 IF (PRESENT(ehartree)) THEN
684 IF (PRESENT(aux_density)) THEN
685 IF (.NOT. PRESENT(h_stress)) CALL pw_pool%create_pw(rhog)
686 CALL pw_transfer(rhor, rhog)
687 ehartree = 0.5_dp*pw_integral_ab(aux_density, rhog)
688 IF (.NOT. PRESENT(h_stress)) CALL pw_pool%give_back_pw(rhog)
689 ELSE
690 IF (.NOT. PRESENT(h_stress)) CALL pw_pool%create_pw(rhog)
691 CALL pw_transfer(rhor, rhog)
692 ehartree = 0.5_dp*pw_integral_ab(density, rhog)
693 IF (.NOT. PRESENT(h_stress)) CALL pw_pool%give_back_pw(rhog)
694 END IF
695 END IF
696 IF (PRESENT(h_stress)) THEN
697 CALL pw_transfer(rhor, rhog)
698 IF (PRESENT(aux_density)) THEN
699 CALL pw_transfer(aux_density, rhor)
700 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
701 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
702 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
703 CALL pw_transfer(rhor, rhog_aux)
704 END IF
705 END IF
706 CALL pw_pool%give_back_pw(rhor)
707
708 END SELECT
709
710 IF (PRESENT(aux_density)) THEN
711 CALL calc_stress_and_gradient_c1d_gs (poisson_env, rhog, ehartree, rhog_aux, h_stress)
712 ELSE
713 CALL calc_stress_and_gradient_c1d_gs (poisson_env, rhog, ehartree, h_stress=h_stress)
714 END IF
715
716 CALL pw_pool%give_back_pw(rhog)
717 IF (PRESENT(aux_density)) THEN
718 CALL pw_pool%give_back_pw(rhog_aux)
719 END IF
720
721 CALL timestop(handle)
722
723 END SUBROUTINE pw_poisson_solve_nov_nodv_c1d_gs
724
725! **************************************************************************************************
726!> \brief Solve Poisson equation in a plane wave basis set
727!> Obtains electrostatic potential and its derivatives with respect to r
728!> from the density
729!> \param poisson_env ...
730!> \param density ...
731!> \param ehartree ...
732!> \param vhartree ...
733!> \param h_stress ...
734!> \param rho_core ...
735!> \param greenfn ...
736!> \param aux_density Hartree energy and stress tensor between 2 different densities
737!> \par History
738!> JGH (13-Mar-2001) : completely revised
739!> \author apsi
740! **************************************************************************************************
741 SUBROUTINE pw_poisson_solve_v_nodv_r3d_rs_r3d_rs (poisson_env, density, ehartree, vhartree, &
742 h_stress, rho_core, greenfn, aux_density)
743
744 TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
745 TYPE(pw_r3d_rs_type), INTENT(IN) :: density
746 REAL(kind=dp), INTENT(out), OPTIONAL :: ehartree
747 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: vhartree
748 REAL(kind=dp), DIMENSION(3, 3), INTENT(OUT), &
749 OPTIONAL :: h_stress
750 TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_core, greenfn
751 TYPE(pw_r3d_rs_type), INTENT(IN), OPTIONAL :: aux_density
752
753 CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_solve'
754
755 INTEGER :: handle
756 LOGICAL :: has_dielectric
757 TYPE(pw_grid_type), POINTER :: pw_grid
758 TYPE(pw_pool_type), POINTER :: pw_pool
759 TYPE(pw_r3d_rs_type) :: &
760 rhor, vhartree_rs
761 TYPE(pw_c1d_gs_type) :: influence_fn, rhog, rhog_aux
762
763 CALL timeset(routinen, handle)
764
765 CALL pw_poisson_rebuild(poisson_env, density)
766
767 has_dielectric = poisson_env%parameters%has_dielectric
768
769 ! point pw
770 pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
771 pw_grid => pw_pool%pw_grid
772 IF (.NOT. pw_grid_compare(pw_pool%pw_grid, vhartree%pw_grid)) THEN
773 cpabort("vhartree has a different grid than the poisson solver")
774 END IF
775 ! density in G space
776 CALL pw_pool%create_pw(rhog)
777 IF (PRESENT(aux_density)) THEN
778 CALL pw_pool%create_pw(rhog_aux)
779 END IF
780
781 SELECT CASE (poisson_env%used_grid)
782 CASE (use_gs_grid)
783
784 SELECT CASE (poisson_env%green_fft%method)
785 CASE (periodic3d, analytic2d, analytic1d, analytic0d, mt2d, mt1d, mt0d, multipole0d)
786
787 CALL pw_transfer(density, rhog)
788 IF (PRESENT(aux_density)) THEN
789 CALL pw_transfer(aux_density, rhog_aux)
790 END IF
791 IF (PRESENT(greenfn)) THEN
792 influence_fn = greenfn
793 ELSE
794 influence_fn = poisson_env%green_fft%influence_fn
795 END IF
796 CALL pw_multiply_with(rhog, influence_fn)
797 IF (PRESENT(aux_density)) THEN
798 CALL pw_multiply_with(rhog_aux, influence_fn)
799 END IF
800 CALL pw_transfer(rhog, vhartree)
801 IF (PRESENT(ehartree)) THEN
802 IF (PRESENT(aux_density)) THEN
803 ehartree = 0.5_dp*pw_integral_ab(aux_density, vhartree)
804 ELSE
805 ehartree = 0.5_dp*pw_integral_ab(density, vhartree)
806 END IF
807 END IF
808
809 CASE (ps_implicit)
810 IF (PRESENT(h_stress)) THEN
811 cpabort("No stress tensor is implemented for the implicit Poisson solver.")
812 END IF
813
814 IF (has_dielectric .AND. PRESENT(rho_core)) THEN
815 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
816 CASE (periodic_bc, mixed_periodic_bc)
817 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
818 poisson_env%diel_rs_grid, &
819 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
820 density, rho_core=rho_core)
821 CASE (neumann_bc, mixed_bc)
822 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
823 poisson_env%diel_rs_grid, &
824 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
825 poisson_env%dct_pw_grid, &
826 poisson_env%parameters%ps_implicit_params%neumann_directions, &
827 poisson_env%implicit_env%dct_env%recv_msgs_bnds, &
828 poisson_env%implicit_env%dct_env%dests_expand, &
829 poisson_env%implicit_env%dct_env%srcs_expand, &
830 poisson_env%implicit_env%dct_env%flipg_stat, &
831 poisson_env%implicit_env%dct_env%bounds_shftd, &
832 density, rho_core=rho_core)
833 END SELECT
834 END IF
835
836 CALL pw_pool%create_pw(rhor)
837 CALL pw_pool%create_pw(vhartree_rs)
838 CALL pw_transfer(density, rhor)
839
840 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
841 CASE (periodic_bc)
842 CALL implicit_poisson_solver_periodic(poisson_env, rhor, vhartree_rs, &
843 ehartree=ehartree)
844 CASE (neumann_bc)
845 CALL implicit_poisson_solver_neumann(poisson_env, rhor, vhartree_rs, &
846 ehartree=ehartree)
847 CASE (mixed_periodic_bc)
848 CALL implicit_poisson_solver_mixed_periodic(poisson_env, rhor, vhartree_rs, &
849 electric_enthalpy=ehartree)
850 CASE (mixed_bc)
851 CALL implicit_poisson_solver_mixed(poisson_env, rhor, vhartree_rs, &
852 electric_enthalpy=ehartree)
853 END SELECT
854
855 IF (PRESENT(aux_density)) THEN
856 CALL pw_transfer(aux_density, rhor)
857 ehartree = 0.5_dp*pw_integral_ab(rhor, vhartree_rs)
858 END IF
859
860 CALL pw_transfer(vhartree_rs, vhartree)
861
862 CALL pw_pool%give_back_pw(rhor)
863 CALL pw_pool%give_back_pw(vhartree_rs)
864
865 CASE DEFAULT
866 CALL cp_abort(__location__, &
867 "unknown poisson method "// &
868 cp_to_string(poisson_env%green_fft%method))
869 END SELECT
870
871 CASE (use_rs_grid)
872
873 CALL pw_pool%create_pw(rhor)
874 CALL pw_transfer(density, rhor)
875 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
876 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
877 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
878 CALL pw_transfer(rhor, vhartree)
879 IF (PRESENT(ehartree)) THEN
880 IF (PRESENT(aux_density)) THEN
881 ehartree = 0.5_dp*pw_integral_ab(aux_density, vhartree)
882 ELSE
883 ehartree = 0.5_dp*pw_integral_ab(density, vhartree)
884 END IF
885 END IF
886 IF (PRESENT(h_stress)) THEN
887 CALL pw_transfer(rhor, rhog)
888 IF (PRESENT(aux_density)) THEN
889 CALL pw_transfer(aux_density, rhor)
890 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
891 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
892 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
893 CALL pw_transfer(rhor, rhog_aux)
894 END IF
895 END IF
896 CALL pw_pool%give_back_pw(rhor)
897
898 END SELECT
899
900 IF (PRESENT(aux_density)) THEN
901 CALL calc_stress_and_gradient_r3d_rs (poisson_env, rhog, ehartree, rhog_aux, h_stress)
902 ELSE
903 CALL calc_stress_and_gradient_r3d_rs (poisson_env, rhog, ehartree, h_stress=h_stress)
904 END IF
905
906 CALL pw_pool%give_back_pw(rhog)
907 IF (PRESENT(aux_density)) THEN
908 CALL pw_pool%give_back_pw(rhog_aux)
909 END IF
910
911 CALL timestop(handle)
912
913 END SUBROUTINE pw_poisson_solve_v_nodv_r3d_rs_r3d_rs
914! **************************************************************************************************
915!> \brief Solve Poisson equation in a plane wave basis set
916!> Obtains electrostatic potential and its derivatives with respect to r
917!> from the density
918!> \param poisson_env ...
919!> \param density ...
920!> \param ehartree ...
921!> \param vhartree ...
922!> \param h_stress ...
923!> \param rho_core ...
924!> \param greenfn ...
925!> \param aux_density Hartree energy and stress tensor between 2 different densities
926!> \par History
927!> JGH (13-Mar-2001) : completely revised
928!> \author apsi
929! **************************************************************************************************
930 SUBROUTINE pw_poisson_solve_v_nodv_r3d_rs_c1d_gs (poisson_env, density, ehartree, vhartree, &
931 h_stress, rho_core, greenfn, aux_density)
932
933 TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
934 TYPE(pw_r3d_rs_type), INTENT(IN) :: density
935 REAL(kind=dp), INTENT(out), OPTIONAL :: ehartree
936 TYPE(pw_c1d_gs_type), INTENT(INOUT) :: vhartree
937 REAL(kind=dp), DIMENSION(3, 3), INTENT(OUT), &
938 OPTIONAL :: h_stress
939 TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_core, greenfn
940 TYPE(pw_r3d_rs_type), INTENT(IN), OPTIONAL :: aux_density
941
942 CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_solve'
943
944 INTEGER :: handle
945 LOGICAL :: has_dielectric
946 TYPE(pw_grid_type), POINTER :: pw_grid
947 TYPE(pw_pool_type), POINTER :: pw_pool
948 TYPE(pw_r3d_rs_type) :: &
949 rhor, vhartree_rs
950 TYPE(pw_c1d_gs_type) :: influence_fn, rhog, rhog_aux
951
952 CALL timeset(routinen, handle)
953
954 CALL pw_poisson_rebuild(poisson_env, density)
955
956 has_dielectric = poisson_env%parameters%has_dielectric
957
958 ! point pw
959 pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
960 pw_grid => pw_pool%pw_grid
961 IF (.NOT. pw_grid_compare(pw_pool%pw_grid, vhartree%pw_grid)) THEN
962 cpabort("vhartree has a different grid than the poisson solver")
963 END IF
964 ! density in G space
965 CALL pw_pool%create_pw(rhog)
966 IF (PRESENT(aux_density)) THEN
967 CALL pw_pool%create_pw(rhog_aux)
968 END IF
969
970 SELECT CASE (poisson_env%used_grid)
971 CASE (use_gs_grid)
972
973 SELECT CASE (poisson_env%green_fft%method)
974 CASE (periodic3d, analytic2d, analytic1d, analytic0d, mt2d, mt1d, mt0d, multipole0d)
975
976 CALL pw_transfer(density, rhog)
977 IF (PRESENT(aux_density)) THEN
978 CALL pw_transfer(aux_density, rhog_aux)
979 END IF
980 IF (PRESENT(greenfn)) THEN
981 influence_fn = greenfn
982 ELSE
983 influence_fn = poisson_env%green_fft%influence_fn
984 END IF
985 CALL pw_multiply_with(rhog, influence_fn)
986 IF (PRESENT(aux_density)) THEN
987 CALL pw_multiply_with(rhog_aux, influence_fn)
988 END IF
989 CALL pw_transfer(rhog, vhartree)
990 IF (PRESENT(ehartree)) THEN
991 IF (PRESENT(aux_density)) THEN
992 CALL pw_transfer(aux_density, rhog)
993 ehartree = 0.5_dp*pw_integral_ab(rhog, vhartree)
994 ELSE
995 CALL pw_transfer(density, rhog)
996 ehartree = 0.5_dp*pw_integral_ab(rhog, vhartree)
997 END IF
998 END IF
999
1000 CASE (ps_implicit)
1001 IF (PRESENT(h_stress)) THEN
1002 cpabort("No stress tensor is implemented for the implicit Poisson solver.")
1003 END IF
1004
1005 IF (has_dielectric .AND. PRESENT(rho_core)) THEN
1006 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
1007 CASE (periodic_bc, mixed_periodic_bc)
1008 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
1009 poisson_env%diel_rs_grid, &
1010 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
1011 density, rho_core=rho_core)
1012 CASE (neumann_bc, mixed_bc)
1013 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
1014 poisson_env%diel_rs_grid, &
1015 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
1016 poisson_env%dct_pw_grid, &
1017 poisson_env%parameters%ps_implicit_params%neumann_directions, &
1018 poisson_env%implicit_env%dct_env%recv_msgs_bnds, &
1019 poisson_env%implicit_env%dct_env%dests_expand, &
1020 poisson_env%implicit_env%dct_env%srcs_expand, &
1021 poisson_env%implicit_env%dct_env%flipg_stat, &
1022 poisson_env%implicit_env%dct_env%bounds_shftd, &
1023 density, rho_core=rho_core)
1024 END SELECT
1025 END IF
1026
1027 CALL pw_pool%create_pw(rhor)
1028 CALL pw_pool%create_pw(vhartree_rs)
1029 CALL pw_transfer(density, rhor)
1030
1031 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
1032 CASE (periodic_bc)
1033 CALL implicit_poisson_solver_periodic(poisson_env, rhor, vhartree_rs, &
1034 ehartree=ehartree)
1035 CASE (neumann_bc)
1036 CALL implicit_poisson_solver_neumann(poisson_env, rhor, vhartree_rs, &
1037 ehartree=ehartree)
1038 CASE (mixed_periodic_bc)
1039 CALL implicit_poisson_solver_mixed_periodic(poisson_env, rhor, vhartree_rs, &
1040 electric_enthalpy=ehartree)
1041 CASE (mixed_bc)
1042 CALL implicit_poisson_solver_mixed(poisson_env, rhor, vhartree_rs, &
1043 electric_enthalpy=ehartree)
1044 END SELECT
1045
1046 IF (PRESENT(aux_density)) THEN
1047 CALL pw_transfer(aux_density, rhor)
1048 ehartree = 0.5_dp*pw_integral_ab(rhor, vhartree_rs)
1049 END IF
1050
1051 CALL pw_transfer(vhartree_rs, vhartree)
1052
1053 CALL pw_pool%give_back_pw(rhor)
1054 CALL pw_pool%give_back_pw(vhartree_rs)
1055
1056 CASE DEFAULT
1057 CALL cp_abort(__location__, &
1058 "unknown poisson method "// &
1059 cp_to_string(poisson_env%green_fft%method))
1060 END SELECT
1061
1062 CASE (use_rs_grid)
1063
1064 CALL pw_pool%create_pw(rhor)
1065 CALL pw_transfer(density, rhor)
1066 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
1067 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
1068 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
1069 CALL pw_transfer(rhor, vhartree)
1070 IF (PRESENT(ehartree)) THEN
1071 IF (PRESENT(aux_density)) THEN
1072 ehartree = 0.5_dp*pw_integral_ab(aux_density, rhor)
1073 ELSE
1074 ehartree = 0.5_dp*pw_integral_ab(density, rhor)
1075 END IF
1076 END IF
1077 IF (PRESENT(h_stress)) THEN
1078 CALL pw_transfer(rhor, rhog)
1079 IF (PRESENT(aux_density)) THEN
1080 CALL pw_transfer(aux_density, rhor)
1081 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
1082 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
1083 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
1084 CALL pw_transfer(rhor, rhog_aux)
1085 END IF
1086 END IF
1087 CALL pw_pool%give_back_pw(rhor)
1088
1089 END SELECT
1090
1091 IF (PRESENT(aux_density)) THEN
1092 CALL calc_stress_and_gradient_r3d_rs (poisson_env, rhog, ehartree, rhog_aux, h_stress)
1093 ELSE
1094 CALL calc_stress_and_gradient_r3d_rs (poisson_env, rhog, ehartree, h_stress=h_stress)
1095 END IF
1096
1097 CALL pw_pool%give_back_pw(rhog)
1098 IF (PRESENT(aux_density)) THEN
1099 CALL pw_pool%give_back_pw(rhog_aux)
1100 END IF
1101
1102 CALL timestop(handle)
1103
1104 END SUBROUTINE pw_poisson_solve_v_nodv_r3d_rs_c1d_gs
1105! **************************************************************************************************
1106!> \brief Solve Poisson equation in a plane wave basis set
1107!> Obtains electrostatic potential and its derivatives with respect to r
1108!> from the density
1109!> \param poisson_env ...
1110!> \param density ...
1111!> \param ehartree ...
1112!> \param vhartree ...
1113!> \param h_stress ...
1114!> \param rho_core ...
1115!> \param greenfn ...
1116!> \param aux_density Hartree energy and stress tensor between 2 different densities
1117!> \par History
1118!> JGH (13-Mar-2001) : completely revised
1119!> \author apsi
1120! **************************************************************************************************
1121 SUBROUTINE pw_poisson_solve_v_nodv_c1d_gs_r3d_rs (poisson_env, density, ehartree, vhartree, &
1122 h_stress, rho_core, greenfn, aux_density)
1123
1124 TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
1125 TYPE(pw_c1d_gs_type), INTENT(IN) :: density
1126 REAL(kind=dp), INTENT(out), OPTIONAL :: ehartree
1127 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: vhartree
1128 REAL(kind=dp), DIMENSION(3, 3), INTENT(OUT), &
1129 OPTIONAL :: h_stress
1130 TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_core, greenfn
1131 TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: aux_density
1132
1133 CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_solve'
1134
1135 INTEGER :: handle
1136 LOGICAL :: has_dielectric
1137 TYPE(pw_grid_type), POINTER :: pw_grid
1138 TYPE(pw_pool_type), POINTER :: pw_pool
1139 TYPE(pw_r3d_rs_type) :: &
1140 rhor, vhartree_rs
1141 TYPE(pw_c1d_gs_type) :: influence_fn, rhog, rhog_aux
1142
1143 CALL timeset(routinen, handle)
1144
1145 CALL pw_poisson_rebuild(poisson_env, density)
1146
1147 has_dielectric = poisson_env%parameters%has_dielectric
1148
1149 ! point pw
1150 pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
1151 pw_grid => pw_pool%pw_grid
1152 IF (.NOT. pw_grid_compare(pw_pool%pw_grid, vhartree%pw_grid)) THEN
1153 cpabort("vhartree has a different grid than the poisson solver")
1154 END IF
1155 ! density in G space
1156 CALL pw_pool%create_pw(rhog)
1157 IF (PRESENT(aux_density)) THEN
1158 CALL pw_pool%create_pw(rhog_aux)
1159 END IF
1160
1161 SELECT CASE (poisson_env%used_grid)
1162 CASE (use_gs_grid)
1163
1164 SELECT CASE (poisson_env%green_fft%method)
1165 CASE (periodic3d, analytic2d, analytic1d, analytic0d, mt2d, mt1d, mt0d, multipole0d)
1166
1167 CALL pw_transfer(density, rhog)
1168 IF (PRESENT(aux_density)) THEN
1169 CALL pw_transfer(aux_density, rhog_aux)
1170 END IF
1171 IF (PRESENT(greenfn)) THEN
1172 influence_fn = greenfn
1173 ELSE
1174 influence_fn = poisson_env%green_fft%influence_fn
1175 END IF
1176 CALL pw_multiply_with(rhog, influence_fn)
1177 IF (PRESENT(aux_density)) THEN
1178 CALL pw_multiply_with(rhog_aux, influence_fn)
1179 END IF
1180 CALL pw_transfer(rhog, vhartree)
1181 IF (PRESENT(ehartree)) THEN
1182 IF (PRESENT(aux_density)) THEN
1183 ehartree = 0.5_dp*pw_integral_ab(aux_density, rhog)
1184 ELSE
1185 ehartree = 0.5_dp*pw_integral_ab(density, rhog)
1186 END IF
1187 END IF
1188
1189 CASE (ps_implicit)
1190 IF (PRESENT(h_stress)) THEN
1191 cpabort("No stress tensor is implemented for the implicit Poisson solver.")
1192 END IF
1193
1194 IF (has_dielectric .AND. PRESENT(rho_core)) THEN
1195 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
1196 CASE (periodic_bc, mixed_periodic_bc)
1197 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
1198 poisson_env%diel_rs_grid, &
1199 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
1200 density, rho_core=rho_core)
1201 CASE (neumann_bc, mixed_bc)
1202 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
1203 poisson_env%diel_rs_grid, &
1204 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
1205 poisson_env%dct_pw_grid, &
1206 poisson_env%parameters%ps_implicit_params%neumann_directions, &
1207 poisson_env%implicit_env%dct_env%recv_msgs_bnds, &
1208 poisson_env%implicit_env%dct_env%dests_expand, &
1209 poisson_env%implicit_env%dct_env%srcs_expand, &
1210 poisson_env%implicit_env%dct_env%flipg_stat, &
1211 poisson_env%implicit_env%dct_env%bounds_shftd, &
1212 density, rho_core=rho_core)
1213 END SELECT
1214 END IF
1215
1216 CALL pw_pool%create_pw(rhor)
1217 CALL pw_pool%create_pw(vhartree_rs)
1218 CALL pw_transfer(density, rhor)
1219
1220 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
1221 CASE (periodic_bc)
1222 CALL implicit_poisson_solver_periodic(poisson_env, rhor, vhartree_rs, &
1223 ehartree=ehartree)
1224 CASE (neumann_bc)
1225 CALL implicit_poisson_solver_neumann(poisson_env, rhor, vhartree_rs, &
1226 ehartree=ehartree)
1227 CASE (mixed_periodic_bc)
1228 CALL implicit_poisson_solver_mixed_periodic(poisson_env, rhor, vhartree_rs, &
1229 electric_enthalpy=ehartree)
1230 CASE (mixed_bc)
1231 CALL implicit_poisson_solver_mixed(poisson_env, rhor, vhartree_rs, &
1232 electric_enthalpy=ehartree)
1233 END SELECT
1234
1235 IF (PRESENT(aux_density)) THEN
1236 CALL pw_transfer(aux_density, rhor)
1237 ehartree = 0.5_dp*pw_integral_ab(rhor, vhartree_rs)
1238 END IF
1239
1240 CALL pw_transfer(vhartree_rs, vhartree)
1241
1242 CALL pw_pool%give_back_pw(rhor)
1243 CALL pw_pool%give_back_pw(vhartree_rs)
1244
1245 CASE DEFAULT
1246 CALL cp_abort(__location__, &
1247 "unknown poisson method "// &
1248 cp_to_string(poisson_env%green_fft%method))
1249 END SELECT
1250
1251 CASE (use_rs_grid)
1252
1253 CALL pw_pool%create_pw(rhor)
1254 CALL pw_transfer(density, rhor)
1255 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
1256 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
1257 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
1258 CALL pw_transfer(rhor, vhartree)
1259 IF (PRESENT(ehartree)) THEN
1260 IF (PRESENT(aux_density)) THEN
1261 CALL pw_transfer(vhartree, rhog)
1262 ehartree = 0.5_dp*pw_integral_ab(aux_density, rhog)
1263 ELSE
1264 CALL pw_transfer(vhartree, rhog)
1265 ehartree = 0.5_dp*pw_integral_ab(density, rhog)
1266 END IF
1267 END IF
1268 IF (PRESENT(h_stress)) THEN
1269 CALL pw_transfer(rhor, rhog)
1270 IF (PRESENT(aux_density)) THEN
1271 CALL pw_transfer(aux_density, rhor)
1272 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
1273 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
1274 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
1275 CALL pw_transfer(rhor, rhog_aux)
1276 END IF
1277 END IF
1278 CALL pw_pool%give_back_pw(rhor)
1279
1280 END SELECT
1281
1282 IF (PRESENT(aux_density)) THEN
1283 CALL calc_stress_and_gradient_c1d_gs (poisson_env, rhog, ehartree, rhog_aux, h_stress)
1284 ELSE
1285 CALL calc_stress_and_gradient_c1d_gs (poisson_env, rhog, ehartree, h_stress=h_stress)
1286 END IF
1287
1288 CALL pw_pool%give_back_pw(rhog)
1289 IF (PRESENT(aux_density)) THEN
1290 CALL pw_pool%give_back_pw(rhog_aux)
1291 END IF
1292
1293 CALL timestop(handle)
1294
1295 END SUBROUTINE pw_poisson_solve_v_nodv_c1d_gs_r3d_rs
1296! **************************************************************************************************
1297!> \brief Solve Poisson equation in a plane wave basis set
1298!> Obtains electrostatic potential and its derivatives with respect to r
1299!> from the density
1300!> \param poisson_env ...
1301!> \param density ...
1302!> \param ehartree ...
1303!> \param vhartree ...
1304!> \param h_stress ...
1305!> \param rho_core ...
1306!> \param greenfn ...
1307!> \param aux_density Hartree energy and stress tensor between 2 different densities
1308!> \par History
1309!> JGH (13-Mar-2001) : completely revised
1310!> \author apsi
1311! **************************************************************************************************
1312 SUBROUTINE pw_poisson_solve_v_nodv_c1d_gs_c1d_gs (poisson_env, density, ehartree, vhartree, &
1313 h_stress, rho_core, greenfn, aux_density)
1314
1315 TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
1316 TYPE(pw_c1d_gs_type), INTENT(IN) :: density
1317 REAL(kind=dp), INTENT(out), OPTIONAL :: ehartree
1318 TYPE(pw_c1d_gs_type), INTENT(INOUT) :: vhartree
1319 REAL(kind=dp), DIMENSION(3, 3), INTENT(OUT), &
1320 OPTIONAL :: h_stress
1321 TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_core, greenfn
1322 TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: aux_density
1323
1324 CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_solve'
1325
1326 INTEGER :: handle
1327 LOGICAL :: has_dielectric
1328 TYPE(pw_grid_type), POINTER :: pw_grid
1329 TYPE(pw_pool_type), POINTER :: pw_pool
1330 TYPE(pw_r3d_rs_type) :: &
1331 rhor, vhartree_rs
1332 TYPE(pw_c1d_gs_type) :: influence_fn, rhog, rhog_aux
1333
1334 CALL timeset(routinen, handle)
1335
1336 CALL pw_poisson_rebuild(poisson_env, density)
1337
1338 has_dielectric = poisson_env%parameters%has_dielectric
1339
1340 ! point pw
1341 pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
1342 pw_grid => pw_pool%pw_grid
1343 IF (.NOT. pw_grid_compare(pw_pool%pw_grid, vhartree%pw_grid)) THEN
1344 cpabort("vhartree has a different grid than the poisson solver")
1345 END IF
1346 ! density in G space
1347 CALL pw_pool%create_pw(rhog)
1348 IF (PRESENT(aux_density)) THEN
1349 CALL pw_pool%create_pw(rhog_aux)
1350 END IF
1351
1352 SELECT CASE (poisson_env%used_grid)
1353 CASE (use_gs_grid)
1354
1355 SELECT CASE (poisson_env%green_fft%method)
1356 CASE (periodic3d, analytic2d, analytic1d, analytic0d, mt2d, mt1d, mt0d, multipole0d)
1357
1358 CALL pw_transfer(density, rhog)
1359 IF (PRESENT(aux_density)) THEN
1360 CALL pw_transfer(aux_density, rhog_aux)
1361 END IF
1362 IF (PRESENT(greenfn)) THEN
1363 influence_fn = greenfn
1364 ELSE
1365 influence_fn = poisson_env%green_fft%influence_fn
1366 END IF
1367 CALL pw_multiply_with(rhog, influence_fn)
1368 IF (PRESENT(aux_density)) THEN
1369 CALL pw_multiply_with(rhog_aux, influence_fn)
1370 END IF
1371 CALL pw_transfer(rhog, vhartree)
1372 IF (PRESENT(ehartree)) THEN
1373 IF (PRESENT(aux_density)) THEN
1374 ehartree = 0.5_dp*pw_integral_ab(aux_density, vhartree)
1375 ELSE
1376 ehartree = 0.5_dp*pw_integral_ab(density, vhartree)
1377 END IF
1378 END IF
1379
1380 CASE (ps_implicit)
1381 IF (PRESENT(h_stress)) THEN
1382 cpabort("No stress tensor is implemented for the implicit Poisson solver.")
1383 END IF
1384
1385 IF (has_dielectric .AND. PRESENT(rho_core)) THEN
1386 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
1387 CASE (periodic_bc, mixed_periodic_bc)
1388 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
1389 poisson_env%diel_rs_grid, &
1390 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
1391 density, rho_core=rho_core)
1392 CASE (neumann_bc, mixed_bc)
1393 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
1394 poisson_env%diel_rs_grid, &
1395 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
1396 poisson_env%dct_pw_grid, &
1397 poisson_env%parameters%ps_implicit_params%neumann_directions, &
1398 poisson_env%implicit_env%dct_env%recv_msgs_bnds, &
1399 poisson_env%implicit_env%dct_env%dests_expand, &
1400 poisson_env%implicit_env%dct_env%srcs_expand, &
1401 poisson_env%implicit_env%dct_env%flipg_stat, &
1402 poisson_env%implicit_env%dct_env%bounds_shftd, &
1403 density, rho_core=rho_core)
1404 END SELECT
1405 END IF
1406
1407 CALL pw_pool%create_pw(rhor)
1408 CALL pw_pool%create_pw(vhartree_rs)
1409 CALL pw_transfer(density, rhor)
1410
1411 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
1412 CASE (periodic_bc)
1413 CALL implicit_poisson_solver_periodic(poisson_env, rhor, vhartree_rs, &
1414 ehartree=ehartree)
1415 CASE (neumann_bc)
1416 CALL implicit_poisson_solver_neumann(poisson_env, rhor, vhartree_rs, &
1417 ehartree=ehartree)
1418 CASE (mixed_periodic_bc)
1419 CALL implicit_poisson_solver_mixed_periodic(poisson_env, rhor, vhartree_rs, &
1420 electric_enthalpy=ehartree)
1421 CASE (mixed_bc)
1422 CALL implicit_poisson_solver_mixed(poisson_env, rhor, vhartree_rs, &
1423 electric_enthalpy=ehartree)
1424 END SELECT
1425
1426 IF (PRESENT(aux_density)) THEN
1427 CALL pw_transfer(aux_density, rhor)
1428 ehartree = 0.5_dp*pw_integral_ab(rhor, vhartree_rs)
1429 END IF
1430
1431 CALL pw_transfer(vhartree_rs, vhartree)
1432
1433 CALL pw_pool%give_back_pw(rhor)
1434 CALL pw_pool%give_back_pw(vhartree_rs)
1435
1436 CASE DEFAULT
1437 CALL cp_abort(__location__, &
1438 "unknown poisson method "// &
1439 cp_to_string(poisson_env%green_fft%method))
1440 END SELECT
1441
1442 CASE (use_rs_grid)
1443
1444 CALL pw_pool%create_pw(rhor)
1445 CALL pw_transfer(density, rhor)
1446 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
1447 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
1448 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
1449 CALL pw_transfer(rhor, vhartree)
1450 IF (PRESENT(ehartree)) THEN
1451 IF (PRESENT(aux_density)) THEN
1452 ehartree = 0.5_dp*pw_integral_ab(aux_density, vhartree)
1453 ELSE
1454 ehartree = 0.5_dp*pw_integral_ab(density, vhartree)
1455 END IF
1456 END IF
1457 IF (PRESENT(h_stress)) THEN
1458 CALL pw_transfer(rhor, rhog)
1459 IF (PRESENT(aux_density)) THEN
1460 CALL pw_transfer(aux_density, rhor)
1461 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
1462 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
1463 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
1464 CALL pw_transfer(rhor, rhog_aux)
1465 END IF
1466 END IF
1467 CALL pw_pool%give_back_pw(rhor)
1468
1469 END SELECT
1470
1471 IF (PRESENT(aux_density)) THEN
1472 CALL calc_stress_and_gradient_c1d_gs (poisson_env, rhog, ehartree, rhog_aux, h_stress)
1473 ELSE
1474 CALL calc_stress_and_gradient_c1d_gs (poisson_env, rhog, ehartree, h_stress=h_stress)
1475 END IF
1476
1477 CALL pw_pool%give_back_pw(rhog)
1478 IF (PRESENT(aux_density)) THEN
1479 CALL pw_pool%give_back_pw(rhog_aux)
1480 END IF
1481
1482 CALL timestop(handle)
1483
1484 END SUBROUTINE pw_poisson_solve_v_nodv_c1d_gs_c1d_gs
1485
1486! **************************************************************************************************
1487!> \brief Solve Poisson equation in a plane wave basis set
1488!> Obtains electrostatic potential and its derivatives with respect to r
1489!> from the density
1490!> \param poisson_env ...
1491!> \param density ...
1492!> \param ehartree ...
1493!> \param dvhartree ...
1494!> \param h_stress ...
1495!> \param rho_core ...
1496!> \param greenfn ...
1497!> \param aux_density Hartree energy and stress tensor between 2 different densities
1498!> \par History
1499!> JGH (13-Mar-2001) : completely revised
1500!> \author apsi
1501! **************************************************************************************************
1502 SUBROUTINE pw_poisson_solve_nov_dv_r3d_rs_r3d_rs (poisson_env, density, ehartree, &
1503 dvhartree, h_stress, rho_core, greenfn, aux_density)
1504
1505 TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
1506 TYPE(pw_r3d_rs_type), INTENT(IN) :: density
1507 REAL(kind=dp), INTENT(out), OPTIONAL :: ehartree
1508 TYPE(pw_r3d_rs_type), DIMENSION(3) :: dvhartree
1509 REAL(kind=dp), DIMENSION(3, 3), INTENT(OUT), &
1510 OPTIONAL :: h_stress
1511 TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_core, greenfn
1512 TYPE(pw_r3d_rs_type), INTENT(IN), OPTIONAL :: aux_density
1513
1514 CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_solve'
1515
1516 INTEGER :: handle
1517 LOGICAL :: has_dielectric
1518 TYPE(pw_grid_type), POINTER :: pw_grid
1519 TYPE(pw_pool_type), POINTER :: pw_pool
1520 TYPE(pw_r3d_rs_type) :: &
1521 rhor, vhartree_rs
1522 TYPE(pw_c1d_gs_type) :: influence_fn, rhog, rhog_aux, tmpg
1523
1524 CALL timeset(routinen, handle)
1525
1526 CALL pw_poisson_rebuild(poisson_env, density)
1527
1528 has_dielectric = poisson_env%parameters%has_dielectric
1529
1530 ! point pw
1531 pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
1532 pw_grid => pw_pool%pw_grid
1533 ! density in G space
1534 CALL pw_pool%create_pw(rhog)
1535 IF (PRESENT(aux_density)) THEN
1536 CALL pw_pool%create_pw(rhog_aux)
1537 END IF
1538
1539 SELECT CASE (poisson_env%used_grid)
1540 CASE (use_gs_grid)
1541
1542 SELECT CASE (poisson_env%green_fft%method)
1543 CASE (periodic3d, analytic2d, analytic1d, analytic0d, mt2d, mt1d, mt0d, multipole0d)
1544
1545 CALL pw_transfer(density, rhog)
1546 IF (PRESENT(aux_density)) THEN
1547 CALL pw_transfer(aux_density, rhog_aux)
1548 END IF
1549 IF (PRESENT(ehartree)) THEN
1550 CALL pw_pool%create_pw(tmpg)
1551 CALL pw_copy(rhog, tmpg)
1552 END IF
1553 IF (PRESENT(greenfn)) THEN
1554 influence_fn = greenfn
1555 ELSE
1556 influence_fn = poisson_env%green_fft%influence_fn
1557 END IF
1558 CALL pw_multiply_with(rhog, influence_fn)
1559 IF (PRESENT(aux_density)) THEN
1560 CALL pw_multiply_with(rhog_aux, influence_fn)
1561 rhog_aux%array(:) = rhog_aux%array(:)*influence_fn%array(:)
1562 END IF
1563 IF (PRESENT(ehartree)) THEN
1564 IF (PRESENT(aux_density)) THEN
1565 ehartree = 0.5_dp*pw_integral_ab(rhog_aux, tmpg)
1566 ELSE
1567 ehartree = 0.5_dp*pw_integral_ab(rhog, tmpg)
1568 END IF
1569 CALL pw_pool%give_back_pw(tmpg)
1570 END IF
1571
1572 CASE (ps_implicit)
1573 IF (PRESENT(h_stress)) THEN
1574 cpabort("No stress tensor is implemented for the implicit Poisson solver.")
1575 END IF
1576
1577 IF (has_dielectric .AND. PRESENT(rho_core)) THEN
1578 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
1579 CASE (periodic_bc, mixed_periodic_bc)
1580 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
1581 poisson_env%diel_rs_grid, &
1582 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
1583 density, rho_core=rho_core)
1584 CASE (neumann_bc, mixed_bc)
1585 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
1586 poisson_env%diel_rs_grid, &
1587 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
1588 poisson_env%dct_pw_grid, &
1589 poisson_env%parameters%ps_implicit_params%neumann_directions, &
1590 poisson_env%implicit_env%dct_env%recv_msgs_bnds, &
1591 poisson_env%implicit_env%dct_env%dests_expand, &
1592 poisson_env%implicit_env%dct_env%srcs_expand, &
1593 poisson_env%implicit_env%dct_env%flipg_stat, &
1594 poisson_env%implicit_env%dct_env%bounds_shftd, &
1595 density, rho_core=rho_core)
1596 END SELECT
1597 END IF
1598
1599 CALL pw_pool%create_pw(rhor)
1600 CALL pw_pool%create_pw(vhartree_rs)
1601 CALL pw_transfer(density, rhor)
1602
1603 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
1604 CASE (periodic_bc)
1605 CALL implicit_poisson_solver_periodic(poisson_env, rhor, vhartree_rs, &
1606 ehartree=ehartree)
1607 CASE (neumann_bc)
1608 CALL implicit_poisson_solver_neumann(poisson_env, rhor, vhartree_rs, &
1609 ehartree=ehartree)
1610 CASE (mixed_periodic_bc)
1611 CALL implicit_poisson_solver_mixed_periodic(poisson_env, rhor, vhartree_rs, &
1612 electric_enthalpy=ehartree)
1613 CASE (mixed_bc)
1614 CALL implicit_poisson_solver_mixed(poisson_env, rhor, vhartree_rs, &
1615 electric_enthalpy=ehartree)
1616 END SELECT
1617
1618 CALL pw_transfer(rhor, rhog)
1619
1620 IF (PRESENT(aux_density)) THEN
1621 CALL pw_transfer(aux_density, rhor)
1622 ehartree = 0.5_dp*pw_integral_ab(rhor, vhartree_rs)
1623 END IF
1624
1625 CALL pw_pool%give_back_pw(rhor)
1626 CALL pw_pool%give_back_pw(vhartree_rs)
1627
1628 CASE DEFAULT
1629 CALL cp_abort(__location__, &
1630 "unknown poisson method "// &
1631 cp_to_string(poisson_env%green_fft%method))
1632 END SELECT
1633
1634 CASE (use_rs_grid)
1635
1636 CALL pw_pool%create_pw(rhor)
1637 CALL pw_transfer(density, rhor)
1638 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
1639 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
1640 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
1641 CALL pw_transfer(rhor, rhog)
1642 IF (PRESENT(ehartree)) THEN
1643 ehartree = 0.5_dp*pw_integral_ab(density, rhor)
1644 END IF
1645 CALL pw_pool%give_back_pw(rhor)
1646
1647 END SELECT
1648
1649 IF (PRESENT(aux_density)) THEN
1650 CALL calc_stress_and_gradient_r3d_rs (poisson_env, rhog, ehartree, rhog_aux, h_stress, dvhartree=dvhartree)
1651 ELSE
1652 CALL calc_stress_and_gradient_r3d_rs (poisson_env, rhog, ehartree, h_stress=h_stress, dvhartree=dvhartree)
1653 END IF
1654
1655 CALL pw_pool%give_back_pw(rhog)
1656 IF (PRESENT(aux_density)) THEN
1657 CALL pw_pool%give_back_pw(rhog_aux)
1658 END IF
1659
1660 CALL timestop(handle)
1661
1662 END SUBROUTINE pw_poisson_solve_nov_dv_r3d_rs_r3d_rs
1663! **************************************************************************************************
1664!> \brief Solve Poisson equation in a plane wave basis set
1665!> Obtains electrostatic potential and its derivatives with respect to r
1666!> from the density
1667!> \param poisson_env ...
1668!> \param density ...
1669!> \param ehartree ...
1670!> \param dvhartree ...
1671!> \param h_stress ...
1672!> \param rho_core ...
1673!> \param greenfn ...
1674!> \param aux_density Hartree energy and stress tensor between 2 different densities
1675!> \par History
1676!> JGH (13-Mar-2001) : completely revised
1677!> \author apsi
1678! **************************************************************************************************
1679 SUBROUTINE pw_poisson_solve_nov_dv_r3d_rs_c1d_gs (poisson_env, density, ehartree, &
1680 dvhartree, h_stress, rho_core, greenfn, aux_density)
1681
1682 TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
1683 TYPE(pw_r3d_rs_type), INTENT(IN) :: density
1684 REAL(kind=dp), INTENT(out), OPTIONAL :: ehartree
1685 TYPE(pw_c1d_gs_type), DIMENSION(3) :: dvhartree
1686 REAL(kind=dp), DIMENSION(3, 3), INTENT(OUT), &
1687 OPTIONAL :: h_stress
1688 TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_core, greenfn
1689 TYPE(pw_r3d_rs_type), INTENT(IN), OPTIONAL :: aux_density
1690
1691 CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_solve'
1692
1693 INTEGER :: handle
1694 LOGICAL :: has_dielectric
1695 TYPE(pw_grid_type), POINTER :: pw_grid
1696 TYPE(pw_pool_type), POINTER :: pw_pool
1697 TYPE(pw_r3d_rs_type) :: &
1698 rhor, vhartree_rs
1699 TYPE(pw_c1d_gs_type) :: influence_fn, rhog, rhog_aux, tmpg
1700
1701 CALL timeset(routinen, handle)
1702
1703 CALL pw_poisson_rebuild(poisson_env, density)
1704
1705 has_dielectric = poisson_env%parameters%has_dielectric
1706
1707 ! point pw
1708 pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
1709 pw_grid => pw_pool%pw_grid
1710 ! density in G space
1711 CALL pw_pool%create_pw(rhog)
1712 IF (PRESENT(aux_density)) THEN
1713 CALL pw_pool%create_pw(rhog_aux)
1714 END IF
1715
1716 SELECT CASE (poisson_env%used_grid)
1717 CASE (use_gs_grid)
1718
1719 SELECT CASE (poisson_env%green_fft%method)
1720 CASE (periodic3d, analytic2d, analytic1d, analytic0d, mt2d, mt1d, mt0d, multipole0d)
1721
1722 CALL pw_transfer(density, rhog)
1723 IF (PRESENT(aux_density)) THEN
1724 CALL pw_transfer(aux_density, rhog_aux)
1725 END IF
1726 IF (PRESENT(ehartree)) THEN
1727 CALL pw_pool%create_pw(tmpg)
1728 CALL pw_copy(rhog, tmpg)
1729 END IF
1730 IF (PRESENT(greenfn)) THEN
1731 influence_fn = greenfn
1732 ELSE
1733 influence_fn = poisson_env%green_fft%influence_fn
1734 END IF
1735 CALL pw_multiply_with(rhog, influence_fn)
1736 IF (PRESENT(aux_density)) THEN
1737 CALL pw_multiply_with(rhog_aux, influence_fn)
1738 rhog_aux%array(:) = rhog_aux%array(:)*influence_fn%array(:)
1739 END IF
1740 IF (PRESENT(ehartree)) THEN
1741 IF (PRESENT(aux_density)) THEN
1742 ehartree = 0.5_dp*pw_integral_ab(rhog_aux, tmpg)
1743 ELSE
1744 ehartree = 0.5_dp*pw_integral_ab(rhog, tmpg)
1745 END IF
1746 CALL pw_pool%give_back_pw(tmpg)
1747 END IF
1748
1749 CASE (ps_implicit)
1750 IF (PRESENT(h_stress)) THEN
1751 cpabort("No stress tensor is implemented for the implicit Poisson solver.")
1752 END IF
1753
1754 IF (has_dielectric .AND. PRESENT(rho_core)) THEN
1755 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
1756 CASE (periodic_bc, mixed_periodic_bc)
1757 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
1758 poisson_env%diel_rs_grid, &
1759 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
1760 density, rho_core=rho_core)
1761 CASE (neumann_bc, mixed_bc)
1762 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
1763 poisson_env%diel_rs_grid, &
1764 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
1765 poisson_env%dct_pw_grid, &
1766 poisson_env%parameters%ps_implicit_params%neumann_directions, &
1767 poisson_env%implicit_env%dct_env%recv_msgs_bnds, &
1768 poisson_env%implicit_env%dct_env%dests_expand, &
1769 poisson_env%implicit_env%dct_env%srcs_expand, &
1770 poisson_env%implicit_env%dct_env%flipg_stat, &
1771 poisson_env%implicit_env%dct_env%bounds_shftd, &
1772 density, rho_core=rho_core)
1773 END SELECT
1774 END IF
1775
1776 CALL pw_pool%create_pw(rhor)
1777 CALL pw_pool%create_pw(vhartree_rs)
1778 CALL pw_transfer(density, rhor)
1779
1780 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
1781 CASE (periodic_bc)
1782 CALL implicit_poisson_solver_periodic(poisson_env, rhor, vhartree_rs, &
1783 ehartree=ehartree)
1784 CASE (neumann_bc)
1785 CALL implicit_poisson_solver_neumann(poisson_env, rhor, vhartree_rs, &
1786 ehartree=ehartree)
1787 CASE (mixed_periodic_bc)
1788 CALL implicit_poisson_solver_mixed_periodic(poisson_env, rhor, vhartree_rs, &
1789 electric_enthalpy=ehartree)
1790 CASE (mixed_bc)
1791 CALL implicit_poisson_solver_mixed(poisson_env, rhor, vhartree_rs, &
1792 electric_enthalpy=ehartree)
1793 END SELECT
1794
1795 CALL pw_transfer(rhor, rhog)
1796
1797 IF (PRESENT(aux_density)) THEN
1798 CALL pw_transfer(aux_density, rhor)
1799 ehartree = 0.5_dp*pw_integral_ab(rhor, vhartree_rs)
1800 END IF
1801
1802 CALL pw_pool%give_back_pw(rhor)
1803 CALL pw_pool%give_back_pw(vhartree_rs)
1804
1805 CASE DEFAULT
1806 CALL cp_abort(__location__, &
1807 "unknown poisson method "// &
1808 cp_to_string(poisson_env%green_fft%method))
1809 END SELECT
1810
1811 CASE (use_rs_grid)
1812
1813 CALL pw_pool%create_pw(rhor)
1814 CALL pw_transfer(density, rhor)
1815 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
1816 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
1817 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
1818 CALL pw_transfer(rhor, rhog)
1819 IF (PRESENT(ehartree)) THEN
1820 ehartree = 0.5_dp*pw_integral_ab(density, rhor)
1821 END IF
1822 CALL pw_pool%give_back_pw(rhor)
1823
1824 END SELECT
1825
1826 IF (PRESENT(aux_density)) THEN
1827 CALL calc_stress_and_gradient_c1d_gs (poisson_env, rhog, ehartree, rhog_aux, h_stress, dvhartree=dvhartree)
1828 ELSE
1829 CALL calc_stress_and_gradient_c1d_gs (poisson_env, rhog, ehartree, h_stress=h_stress, dvhartree=dvhartree)
1830 END IF
1831
1832 CALL pw_pool%give_back_pw(rhog)
1833 IF (PRESENT(aux_density)) THEN
1834 CALL pw_pool%give_back_pw(rhog_aux)
1835 END IF
1836
1837 CALL timestop(handle)
1838
1839 END SUBROUTINE pw_poisson_solve_nov_dv_r3d_rs_c1d_gs
1840! **************************************************************************************************
1841!> \brief Solve Poisson equation in a plane wave basis set
1842!> Obtains electrostatic potential and its derivatives with respect to r
1843!> from the density
1844!> \param poisson_env ...
1845!> \param density ...
1846!> \param ehartree ...
1847!> \param dvhartree ...
1848!> \param h_stress ...
1849!> \param rho_core ...
1850!> \param greenfn ...
1851!> \param aux_density Hartree energy and stress tensor between 2 different densities
1852!> \par History
1853!> JGH (13-Mar-2001) : completely revised
1854!> \author apsi
1855! **************************************************************************************************
1856 SUBROUTINE pw_poisson_solve_nov_dv_c1d_gs_r3d_rs (poisson_env, density, ehartree, &
1857 dvhartree, h_stress, rho_core, greenfn, aux_density)
1858
1859 TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
1860 TYPE(pw_c1d_gs_type), INTENT(IN) :: density
1861 REAL(kind=dp), INTENT(out), OPTIONAL :: ehartree
1862 TYPE(pw_r3d_rs_type), DIMENSION(3) :: dvhartree
1863 REAL(kind=dp), DIMENSION(3, 3), INTENT(OUT), &
1864 OPTIONAL :: h_stress
1865 TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_core, greenfn
1866 TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: aux_density
1867
1868 CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_solve'
1869
1870 INTEGER :: handle
1871 LOGICAL :: has_dielectric
1872 TYPE(pw_grid_type), POINTER :: pw_grid
1873 TYPE(pw_pool_type), POINTER :: pw_pool
1874 TYPE(pw_r3d_rs_type) :: &
1875 rhor, vhartree_rs
1876 TYPE(pw_c1d_gs_type) :: influence_fn, rhog, rhog_aux, tmpg
1877
1878 CALL timeset(routinen, handle)
1879
1880 CALL pw_poisson_rebuild(poisson_env, density)
1881
1882 has_dielectric = poisson_env%parameters%has_dielectric
1883
1884 ! point pw
1885 pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
1886 pw_grid => pw_pool%pw_grid
1887 ! density in G space
1888 CALL pw_pool%create_pw(rhog)
1889 IF (PRESENT(aux_density)) THEN
1890 CALL pw_pool%create_pw(rhog_aux)
1891 END IF
1892
1893 SELECT CASE (poisson_env%used_grid)
1894 CASE (use_gs_grid)
1895
1896 SELECT CASE (poisson_env%green_fft%method)
1897 CASE (periodic3d, analytic2d, analytic1d, analytic0d, mt2d, mt1d, mt0d, multipole0d)
1898
1899 CALL pw_transfer(density, rhog)
1900 IF (PRESENT(aux_density)) THEN
1901 CALL pw_transfer(aux_density, rhog_aux)
1902 END IF
1903 IF (PRESENT(ehartree)) THEN
1904 CALL pw_pool%create_pw(tmpg)
1905 CALL pw_copy(rhog, tmpg)
1906 END IF
1907 IF (PRESENT(greenfn)) THEN
1908 influence_fn = greenfn
1909 ELSE
1910 influence_fn = poisson_env%green_fft%influence_fn
1911 END IF
1912 CALL pw_multiply_with(rhog, influence_fn)
1913 IF (PRESENT(aux_density)) THEN
1914 CALL pw_multiply_with(rhog_aux, influence_fn)
1915 rhog_aux%array(:) = rhog_aux%array(:)*influence_fn%array(:)
1916 END IF
1917 IF (PRESENT(ehartree)) THEN
1918 IF (PRESENT(aux_density)) THEN
1919 ehartree = 0.5_dp*pw_integral_ab(rhog_aux, tmpg)
1920 ELSE
1921 ehartree = 0.5_dp*pw_integral_ab(rhog, tmpg)
1922 END IF
1923 CALL pw_pool%give_back_pw(tmpg)
1924 END IF
1925
1926 CASE (ps_implicit)
1927 IF (PRESENT(h_stress)) THEN
1928 cpabort("No stress tensor is implemented for the implicit Poisson solver.")
1929 END IF
1930
1931 IF (has_dielectric .AND. PRESENT(rho_core)) THEN
1932 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
1933 CASE (periodic_bc, mixed_periodic_bc)
1934 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
1935 poisson_env%diel_rs_grid, &
1936 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
1937 density, rho_core=rho_core)
1938 CASE (neumann_bc, mixed_bc)
1939 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
1940 poisson_env%diel_rs_grid, &
1941 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
1942 poisson_env%dct_pw_grid, &
1943 poisson_env%parameters%ps_implicit_params%neumann_directions, &
1944 poisson_env%implicit_env%dct_env%recv_msgs_bnds, &
1945 poisson_env%implicit_env%dct_env%dests_expand, &
1946 poisson_env%implicit_env%dct_env%srcs_expand, &
1947 poisson_env%implicit_env%dct_env%flipg_stat, &
1948 poisson_env%implicit_env%dct_env%bounds_shftd, &
1949 density, rho_core=rho_core)
1950 END SELECT
1951 END IF
1952
1953 CALL pw_pool%create_pw(rhor)
1954 CALL pw_pool%create_pw(vhartree_rs)
1955 CALL pw_transfer(density, rhor)
1956
1957 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
1958 CASE (periodic_bc)
1959 CALL implicit_poisson_solver_periodic(poisson_env, rhor, vhartree_rs, &
1960 ehartree=ehartree)
1961 CASE (neumann_bc)
1962 CALL implicit_poisson_solver_neumann(poisson_env, rhor, vhartree_rs, &
1963 ehartree=ehartree)
1964 CASE (mixed_periodic_bc)
1965 CALL implicit_poisson_solver_mixed_periodic(poisson_env, rhor, vhartree_rs, &
1966 electric_enthalpy=ehartree)
1967 CASE (mixed_bc)
1968 CALL implicit_poisson_solver_mixed(poisson_env, rhor, vhartree_rs, &
1969 electric_enthalpy=ehartree)
1970 END SELECT
1971
1972 CALL pw_transfer(rhor, rhog)
1973
1974 IF (PRESENT(aux_density)) THEN
1975 CALL pw_transfer(aux_density, rhor)
1976 ehartree = 0.5_dp*pw_integral_ab(rhor, vhartree_rs)
1977 END IF
1978
1979 CALL pw_pool%give_back_pw(rhor)
1980 CALL pw_pool%give_back_pw(vhartree_rs)
1981
1982 CASE DEFAULT
1983 CALL cp_abort(__location__, &
1984 "unknown poisson method "// &
1985 cp_to_string(poisson_env%green_fft%method))
1986 END SELECT
1987
1988 CASE (use_rs_grid)
1989
1990 CALL pw_pool%create_pw(rhor)
1991 CALL pw_transfer(density, rhor)
1992 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
1993 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
1994 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
1995 CALL pw_transfer(rhor, rhog)
1996 IF (PRESENT(ehartree)) THEN
1997 ehartree = 0.5_dp*pw_integral_ab(density, rhog)
1998 END IF
1999 CALL pw_pool%give_back_pw(rhor)
2000
2001 END SELECT
2002
2003 IF (PRESENT(aux_density)) THEN
2004 CALL calc_stress_and_gradient_r3d_rs (poisson_env, rhog, ehartree, rhog_aux, h_stress, dvhartree=dvhartree)
2005 ELSE
2006 CALL calc_stress_and_gradient_r3d_rs (poisson_env, rhog, ehartree, h_stress=h_stress, dvhartree=dvhartree)
2007 END IF
2008
2009 CALL pw_pool%give_back_pw(rhog)
2010 IF (PRESENT(aux_density)) THEN
2011 CALL pw_pool%give_back_pw(rhog_aux)
2012 END IF
2013
2014 CALL timestop(handle)
2015
2016 END SUBROUTINE pw_poisson_solve_nov_dv_c1d_gs_r3d_rs
2017! **************************************************************************************************
2018!> \brief Solve Poisson equation in a plane wave basis set
2019!> Obtains electrostatic potential and its derivatives with respect to r
2020!> from the density
2021!> \param poisson_env ...
2022!> \param density ...
2023!> \param ehartree ...
2024!> \param dvhartree ...
2025!> \param h_stress ...
2026!> \param rho_core ...
2027!> \param greenfn ...
2028!> \param aux_density Hartree energy and stress tensor between 2 different densities
2029!> \par History
2030!> JGH (13-Mar-2001) : completely revised
2031!> \author apsi
2032! **************************************************************************************************
2033 SUBROUTINE pw_poisson_solve_nov_dv_c1d_gs_c1d_gs (poisson_env, density, ehartree, &
2034 dvhartree, h_stress, rho_core, greenfn, aux_density)
2035
2036 TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
2037 TYPE(pw_c1d_gs_type), INTENT(IN) :: density
2038 REAL(kind=dp), INTENT(out), OPTIONAL :: ehartree
2039 TYPE(pw_c1d_gs_type), DIMENSION(3) :: dvhartree
2040 REAL(kind=dp), DIMENSION(3, 3), INTENT(OUT), &
2041 OPTIONAL :: h_stress
2042 TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_core, greenfn
2043 TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: aux_density
2044
2045 CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_solve'
2046
2047 INTEGER :: handle
2048 LOGICAL :: has_dielectric
2049 TYPE(pw_grid_type), POINTER :: pw_grid
2050 TYPE(pw_pool_type), POINTER :: pw_pool
2051 TYPE(pw_r3d_rs_type) :: &
2052 rhor, vhartree_rs
2053 TYPE(pw_c1d_gs_type) :: influence_fn, rhog, rhog_aux, tmpg
2054
2055 CALL timeset(routinen, handle)
2056
2057 CALL pw_poisson_rebuild(poisson_env, density)
2058
2059 has_dielectric = poisson_env%parameters%has_dielectric
2060
2061 ! point pw
2062 pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
2063 pw_grid => pw_pool%pw_grid
2064 ! density in G space
2065 CALL pw_pool%create_pw(rhog)
2066 IF (PRESENT(aux_density)) THEN
2067 CALL pw_pool%create_pw(rhog_aux)
2068 END IF
2069
2070 SELECT CASE (poisson_env%used_grid)
2071 CASE (use_gs_grid)
2072
2073 SELECT CASE (poisson_env%green_fft%method)
2074 CASE (periodic3d, analytic2d, analytic1d, analytic0d, mt2d, mt1d, mt0d, multipole0d)
2075
2076 CALL pw_transfer(density, rhog)
2077 IF (PRESENT(aux_density)) THEN
2078 CALL pw_transfer(aux_density, rhog_aux)
2079 END IF
2080 IF (PRESENT(ehartree)) THEN
2081 CALL pw_pool%create_pw(tmpg)
2082 CALL pw_copy(rhog, tmpg)
2083 END IF
2084 IF (PRESENT(greenfn)) THEN
2085 influence_fn = greenfn
2086 ELSE
2087 influence_fn = poisson_env%green_fft%influence_fn
2088 END IF
2089 CALL pw_multiply_with(rhog, influence_fn)
2090 IF (PRESENT(aux_density)) THEN
2091 CALL pw_multiply_with(rhog_aux, influence_fn)
2092 rhog_aux%array(:) = rhog_aux%array(:)*influence_fn%array(:)
2093 END IF
2094 IF (PRESENT(ehartree)) THEN
2095 IF (PRESENT(aux_density)) THEN
2096 ehartree = 0.5_dp*pw_integral_ab(rhog_aux, tmpg)
2097 ELSE
2098 ehartree = 0.5_dp*pw_integral_ab(rhog, tmpg)
2099 END IF
2100 CALL pw_pool%give_back_pw(tmpg)
2101 END IF
2102
2103 CASE (ps_implicit)
2104 IF (PRESENT(h_stress)) THEN
2105 cpabort("No stress tensor is implemented for the implicit Poisson solver.")
2106 END IF
2107
2108 IF (has_dielectric .AND. PRESENT(rho_core)) THEN
2109 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
2110 CASE (periodic_bc, mixed_periodic_bc)
2111 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
2112 poisson_env%diel_rs_grid, &
2113 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
2114 density, rho_core=rho_core)
2115 CASE (neumann_bc, mixed_bc)
2116 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
2117 poisson_env%diel_rs_grid, &
2118 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
2119 poisson_env%dct_pw_grid, &
2120 poisson_env%parameters%ps_implicit_params%neumann_directions, &
2121 poisson_env%implicit_env%dct_env%recv_msgs_bnds, &
2122 poisson_env%implicit_env%dct_env%dests_expand, &
2123 poisson_env%implicit_env%dct_env%srcs_expand, &
2124 poisson_env%implicit_env%dct_env%flipg_stat, &
2125 poisson_env%implicit_env%dct_env%bounds_shftd, &
2126 density, rho_core=rho_core)
2127 END SELECT
2128 END IF
2129
2130 CALL pw_pool%create_pw(rhor)
2131 CALL pw_pool%create_pw(vhartree_rs)
2132 CALL pw_transfer(density, rhor)
2133
2134 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
2135 CASE (periodic_bc)
2136 CALL implicit_poisson_solver_periodic(poisson_env, rhor, vhartree_rs, &
2137 ehartree=ehartree)
2138 CASE (neumann_bc)
2139 CALL implicit_poisson_solver_neumann(poisson_env, rhor, vhartree_rs, &
2140 ehartree=ehartree)
2141 CASE (mixed_periodic_bc)
2142 CALL implicit_poisson_solver_mixed_periodic(poisson_env, rhor, vhartree_rs, &
2143 electric_enthalpy=ehartree)
2144 CASE (mixed_bc)
2145 CALL implicit_poisson_solver_mixed(poisson_env, rhor, vhartree_rs, &
2146 electric_enthalpy=ehartree)
2147 END SELECT
2148
2149 CALL pw_transfer(rhor, rhog)
2150
2151 IF (PRESENT(aux_density)) THEN
2152 CALL pw_transfer(aux_density, rhor)
2153 ehartree = 0.5_dp*pw_integral_ab(rhor, vhartree_rs)
2154 END IF
2155
2156 CALL pw_pool%give_back_pw(rhor)
2157 CALL pw_pool%give_back_pw(vhartree_rs)
2158
2159 CASE DEFAULT
2160 CALL cp_abort(__location__, &
2161 "unknown poisson method "// &
2162 cp_to_string(poisson_env%green_fft%method))
2163 END SELECT
2164
2165 CASE (use_rs_grid)
2166
2167 CALL pw_pool%create_pw(rhor)
2168 CALL pw_transfer(density, rhor)
2169 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
2170 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
2171 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
2172 CALL pw_transfer(rhor, rhog)
2173 IF (PRESENT(ehartree)) THEN
2174 ehartree = 0.5_dp*pw_integral_ab(density, rhog)
2175 END IF
2176 CALL pw_pool%give_back_pw(rhor)
2177
2178 END SELECT
2179
2180 IF (PRESENT(aux_density)) THEN
2181 CALL calc_stress_and_gradient_c1d_gs (poisson_env, rhog, ehartree, rhog_aux, h_stress, dvhartree=dvhartree)
2182 ELSE
2183 CALL calc_stress_and_gradient_c1d_gs (poisson_env, rhog, ehartree, h_stress=h_stress, dvhartree=dvhartree)
2184 END IF
2185
2186 CALL pw_pool%give_back_pw(rhog)
2187 IF (PRESENT(aux_density)) THEN
2188 CALL pw_pool%give_back_pw(rhog_aux)
2189 END IF
2190
2191 CALL timestop(handle)
2192
2193 END SUBROUTINE pw_poisson_solve_nov_dv_c1d_gs_c1d_gs
2194
2195! **************************************************************************************************
2196!> \brief Solve Poisson equation in a plane wave basis set
2197!> Obtains electrostatic potential and its derivatives with respect to r
2198!> from the density
2199!> \param poisson_env ...
2200!> \param density ...
2201!> \param ehartree ...
2202!> \param vhartree ...
2203!> \param dvhartree ...
2204!> \param h_stress ...
2205!> \param rho_core ...
2206!> \param greenfn ...
2207!> \param aux_density Hartree energy and stress tensor between 2 different densities
2208!> \par History
2209!> JGH (13-Mar-2001) : completely revised
2210!> \author apsi
2211! **************************************************************************************************
2212 SUBROUTINE pw_poisson_solve_v_dv_r3d_rs_r3d_rs_r3d_rs (poisson_env, density, ehartree, vhartree, &
2213 dvhartree, h_stress, rho_core, greenfn, aux_density)
2214
2215 TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
2216 TYPE(pw_r3d_rs_type), INTENT(IN) :: density
2217 REAL(kind=dp), INTENT(out), OPTIONAL :: ehartree
2218 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: vhartree
2219 TYPE(pw_r3d_rs_type), DIMENSION(3) :: dvhartree
2220 REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT), &
2221 OPTIONAL :: h_stress
2222 TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_core, greenfn
2223 TYPE(pw_r3d_rs_type), INTENT(IN), OPTIONAL :: aux_density
2224
2225 CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_solve'
2226
2227 INTEGER :: handle
2228 LOGICAL :: has_dielectric
2229 TYPE(pw_grid_type), POINTER :: pw_grid
2230 TYPE(pw_pool_type), POINTER :: pw_pool
2231 TYPE(pw_r3d_rs_type) :: rhor, vhartree_rs
2232 TYPE(pw_c1d_gs_type) :: influence_fn, rhog, rhog_aux
2233
2234 CALL timeset(routinen, handle)
2235
2236 CALL pw_poisson_rebuild(poisson_env, density)
2237
2238 has_dielectric = poisson_env%parameters%has_dielectric
2239
2240 ! point pw
2241 pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
2242 pw_grid => pw_pool%pw_grid
2243 IF (.NOT. pw_grid_compare(pw_pool%pw_grid, vhartree%pw_grid)) THEN
2244 cpabort("vhartree has a different grid than the poisson solver")
2245 END IF
2246 ! density in G space
2247 CALL pw_pool%create_pw(rhog)
2248 IF (PRESENT(aux_density)) THEN
2249 CALL pw_pool%create_pw(rhog_aux)
2250 END IF
2251
2252 SELECT CASE (poisson_env%used_grid)
2253 CASE (use_gs_grid)
2254
2255 SELECT CASE (poisson_env%green_fft%method)
2256 CASE (periodic3d, analytic2d, analytic1d, analytic0d, mt2d, mt1d, mt0d, multipole0d)
2257
2258 CALL pw_transfer(density, rhog)
2259 IF (PRESENT(aux_density)) THEN
2260 CALL pw_transfer(aux_density, rhog_aux)
2261 END IF
2262 IF (PRESENT(greenfn)) THEN
2263 influence_fn = greenfn
2264 ELSE
2265 influence_fn = poisson_env%green_fft%influence_fn
2266 END IF
2267 CALL pw_multiply_with(rhog, influence_fn)
2268 IF (PRESENT(aux_density)) THEN
2269 CALL pw_multiply_with(rhog_aux, influence_fn)
2270 END IF
2271 CALL pw_transfer(rhog, vhartree)
2272 IF (PRESENT(ehartree)) THEN
2273 IF (PRESENT(aux_density)) THEN
2274 ehartree = 0.5_dp*pw_integral_ab(aux_density, vhartree)
2275 ELSE
2276 ehartree = 0.5_dp*pw_integral_ab(density, vhartree)
2277 END IF
2278 END IF
2279
2280 CASE (ps_implicit)
2281 IF (PRESENT(h_stress)) THEN
2282 CALL cp_abort(__location__, &
2283 "No stress tensor is implemented for the implicit Poisson solver.")
2284 END IF
2285
2286 IF (has_dielectric .AND. PRESENT(rho_core)) THEN
2287 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
2288 CASE (periodic_bc, mixed_periodic_bc)
2289 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
2290 poisson_env%diel_rs_grid, &
2291 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
2292 density, rho_core=rho_core)
2293 CASE (neumann_bc, mixed_bc)
2294 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
2295 poisson_env%diel_rs_grid, &
2296 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
2297 poisson_env%dct_pw_grid, &
2298 poisson_env%parameters%ps_implicit_params%neumann_directions, &
2299 poisson_env%implicit_env%dct_env%recv_msgs_bnds, &
2300 poisson_env%implicit_env%dct_env%dests_expand, &
2301 poisson_env%implicit_env%dct_env%srcs_expand, &
2302 poisson_env%implicit_env%dct_env%flipg_stat, &
2303 poisson_env%implicit_env%dct_env%bounds_shftd, &
2304 density, rho_core=rho_core)
2305 END SELECT
2306 END IF
2307
2308 CALL pw_pool%create_pw(rhor)
2309 CALL pw_pool%create_pw(vhartree_rs)
2310 CALL pw_transfer(density, rhor)
2311
2312 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
2313 CASE (periodic_bc)
2314 CALL implicit_poisson_solver_periodic(poisson_env, rhor, vhartree_rs, &
2315 ehartree=ehartree)
2316 CASE (neumann_bc)
2317 CALL implicit_poisson_solver_neumann(poisson_env, rhor, vhartree_rs, &
2318 ehartree=ehartree)
2319 CASE (mixed_periodic_bc)
2320 CALL implicit_poisson_solver_mixed_periodic(poisson_env, rhor, vhartree_rs, &
2321 electric_enthalpy=ehartree)
2322 CASE (mixed_bc)
2323 CALL implicit_poisson_solver_mixed(poisson_env, rhor, vhartree_rs, &
2324 electric_enthalpy=ehartree)
2325 END SELECT
2326
2327 CALL pw_transfer(vhartree_rs, vhartree)
2328 CALL pw_transfer(rhor, rhog)
2329
2330 IF (PRESENT(aux_density)) THEN
2331 CALL pw_transfer(aux_density, rhor)
2332 ehartree = 0.5_dp*pw_integral_ab(rhor, vhartree_rs)
2333 END IF
2334
2335 CALL pw_pool%give_back_pw(rhor)
2336 CALL pw_pool%give_back_pw(vhartree_rs)
2337
2338 CASE DEFAULT
2339 CALL cp_abort(__location__, &
2340 "unknown poisson method "// &
2341 cp_to_string(poisson_env%green_fft%method))
2342 END SELECT
2343
2344 CASE (use_rs_grid)
2345
2346 CALL pw_pool%create_pw(rhor)
2347 CALL pw_transfer(density, rhor)
2348 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
2349 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
2350 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
2351 CALL pw_transfer(rhor, vhartree)
2352 CALL pw_transfer(rhor, rhog)
2353 IF (PRESENT(ehartree)) THEN
2354 IF (PRESENT(aux_density)) THEN
2355 ehartree = 0.5_dp*pw_integral_ab(aux_density, vhartree)
2356 ELSE
2357 ehartree = 0.5_dp*pw_integral_ab(density, vhartree)
2358 END IF
2359 END IF
2360 CALL pw_transfer(rhor, rhog)
2361 IF (PRESENT(aux_density)) THEN
2362 CALL pw_transfer(aux_density, rhor)
2363 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
2364 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
2365 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
2366 CALL pw_transfer(rhor, rhog_aux)
2367 END IF
2368 CALL pw_pool%give_back_pw(rhor)
2369
2370 END SELECT
2371
2372 IF (PRESENT(aux_density)) THEN
2373 CALL calc_stress_and_gradient_r3d_rs (poisson_env, rhog, ehartree, rhog_aux, h_stress, dvhartree=dvhartree)
2374 ELSE
2375 CALL calc_stress_and_gradient_r3d_rs (poisson_env, rhog, ehartree, h_stress=h_stress, dvhartree=dvhartree)
2376 END IF
2377
2378 CALL pw_pool%give_back_pw(rhog)
2379 IF (PRESENT(aux_density)) THEN
2380 CALL pw_pool%give_back_pw(rhog_aux)
2381 END IF
2382
2383 CALL timestop(handle)
2384
2385 END SUBROUTINE pw_poisson_solve_v_dv_r3d_rs_r3d_rs_r3d_rs
2386! **************************************************************************************************
2387!> \brief Solve Poisson equation in a plane wave basis set
2388!> Obtains electrostatic potential and its derivatives with respect to r
2389!> from the density
2390!> \param poisson_env ...
2391!> \param density ...
2392!> \param ehartree ...
2393!> \param vhartree ...
2394!> \param dvhartree ...
2395!> \param h_stress ...
2396!> \param rho_core ...
2397!> \param greenfn ...
2398!> \param aux_density Hartree energy and stress tensor between 2 different densities
2399!> \par History
2400!> JGH (13-Mar-2001) : completely revised
2401!> \author apsi
2402! **************************************************************************************************
2403 SUBROUTINE pw_poisson_solve_v_dv_r3d_rs_c1d_gs_r3d_rs (poisson_env, density, ehartree, vhartree, &
2404 dvhartree, h_stress, rho_core, greenfn, aux_density)
2405
2406 TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
2407 TYPE(pw_r3d_rs_type), INTENT(IN) :: density
2408 REAL(kind=dp), INTENT(out), OPTIONAL :: ehartree
2409 TYPE(pw_c1d_gs_type), INTENT(INOUT) :: vhartree
2410 TYPE(pw_r3d_rs_type), DIMENSION(3) :: dvhartree
2411 REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT), &
2412 OPTIONAL :: h_stress
2413 TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_core, greenfn
2414 TYPE(pw_r3d_rs_type), INTENT(IN), OPTIONAL :: aux_density
2415
2416 CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_solve'
2417
2418 INTEGER :: handle
2419 LOGICAL :: has_dielectric
2420 TYPE(pw_grid_type), POINTER :: pw_grid
2421 TYPE(pw_pool_type), POINTER :: pw_pool
2422 TYPE(pw_r3d_rs_type) :: rhor, vhartree_rs
2423 TYPE(pw_c1d_gs_type) :: influence_fn, rhog, rhog_aux
2424
2425 CALL timeset(routinen, handle)
2426
2427 CALL pw_poisson_rebuild(poisson_env, density)
2428
2429 has_dielectric = poisson_env%parameters%has_dielectric
2430
2431 ! point pw
2432 pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
2433 pw_grid => pw_pool%pw_grid
2434 IF (.NOT. pw_grid_compare(pw_pool%pw_grid, vhartree%pw_grid)) THEN
2435 cpabort("vhartree has a different grid than the poisson solver")
2436 END IF
2437 ! density in G space
2438 CALL pw_pool%create_pw(rhog)
2439 IF (PRESENT(aux_density)) THEN
2440 CALL pw_pool%create_pw(rhog_aux)
2441 END IF
2442
2443 SELECT CASE (poisson_env%used_grid)
2444 CASE (use_gs_grid)
2445
2446 SELECT CASE (poisson_env%green_fft%method)
2447 CASE (periodic3d, analytic2d, analytic1d, analytic0d, mt2d, mt1d, mt0d, multipole0d)
2448
2449 CALL pw_transfer(density, rhog)
2450 IF (PRESENT(aux_density)) THEN
2451 CALL pw_transfer(aux_density, rhog_aux)
2452 END IF
2453 IF (PRESENT(greenfn)) THEN
2454 influence_fn = greenfn
2455 ELSE
2456 influence_fn = poisson_env%green_fft%influence_fn
2457 END IF
2458 CALL pw_multiply_with(rhog, influence_fn)
2459 IF (PRESENT(aux_density)) THEN
2460 CALL pw_multiply_with(rhog_aux, influence_fn)
2461 END IF
2462 CALL pw_transfer(rhog, vhartree)
2463 IF (PRESENT(ehartree)) THEN
2464 IF (PRESENT(aux_density)) THEN
2465 CALL pw_transfer(aux_density, rhog)
2466 ehartree = 0.5_dp*pw_integral_ab(rhog, vhartree)
2467 ELSE
2468 CALL pw_transfer(density, rhog)
2469 ehartree = 0.5_dp*pw_integral_ab(rhog, vhartree)
2470 END IF
2471 END IF
2472
2473 CASE (ps_implicit)
2474 IF (PRESENT(h_stress)) THEN
2475 CALL cp_abort(__location__, &
2476 "No stress tensor is implemented for the implicit Poisson solver.")
2477 END IF
2478
2479 IF (has_dielectric .AND. PRESENT(rho_core)) THEN
2480 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
2481 CASE (periodic_bc, mixed_periodic_bc)
2482 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
2483 poisson_env%diel_rs_grid, &
2484 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
2485 density, rho_core=rho_core)
2486 CASE (neumann_bc, mixed_bc)
2487 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
2488 poisson_env%diel_rs_grid, &
2489 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
2490 poisson_env%dct_pw_grid, &
2491 poisson_env%parameters%ps_implicit_params%neumann_directions, &
2492 poisson_env%implicit_env%dct_env%recv_msgs_bnds, &
2493 poisson_env%implicit_env%dct_env%dests_expand, &
2494 poisson_env%implicit_env%dct_env%srcs_expand, &
2495 poisson_env%implicit_env%dct_env%flipg_stat, &
2496 poisson_env%implicit_env%dct_env%bounds_shftd, &
2497 density, rho_core=rho_core)
2498 END SELECT
2499 END IF
2500
2501 CALL pw_pool%create_pw(rhor)
2502 CALL pw_pool%create_pw(vhartree_rs)
2503 CALL pw_transfer(density, rhor)
2504
2505 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
2506 CASE (periodic_bc)
2507 CALL implicit_poisson_solver_periodic(poisson_env, rhor, vhartree_rs, &
2508 ehartree=ehartree)
2509 CASE (neumann_bc)
2510 CALL implicit_poisson_solver_neumann(poisson_env, rhor, vhartree_rs, &
2511 ehartree=ehartree)
2512 CASE (mixed_periodic_bc)
2513 CALL implicit_poisson_solver_mixed_periodic(poisson_env, rhor, vhartree_rs, &
2514 electric_enthalpy=ehartree)
2515 CASE (mixed_bc)
2516 CALL implicit_poisson_solver_mixed(poisson_env, rhor, vhartree_rs, &
2517 electric_enthalpy=ehartree)
2518 END SELECT
2519
2520 CALL pw_transfer(vhartree_rs, vhartree)
2521 CALL pw_transfer(rhor, rhog)
2522
2523 IF (PRESENT(aux_density)) THEN
2524 CALL pw_transfer(aux_density, rhor)
2525 ehartree = 0.5_dp*pw_integral_ab(rhor, vhartree_rs)
2526 END IF
2527
2528 CALL pw_pool%give_back_pw(rhor)
2529 CALL pw_pool%give_back_pw(vhartree_rs)
2530
2531 CASE DEFAULT
2532 CALL cp_abort(__location__, &
2533 "unknown poisson method "// &
2534 cp_to_string(poisson_env%green_fft%method))
2535 END SELECT
2536
2537 CASE (use_rs_grid)
2538
2539 CALL pw_pool%create_pw(rhor)
2540 CALL pw_transfer(density, rhor)
2541 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
2542 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
2543 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
2544 CALL pw_transfer(rhor, vhartree)
2545 CALL pw_transfer(rhor, rhog)
2546 IF (PRESENT(ehartree)) THEN
2547 IF (PRESENT(aux_density)) THEN
2548 ehartree = 0.5_dp*pw_integral_ab(aux_density, rhor)
2549 ELSE
2550 ehartree = 0.5_dp*pw_integral_ab(density, rhor)
2551 END IF
2552 END IF
2553 CALL pw_transfer(rhor, rhog)
2554 IF (PRESENT(aux_density)) THEN
2555 CALL pw_transfer(aux_density, rhor)
2556 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
2557 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
2558 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
2559 CALL pw_transfer(rhor, rhog_aux)
2560 END IF
2561 CALL pw_pool%give_back_pw(rhor)
2562
2563 END SELECT
2564
2565 IF (PRESENT(aux_density)) THEN
2566 CALL calc_stress_and_gradient_r3d_rs (poisson_env, rhog, ehartree, rhog_aux, h_stress, dvhartree=dvhartree)
2567 ELSE
2568 CALL calc_stress_and_gradient_r3d_rs (poisson_env, rhog, ehartree, h_stress=h_stress, dvhartree=dvhartree)
2569 END IF
2570
2571 CALL pw_pool%give_back_pw(rhog)
2572 IF (PRESENT(aux_density)) THEN
2573 CALL pw_pool%give_back_pw(rhog_aux)
2574 END IF
2575
2576 CALL timestop(handle)
2577
2578 END SUBROUTINE pw_poisson_solve_v_dv_r3d_rs_c1d_gs_r3d_rs
2579! **************************************************************************************************
2580!> \brief Solve Poisson equation in a plane wave basis set
2581!> Obtains electrostatic potential and its derivatives with respect to r
2582!> from the density
2583!> \param poisson_env ...
2584!> \param density ...
2585!> \param ehartree ...
2586!> \param vhartree ...
2587!> \param dvhartree ...
2588!> \param h_stress ...
2589!> \param rho_core ...
2590!> \param greenfn ...
2591!> \param aux_density Hartree energy and stress tensor between 2 different densities
2592!> \par History
2593!> JGH (13-Mar-2001) : completely revised
2594!> \author apsi
2595! **************************************************************************************************
2596 SUBROUTINE pw_poisson_solve_v_dv_r3d_rs_r3d_rs_c1d_gs (poisson_env, density, ehartree, vhartree, &
2597 dvhartree, h_stress, rho_core, greenfn, aux_density)
2598
2599 TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
2600 TYPE(pw_r3d_rs_type), INTENT(IN) :: density
2601 REAL(kind=dp), INTENT(out), OPTIONAL :: ehartree
2602 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: vhartree
2603 TYPE(pw_c1d_gs_type), DIMENSION(3) :: dvhartree
2604 REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT), &
2605 OPTIONAL :: h_stress
2606 TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_core, greenfn
2607 TYPE(pw_r3d_rs_type), INTENT(IN), OPTIONAL :: aux_density
2608
2609 CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_solve'
2610
2611 INTEGER :: handle
2612 LOGICAL :: has_dielectric
2613 TYPE(pw_grid_type), POINTER :: pw_grid
2614 TYPE(pw_pool_type), POINTER :: pw_pool
2615 TYPE(pw_r3d_rs_type) :: rhor, vhartree_rs
2616 TYPE(pw_c1d_gs_type) :: influence_fn, rhog, rhog_aux
2617
2618 CALL timeset(routinen, handle)
2619
2620 CALL pw_poisson_rebuild(poisson_env, density)
2621
2622 has_dielectric = poisson_env%parameters%has_dielectric
2623
2624 ! point pw
2625 pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
2626 pw_grid => pw_pool%pw_grid
2627 IF (.NOT. pw_grid_compare(pw_pool%pw_grid, vhartree%pw_grid)) THEN
2628 cpabort("vhartree has a different grid than the poisson solver")
2629 END IF
2630 ! density in G space
2631 CALL pw_pool%create_pw(rhog)
2632 IF (PRESENT(aux_density)) THEN
2633 CALL pw_pool%create_pw(rhog_aux)
2634 END IF
2635
2636 SELECT CASE (poisson_env%used_grid)
2637 CASE (use_gs_grid)
2638
2639 SELECT CASE (poisson_env%green_fft%method)
2640 CASE (periodic3d, analytic2d, analytic1d, analytic0d, mt2d, mt1d, mt0d, multipole0d)
2641
2642 CALL pw_transfer(density, rhog)
2643 IF (PRESENT(aux_density)) THEN
2644 CALL pw_transfer(aux_density, rhog_aux)
2645 END IF
2646 IF (PRESENT(greenfn)) THEN
2647 influence_fn = greenfn
2648 ELSE
2649 influence_fn = poisson_env%green_fft%influence_fn
2650 END IF
2651 CALL pw_multiply_with(rhog, influence_fn)
2652 IF (PRESENT(aux_density)) THEN
2653 CALL pw_multiply_with(rhog_aux, influence_fn)
2654 END IF
2655 CALL pw_transfer(rhog, vhartree)
2656 IF (PRESENT(ehartree)) THEN
2657 IF (PRESENT(aux_density)) THEN
2658 ehartree = 0.5_dp*pw_integral_ab(aux_density, vhartree)
2659 ELSE
2660 ehartree = 0.5_dp*pw_integral_ab(density, vhartree)
2661 END IF
2662 END IF
2663
2664 CASE (ps_implicit)
2665 IF (PRESENT(h_stress)) THEN
2666 CALL cp_abort(__location__, &
2667 "No stress tensor is implemented for the implicit Poisson solver.")
2668 END IF
2669
2670 IF (has_dielectric .AND. PRESENT(rho_core)) THEN
2671 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
2672 CASE (periodic_bc, mixed_periodic_bc)
2673 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
2674 poisson_env%diel_rs_grid, &
2675 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
2676 density, rho_core=rho_core)
2677 CASE (neumann_bc, mixed_bc)
2678 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
2679 poisson_env%diel_rs_grid, &
2680 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
2681 poisson_env%dct_pw_grid, &
2682 poisson_env%parameters%ps_implicit_params%neumann_directions, &
2683 poisson_env%implicit_env%dct_env%recv_msgs_bnds, &
2684 poisson_env%implicit_env%dct_env%dests_expand, &
2685 poisson_env%implicit_env%dct_env%srcs_expand, &
2686 poisson_env%implicit_env%dct_env%flipg_stat, &
2687 poisson_env%implicit_env%dct_env%bounds_shftd, &
2688 density, rho_core=rho_core)
2689 END SELECT
2690 END IF
2691
2692 CALL pw_pool%create_pw(rhor)
2693 CALL pw_pool%create_pw(vhartree_rs)
2694 CALL pw_transfer(density, rhor)
2695
2696 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
2697 CASE (periodic_bc)
2698 CALL implicit_poisson_solver_periodic(poisson_env, rhor, vhartree_rs, &
2699 ehartree=ehartree)
2700 CASE (neumann_bc)
2701 CALL implicit_poisson_solver_neumann(poisson_env, rhor, vhartree_rs, &
2702 ehartree=ehartree)
2703 CASE (mixed_periodic_bc)
2704 CALL implicit_poisson_solver_mixed_periodic(poisson_env, rhor, vhartree_rs, &
2705 electric_enthalpy=ehartree)
2706 CASE (mixed_bc)
2707 CALL implicit_poisson_solver_mixed(poisson_env, rhor, vhartree_rs, &
2708 electric_enthalpy=ehartree)
2709 END SELECT
2710
2711 CALL pw_transfer(vhartree_rs, vhartree)
2712 CALL pw_transfer(rhor, rhog)
2713
2714 IF (PRESENT(aux_density)) THEN
2715 CALL pw_transfer(aux_density, rhor)
2716 ehartree = 0.5_dp*pw_integral_ab(rhor, vhartree_rs)
2717 END IF
2718
2719 CALL pw_pool%give_back_pw(rhor)
2720 CALL pw_pool%give_back_pw(vhartree_rs)
2721
2722 CASE DEFAULT
2723 CALL cp_abort(__location__, &
2724 "unknown poisson method "// &
2725 cp_to_string(poisson_env%green_fft%method))
2726 END SELECT
2727
2728 CASE (use_rs_grid)
2729
2730 CALL pw_pool%create_pw(rhor)
2731 CALL pw_transfer(density, rhor)
2732 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
2733 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
2734 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
2735 CALL pw_transfer(rhor, vhartree)
2736 CALL pw_transfer(rhor, rhog)
2737 IF (PRESENT(ehartree)) THEN
2738 IF (PRESENT(aux_density)) THEN
2739 ehartree = 0.5_dp*pw_integral_ab(aux_density, vhartree)
2740 ELSE
2741 ehartree = 0.5_dp*pw_integral_ab(density, vhartree)
2742 END IF
2743 END IF
2744 CALL pw_transfer(rhor, rhog)
2745 IF (PRESENT(aux_density)) THEN
2746 CALL pw_transfer(aux_density, rhor)
2747 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
2748 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
2749 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
2750 CALL pw_transfer(rhor, rhog_aux)
2751 END IF
2752 CALL pw_pool%give_back_pw(rhor)
2753
2754 END SELECT
2755
2756 IF (PRESENT(aux_density)) THEN
2757 CALL calc_stress_and_gradient_c1d_gs (poisson_env, rhog, ehartree, rhog_aux, h_stress, dvhartree=dvhartree)
2758 ELSE
2759 CALL calc_stress_and_gradient_c1d_gs (poisson_env, rhog, ehartree, h_stress=h_stress, dvhartree=dvhartree)
2760 END IF
2761
2762 CALL pw_pool%give_back_pw(rhog)
2763 IF (PRESENT(aux_density)) THEN
2764 CALL pw_pool%give_back_pw(rhog_aux)
2765 END IF
2766
2767 CALL timestop(handle)
2768
2769 END SUBROUTINE pw_poisson_solve_v_dv_r3d_rs_r3d_rs_c1d_gs
2770! **************************************************************************************************
2771!> \brief Solve Poisson equation in a plane wave basis set
2772!> Obtains electrostatic potential and its derivatives with respect to r
2773!> from the density
2774!> \param poisson_env ...
2775!> \param density ...
2776!> \param ehartree ...
2777!> \param vhartree ...
2778!> \param dvhartree ...
2779!> \param h_stress ...
2780!> \param rho_core ...
2781!> \param greenfn ...
2782!> \param aux_density Hartree energy and stress tensor between 2 different densities
2783!> \par History
2784!> JGH (13-Mar-2001) : completely revised
2785!> \author apsi
2786! **************************************************************************************************
2787 SUBROUTINE pw_poisson_solve_v_dv_r3d_rs_c1d_gs_c1d_gs (poisson_env, density, ehartree, vhartree, &
2788 dvhartree, h_stress, rho_core, greenfn, aux_density)
2789
2790 TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
2791 TYPE(pw_r3d_rs_type), INTENT(IN) :: density
2792 REAL(kind=dp), INTENT(out), OPTIONAL :: ehartree
2793 TYPE(pw_c1d_gs_type), INTENT(INOUT) :: vhartree
2794 TYPE(pw_c1d_gs_type), DIMENSION(3) :: dvhartree
2795 REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT), &
2796 OPTIONAL :: h_stress
2797 TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_core, greenfn
2798 TYPE(pw_r3d_rs_type), INTENT(IN), OPTIONAL :: aux_density
2799
2800 CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_solve'
2801
2802 INTEGER :: handle
2803 LOGICAL :: has_dielectric
2804 TYPE(pw_grid_type), POINTER :: pw_grid
2805 TYPE(pw_pool_type), POINTER :: pw_pool
2806 TYPE(pw_r3d_rs_type) :: rhor, vhartree_rs
2807 TYPE(pw_c1d_gs_type) :: influence_fn, rhog, rhog_aux
2808
2809 CALL timeset(routinen, handle)
2810
2811 CALL pw_poisson_rebuild(poisson_env, density)
2812
2813 has_dielectric = poisson_env%parameters%has_dielectric
2814
2815 ! point pw
2816 pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
2817 pw_grid => pw_pool%pw_grid
2818 IF (.NOT. pw_grid_compare(pw_pool%pw_grid, vhartree%pw_grid)) THEN
2819 cpabort("vhartree has a different grid than the poisson solver")
2820 END IF
2821 ! density in G space
2822 CALL pw_pool%create_pw(rhog)
2823 IF (PRESENT(aux_density)) THEN
2824 CALL pw_pool%create_pw(rhog_aux)
2825 END IF
2826
2827 SELECT CASE (poisson_env%used_grid)
2828 CASE (use_gs_grid)
2829
2830 SELECT CASE (poisson_env%green_fft%method)
2831 CASE (periodic3d, analytic2d, analytic1d, analytic0d, mt2d, mt1d, mt0d, multipole0d)
2832
2833 CALL pw_transfer(density, rhog)
2834 IF (PRESENT(aux_density)) THEN
2835 CALL pw_transfer(aux_density, rhog_aux)
2836 END IF
2837 IF (PRESENT(greenfn)) THEN
2838 influence_fn = greenfn
2839 ELSE
2840 influence_fn = poisson_env%green_fft%influence_fn
2841 END IF
2842 CALL pw_multiply_with(rhog, influence_fn)
2843 IF (PRESENT(aux_density)) THEN
2844 CALL pw_multiply_with(rhog_aux, influence_fn)
2845 END IF
2846 CALL pw_transfer(rhog, vhartree)
2847 IF (PRESENT(ehartree)) THEN
2848 IF (PRESENT(aux_density)) THEN
2849 CALL pw_transfer(aux_density, rhog)
2850 ehartree = 0.5_dp*pw_integral_ab(rhog, vhartree)
2851 ELSE
2852 CALL pw_transfer(density, rhog)
2853 ehartree = 0.5_dp*pw_integral_ab(rhog, vhartree)
2854 END IF
2855 END IF
2856
2857 CASE (ps_implicit)
2858 IF (PRESENT(h_stress)) THEN
2859 CALL cp_abort(__location__, &
2860 "No stress tensor is implemented for the implicit Poisson solver.")
2861 END IF
2862
2863 IF (has_dielectric .AND. PRESENT(rho_core)) THEN
2864 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
2865 CASE (periodic_bc, mixed_periodic_bc)
2866 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
2867 poisson_env%diel_rs_grid, &
2868 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
2869 density, rho_core=rho_core)
2870 CASE (neumann_bc, mixed_bc)
2871 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
2872 poisson_env%diel_rs_grid, &
2873 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
2874 poisson_env%dct_pw_grid, &
2875 poisson_env%parameters%ps_implicit_params%neumann_directions, &
2876 poisson_env%implicit_env%dct_env%recv_msgs_bnds, &
2877 poisson_env%implicit_env%dct_env%dests_expand, &
2878 poisson_env%implicit_env%dct_env%srcs_expand, &
2879 poisson_env%implicit_env%dct_env%flipg_stat, &
2880 poisson_env%implicit_env%dct_env%bounds_shftd, &
2881 density, rho_core=rho_core)
2882 END SELECT
2883 END IF
2884
2885 CALL pw_pool%create_pw(rhor)
2886 CALL pw_pool%create_pw(vhartree_rs)
2887 CALL pw_transfer(density, rhor)
2888
2889 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
2890 CASE (periodic_bc)
2891 CALL implicit_poisson_solver_periodic(poisson_env, rhor, vhartree_rs, &
2892 ehartree=ehartree)
2893 CASE (neumann_bc)
2894 CALL implicit_poisson_solver_neumann(poisson_env, rhor, vhartree_rs, &
2895 ehartree=ehartree)
2896 CASE (mixed_periodic_bc)
2897 CALL implicit_poisson_solver_mixed_periodic(poisson_env, rhor, vhartree_rs, &
2898 electric_enthalpy=ehartree)
2899 CASE (mixed_bc)
2900 CALL implicit_poisson_solver_mixed(poisson_env, rhor, vhartree_rs, &
2901 electric_enthalpy=ehartree)
2902 END SELECT
2903
2904 CALL pw_transfer(vhartree_rs, vhartree)
2905 CALL pw_transfer(rhor, rhog)
2906
2907 IF (PRESENT(aux_density)) THEN
2908 CALL pw_transfer(aux_density, rhor)
2909 ehartree = 0.5_dp*pw_integral_ab(rhor, vhartree_rs)
2910 END IF
2911
2912 CALL pw_pool%give_back_pw(rhor)
2913 CALL pw_pool%give_back_pw(vhartree_rs)
2914
2915 CASE DEFAULT
2916 CALL cp_abort(__location__, &
2917 "unknown poisson method "// &
2918 cp_to_string(poisson_env%green_fft%method))
2919 END SELECT
2920
2921 CASE (use_rs_grid)
2922
2923 CALL pw_pool%create_pw(rhor)
2924 CALL pw_transfer(density, rhor)
2925 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
2926 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
2927 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
2928 CALL pw_transfer(rhor, vhartree)
2929 CALL pw_transfer(rhor, rhog)
2930 IF (PRESENT(ehartree)) THEN
2931 IF (PRESENT(aux_density)) THEN
2932 ehartree = 0.5_dp*pw_integral_ab(aux_density, rhor)
2933 ELSE
2934 ehartree = 0.5_dp*pw_integral_ab(density, rhor)
2935 END IF
2936 END IF
2937 CALL pw_transfer(rhor, rhog)
2938 IF (PRESENT(aux_density)) THEN
2939 CALL pw_transfer(aux_density, rhor)
2940 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
2941 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
2942 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
2943 CALL pw_transfer(rhor, rhog_aux)
2944 END IF
2945 CALL pw_pool%give_back_pw(rhor)
2946
2947 END SELECT
2948
2949 IF (PRESENT(aux_density)) THEN
2950 CALL calc_stress_and_gradient_c1d_gs (poisson_env, rhog, ehartree, rhog_aux, h_stress, dvhartree=dvhartree)
2951 ELSE
2952 CALL calc_stress_and_gradient_c1d_gs (poisson_env, rhog, ehartree, h_stress=h_stress, dvhartree=dvhartree)
2953 END IF
2954
2955 CALL pw_pool%give_back_pw(rhog)
2956 IF (PRESENT(aux_density)) THEN
2957 CALL pw_pool%give_back_pw(rhog_aux)
2958 END IF
2959
2960 CALL timestop(handle)
2961
2962 END SUBROUTINE pw_poisson_solve_v_dv_r3d_rs_c1d_gs_c1d_gs
2963! **************************************************************************************************
2964!> \brief Solve Poisson equation in a plane wave basis set
2965!> Obtains electrostatic potential and its derivatives with respect to r
2966!> from the density
2967!> \param poisson_env ...
2968!> \param density ...
2969!> \param ehartree ...
2970!> \param vhartree ...
2971!> \param dvhartree ...
2972!> \param h_stress ...
2973!> \param rho_core ...
2974!> \param greenfn ...
2975!> \param aux_density Hartree energy and stress tensor between 2 different densities
2976!> \par History
2977!> JGH (13-Mar-2001) : completely revised
2978!> \author apsi
2979! **************************************************************************************************
2980 SUBROUTINE pw_poisson_solve_v_dv_c1d_gs_r3d_rs_r3d_rs (poisson_env, density, ehartree, vhartree, &
2981 dvhartree, h_stress, rho_core, greenfn, aux_density)
2982
2983 TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
2984 TYPE(pw_c1d_gs_type), INTENT(IN) :: density
2985 REAL(kind=dp), INTENT(out), OPTIONAL :: ehartree
2986 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: vhartree
2987 TYPE(pw_r3d_rs_type), DIMENSION(3) :: dvhartree
2988 REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT), &
2989 OPTIONAL :: h_stress
2990 TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_core, greenfn
2991 TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: aux_density
2992
2993 CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_solve'
2994
2995 INTEGER :: handle
2996 LOGICAL :: has_dielectric
2997 TYPE(pw_grid_type), POINTER :: pw_grid
2998 TYPE(pw_pool_type), POINTER :: pw_pool
2999 TYPE(pw_r3d_rs_type) :: rhor, vhartree_rs
3000 TYPE(pw_c1d_gs_type) :: influence_fn, rhog, rhog_aux
3001
3002 CALL timeset(routinen, handle)
3003
3004 CALL pw_poisson_rebuild(poisson_env, density)
3005
3006 has_dielectric = poisson_env%parameters%has_dielectric
3007
3008 ! point pw
3009 pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
3010 pw_grid => pw_pool%pw_grid
3011 IF (.NOT. pw_grid_compare(pw_pool%pw_grid, vhartree%pw_grid)) THEN
3012 cpabort("vhartree has a different grid than the poisson solver")
3013 END IF
3014 ! density in G space
3015 CALL pw_pool%create_pw(rhog)
3016 IF (PRESENT(aux_density)) THEN
3017 CALL pw_pool%create_pw(rhog_aux)
3018 END IF
3019
3020 SELECT CASE (poisson_env%used_grid)
3021 CASE (use_gs_grid)
3022
3023 SELECT CASE (poisson_env%green_fft%method)
3024 CASE (periodic3d, analytic2d, analytic1d, analytic0d, mt2d, mt1d, mt0d, multipole0d)
3025
3026 CALL pw_transfer(density, rhog)
3027 IF (PRESENT(aux_density)) THEN
3028 CALL pw_transfer(aux_density, rhog_aux)
3029 END IF
3030 IF (PRESENT(greenfn)) THEN
3031 influence_fn = greenfn
3032 ELSE
3033 influence_fn = poisson_env%green_fft%influence_fn
3034 END IF
3035 CALL pw_multiply_with(rhog, influence_fn)
3036 IF (PRESENT(aux_density)) THEN
3037 CALL pw_multiply_with(rhog_aux, influence_fn)
3038 END IF
3039 CALL pw_transfer(rhog, vhartree)
3040 IF (PRESENT(ehartree)) THEN
3041 IF (PRESENT(aux_density)) THEN
3042 ehartree = 0.5_dp*pw_integral_ab(aux_density, rhog)
3043 ELSE
3044 ehartree = 0.5_dp*pw_integral_ab(density, rhog)
3045 END IF
3046 END IF
3047
3048 CASE (ps_implicit)
3049 IF (PRESENT(h_stress)) THEN
3050 CALL cp_abort(__location__, &
3051 "No stress tensor is implemented for the implicit Poisson solver.")
3052 END IF
3053
3054 IF (has_dielectric .AND. PRESENT(rho_core)) THEN
3055 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
3056 CASE (periodic_bc, mixed_periodic_bc)
3057 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
3058 poisson_env%diel_rs_grid, &
3059 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
3060 density, rho_core=rho_core)
3061 CASE (neumann_bc, mixed_bc)
3062 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
3063 poisson_env%diel_rs_grid, &
3064 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
3065 poisson_env%dct_pw_grid, &
3066 poisson_env%parameters%ps_implicit_params%neumann_directions, &
3067 poisson_env%implicit_env%dct_env%recv_msgs_bnds, &
3068 poisson_env%implicit_env%dct_env%dests_expand, &
3069 poisson_env%implicit_env%dct_env%srcs_expand, &
3070 poisson_env%implicit_env%dct_env%flipg_stat, &
3071 poisson_env%implicit_env%dct_env%bounds_shftd, &
3072 density, rho_core=rho_core)
3073 END SELECT
3074 END IF
3075
3076 CALL pw_pool%create_pw(rhor)
3077 CALL pw_pool%create_pw(vhartree_rs)
3078 CALL pw_transfer(density, rhor)
3079
3080 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
3081 CASE (periodic_bc)
3082 CALL implicit_poisson_solver_periodic(poisson_env, rhor, vhartree_rs, &
3083 ehartree=ehartree)
3084 CASE (neumann_bc)
3085 CALL implicit_poisson_solver_neumann(poisson_env, rhor, vhartree_rs, &
3086 ehartree=ehartree)
3087 CASE (mixed_periodic_bc)
3088 CALL implicit_poisson_solver_mixed_periodic(poisson_env, rhor, vhartree_rs, &
3089 electric_enthalpy=ehartree)
3090 CASE (mixed_bc)
3091 CALL implicit_poisson_solver_mixed(poisson_env, rhor, vhartree_rs, &
3092 electric_enthalpy=ehartree)
3093 END SELECT
3094
3095 CALL pw_transfer(vhartree_rs, vhartree)
3096 CALL pw_transfer(rhor, rhog)
3097
3098 IF (PRESENT(aux_density)) THEN
3099 CALL pw_transfer(aux_density, rhor)
3100 ehartree = 0.5_dp*pw_integral_ab(rhor, vhartree_rs)
3101 END IF
3102
3103 CALL pw_pool%give_back_pw(rhor)
3104 CALL pw_pool%give_back_pw(vhartree_rs)
3105
3106 CASE DEFAULT
3107 CALL cp_abort(__location__, &
3108 "unknown poisson method "// &
3109 cp_to_string(poisson_env%green_fft%method))
3110 END SELECT
3111
3112 CASE (use_rs_grid)
3113
3114 CALL pw_pool%create_pw(rhor)
3115 CALL pw_transfer(density, rhor)
3116 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
3117 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
3118 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
3119 CALL pw_transfer(rhor, vhartree)
3120 CALL pw_transfer(rhor, rhog)
3121 IF (PRESENT(ehartree)) THEN
3122 IF (PRESENT(aux_density)) THEN
3123 ehartree = 0.5_dp*pw_integral_ab(aux_density, rhog)
3124 ELSE
3125 ehartree = 0.5_dp*pw_integral_ab(density, rhog)
3126 END IF
3127 END IF
3128 CALL pw_transfer(rhor, rhog)
3129 IF (PRESENT(aux_density)) THEN
3130 CALL pw_transfer(aux_density, rhor)
3131 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
3132 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
3133 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
3134 CALL pw_transfer(rhor, rhog_aux)
3135 END IF
3136 CALL pw_pool%give_back_pw(rhor)
3137
3138 END SELECT
3139
3140 IF (PRESENT(aux_density)) THEN
3141 CALL calc_stress_and_gradient_r3d_rs (poisson_env, rhog, ehartree, rhog_aux, h_stress, dvhartree=dvhartree)
3142 ELSE
3143 CALL calc_stress_and_gradient_r3d_rs (poisson_env, rhog, ehartree, h_stress=h_stress, dvhartree=dvhartree)
3144 END IF
3145
3146 CALL pw_pool%give_back_pw(rhog)
3147 IF (PRESENT(aux_density)) THEN
3148 CALL pw_pool%give_back_pw(rhog_aux)
3149 END IF
3150
3151 CALL timestop(handle)
3152
3153 END SUBROUTINE pw_poisson_solve_v_dv_c1d_gs_r3d_rs_r3d_rs
3154! **************************************************************************************************
3155!> \brief Solve Poisson equation in a plane wave basis set
3156!> Obtains electrostatic potential and its derivatives with respect to r
3157!> from the density
3158!> \param poisson_env ...
3159!> \param density ...
3160!> \param ehartree ...
3161!> \param vhartree ...
3162!> \param dvhartree ...
3163!> \param h_stress ...
3164!> \param rho_core ...
3165!> \param greenfn ...
3166!> \param aux_density Hartree energy and stress tensor between 2 different densities
3167!> \par History
3168!> JGH (13-Mar-2001) : completely revised
3169!> \author apsi
3170! **************************************************************************************************
3171 SUBROUTINE pw_poisson_solve_v_dv_c1d_gs_c1d_gs_r3d_rs (poisson_env, density, ehartree, vhartree, &
3172 dvhartree, h_stress, rho_core, greenfn, aux_density)
3173
3174 TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
3175 TYPE(pw_c1d_gs_type), INTENT(IN) :: density
3176 REAL(kind=dp), INTENT(out), OPTIONAL :: ehartree
3177 TYPE(pw_c1d_gs_type), INTENT(INOUT) :: vhartree
3178 TYPE(pw_r3d_rs_type), DIMENSION(3) :: dvhartree
3179 REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT), &
3180 OPTIONAL :: h_stress
3181 TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_core, greenfn
3182 TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: aux_density
3183
3184 CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_solve'
3185
3186 INTEGER :: handle
3187 LOGICAL :: has_dielectric
3188 TYPE(pw_grid_type), POINTER :: pw_grid
3189 TYPE(pw_pool_type), POINTER :: pw_pool
3190 TYPE(pw_r3d_rs_type) :: rhor, vhartree_rs
3191 TYPE(pw_c1d_gs_type) :: influence_fn, rhog, rhog_aux
3192
3193 CALL timeset(routinen, handle)
3194
3195 CALL pw_poisson_rebuild(poisson_env, density)
3196
3197 has_dielectric = poisson_env%parameters%has_dielectric
3198
3199 ! point pw
3200 pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
3201 pw_grid => pw_pool%pw_grid
3202 IF (.NOT. pw_grid_compare(pw_pool%pw_grid, vhartree%pw_grid)) THEN
3203 cpabort("vhartree has a different grid than the poisson solver")
3204 END IF
3205 ! density in G space
3206 CALL pw_pool%create_pw(rhog)
3207 IF (PRESENT(aux_density)) THEN
3208 CALL pw_pool%create_pw(rhog_aux)
3209 END IF
3210
3211 SELECT CASE (poisson_env%used_grid)
3212 CASE (use_gs_grid)
3213
3214 SELECT CASE (poisson_env%green_fft%method)
3215 CASE (periodic3d, analytic2d, analytic1d, analytic0d, mt2d, mt1d, mt0d, multipole0d)
3216
3217 CALL pw_transfer(density, rhog)
3218 IF (PRESENT(aux_density)) THEN
3219 CALL pw_transfer(aux_density, rhog_aux)
3220 END IF
3221 IF (PRESENT(greenfn)) THEN
3222 influence_fn = greenfn
3223 ELSE
3224 influence_fn = poisson_env%green_fft%influence_fn
3225 END IF
3226 CALL pw_multiply_with(rhog, influence_fn)
3227 IF (PRESENT(aux_density)) THEN
3228 CALL pw_multiply_with(rhog_aux, influence_fn)
3229 END IF
3230 CALL pw_transfer(rhog, vhartree)
3231 IF (PRESENT(ehartree)) THEN
3232 IF (PRESENT(aux_density)) THEN
3233 ehartree = 0.5_dp*pw_integral_ab(aux_density, vhartree)
3234 ELSE
3235 ehartree = 0.5_dp*pw_integral_ab(density, vhartree)
3236 END IF
3237 END IF
3238
3239 CASE (ps_implicit)
3240 IF (PRESENT(h_stress)) THEN
3241 CALL cp_abort(__location__, &
3242 "No stress tensor is implemented for the implicit Poisson solver.")
3243 END IF
3244
3245 IF (has_dielectric .AND. PRESENT(rho_core)) THEN
3246 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
3247 CASE (periodic_bc, mixed_periodic_bc)
3248 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
3249 poisson_env%diel_rs_grid, &
3250 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
3251 density, rho_core=rho_core)
3252 CASE (neumann_bc, mixed_bc)
3253 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
3254 poisson_env%diel_rs_grid, &
3255 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
3256 poisson_env%dct_pw_grid, &
3257 poisson_env%parameters%ps_implicit_params%neumann_directions, &
3258 poisson_env%implicit_env%dct_env%recv_msgs_bnds, &
3259 poisson_env%implicit_env%dct_env%dests_expand, &
3260 poisson_env%implicit_env%dct_env%srcs_expand, &
3261 poisson_env%implicit_env%dct_env%flipg_stat, &
3262 poisson_env%implicit_env%dct_env%bounds_shftd, &
3263 density, rho_core=rho_core)
3264 END SELECT
3265 END IF
3266
3267 CALL pw_pool%create_pw(rhor)
3268 CALL pw_pool%create_pw(vhartree_rs)
3269 CALL pw_transfer(density, rhor)
3270
3271 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
3272 CASE (periodic_bc)
3273 CALL implicit_poisson_solver_periodic(poisson_env, rhor, vhartree_rs, &
3274 ehartree=ehartree)
3275 CASE (neumann_bc)
3276 CALL implicit_poisson_solver_neumann(poisson_env, rhor, vhartree_rs, &
3277 ehartree=ehartree)
3278 CASE (mixed_periodic_bc)
3279 CALL implicit_poisson_solver_mixed_periodic(poisson_env, rhor, vhartree_rs, &
3280 electric_enthalpy=ehartree)
3281 CASE (mixed_bc)
3282 CALL implicit_poisson_solver_mixed(poisson_env, rhor, vhartree_rs, &
3283 electric_enthalpy=ehartree)
3284 END SELECT
3285
3286 CALL pw_transfer(vhartree_rs, vhartree)
3287 CALL pw_transfer(rhor, rhog)
3288
3289 IF (PRESENT(aux_density)) THEN
3290 CALL pw_transfer(aux_density, rhor)
3291 ehartree = 0.5_dp*pw_integral_ab(rhor, vhartree_rs)
3292 END IF
3293
3294 CALL pw_pool%give_back_pw(rhor)
3295 CALL pw_pool%give_back_pw(vhartree_rs)
3296
3297 CASE DEFAULT
3298 CALL cp_abort(__location__, &
3299 "unknown poisson method "// &
3300 cp_to_string(poisson_env%green_fft%method))
3301 END SELECT
3302
3303 CASE (use_rs_grid)
3304
3305 CALL pw_pool%create_pw(rhor)
3306 CALL pw_transfer(density, rhor)
3307 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
3308 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
3309 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
3310 CALL pw_transfer(rhor, vhartree)
3311 CALL pw_transfer(rhor, rhog)
3312 IF (PRESENT(ehartree)) THEN
3313 IF (PRESENT(aux_density)) THEN
3314 ehartree = 0.5_dp*pw_integral_ab(aux_density, vhartree)
3315 ELSE
3316 ehartree = 0.5_dp*pw_integral_ab(density, vhartree)
3317 END IF
3318 END IF
3319 CALL pw_transfer(rhor, rhog)
3320 IF (PRESENT(aux_density)) THEN
3321 CALL pw_transfer(aux_density, rhor)
3322 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
3323 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
3324 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
3325 CALL pw_transfer(rhor, rhog_aux)
3326 END IF
3327 CALL pw_pool%give_back_pw(rhor)
3328
3329 END SELECT
3330
3331 IF (PRESENT(aux_density)) THEN
3332 CALL calc_stress_and_gradient_r3d_rs (poisson_env, rhog, ehartree, rhog_aux, h_stress, dvhartree=dvhartree)
3333 ELSE
3334 CALL calc_stress_and_gradient_r3d_rs (poisson_env, rhog, ehartree, h_stress=h_stress, dvhartree=dvhartree)
3335 END IF
3336
3337 CALL pw_pool%give_back_pw(rhog)
3338 IF (PRESENT(aux_density)) THEN
3339 CALL pw_pool%give_back_pw(rhog_aux)
3340 END IF
3341
3342 CALL timestop(handle)
3343
3344 END SUBROUTINE pw_poisson_solve_v_dv_c1d_gs_c1d_gs_r3d_rs
3345! **************************************************************************************************
3346!> \brief Solve Poisson equation in a plane wave basis set
3347!> Obtains electrostatic potential and its derivatives with respect to r
3348!> from the density
3349!> \param poisson_env ...
3350!> \param density ...
3351!> \param ehartree ...
3352!> \param vhartree ...
3353!> \param dvhartree ...
3354!> \param h_stress ...
3355!> \param rho_core ...
3356!> \param greenfn ...
3357!> \param aux_density Hartree energy and stress tensor between 2 different densities
3358!> \par History
3359!> JGH (13-Mar-2001) : completely revised
3360!> \author apsi
3361! **************************************************************************************************
3362 SUBROUTINE pw_poisson_solve_v_dv_c1d_gs_r3d_rs_c1d_gs (poisson_env, density, ehartree, vhartree, &
3363 dvhartree, h_stress, rho_core, greenfn, aux_density)
3364
3365 TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
3366 TYPE(pw_c1d_gs_type), INTENT(IN) :: density
3367 REAL(kind=dp), INTENT(out), OPTIONAL :: ehartree
3368 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: vhartree
3369 TYPE(pw_c1d_gs_type), DIMENSION(3) :: dvhartree
3370 REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT), &
3371 OPTIONAL :: h_stress
3372 TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_core, greenfn
3373 TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: aux_density
3374
3375 CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_solve'
3376
3377 INTEGER :: handle
3378 LOGICAL :: has_dielectric
3379 TYPE(pw_grid_type), POINTER :: pw_grid
3380 TYPE(pw_pool_type), POINTER :: pw_pool
3381 TYPE(pw_r3d_rs_type) :: rhor, vhartree_rs
3382 TYPE(pw_c1d_gs_type) :: influence_fn, rhog, rhog_aux
3383
3384 CALL timeset(routinen, handle)
3385
3386 CALL pw_poisson_rebuild(poisson_env, density)
3387
3388 has_dielectric = poisson_env%parameters%has_dielectric
3389
3390 ! point pw
3391 pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
3392 pw_grid => pw_pool%pw_grid
3393 IF (.NOT. pw_grid_compare(pw_pool%pw_grid, vhartree%pw_grid)) THEN
3394 cpabort("vhartree has a different grid than the poisson solver")
3395 END IF
3396 ! density in G space
3397 CALL pw_pool%create_pw(rhog)
3398 IF (PRESENT(aux_density)) THEN
3399 CALL pw_pool%create_pw(rhog_aux)
3400 END IF
3401
3402 SELECT CASE (poisson_env%used_grid)
3403 CASE (use_gs_grid)
3404
3405 SELECT CASE (poisson_env%green_fft%method)
3406 CASE (periodic3d, analytic2d, analytic1d, analytic0d, mt2d, mt1d, mt0d, multipole0d)
3407
3408 CALL pw_transfer(density, rhog)
3409 IF (PRESENT(aux_density)) THEN
3410 CALL pw_transfer(aux_density, rhog_aux)
3411 END IF
3412 IF (PRESENT(greenfn)) THEN
3413 influence_fn = greenfn
3414 ELSE
3415 influence_fn = poisson_env%green_fft%influence_fn
3416 END IF
3417 CALL pw_multiply_with(rhog, influence_fn)
3418 IF (PRESENT(aux_density)) THEN
3419 CALL pw_multiply_with(rhog_aux, influence_fn)
3420 END IF
3421 CALL pw_transfer(rhog, vhartree)
3422 IF (PRESENT(ehartree)) THEN
3423 IF (PRESENT(aux_density)) THEN
3424 ehartree = 0.5_dp*pw_integral_ab(aux_density, rhog)
3425 ELSE
3426 ehartree = 0.5_dp*pw_integral_ab(density, rhog)
3427 END IF
3428 END IF
3429
3430 CASE (ps_implicit)
3431 IF (PRESENT(h_stress)) THEN
3432 CALL cp_abort(__location__, &
3433 "No stress tensor is implemented for the implicit Poisson solver.")
3434 END IF
3435
3436 IF (has_dielectric .AND. PRESENT(rho_core)) THEN
3437 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
3438 CASE (periodic_bc, mixed_periodic_bc)
3439 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
3440 poisson_env%diel_rs_grid, &
3441 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
3442 density, rho_core=rho_core)
3443 CASE (neumann_bc, mixed_bc)
3444 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
3445 poisson_env%diel_rs_grid, &
3446 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
3447 poisson_env%dct_pw_grid, &
3448 poisson_env%parameters%ps_implicit_params%neumann_directions, &
3449 poisson_env%implicit_env%dct_env%recv_msgs_bnds, &
3450 poisson_env%implicit_env%dct_env%dests_expand, &
3451 poisson_env%implicit_env%dct_env%srcs_expand, &
3452 poisson_env%implicit_env%dct_env%flipg_stat, &
3453 poisson_env%implicit_env%dct_env%bounds_shftd, &
3454 density, rho_core=rho_core)
3455 END SELECT
3456 END IF
3457
3458 CALL pw_pool%create_pw(rhor)
3459 CALL pw_pool%create_pw(vhartree_rs)
3460 CALL pw_transfer(density, rhor)
3461
3462 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
3463 CASE (periodic_bc)
3464 CALL implicit_poisson_solver_periodic(poisson_env, rhor, vhartree_rs, &
3465 ehartree=ehartree)
3466 CASE (neumann_bc)
3467 CALL implicit_poisson_solver_neumann(poisson_env, rhor, vhartree_rs, &
3468 ehartree=ehartree)
3469 CASE (mixed_periodic_bc)
3470 CALL implicit_poisson_solver_mixed_periodic(poisson_env, rhor, vhartree_rs, &
3471 electric_enthalpy=ehartree)
3472 CASE (mixed_bc)
3473 CALL implicit_poisson_solver_mixed(poisson_env, rhor, vhartree_rs, &
3474 electric_enthalpy=ehartree)
3475 END SELECT
3476
3477 CALL pw_transfer(vhartree_rs, vhartree)
3478 CALL pw_transfer(rhor, rhog)
3479
3480 IF (PRESENT(aux_density)) THEN
3481 CALL pw_transfer(aux_density, rhor)
3482 ehartree = 0.5_dp*pw_integral_ab(rhor, vhartree_rs)
3483 END IF
3484
3485 CALL pw_pool%give_back_pw(rhor)
3486 CALL pw_pool%give_back_pw(vhartree_rs)
3487
3488 CASE DEFAULT
3489 CALL cp_abort(__location__, &
3490 "unknown poisson method "// &
3491 cp_to_string(poisson_env%green_fft%method))
3492 END SELECT
3493
3494 CASE (use_rs_grid)
3495
3496 CALL pw_pool%create_pw(rhor)
3497 CALL pw_transfer(density, rhor)
3498 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
3499 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
3500 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
3501 CALL pw_transfer(rhor, vhartree)
3502 CALL pw_transfer(rhor, rhog)
3503 IF (PRESENT(ehartree)) THEN
3504 IF (PRESENT(aux_density)) THEN
3505 ehartree = 0.5_dp*pw_integral_ab(aux_density, rhog)
3506 ELSE
3507 ehartree = 0.5_dp*pw_integral_ab(density, rhog)
3508 END IF
3509 END IF
3510 CALL pw_transfer(rhor, rhog)
3511 IF (PRESENT(aux_density)) THEN
3512 CALL pw_transfer(aux_density, rhor)
3513 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
3514 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
3515 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
3516 CALL pw_transfer(rhor, rhog_aux)
3517 END IF
3518 CALL pw_pool%give_back_pw(rhor)
3519
3520 END SELECT
3521
3522 IF (PRESENT(aux_density)) THEN
3523 CALL calc_stress_and_gradient_c1d_gs (poisson_env, rhog, ehartree, rhog_aux, h_stress, dvhartree=dvhartree)
3524 ELSE
3525 CALL calc_stress_and_gradient_c1d_gs (poisson_env, rhog, ehartree, h_stress=h_stress, dvhartree=dvhartree)
3526 END IF
3527
3528 CALL pw_pool%give_back_pw(rhog)
3529 IF (PRESENT(aux_density)) THEN
3530 CALL pw_pool%give_back_pw(rhog_aux)
3531 END IF
3532
3533 CALL timestop(handle)
3534
3535 END SUBROUTINE pw_poisson_solve_v_dv_c1d_gs_r3d_rs_c1d_gs
3536! **************************************************************************************************
3537!> \brief Solve Poisson equation in a plane wave basis set
3538!> Obtains electrostatic potential and its derivatives with respect to r
3539!> from the density
3540!> \param poisson_env ...
3541!> \param density ...
3542!> \param ehartree ...
3543!> \param vhartree ...
3544!> \param dvhartree ...
3545!> \param h_stress ...
3546!> \param rho_core ...
3547!> \param greenfn ...
3548!> \param aux_density Hartree energy and stress tensor between 2 different densities
3549!> \par History
3550!> JGH (13-Mar-2001) : completely revised
3551!> \author apsi
3552! **************************************************************************************************
3553 SUBROUTINE pw_poisson_solve_v_dv_c1d_gs_c1d_gs_c1d_gs (poisson_env, density, ehartree, vhartree, &
3554 dvhartree, h_stress, rho_core, greenfn, aux_density)
3555
3556 TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
3557 TYPE(pw_c1d_gs_type), INTENT(IN) :: density
3558 REAL(kind=dp), INTENT(out), OPTIONAL :: ehartree
3559 TYPE(pw_c1d_gs_type), INTENT(INOUT) :: vhartree
3560 TYPE(pw_c1d_gs_type), DIMENSION(3) :: dvhartree
3561 REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT), &
3562 OPTIONAL :: h_stress
3563 TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_core, greenfn
3564 TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: aux_density
3565
3566 CHARACTER(len=*), PARAMETER :: routineN = 'pw_poisson_solve'
3567
3568 INTEGER :: handle
3569 LOGICAL :: has_dielectric
3570 TYPE(pw_grid_type), POINTER :: pw_grid
3571 TYPE(pw_pool_type), POINTER :: pw_pool
3572 TYPE(pw_r3d_rs_type) :: rhor, vhartree_rs
3573 TYPE(pw_c1d_gs_type) :: influence_fn, rhog, rhog_aux
3574
3575 CALL timeset(routinen, handle)
3576
3577 CALL pw_poisson_rebuild(poisson_env, density)
3578
3579 has_dielectric = poisson_env%parameters%has_dielectric
3580
3581 ! point pw
3582 pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
3583 pw_grid => pw_pool%pw_grid
3584 IF (.NOT. pw_grid_compare(pw_pool%pw_grid, vhartree%pw_grid)) THEN
3585 cpabort("vhartree has a different grid than the poisson solver")
3586 END IF
3587 ! density in G space
3588 CALL pw_pool%create_pw(rhog)
3589 IF (PRESENT(aux_density)) THEN
3590 CALL pw_pool%create_pw(rhog_aux)
3591 END IF
3592
3593 SELECT CASE (poisson_env%used_grid)
3594 CASE (use_gs_grid)
3595
3596 SELECT CASE (poisson_env%green_fft%method)
3597 CASE (periodic3d, analytic2d, analytic1d, analytic0d, mt2d, mt1d, mt0d, multipole0d)
3598
3599 CALL pw_transfer(density, rhog)
3600 IF (PRESENT(aux_density)) THEN
3601 CALL pw_transfer(aux_density, rhog_aux)
3602 END IF
3603 IF (PRESENT(greenfn)) THEN
3604 influence_fn = greenfn
3605 ELSE
3606 influence_fn = poisson_env%green_fft%influence_fn
3607 END IF
3608 CALL pw_multiply_with(rhog, influence_fn)
3609 IF (PRESENT(aux_density)) THEN
3610 CALL pw_multiply_with(rhog_aux, influence_fn)
3611 END IF
3612 CALL pw_transfer(rhog, vhartree)
3613 IF (PRESENT(ehartree)) THEN
3614 IF (PRESENT(aux_density)) THEN
3615 ehartree = 0.5_dp*pw_integral_ab(aux_density, vhartree)
3616 ELSE
3617 ehartree = 0.5_dp*pw_integral_ab(density, vhartree)
3618 END IF
3619 END IF
3620
3621 CASE (ps_implicit)
3622 IF (PRESENT(h_stress)) THEN
3623 CALL cp_abort(__location__, &
3624 "No stress tensor is implemented for the implicit Poisson solver.")
3625 END IF
3626
3627 IF (has_dielectric .AND. PRESENT(rho_core)) THEN
3628 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
3629 CASE (periodic_bc, mixed_periodic_bc)
3630 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
3631 poisson_env%diel_rs_grid, &
3632 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
3633 density, rho_core=rho_core)
3634 CASE (neumann_bc, mixed_bc)
3635 CALL dielectric_compute(poisson_env%implicit_env%dielectric, &
3636 poisson_env%diel_rs_grid, &
3637 poisson_env%pw_pools(poisson_env%pw_level)%pool, &
3638 poisson_env%dct_pw_grid, &
3639 poisson_env%parameters%ps_implicit_params%neumann_directions, &
3640 poisson_env%implicit_env%dct_env%recv_msgs_bnds, &
3641 poisson_env%implicit_env%dct_env%dests_expand, &
3642 poisson_env%implicit_env%dct_env%srcs_expand, &
3643 poisson_env%implicit_env%dct_env%flipg_stat, &
3644 poisson_env%implicit_env%dct_env%bounds_shftd, &
3645 density, rho_core=rho_core)
3646 END SELECT
3647 END IF
3648
3649 CALL pw_pool%create_pw(rhor)
3650 CALL pw_pool%create_pw(vhartree_rs)
3651 CALL pw_transfer(density, rhor)
3652
3653 SELECT CASE (poisson_env%parameters%ps_implicit_params%boundary_condition)
3654 CASE (periodic_bc)
3655 CALL implicit_poisson_solver_periodic(poisson_env, rhor, vhartree_rs, &
3656 ehartree=ehartree)
3657 CASE (neumann_bc)
3658 CALL implicit_poisson_solver_neumann(poisson_env, rhor, vhartree_rs, &
3659 ehartree=ehartree)
3660 CASE (mixed_periodic_bc)
3661 CALL implicit_poisson_solver_mixed_periodic(poisson_env, rhor, vhartree_rs, &
3662 electric_enthalpy=ehartree)
3663 CASE (mixed_bc)
3664 CALL implicit_poisson_solver_mixed(poisson_env, rhor, vhartree_rs, &
3665 electric_enthalpy=ehartree)
3666 END SELECT
3667
3668 CALL pw_transfer(vhartree_rs, vhartree)
3669 CALL pw_transfer(rhor, rhog)
3670
3671 IF (PRESENT(aux_density)) THEN
3672 CALL pw_transfer(aux_density, rhor)
3673 ehartree = 0.5_dp*pw_integral_ab(rhor, vhartree_rs)
3674 END IF
3675
3676 CALL pw_pool%give_back_pw(rhor)
3677 CALL pw_pool%give_back_pw(vhartree_rs)
3678
3679 CASE DEFAULT
3680 CALL cp_abort(__location__, &
3681 "unknown poisson method "// &
3682 cp_to_string(poisson_env%green_fft%method))
3683 END SELECT
3684
3685 CASE (use_rs_grid)
3686
3687 CALL pw_pool%create_pw(rhor)
3688 CALL pw_transfer(density, rhor)
3689 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
3690 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
3691 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
3692 CALL pw_transfer(rhor, vhartree)
3693 CALL pw_transfer(rhor, rhog)
3694 IF (PRESENT(ehartree)) THEN
3695 IF (PRESENT(aux_density)) THEN
3696 ehartree = 0.5_dp*pw_integral_ab(aux_density, vhartree)
3697 ELSE
3698 ehartree = 0.5_dp*pw_integral_ab(density, vhartree)
3699 END IF
3700 END IF
3701 CALL pw_transfer(rhor, rhog)
3702 IF (PRESENT(aux_density)) THEN
3703 CALL pw_transfer(aux_density, rhor)
3704 CALL cp2k_distribution_to_z_slices(rhor, poisson_env%wavelet, rhor%pw_grid)
3705 CALL ps_wavelet_solve(poisson_env%wavelet, rhor%pw_grid)
3706 CALL z_slices_to_cp2k_distribution(rhor, poisson_env%wavelet, rhor%pw_grid)
3707 CALL pw_transfer(rhor, rhog_aux)
3708 END IF
3709 CALL pw_pool%give_back_pw(rhor)
3710
3711 END SELECT
3712
3713 IF (PRESENT(aux_density)) THEN
3714 CALL calc_stress_and_gradient_c1d_gs (poisson_env, rhog, ehartree, rhog_aux, h_stress, dvhartree=dvhartree)
3715 ELSE
3716 CALL calc_stress_and_gradient_c1d_gs (poisson_env, rhog, ehartree, h_stress=h_stress, dvhartree=dvhartree)
3717 END IF
3718
3719 CALL pw_pool%give_back_pw(rhog)
3720 IF (PRESENT(aux_density)) THEN
3721 CALL pw_pool%give_back_pw(rhog_aux)
3722 END IF
3723
3724 CALL timestop(handle)
3725
3726 END SUBROUTINE pw_poisson_solve_v_dv_c1d_gs_c1d_gs_c1d_gs
3727
3728 SUBROUTINE calc_stress_and_gradient_c1d_gs (poisson_env, rhog, ehartree, rhog_aux, h_stress, dvhartree)
3729 TYPE(pw_poisson_type), INTENT(IN) :: poisson_env
3730 TYPE(pw_c1d_gs_type), INTENT(IN) :: rhog
3731 REAL(KIND=dp) :: ehartree
3732 TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rhog_aux
3733 REAL(kind=dp), DIMENSION(3, 3), INTENT(OUT), OPTIONAL :: h_stress
3734 TYPE(pw_c1d_gs_type), DIMENSION(3), INTENT(INOUT), OPTIONAL :: dvhartree
3735
3736 CHARACTER(len=*), PARAMETER :: routinen = 'pw_poisson_set'
3737
3738 REAL(kind=dp) :: ffa
3739 INTEGER :: alpha, beta, n(3), handle, i
3740 TYPE(pw_c1d_gs_type) :: dvg(3), dvg_aux(3)
3741 TYPE(pw_pool_type), POINTER :: pw_pool
3742
3743 CALL timeset(routinen, handle)
3744
3745 pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
3746
3747 DO i = 1, 3
3748 CALL pw_pool%create_pw(dvg(i))
3749 n = 0
3750 n(i) = 1
3751 CALL pw_copy(rhog, dvg(i))
3752 CALL pw_derive(dvg(i), n)
3753 IF (PRESENT(rhog_aux)) THEN
3754 CALL pw_pool%create_pw(dvg_aux(i))
3755 CALL pw_copy(rhog_aux, dvg_aux(i))
3756 CALL pw_derive(dvg_aux(i), n)
3757 END IF
3758 END DO
3759 ! save the derivatives
3760 IF (PRESENT(dvhartree)) THEN
3761 DO i = 1, 3
3762 CALL pw_transfer(dvg(i), dvhartree(i))
3763 END DO
3764 END IF
3765 ! Calculate the contribution to the stress tensor this is only the contribution from
3766 ! the Greens FUNCTION and the volume factor of the plane waves
3767 IF (PRESENT(h_stress)) THEN
3768 ffa = -1.0_dp/fourpi
3769 h_stress = 0.0_dp
3770 DO alpha = 1, 3
3771 h_stress(alpha, alpha) = ehartree
3772 IF (PRESENT(rhog_aux)) THEN
3773 DO beta = alpha, 3
3774 h_stress(alpha, beta) = h_stress(alpha, beta) &
3775 + ffa*pw_integral_ab(dvg_aux(alpha), dvg(beta))
3776 h_stress(beta, alpha) = h_stress(alpha, beta)
3777 END DO
3778 ELSE
3779 DO beta = alpha, 3
3780 h_stress(alpha, beta) = h_stress(alpha, beta) &
3781 + ffa*pw_integral_ab(dvg(alpha), dvg(beta))
3782 h_stress(beta, alpha) = h_stress(alpha, beta)
3783 END DO
3784 END IF
3785 END DO
3786
3787 ! Handle the periodicity cases for the Stress Tensor
3788 SELECT CASE (poisson_env%used_grid)
3789 CASE (use_gs_grid)
3790
3791 ! FFT based Poisson-Solver
3792 SELECT CASE (poisson_env%green_fft%method)
3793 CASE (periodic3d, ps_implicit)
3794 ! Do Nothing
3795 CASE (analytic2d, mt2d)
3796 ! Zero the 1 non-periodic component. Only the periodic-subspace
3797 ! stress components are defined for 2D systems.
3798 alpha = poisson_env%green_fft%special_dimension
3799 h_stress(:, alpha) = 0.0_dp
3800 h_stress(alpha, :) = 0.0_dp
3801 CASE (analytic1d, mt1d)
3802 ! Zero the 2 non-periodic components
3803 DO alpha = 1, 3
3804 DO beta = alpha, 3
3805 IF ((alpha /= poisson_env%green_fft%special_dimension) .OR. &
3806 (beta /= poisson_env%green_fft%special_dimension)) THEN
3807 h_stress(alpha, beta) = 0.0_dp
3808 h_stress(beta, alpha) = 0.0_dp
3809 END IF
3810 END DO
3811 END DO
3812 cpabort("Stress Tensor not tested for 1D systems.")
3813 CASE (analytic0d, mt0d, multipole0d)
3814 ! Zero the full stress tensor
3815 h_stress = 0.0_dp
3816 CASE DEFAULT
3817 CALL cp_abort(__location__, &
3818 "unknown poisson method"// &
3819 cp_to_string(poisson_env%green_fft%method))
3820 END SELECT
3821
3822 CASE (use_rs_grid)
3823
3824 ! Wavelet based Poisson-Solver
3825 SELECT CASE (poisson_env%wavelet%method)
3826 CASE (wavelet3d)
3827 ! Do Nothing
3828 CASE (wavelet2d)
3829 ! Zero the 1 non-periodic component
3830 alpha = poisson_env%wavelet%special_dimension
3831 h_stress(:, alpha) = 0.0_dp
3832 h_stress(alpha, :) = 0.0_dp
3833 cpabort("Stress Tensor not tested for WAVELET 2D.")
3834 CASE (wavelet1d)
3835 ! Zero the 2 non-periodic components
3836 cpabort("WAVELET 1D not implemented!")
3837 CASE (wavelet0d)
3838 ! Zero the full stress tensor
3839 h_stress = 0.0_dp
3840 END SELECT
3841
3842 END SELECT
3843 END IF
3844
3845 DO i = 1, 3
3846 CALL pw_pool%give_back_pw(dvg(i))
3847 IF (PRESENT(rhog_aux)) THEN
3848 CALL pw_pool%give_back_pw(dvg_aux(i))
3849 END IF
3850 END DO
3851
3852 CALL timestop(handle)
3853
3854 END SUBROUTINE calc_stress_and_gradient_c1d_gs
3855 SUBROUTINE calc_stress_and_gradient_r3d_rs (poisson_env, rhog, ehartree, rhog_aux, h_stress, dvhartree)
3856 TYPE(pw_poisson_type), INTENT(IN) :: poisson_env
3857 TYPE(pw_c1d_gs_type), INTENT(IN) :: rhog
3858 REAL(KIND=dp) :: ehartree
3859 TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rhog_aux
3860 REAL(kind=dp), DIMENSION(3, 3), INTENT(OUT), OPTIONAL :: h_stress
3861 TYPE(pw_r3d_rs_type), DIMENSION(3), INTENT(INOUT), OPTIONAL :: dvhartree
3862
3863 CHARACTER(len=*), PARAMETER :: routinen = 'pw_poisson_set'
3864
3865 REAL(kind=dp) :: ffa
3866 INTEGER :: alpha, beta, n(3), handle, i
3867 TYPE(pw_c1d_gs_type) :: dvg(3), dvg_aux(3)
3868 TYPE(pw_pool_type), POINTER :: pw_pool
3869
3870 CALL timeset(routinen, handle)
3871
3872 pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
3873
3874 DO i = 1, 3
3875 CALL pw_pool%create_pw(dvg(i))
3876 n = 0
3877 n(i) = 1
3878 CALL pw_copy(rhog, dvg(i))
3879 CALL pw_derive(dvg(i), n)
3880 IF (PRESENT(rhog_aux)) THEN
3881 CALL pw_pool%create_pw(dvg_aux(i))
3882 CALL pw_copy(rhog_aux, dvg_aux(i))
3883 CALL pw_derive(dvg_aux(i), n)
3884 END IF
3885 END DO
3886 ! save the derivatives
3887 IF (PRESENT(dvhartree)) THEN
3888 DO i = 1, 3
3889 CALL pw_transfer(dvg(i), dvhartree(i))
3890 END DO
3891 END IF
3892 ! Calculate the contribution to the stress tensor this is only the contribution from
3893 ! the Greens FUNCTION and the volume factor of the plane waves
3894 IF (PRESENT(h_stress)) THEN
3895 ffa = -1.0_dp/fourpi
3896 h_stress = 0.0_dp
3897 DO alpha = 1, 3
3898 h_stress(alpha, alpha) = ehartree
3899 IF (PRESENT(rhog_aux)) THEN
3900 DO beta = alpha, 3
3901 h_stress(alpha, beta) = h_stress(alpha, beta) &
3902 + ffa*pw_integral_ab(dvg_aux(alpha), dvg(beta))
3903 h_stress(beta, alpha) = h_stress(alpha, beta)
3904 END DO
3905 ELSE
3906 DO beta = alpha, 3
3907 h_stress(alpha, beta) = h_stress(alpha, beta) &
3908 + ffa*pw_integral_ab(dvg(alpha), dvg(beta))
3909 h_stress(beta, alpha) = h_stress(alpha, beta)
3910 END DO
3911 END IF
3912 END DO
3913
3914 ! Handle the periodicity cases for the Stress Tensor
3915 SELECT CASE (poisson_env%used_grid)
3916 CASE (use_gs_grid)
3917
3918 ! FFT based Poisson-Solver
3919 SELECT CASE (poisson_env%green_fft%method)
3920 CASE (periodic3d, ps_implicit)
3921 ! Do Nothing
3922 CASE (analytic2d, mt2d)
3923 ! Zero the 1 non-periodic component. Only the periodic-subspace
3924 ! stress components are defined for 2D systems.
3925 alpha = poisson_env%green_fft%special_dimension
3926 h_stress(:, alpha) = 0.0_dp
3927 h_stress(alpha, :) = 0.0_dp
3928 CASE (analytic1d, mt1d)
3929 ! Zero the 2 non-periodic components
3930 DO alpha = 1, 3
3931 DO beta = alpha, 3
3932 IF ((alpha /= poisson_env%green_fft%special_dimension) .OR. &
3933 (beta /= poisson_env%green_fft%special_dimension)) THEN
3934 h_stress(alpha, beta) = 0.0_dp
3935 h_stress(beta, alpha) = 0.0_dp
3936 END IF
3937 END DO
3938 END DO
3939 cpabort("Stress Tensor not tested for 1D systems.")
3940 CASE (analytic0d, mt0d, multipole0d)
3941 ! Zero the full stress tensor
3942 h_stress = 0.0_dp
3943 CASE DEFAULT
3944 CALL cp_abort(__location__, &
3945 "unknown poisson method"// &
3946 cp_to_string(poisson_env%green_fft%method))
3947 END SELECT
3948
3949 CASE (use_rs_grid)
3950
3951 ! Wavelet based Poisson-Solver
3952 SELECT CASE (poisson_env%wavelet%method)
3953 CASE (wavelet3d)
3954 ! Do Nothing
3955 CASE (wavelet2d)
3956 ! Zero the 1 non-periodic component
3957 alpha = poisson_env%wavelet%special_dimension
3958 h_stress(:, alpha) = 0.0_dp
3959 h_stress(alpha, :) = 0.0_dp
3960 cpabort("Stress Tensor not tested for WAVELET 2D.")
3961 CASE (wavelet1d)
3962 ! Zero the 2 non-periodic components
3963 cpabort("WAVELET 1D not implemented!")
3964 CASE (wavelet0d)
3965 ! Zero the full stress tensor
3966 h_stress = 0.0_dp
3967 END SELECT
3968
3969 END SELECT
3970 END IF
3971
3972 DO i = 1, 3
3973 CALL pw_pool%give_back_pw(dvg(i))
3974 IF (PRESENT(rhog_aux)) THEN
3975 CALL pw_pool%give_back_pw(dvg_aux(i))
3976 END IF
3977 END DO
3978
3979 CALL timestop(handle)
3980
3981 END SUBROUTINE calc_stress_and_gradient_r3d_rs
3982
3983! **************************************************************************************************
3984!> \brief sets cell, grids and parameters used by the poisson solver
3985!> You should call this at least once (and set everything)
3986!> before using the poisson solver.
3987!> Smart, doesn't set the thing twice to the same value
3988!> Keeps track of the need to rebuild the poisson_env
3989!> \param poisson_env ...
3990!> \param cell_hmat ...
3991!> \param parameters ...
3992!> \param pw_pools ...
3993!> \param use_level ...
3994!> \param mt_super_ref_pw_grid ...
3995!> \param dct_pw_grid ...
3996!> \param force_rebuild ...
3997!> \author fawzi
3998!> \note
3999!> Checks everything at the end. This means that after *each* call to
4000!> this method the poisson env must be fully ready, so the first time
4001!> you have to set everything at once. Change this behaviour?
4002! **************************************************************************************************
4003 SUBROUTINE pw_poisson_set(poisson_env, cell_hmat, parameters, pw_pools, use_level, &
4004 mt_super_ref_pw_grid, dct_pw_grid, force_rebuild)
4005
4006 TYPE(pw_poisson_type), INTENT(INOUT) :: poisson_env
4007 REAL(kind=dp), DIMENSION(3, 3), INTENT(IN), &
4008 OPTIONAL :: cell_hmat
4009 TYPE(pw_poisson_parameter_type), INTENT(IN), &
4010 OPTIONAL :: parameters
4011 TYPE(pw_pool_p_type), DIMENSION(:), OPTIONAL, &
4012 POINTER :: pw_pools
4013 INTEGER, INTENT(in), OPTIONAL :: use_level
4014 TYPE(pw_grid_type), OPTIONAL, POINTER :: mt_super_ref_pw_grid, dct_pw_grid
4015 LOGICAL, INTENT(in), OPTIONAL :: force_rebuild
4016
4017 CHARACTER(len=*), PARAMETER :: routinen = 'pw_poisson_set'
4018
4019 INTEGER :: handle, i
4020 LOGICAL :: same
4021 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: tmp_pools
4022
4023 CALL timeset(routinen, handle)
4024
4025 IF (PRESENT(parameters)) poisson_env%parameters = parameters
4026
4027 IF (PRESENT(cell_hmat)) THEN
4028 IF (any(poisson_env%cell_hmat /= cell_hmat)) THEN
4029 CALL pw_poisson_cleanup(poisson_env)
4030 END IF
4031 poisson_env%cell_hmat(:, :) = cell_hmat(:, :)
4032 poisson_env%rebuild = .true.
4033 END IF
4034
4035 IF (PRESENT(pw_pools)) THEN
4036 cpassert(ASSOCIATED(pw_pools))
4037 same = .false.
4038 IF (ASSOCIATED(poisson_env%pw_pools)) THEN
4039 same = SIZE(poisson_env%pw_pools) == SIZE(pw_pools)
4040 IF (same) THEN
4041 DO i = 1, SIZE(pw_pools)
4042 IF (.NOT. ASSOCIATED(poisson_env%pw_pools(i)%pool, &
4043 pw_pools(i)%pool)) same = .false.
4044 END DO
4045 END IF
4046 END IF
4047 IF (.NOT. same) THEN
4048 poisson_env%rebuild = .true.
4049 CALL pw_pools_copy(pw_pools, tmp_pools)
4050 CALL pw_pools_dealloc(poisson_env%pw_pools)
4051 poisson_env%pw_pools => tmp_pools
4052 END IF
4053 END IF
4054
4055 IF (PRESENT(use_level)) poisson_env%pw_level = use_level
4056
4057 IF (PRESENT(dct_pw_grid)) THEN
4058 IF (ASSOCIATED(dct_pw_grid)) THEN
4059 CALL pw_grid_retain(dct_pw_grid)
4060 END IF
4061 CALL pw_grid_release(poisson_env%dct_pw_grid)
4062 poisson_env%dct_pw_grid => dct_pw_grid
4063 END IF
4064
4065 IF (PRESENT(mt_super_ref_pw_grid)) THEN
4066 IF (ASSOCIATED(mt_super_ref_pw_grid)) THEN
4067 CALL pw_grid_retain(mt_super_ref_pw_grid)
4068 END IF
4069 CALL pw_grid_release(poisson_env%mt_super_ref_pw_grid)
4070 poisson_env%mt_super_ref_pw_grid => mt_super_ref_pw_grid
4071 END IF
4072
4073 IF (PRESENT(force_rebuild)) THEN
4074 IF (force_rebuild) poisson_env%rebuild = .true.
4075 END IF
4076
4077 CALL pw_poisson_check(poisson_env)
4078
4079 CALL timestop(handle)
4080
4081 END SUBROUTINE pw_poisson_set
4082
4083! **************************************************************************************************
4084!> \brief computes the convolution f * u using reciprocal-space multiplication
4085!> \par
4086!> Given two real-space functions f and u, this routine evaluates their
4087!> convolution by transforming both to G space, multiplying their Fourier
4088!> components, and transforming the product back to real space.
4089!> \param poisson_env poisson environment
4090!> \param func input function f in real space
4091!> \param convolution convolution result (f * u) in real space
4092!> \param u input function/kernel u in real space
4093!> \author Ziwei Chai
4094!> \date 12.03.2026
4095! **************************************************************************************************
4096 SUBROUTINE pw_func_u_convolution(poisson_env, func, convolution, u)
4097
4098 TYPE(pw_poisson_type), POINTER :: poisson_env
4099 TYPE(pw_r3d_rs_type), INTENT(IN) :: func, u
4100 TYPE(pw_r3d_rs_type), INTENT(INOUT), OPTIONAL :: convolution
4101 CHARACTER(len=*), PARAMETER :: routinen = 'pw_func_u_convolution'
4102 INTEGER :: handle
4103 TYPE(pw_grid_type), POINTER :: pw_grid
4104 TYPE(pw_pool_type), POINTER :: pw_pool
4105 TYPE(pw_c1d_gs_type) :: funcg, ug
4106
4107 CALL timeset(routinen, handle)
4108
4109 cpassert(ASSOCIATED(poisson_env))
4110 CALL pw_poisson_rebuild(poisson_env, func)
4111
4112 ! point pw
4113 pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
4114 pw_grid => pw_pool%pw_grid
4115 ! func and u in G space
4116 CALL pw_pool%create_pw(funcg)
4117 CALL pw_pool%create_pw(ug)
4118
4119 CALL pw_transfer(func, funcg)
4120 CALL pw_transfer(u, ug)
4121
4122 funcg%array(:) = funcg%array(:)*ug%array(:)
4123
4124 CALL pw_transfer(funcg, convolution)
4125
4126 CALL pw_pool%give_back_pw(funcg)
4127 CALL pw_pool%give_back_pw(ug)
4128
4129 CALL timestop(handle)
4130
4131 END SUBROUTINE pw_func_u_convolution
4132
4133END MODULE pw_poisson_methods
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
dielectric_methods::dielectric_compute
Definition dielectric_methods.F:49
pw_methods::pw_copy
Definition pw_methods.F:146
pw_methods::pw_integral_ab
Definition pw_methods.F:222
pw_methods::pw_multiply_with
Definition pw_methods.F:203
pw_methods::pw_transfer
Definition pw_methods.F:304
pw_poisson_methods::pw_poisson_rebuild
Definition pw_poisson_methods.F:73
pw_poisson_methods::pw_poisson_solve
Definition pw_poisson_methods.F:78
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
dielectric_methods
methods for evaluating the dielectric constant
Definition dielectric_methods.F:14
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::fourpi
real(kind=dp), parameter, public fourpi
Definition mathconstants.F:113
mt_util
Definition mt_util.F:9
mt_util::mt0d
integer, parameter, public mt0d
Definition mt_util.F:33
mt_util::mt2d
integer, parameter, public mt2d
Definition mt_util.F:33
mt_util::mt1d
integer, parameter, public mt1d
Definition mt_util.F:33
ps_implicit_methods
The implicit (generalized) Poisson solver.
Definition ps_implicit_methods.F:16
ps_implicit_methods::implicit_poisson_solver_periodic
subroutine, public implicit_poisson_solver_periodic(poisson_env, density, v_new, ehartree)
implicit Poisson solver for periodic boundary conditions
Definition ps_implicit_methods.F:187
ps_implicit_methods::implicit_poisson_solver_mixed
subroutine, public implicit_poisson_solver_mixed(poisson_env, density, v_new, electric_enthalpy)
implicit Poisson solver for mixed boundary conditions (Neumann + Dirichlet)
Definition ps_implicit_methods.F:765
ps_implicit_methods::implicit_poisson_solver_neumann
subroutine, public implicit_poisson_solver_neumann(poisson_env, density, v_new, ehartree)
implicit Poisson solver: zero-average solution of the Poisson equation subject to homogeneous Neumann...
Definition ps_implicit_methods.F:327
ps_implicit_methods::implicit_poisson_solver_mixed_periodic
subroutine, public implicit_poisson_solver_mixed_periodic(poisson_env, density, v_new, electric_enthalpy)
implicit Poisson solver for mixed-periodic boundary conditions (periodic + Dirichlet)
Definition ps_implicit_methods.F:502
ps_implicit_methods::ps_implicit_create
subroutine, public ps_implicit_create(pw_pool, poisson_params, dct_pw_grid, green, ps_implicit_env)
Creates implicit Poisson solver environment.
Definition ps_implicit_methods.F:88
ps_implicit_types
Types containing essential information for running implicit (iterative) Poisson solver.
Definition ps_implicit_types.F:15
ps_implicit_types::neumann_bc
integer, parameter, public neumann_bc
Definition ps_implicit_types.F:32
ps_implicit_types::mixed_bc
integer, parameter, public mixed_bc
Definition ps_implicit_types.F:32
ps_implicit_types::mixed_periodic_bc
integer, parameter, public mixed_periodic_bc
Definition ps_implicit_types.F:32
ps_implicit_types::periodic_bc
integer, parameter, public periodic_bc
Definition ps_implicit_types.F:32
ps_wavelet_methods
Definition and initialisation of the ps_wavelet data type. \history 01.2014 Renamed from ps_wavelet_t...
Definition ps_wavelet_methods.F:13
ps_wavelet_methods::ps_wavelet_create
subroutine, public ps_wavelet_create(poisson_params, wavelet, pw_grid)
creates the ps_wavelet_type which is needed for the link to the Poisson Solver of Luigi Genovese
Definition ps_wavelet_methods.F:58
ps_wavelet_methods::ps_wavelet_solve
subroutine, public ps_wavelet_solve(wavelet, pw_grid)
...
Definition ps_wavelet_methods.F:862
ps_wavelet_methods::z_slices_to_cp2k_distribution
subroutine, public z_slices_to_cp2k_distribution(density, wavelet, pw_grid)
...
Definition ps_wavelet_methods.F:323
ps_wavelet_methods::cp2k_distribution_to_z_slices
subroutine, public cp2k_distribution_to_z_slices(density, wavelet, pw_grid)
...
Definition ps_wavelet_methods.F:160
ps_wavelet_types
Definition and initialisation of the ps_wavelet data type.
Definition ps_wavelet_types.F:12
ps_wavelet_types::wavelet3d
integer, parameter, public wavelet3d
Definition ps_wavelet_types.F:23
ps_wavelet_types::wavelet0d
integer, parameter, public wavelet0d
Definition ps_wavelet_types.F:23
ps_wavelet_types::wavelet1d
integer, parameter, public wavelet1d
Definition ps_wavelet_types.F:23
ps_wavelet_types::wavelet2d
integer, parameter, public wavelet2d
Definition ps_wavelet_types.F:23
pw_grid_types
Definition pw_grid_types.F:13
pw_grids
This module defines the grid data type and some basic operations on it.
Definition pw_grids.F:36
pw_grids::pw_grid_compare
logical function, public pw_grid_compare(grida, gridb)
Check if two pw_grids are equal.
Definition pw_grids.F:148
pw_grids::pw_grid_release
subroutine, public pw_grid_release(pw_grid)
releases the given pw grid
Definition pw_grids.F:2163
pw_grids::pw_grid_retain
subroutine, public pw_grid_retain(pw_grid)
retains the given pw grid
Definition pw_grids.F:2147
pw_methods
Definition pw_methods.F:25
pw_methods::pw_derive
subroutine, public pw_derive(pw, n)
Calculate the derivative of a plane wave vector.
Definition pw_methods.F:10257
pw_poisson_methods
Definition pw_poisson_methods.F:13
pw_poisson_methods::pw_poisson_set
subroutine, public pw_poisson_set(poisson_env, cell_hmat, parameters, pw_pools, use_level, mt_super_ref_pw_grid, dct_pw_grid, force_rebuild)
sets cell, grids and parameters used by the poisson solver You should call this at least once (and se...
Definition pw_poisson_methods.F:4005
pw_poisson_methods::pw_func_u_convolution
subroutine, public pw_func_u_convolution(poisson_env, func, convolution, u)
computes the convolution f * u using reciprocal-space multiplication
Definition pw_poisson_methods.F:4097
pw_poisson_types
functions related to the poisson solver on regular grids
Definition pw_poisson_types.F:22
pw_poisson_types::pw_poisson_wavelet
integer, parameter, public pw_poisson_wavelet
Definition pw_poisson_types.F:67
pw_poisson_types::pw_green_create
subroutine, public pw_green_create(green, poisson_params, cell_hmat, pw_pool, mt_super_ref_pw_grid, dct_pw_grid)
Allocates and sets up the green functions for the fft based poisson solvers.
Definition pw_poisson_types.F:183
pw_poisson_types::pw_poisson_periodic
integer, parameter, public pw_poisson_periodic
Definition pw_poisson_types.F:67
pw_poisson_types::pw_green_release
subroutine, public pw_green_release(gftype, pw_pool)
destroys the type (deallocates data)
Definition pw_poisson_types.F:421
pw_poisson_types::pw_poisson_none
integer, parameter, public pw_poisson_none
Definition pw_poisson_types.F:67
pw_poisson_types::pw_poisson_mt
integer, parameter, public pw_poisson_mt
Definition pw_poisson_types.F:67
pw_poisson_types::analytic2d
integer, parameter, public analytic2d
Definition pw_poisson_types.F:81
pw_poisson_types::analytic1d
integer, parameter, public analytic1d
Definition pw_poisson_types.F:81
pw_poisson_types::periodic3d
integer, parameter, public periodic3d
Definition pw_poisson_types.F:81
pw_poisson_types::ps_implicit
integer, parameter, public ps_implicit
Definition pw_poisson_types.F:81
pw_poisson_types::analytic0d
integer, parameter, public analytic0d
Definition pw_poisson_types.F:81
pw_poisson_types::multipole0d
integer, parameter, public multipole0d
Definition pw_poisson_types.F:81
pw_poisson_types::pw_poisson_implicit
integer, parameter, public pw_poisson_implicit
Definition pw_poisson_types.F:67
pw_poisson_types::pw_poisson_analytic
integer, parameter, public pw_poisson_analytic
Definition pw_poisson_types.F:67
pw_poisson_types::do_ewald_spme
integer, parameter, public do_ewald_spme
Definition pw_poisson_types.F:76
pw_poisson_types::pw_poisson_multipole
integer, parameter, public pw_poisson_multipole
Definition pw_poisson_types.F:67
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
pw_pool_types::pw_pools_dealloc
subroutine, public pw_pools_dealloc(pools)
deallocates the given pools (releasing each of the underlying pools)
Definition pw_pool_types.F:1268
pw_pool_types::pw_pools_copy
subroutine, public pw_pools_copy(source_pools, target_pools)
copies a multigrid pool, the underlying pools are shared
Definition pw_pool_types.F:1247
pw_types
Definition pw_types.F:24
ps_wavelet_types::ps_wavelet_type
Definition ps_wavelet_types.F:36
pw_grid_types::pw_grid_type
Definition pw_grid_types.F:53
pw_poisson_types::greens_fn_type
contains all the informations needed by the fft based poisson solvers
Definition pw_poisson_types.F:146
pw_poisson_types::pw_poisson_parameter_type
parameters for the poisson solver independet of input_section
Definition pw_poisson_types.F:97
pw_poisson_types::pw_poisson_type
environment for the poisson solver
Definition pw_poisson_types.F:123
pw_pool_types::pw_pool_p_type
to create arrays of pools
Definition pw_pool_types.F:135
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62