(git:b6ef100)
Loading...
Searching...
No Matches
ps_implicit_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!> \brief The implicit (generalized) Poisson solver
10!> \par History
11!> 06.2014 created [Hossein Bani-Hashemian]
12!> 11.2015 - dealt with missing grid points of periodic grids while performing dct;
13!> - revised solver for Neumann and mixed boundary setups.
14!> \author Hossein Bani-Hashemian
15! **************************************************************************************************
16MODULE ps_implicit_methods
17 USE bibliography, ONLY: banihashemian2016,&
18 cite_reference
19 USE cp_log_handling, ONLY: cp_get_default_logger,&
20 cp_logger_get_default_unit_nr,&
21 cp_logger_type
22 USE dct, ONLY: &
23 dct_type, dct_type_init, neumannx, neumannxy, neumannxyz, neumannxz, neumanny, neumannyz, &
24 neumannz, pw_expand, pw_shrink
25 USE dielectric_methods, ONLY: derive_fft,&
26 dielectric_create
27 USE dielectric_types, ONLY: dielectric_type
28 USE dirichlet_bc_methods, ONLY: dirichlet_boundary_region_setup
29 USE dirichlet_bc_types, ONLY: dbc_tile_release
30 USE kahan_sum, ONLY: accurate_sum
31 USE kinds, ONLY: dp,&
32 int_8
33 USE mathconstants, ONLY: fourpi,&
34 pi
35 USE ps_implicit_types, ONLY: mixed_bc,&
36 mixed_periodic_bc,&
37 neumann_bc,&
38 periodic_bc,&
39 ps_implicit_type
40 USE pw_grid_types, ONLY: pw_grid_type
41 USE pw_methods, ONLY: pw_axpy,&
42 pw_copy,&
43 pw_integral_ab,&
44 pw_scale,&
45 pw_transfer,&
46 pw_zero
47 USE pw_poisson_types, ONLY: greens_fn_type,&
48 pw_poisson_parameter_type,&
49 pw_poisson_type
50 USE pw_pool_types, ONLY: pw_pool_create,&
51 pw_pool_release,&
52 pw_pool_type
53 USE pw_types, ONLY: pw_c1d_gs_type,&
54 pw_r3d_rs_type
55#include "../base/base_uses.f90"
56
57 IMPLICIT NONE
58 PRIVATE
59 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ps_implicit_methods'
60
61 PUBLIC ps_implicit_create, &
62 implicit_poisson_solver_periodic, &
63 implicit_poisson_solver_neumann, &
64 implicit_poisson_solver_mixed_periodic, &
65 implicit_poisson_solver_mixed
66
67 INTERFACE ps_implicit_compute_ehartree
68 MODULE PROCEDURE compute_ehartree_periodic_bc, &
69 compute_ehartree_mixed_bc
70 END INTERFACE ps_implicit_compute_ehartree
71
72 REAL(dp), PRIVATE, PARAMETER :: large_error = 1.0e4_dp
73
74CONTAINS
75
76! **************************************************************************************************
77!> \brief Creates implicit Poisson solver environment
78!> \param pw_pool pool of pw grid
79!> \param poisson_params poisson_env parameters
80!> \param dct_pw_grid discrete cosine transform (extended) grid
81!> \param green green function for FFT based inverse Laplacian
82!> \param ps_implicit_env implicit env to be created
83!> \par History
84!> 06.2014 created [Hossein Bani-Hashemian]
85!> \author Mohammad Hossein Bani-Hashemian
86! **************************************************************************************************
87 SUBROUTINE ps_implicit_create(pw_pool, poisson_params, dct_pw_grid, green, ps_implicit_env)
88
89 TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
90 TYPE(pw_poisson_parameter_type), INTENT(INOUT) :: poisson_params
91 TYPE(pw_grid_type), INTENT(IN), POINTER :: dct_pw_grid
92 TYPE(greens_fn_type), INTENT(IN), POINTER :: green
93 TYPE(ps_implicit_type), INTENT(INOUT), POINTER :: ps_implicit_env
94
95 CHARACTER(LEN=*), PARAMETER :: routinen = 'ps_implicit_create'
96
97 INTEGER :: boundary_condition, handle, j, &
98 n_contacts, neumann_directions
99 TYPE(pw_pool_type), POINTER :: pw_pool_xpndd
100
101 CALL timeset(routinen, handle)
102
103 CALL cite_reference(banihashemian2016)
104
105 IF (.NOT. ASSOCIATED(ps_implicit_env)) THEN
106 ALLOCATE (ps_implicit_env)
107
108 ps_implicit_env%do_dbc_cube = poisson_params%dbc_params%do_dbc_cube
109 boundary_condition = poisson_params%ps_implicit_params%boundary_condition
110 neumann_directions = poisson_params%ps_implicit_params%neumann_directions
111
112! create dielectric
113 NULLIFY (ps_implicit_env%dielectric)
114 SELECT CASE (boundary_condition)
115 CASE (periodic_bc, mixed_periodic_bc)
116 CALL dielectric_create(ps_implicit_env%dielectric, pw_pool, poisson_params%dielectric_params)
117 CASE (neumann_bc, mixed_bc)
118 CALL pw_pool_create(pw_pool_xpndd, pw_grid=dct_pw_grid)
119 CALL dielectric_create(ps_implicit_env%dielectric, pw_pool_xpndd, poisson_params%dielectric_params)
120 CALL pw_pool_release(pw_pool_xpndd)
121 END SELECT
122
123! initial guess
124 NULLIFY (ps_implicit_env%initial_guess)
125
126! v_eps
127 NULLIFY (ps_implicit_env%v_eps)
128 ALLOCATE (ps_implicit_env%v_eps)
129 CALL pw_pool%create_pw(ps_implicit_env%v_eps)
130 CALL pw_zero(ps_implicit_env%v_eps)
131
132! constraint charge
133 NULLIFY (ps_implicit_env%cstr_charge)
134 SELECT CASE (boundary_condition)
135 CASE (mixed_periodic_bc)
136 ALLOCATE (ps_implicit_env%cstr_charge)
137 CALL pw_pool%create_pw(ps_implicit_env%cstr_charge)
138 CALL pw_zero(ps_implicit_env%cstr_charge)
139 CASE (mixed_bc)
140 CALL pw_pool_create(pw_pool_xpndd, pw_grid=dct_pw_grid)
141 ALLOCATE (ps_implicit_env%cstr_charge)
142 CALL pw_pool_xpndd%create_pw(ps_implicit_env%cstr_charge)
143 CALL pw_zero(ps_implicit_env%cstr_charge)
144 CALL pw_pool_release(pw_pool_xpndd)
145 END SELECT
146
147! initialize energies
148 ps_implicit_env%ehartree = 0.0_dp
149 ps_implicit_env%electric_enthalpy = 0.0_dp
150! times called
151 ps_implicit_env%times_called = 0
152
153! dct env
154 IF (boundary_condition == mixed_bc .OR. boundary_condition == neumann_bc) THEN
155 CALL dct_type_init(pw_pool%pw_grid, neumann_directions, ps_implicit_env%dct_env)
156 END IF
157
158! prepare dirichlet bc
159 CALL dirichlet_boundary_region_setup(pw_pool, poisson_params, ps_implicit_env%contacts)
160 CALL ps_implicit_prepare_blocks(pw_pool, dct_pw_grid, green, poisson_params, ps_implicit_env)
161 ! release tiles if they are not supposed to be written into cube files
162 IF ((boundary_condition == mixed_periodic_bc .OR. boundary_condition == mixed_bc) .AND. &
163 (.NOT. poisson_params%dbc_params%do_dbc_cube)) THEN
164 n_contacts = SIZE(ps_implicit_env%contacts)
165 DO j = 1, n_contacts
166 CALL dbc_tile_release(ps_implicit_env%contacts(j)%dirichlet_bc, pw_pool)
167 END DO
168 END IF
169
170 END IF
171
172 CALL timestop(handle)
173
174 END SUBROUTINE ps_implicit_create
175
176! **************************************************************************************************
177!> \brief implicit Poisson solver for periodic boundary conditions
178!> \param poisson_env poisson environment
179!> \param density electron density
180!> \param v_new electrostatic potential
181!> \param ehartree Hartree energy
182!> \par History
183!> 07.2014 created [Hossein Bani-Hashemian]
184!> \author Mohammad Hossein Bani-Hashemian
185! **************************************************************************************************
186 SUBROUTINE implicit_poisson_solver_periodic(poisson_env, density, v_new, ehartree)
187
188 TYPE(pw_poisson_type), INTENT(IN) :: poisson_env
189 TYPE(pw_r3d_rs_type), INTENT(IN) :: density
190 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: v_new
191 REAL(dp), INTENT(OUT), OPTIONAL :: ehartree
192
193 CHARACTER(LEN=*), PARAMETER :: routinen = 'implicit_poisson_solver_periodic'
194
195 INTEGER :: handle, iter, max_iter, outp_unit, &
196 times_called
197 LOGICAL :: reached_max_iter, reached_tol, &
198 use_zero_initial_guess
199 REAL(dp) :: nabs_error, omega, pres_error, tol
200 TYPE(dielectric_type), POINTER :: dielectric
201 TYPE(greens_fn_type), POINTER :: green
202 TYPE(ps_implicit_type), POINTER :: ps_implicit_env
203 TYPE(pw_pool_type), POINTER :: pw_pool
204 TYPE(pw_r3d_rs_type) :: g, pxqainvxres, qainvxres, res_new, &
205 res_old, v_old
206
207 CALL timeset(routinen, handle)
208
209 pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
210 dielectric => poisson_env%implicit_env%dielectric
211 green => poisson_env%green_fft
212 ps_implicit_env => poisson_env%implicit_env
213
214 tol = poisson_env%parameters%ps_implicit_params%tol
215 omega = poisson_env%parameters%ps_implicit_params%omega
216 max_iter = poisson_env%parameters%ps_implicit_params%max_iter
217 use_zero_initial_guess = poisson_env%parameters%ps_implicit_params%zero_initial_guess
218 times_called = ps_implicit_env%times_called
219
220! check if this is the first scf iteration
221 IF (times_called == 0) CALL ps_implicit_initial_guess_create(ps_implicit_env, pw_pool)
222
223 CALL pw_pool%create_pw(g)
224 CALL pw_pool%create_pw(v_old)
225 CALL pw_pool%create_pw(res_old)
226 CALL pw_pool%create_pw(res_new)
227 CALL pw_pool%create_pw(qainvxres)
228 CALL pw_pool%create_pw(pxqainvxres)
229
230 IF (use_zero_initial_guess) THEN
231 CALL pw_zero(v_old)
232 ELSE
233 CALL pw_copy(ps_implicit_env%initial_guess, v_old)
234 END IF
235
236 g%array = fourpi*density%array/dielectric%eps%array
237
238! res_old = g - \Delta(v_old) - P(v_old)
239 CALL apply_poisson_operator_fft(pw_pool, green, dielectric, v_old, res_old)
240 CALL pw_scale(res_old, -1.0_dp)
241 CALL pw_axpy(g, res_old)
242
243! evaluate \Delta^-1(res_old)
244 CALL apply_inv_laplace_operator_fft(pw_pool, green, res_old, qainvxres)
245
246 iter = 1
247 DO
248
249! v_new = v_old + \omega * QAinvxres_old
250 CALL pw_scale(qainvxres, omega)
251 CALL pw_copy(qainvxres, v_new)
252 CALL pw_axpy(v_old, v_new)
253
254! res_new = res_old - \omega * ( \Delta(QAinvxres_old) + P(QAinvxres_old) )
255! = (1 - \omega) * res_old - \omega * PxQAinvxres
256 CALL apply_p_operator(pw_pool, dielectric, qainvxres, pxqainvxres)
257 CALL pw_copy(pxqainvxres, res_new)
258 CALL pw_scale(res_new, -1.0_dp)
259 CALL pw_axpy(res_old, res_new, 1.0_dp - omega)
260
261! compute the error
262 CALL ps_implicit_compute_error_fft(pw_pool, green, res_new, v_old, v_new, qainvxres, &
263 pres_error, nabs_error)
264! output
265 CALL ps_implicit_output(iter, pres_error, nabs_error, outp_unit)
266 IF (PRESENT(ehartree)) THEN
267 CALL ps_implicit_compute_ehartree(density, v_new, ehartree)
268 CALL ps_implicit_report_ehartree(ps_implicit_env, outp_unit, ehartree)
269 ps_implicit_env%ehartree = ehartree
270 ELSE
271 IF (outp_unit > 0) WRITE (outp_unit, '(A1,/)')
272 END IF
273
274 iter = iter + 1
275 reached_max_iter = iter > max_iter
276 reached_tol = pres_error <= tol
277 IF (pres_error > large_error) THEN
278 cpabort("Poisson solver did not converge.")
279 END IF
280 ps_implicit_env%times_called = ps_implicit_env%times_called + 1
281 IF (reached_max_iter .OR. reached_tol) EXIT
282
283! v_old = v_new, res_old = res_new
284 CALL pw_copy(v_new, v_old)
285 CALL pw_copy(res_new, res_old)
286
287 END DO
288 CALL ps_implicit_print_convergence_msg(iter, max_iter, outp_unit)
289
290 IF ((times_called /= 0) .AND. (.NOT. use_zero_initial_guess)) THEN
291 CALL pw_copy(v_new, ps_implicit_env%initial_guess)
292 END IF
293
294 IF (PRESENT(ehartree)) ehartree = ps_implicit_env%ehartree
295! compute the extra contribution to the Hamiltonian due to the presence of dielectric
296 block
297 TYPE(pw_r3d_rs_type) :: v_eps
298 v_eps%pw_grid => ps_implicit_env%v_eps%pw_grid
299 v_eps%array => ps_implicit_env%v_eps%array
300 CALL ps_implicit_compute_veps(pw_pool, dielectric, v_new, v_eps)
301 END block
302
303 CALL pw_pool%give_back_pw(g)
304 CALL pw_pool%give_back_pw(v_old)
305 CALL pw_pool%give_back_pw(res_old)
306 CALL pw_pool%give_back_pw(res_new)
307 CALL pw_pool%give_back_pw(qainvxres)
308 CALL pw_pool%give_back_pw(pxqainvxres)
309
310 CALL timestop(handle)
311
312 END SUBROUTINE implicit_poisson_solver_periodic
313
314! **************************************************************************************************
315!> \brief implicit Poisson solver: zero-average solution of the Poisson equation
316!> subject to homogeneous Neumann boundary conditions
317!> \param poisson_env poisson environment
318!> \param density electron density
319!> \param v_new electrostatic potential
320!> \param ehartree Hartree energy
321!> \par History
322!> 02.2015 created [Hossein Bani-Hashemian]
323!> 11.2015 revised [Hossein Bani-Hashemian]
324!> \author Mohammad Hossein Bani-Hashemian
325! **************************************************************************************************
326 SUBROUTINE implicit_poisson_solver_neumann(poisson_env, density, v_new, ehartree)
327
328 TYPE(pw_poisson_type), INTENT(IN) :: poisson_env
329 TYPE(pw_r3d_rs_type), INTENT(IN) :: density
330 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: v_new
331 REAL(dp), INTENT(OUT), OPTIONAL :: ehartree
332
333 CHARACTER(LEN=*), PARAMETER :: routinen = 'implicit_poisson_solver_neumann'
334
335 INTEGER :: handle, iter, max_iter, &
336 neumann_directions, outp_unit, &
337 times_called
338 LOGICAL :: reached_max_iter, reached_tol, &
339 use_zero_initial_guess
340 REAL(dp) :: nabs_error, omega, pres_error, tol, &
341 vol_scfac
342 TYPE(dct_type), POINTER :: dct_env
343 TYPE(dielectric_type), POINTER :: dielectric
344 TYPE(greens_fn_type), POINTER :: green
345 TYPE(ps_implicit_type), POINTER :: ps_implicit_env
346 TYPE(pw_pool_type), POINTER :: pw_pool, pw_pool_xpndd
347 TYPE(pw_r3d_rs_type) :: density_xpndd, g, pxqainvxres, &
348 qainvxres, res_new, res_old, &
349 v_eps_xpndd, v_new_xpndd, v_old
350
351 CALL timeset(routinen, handle)
352
353 pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
354 dielectric => poisson_env%implicit_env%dielectric
355 green => poisson_env%green_fft
356 ps_implicit_env => poisson_env%implicit_env
357 dct_env => ps_implicit_env%dct_env
358
359 tol = poisson_env%parameters%ps_implicit_params%tol
360 omega = poisson_env%parameters%ps_implicit_params%omega
361 max_iter = poisson_env%parameters%ps_implicit_params%max_iter
362 use_zero_initial_guess = poisson_env%parameters%ps_implicit_params%zero_initial_guess
363 neumann_directions = poisson_env%parameters%ps_implicit_params%neumann_directions
364 times_called = ps_implicit_env%times_called
365
366 SELECT CASE (neumann_directions)
367 CASE (neumannxyz)
368 vol_scfac = 8.0_dp
369 CASE (neumannxy, neumannxz, neumannyz)
370 vol_scfac = 4.0_dp
371 CASE (neumannx, neumanny, neumannz)
372 vol_scfac = 2.0_dp
373 END SELECT
374
375 CALL pw_pool_create(pw_pool_xpndd, pw_grid=poisson_env%dct_pw_grid)
376
377! check if this is the first scf iteration
378 IF (times_called == 0) CALL ps_implicit_initial_guess_create(ps_implicit_env, pw_pool_xpndd)
379
380 CALL pw_pool_xpndd%create_pw(g)
381 CALL pw_pool_xpndd%create_pw(v_old)
382 CALL pw_pool_xpndd%create_pw(res_old)
383 CALL pw_pool_xpndd%create_pw(res_new)
384 CALL pw_pool_xpndd%create_pw(qainvxres)
385 CALL pw_pool_xpndd%create_pw(pxqainvxres)
386 CALL pw_pool_xpndd%create_pw(density_xpndd)
387 CALL pw_pool_xpndd%create_pw(v_new_xpndd)
388 CALL pw_pool_xpndd%create_pw(v_eps_xpndd)
389
390 IF (use_zero_initial_guess) THEN
391 CALL pw_zero(v_old)
392 ELSE
393 CALL pw_copy(ps_implicit_env%initial_guess, v_old)
394 END IF
395
396 CALL pw_expand(neumann_directions, &
397 dct_env%recv_msgs_bnds, dct_env%dests_expand, dct_env%srcs_expand, &
398 dct_env%flipg_stat, dct_env%bounds_shftd, density, density_xpndd)
399 CALL pw_expand(neumann_directions, &
400 dct_env%recv_msgs_bnds, dct_env%dests_expand, dct_env%srcs_expand, &
401 dct_env%flipg_stat, dct_env%bounds_shftd, v_new, v_new_xpndd)
402
403 g%array = fourpi*density_xpndd%array/dielectric%eps%array
404
405! res_old = g - \Delta(v_old) - P(v_old)
406 CALL apply_poisson_operator_dct(pw_pool_xpndd, green, dielectric, v_old, res_old)
407 CALL pw_scale(res_old, -1.0_dp)
408 CALL pw_axpy(g, res_old)
409
410! evaluate \Delta^-1(res_old)
411 CALL apply_inv_laplace_operator_dct(pw_pool_xpndd, green, res_old, qainvxres)
412
413 iter = 1
414 DO
415
416! v_new = v_old + \omega * QAinvxres_old
417 CALL pw_scale(qainvxres, omega)
418 CALL pw_copy(qainvxres, v_new_xpndd)
419 CALL pw_axpy(v_old, v_new_xpndd)
420
421! res_new = res_old - \omega * ( \Delta(QAinvxres_old) + P(QAinvxres_old) )
422! = (1 - \omega) * res_old - \omega * PxQAinvxres
423 CALL apply_p_operator(pw_pool_xpndd, dielectric, qainvxres, pxqainvxres)
424 CALL pw_copy(pxqainvxres, res_new)
425 CALL pw_scale(res_new, -1.0_dp)
426 CALL pw_axpy(res_old, res_new, 1.0_dp - omega)
427
428! compute the error
429 CALL ps_implicit_compute_error_dct(pw_pool_xpndd, green, res_new, v_old, v_new_xpndd, qainvxres, &
430 pres_error, nabs_error)
431! output
432 CALL ps_implicit_output(iter, pres_error, nabs_error, outp_unit)
433 IF (PRESENT(ehartree)) THEN
434 CALL ps_implicit_compute_ehartree(density_xpndd, v_new_xpndd, ehartree)
435 CALL ps_implicit_report_ehartree(ps_implicit_env, outp_unit, ehartree/vol_scfac)
436 ps_implicit_env%ehartree = ehartree/vol_scfac
437 ELSE
438 IF (outp_unit > 0) WRITE (outp_unit, '(A1,/)')
439 END IF
440
441 iter = iter + 1
442 reached_max_iter = iter > max_iter
443 reached_tol = pres_error <= tol
444 IF (pres_error > large_error) THEN
445 cpabort("Poisson solver did not converge.")
446 END IF
447 ps_implicit_env%times_called = ps_implicit_env%times_called + 1
448 IF (reached_max_iter .OR. reached_tol) EXIT
449
450! v_old = v_new, res_old = res_new
451 CALL pw_copy(v_new_xpndd, v_old)
452 CALL pw_copy(res_new, res_old)
453
454 END DO
455 CALL ps_implicit_print_convergence_msg(iter, max_iter, outp_unit)
456
457 CALL pw_shrink(neumann_directions, dct_env%dests_shrink, dct_env%srcs_shrink, &
458 dct_env%bounds_local_shftd, v_new_xpndd, v_new)
459
460 IF ((times_called /= 0) .AND. (.NOT. use_zero_initial_guess)) THEN
461 CALL pw_copy(v_new_xpndd, ps_implicit_env%initial_guess)
462 END IF
463
464 IF (PRESENT(ehartree)) ehartree = ps_implicit_env%ehartree
465! compute the extra contribution to the Hamiltonian due to the presence of dielectric
466! veps has to be computed for the expanded data and then shrunk otherwise we loose accuracy
467 CALL ps_implicit_compute_veps(pw_pool_xpndd, dielectric, v_new_xpndd, v_eps_xpndd)
468 block
469 TYPE(pw_r3d_rs_type) :: v_eps
470 v_eps%pw_grid => ps_implicit_env%v_eps%pw_grid
471 v_eps%array => ps_implicit_env%v_eps%array
472 CALL pw_shrink(neumann_directions, dct_env%dests_shrink, dct_env%srcs_shrink, &
473 dct_env%bounds_local_shftd, v_eps_xpndd, v_eps)
474 END block
475
476 CALL pw_pool_xpndd%give_back_pw(g)
477 CALL pw_pool_xpndd%give_back_pw(v_old)
478 CALL pw_pool_xpndd%give_back_pw(res_old)
479 CALL pw_pool_xpndd%give_back_pw(res_new)
480 CALL pw_pool_xpndd%give_back_pw(qainvxres)
481 CALL pw_pool_xpndd%give_back_pw(pxqainvxres)
482 CALL pw_pool_xpndd%give_back_pw(density_xpndd)
483 CALL pw_pool_xpndd%give_back_pw(v_new_xpndd)
484 CALL pw_pool_xpndd%give_back_pw(v_eps_xpndd)
485 CALL pw_pool_release(pw_pool_xpndd)
486
487 CALL timestop(handle)
488
489 END SUBROUTINE implicit_poisson_solver_neumann
490
491! **************************************************************************************************
492!> \brief implicit Poisson solver for mixed-periodic boundary conditions (periodic + Dirichlet)
493!> \param poisson_env poisson environment
494!> \param density electron density
495!> \param v_new electrostatic potential
496!> \param electric_enthalpy electric enthalpy
497!> \par History
498!> 07.2014 created [Hossein Bani-Hashemian]
499!> \author Mohammad Hossein Bani-Hashemian
500! **************************************************************************************************
501 SUBROUTINE implicit_poisson_solver_mixed_periodic(poisson_env, density, v_new, electric_enthalpy)
502
503 TYPE(pw_poisson_type), INTENT(IN) :: poisson_env
504 TYPE(pw_r3d_rs_type), INTENT(IN) :: density
505 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: v_new
506 REAL(dp), INTENT(OUT), OPTIONAL :: electric_enthalpy
507
508 CHARACTER(LEN=*), PARAMETER :: routinen = 'implicit_poisson_solver_mixed_periodic'
509
510 INTEGER :: data_size, handle, iter, j, lb1, lb2, lb3, max_iter, n_contacts, n_tiles_tot, ng, &
511 ngpts_local, nt, nt_tot, outp_unit, times_called, ub1, ub2, ub3
512 INTEGER(KIND=int_8) :: ngpts
513 INTEGER, DIMENSION(2, 3) :: bounds_local
514 INTEGER, DIMENSION(3) :: npts_local
515 LOGICAL :: reached_max_iter, reached_tol, &
516 use_zero_initial_guess
517 REAL(dp) :: axvbar_avg, ehartree, eta, g_avg, &
518 gminusaxvbar_avg, nabs_error, omega, &
519 pres_error, tol
520 REAL(dp), ALLOCATABLE, DIMENSION(:) :: btxlambda_new, btxlambda_old, bxv_bar, bxv_new, &
521 lambda0, lambda_new, lambda_newneta, lambda_old, qsxlambda, v_bar1d, v_d, v_new1d, w
522 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: b, bt, qs, rinv
523 REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: btxlambda_new3d, btxlambda_old3d
524 TYPE(dielectric_type), POINTER :: dielectric
525 TYPE(greens_fn_type), POINTER :: green
526 TYPE(ps_implicit_type), POINTER :: ps_implicit_env
527 TYPE(pw_grid_type), POINTER :: pw_grid
528 TYPE(pw_pool_type), POINTER :: pw_pool
529 TYPE(pw_r3d_rs_type) :: axvbar, g, pxqainvxres, qainvxres, &
530 res_new, res_old, v_old
531
532 CALL timeset(routinen, handle)
533
534 pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
535 pw_grid => pw_pool%pw_grid
536 dielectric => poisson_env%implicit_env%dielectric
537 green => poisson_env%green_fft
538 ps_implicit_env => poisson_env%implicit_env
539
540 ngpts_local = pw_grid%ngpts_local
541 ngpts = pw_grid%ngpts
542 npts_local = pw_grid%npts_local
543 bounds_local = pw_grid%bounds_local
544 tol = poisson_env%parameters%ps_implicit_params%tol
545 omega = poisson_env%parameters%ps_implicit_params%omega
546 max_iter = poisson_env%parameters%ps_implicit_params%max_iter
547 use_zero_initial_guess = poisson_env%parameters%ps_implicit_params%zero_initial_guess
548 times_called = ps_implicit_env%times_called
549
550 n_contacts = SIZE(ps_implicit_env%contacts)
551 n_tiles_tot = 0
552 DO j = 1, n_contacts
553 n_tiles_tot = n_tiles_tot + ps_implicit_env%contacts(j)%dirichlet_bc%n_tiles
554 END DO
555
556 IF (pw_grid%para%blocked) THEN
557 data_size = product(npts_local)
558 ELSE IF (pw_grid%para%ray_distribution) THEN
559 data_size = ngpts_local
560 ELSE ! parallel run with np = 1
561 data_size = product(npts_local)
562 END IF
563
564! check if this is the first scf iteration
565 IF (times_called == 0) CALL ps_implicit_initial_guess_create(ps_implicit_env, pw_pool)
566
567 ALLOCATE (b(n_tiles_tot, data_size))
568 ALLOCATE (bt(data_size, n_tiles_tot))
569 ALLOCATE (qs(n_tiles_tot, n_tiles_tot))
570 ALLOCATE (rinv(n_tiles_tot + 1, n_tiles_tot + 1))
571
572 b(:, :) = ps_implicit_env%B
573 bt(:, :) = ps_implicit_env%Bt
574 qs(:, :) = ps_implicit_env%QS
575 rinv(:, :) = ps_implicit_env%Rinv
576 CALL get_voltage(poisson_env%parameters%dbc_params%time, ps_implicit_env%v_D, ps_implicit_env%osc_frac, &
577 ps_implicit_env%frequency, ps_implicit_env%phase, v_d)
578
579 lb1 = bounds_local(1, 1); ub1 = bounds_local(2, 1)
580 lb2 = bounds_local(1, 2); ub2 = bounds_local(2, 2)
581 lb3 = bounds_local(1, 3); ub3 = bounds_local(2, 3)
582
583 ALLOCATE (lambda0(n_tiles_tot), lambda_old(n_tiles_tot), lambda_new(n_tiles_tot))
584 ALLOCATE (btxlambda_old(data_size), btxlambda_new(data_size))
585 ALLOCATE (btxlambda_old3d(lb1:ub1, lb2:ub2, lb3:ub3), btxlambda_new3d(lb1:ub1, lb2:ub2, lb3:ub3))
586 ALLOCATE (qsxlambda(n_tiles_tot))
587 ALLOCATE (w(n_tiles_tot + 1))
588 ALLOCATE (lambda_newneta(n_tiles_tot + 1))
589 ALLOCATE (v_bar1d(data_size))
590 ALLOCATE (bxv_bar(n_tiles_tot))
591
592 ALLOCATE (v_new1d(data_size))
593 ALLOCATE (bxv_new(n_tiles_tot))
594
595 CALL pw_pool%create_pw(g)
596 CALL pw_pool%create_pw(v_old)
597 CALL pw_pool%create_pw(res_old)
598 CALL pw_pool%create_pw(res_new)
599 CALL pw_pool%create_pw(qainvxres)
600 CALL pw_pool%create_pw(pxqainvxres)
601 CALL pw_pool%create_pw(axvbar)
602
603 IF (use_zero_initial_guess) THEN
604 CALL pw_zero(v_old)
605 lambda0 = 0.0_dp
606 ELSE
607 CALL pw_copy(ps_implicit_env%initial_guess, v_old)
608 lambda0(:) = ps_implicit_env%initial_lambda
609 END IF
610
611 g%array = fourpi*density%array/dielectric%eps%array
612 g_avg = accurate_sum(g%array)/ngpts
613
614 lambda_old(:) = lambda0
615
616! res_old = g - \Delta(v_old) - P(v_old) - B^t * \lambda_old
617 CALL apply_poisson_operator_fft(pw_pool, green, dielectric, v_old, res_old)
618 CALL pw_scale(res_old, -1.0_dp)
619 CALL pw_axpy(g, res_old)
620 IF (data_size /= 0) THEN
621 CALL dgemv('N', data_size, n_tiles_tot, 1.0_dp, bt, data_size, lambda_old, 1, 0.0_dp, btxlambda_old, 1)
622 END IF
623 CALL convert_1dto3d(ps_implicit_env%idx_1dto3d, btxlambda_old, btxlambda_old3d)
624 res_old%array = res_old%array - btxlambda_old3d
625
626! evaluate \Delta^-1(res_old)
627 CALL apply_inv_laplace_operator_fft(pw_pool, green, res_old, qainvxres)
628
629 iter = 1
630 DO
631
632! v_new (v_bar) = v_old + \omega * QAinvxres_old
633 CALL pw_scale(qainvxres, omega)
634 CALL pw_copy(qainvxres, v_new)
635 CALL pw_axpy(v_old, v_new)
636
637! evaluate 1^t * (g - \Delta(\bar{v}) - P(\bar{v}))
638! = 1^t * (g - P(\bar{v}))
639 CALL apply_p_operator(pw_pool, dielectric, v_new, axvbar)
640 axvbar_avg = accurate_sum(axvbar%array)/ngpts
641 gminusaxvbar_avg = g_avg - axvbar_avg
642 CALL pw_grid%para%group%sum(gminusaxvbar_avg)
643
644! evaluate Q_S * \lambda + v_D - B * \bar{v}
645 CALL dgemv('N', n_tiles_tot, n_tiles_tot, 1.0_dp, qs, n_tiles_tot, lambda_old, 1, 0.0_dp, qsxlambda, 1)
646 v_bar1d(ps_implicit_env%idx_1dto3d) = reshape(v_new%array, [data_size])
647 CALL dgemv('N', n_tiles_tot, data_size, 1.0_dp, b, n_tiles_tot, v_bar1d, 1, 0.0_dp, bxv_bar, 1)
648 CALL pw_grid%para%group%sum(bxv_bar)
649! solve R [\lambda; \eta] = [Q_S * \lambda + v_D - B * \bar{v}; 1^t * (g - \Delta(\bar{v}) - P(\bar{v}))]
650 w = 0.0_dp
651 w(:) = [qsxlambda + v_d - bxv_bar, gminusaxvbar_avg]
652 CALL dgemv('N', n_tiles_tot + 1, n_tiles_tot + 1, 1.0_dp, rinv, n_tiles_tot + 1, w, 1, 0.0_dp, lambda_newneta, 1)
653 lambda_new(:) = lambda_newneta(1:n_tiles_tot)
654 eta = lambda_newneta(n_tiles_tot + 1)
655
656! v_new = v_bar + 1 * \eta
657 v_new%array = v_new%array + eta/ngpts
658
659! evaluate B^t * \lambda_new
660 IF (data_size /= 0) THEN
661 CALL dgemv('N', data_size, n_tiles_tot, 1.0_dp, bt, data_size, lambda_new, 1, 0.0_dp, btxlambda_new, 1)
662 END IF
663 CALL convert_1dto3d(ps_implicit_env%idx_1dto3d, btxlambda_new, btxlambda_new3d)
664
665! res_new = res_old - \omega * ( \Delta(QAinvxres_old) + P(QAinvxres_old) ) - B^t * ( \lambda_new - \lambda_old )
666! = (1 - \omega) * res_old - \omega * P(QAinvxres_old) - B^t * ( \lambda_new - \lambda_old )
667 CALL pw_zero(res_new)
668 CALL apply_p_operator(pw_pool, dielectric, qainvxres, pxqainvxres)
669 CALL pw_axpy(pxqainvxres, res_new, -1.0_dp)
670 CALL pw_axpy(res_old, res_new, 1.0_dp - omega)
671 res_new%array = res_new%array + btxlambda_old3d - btxlambda_new3d
672
673! compute the error
674 CALL ps_implicit_compute_error_fft(pw_pool, green, res_new, v_old, v_new, qainvxres, &
675 pres_error, nabs_error)
676! output
677 CALL ps_implicit_output(iter, pres_error, nabs_error, outp_unit)
678 IF (PRESENT(electric_enthalpy)) THEN
679 CALL ps_implicit_compute_ehartree(dielectric, density, btxlambda_new3d, v_new, ehartree, electric_enthalpy)
680 CALL ps_implicit_report_ehartree(ps_implicit_env, outp_unit, ehartree)
681 ps_implicit_env%ehartree = ehartree
682 ps_implicit_env%electric_enthalpy = electric_enthalpy
683 ELSE
684 IF (outp_unit > 0) WRITE (outp_unit, '(A1,/)')
685 END IF
686
687! verbose output
688 IF (poisson_env%parameters%dbc_params%verbose_output) THEN
689 v_new1d(ps_implicit_env%idx_1dto3d) = reshape(v_new%array, [data_size])
690 CALL dgemv('N', n_tiles_tot, data_size, 1.0_dp, b, n_tiles_tot, v_new1d, 1, 0.0_dp, bxv_new, 1)
691 CALL pw_grid%para%group%sum(bxv_new)
692 IF (outp_unit > 0) THEN
693 WRITE (outp_unit, '(T3,A,A)') "======== verbose ", repeat('=', 61)
694 WRITE (outp_unit, '(T20,A)') "Drgn tile vhartree lambda "
695 WRITE (outp_unit, '(T19,A)') repeat('-', 46)
696 nt_tot = 1
697 DO ng = 1, n_contacts
698 DO nt = 1, ps_implicit_env%contacts(ng)%dirichlet_bc%n_tiles
699 WRITE (outp_unit, '(T17,I6,5X,I6,3X,E13.4,E13.4)') ng, nt, bxv_new(nt_tot), lambda_new(nt_tot)
700 nt_tot = nt_tot + 1
701 END DO
702 END DO
703 WRITE (outp_unit, '(T3,A)') repeat('=', 78)
704 END IF
705 END IF
706
707! check the convergence
708 iter = iter + 1
709 reached_max_iter = iter > max_iter
710 reached_tol = pres_error <= tol
711 ps_implicit_env%times_called = ps_implicit_env%times_called + 1
712 IF (pres_error > large_error) THEN
713 cpabort("Poisson solver did not converge.")
714 END IF
715 IF (reached_max_iter .OR. reached_tol) EXIT
716
717! update
718 CALL pw_copy(v_new, v_old)
719 lambda_old(:) = lambda_new
720 CALL pw_copy(res_new, res_old)
721 btxlambda_old3d(:, :, :) = btxlambda_new3d
722
723 END DO
724 CALL ps_implicit_print_convergence_msg(iter, max_iter, outp_unit)
725
726 IF ((times_called /= 0) .AND. (.NOT. use_zero_initial_guess)) THEN
727 CALL pw_copy(v_new, ps_implicit_env%initial_guess)
728 ps_implicit_env%initial_lambda(:) = lambda_new
729 END IF
730
731 ps_implicit_env%cstr_charge%array = btxlambda_new3d
732 IF (PRESENT(electric_enthalpy)) electric_enthalpy = ps_implicit_env%electric_enthalpy
733! compute the extra contribution to the Hamiltonian due to the presence of dielectric
734 block
735 TYPE(pw_r3d_rs_type) :: tmp
736 tmp%pw_grid => ps_implicit_env%v_eps%pw_grid
737 tmp%array => ps_implicit_env%v_eps%array
738 CALL ps_implicit_compute_veps(pw_pool, dielectric, v_new, tmp)
739 END block
740
741 CALL pw_pool%give_back_pw(g)
742 CALL pw_pool%give_back_pw(v_old)
743 CALL pw_pool%give_back_pw(res_old)
744 CALL pw_pool%give_back_pw(res_new)
745 CALL pw_pool%give_back_pw(qainvxres)
746 CALL pw_pool%give_back_pw(pxqainvxres)
747 CALL pw_pool%give_back_pw(axvbar)
748
749 CALL timestop(handle)
750
751 END SUBROUTINE implicit_poisson_solver_mixed_periodic
752
753! **************************************************************************************************
754!> \brief implicit Poisson solver for mixed boundary conditions (Neumann + Dirichlet)
755!> \param poisson_env poisson environment
756!> \param density electron density
757!> \param v_new electrostatic potential
758!> \param electric_enthalpy electric enthalpy
759!> \par History
760!> 10.2014 created [Hossein Bani-Hashemian]
761!> 11.2015 revised [Hossein Bani-Hashemian]
762!> \author Mohammad Hossein Bani-Hashemian
763! **************************************************************************************************
764 SUBROUTINE implicit_poisson_solver_mixed(poisson_env, density, v_new, electric_enthalpy)
765
766 TYPE(pw_poisson_type), INTENT(IN) :: poisson_env
767 TYPE(pw_r3d_rs_type), INTENT(IN) :: density
768 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: v_new
769 REAL(dp), INTENT(OUT), OPTIONAL :: electric_enthalpy
770
771 CHARACTER(LEN=*), PARAMETER :: routinen = 'implicit_poisson_solver_mixed'
772
773 INTEGER :: data_size, handle, iter, j, lb1, lb2, lb3, max_iter, n_contacts, n_tiles_tot, &
774 neumann_directions, ng, ngpts_local, nt, nt_tot, outp_unit, times_called, ub1, ub2, ub3
775 INTEGER(KIND=int_8) :: ngpts
776 INTEGER, DIMENSION(2, 3) :: bounds_local
777 INTEGER, DIMENSION(3) :: npts_local
778 LOGICAL :: reached_max_iter, reached_tol, &
779 use_zero_initial_guess
780 REAL(dp) :: axvbar_avg, ehartree, eta, g_avg, &
781 gminusaxvbar_avg, nabs_error, omega, &
782 pres_error, tol, vol_scfac
783 REAL(dp), ALLOCATABLE, DIMENSION(:) :: btxlambda_new, btxlambda_old, bxv_bar, bxv_new, &
784 lambda0, lambda_new, lambda_newneta, lambda_old, qsxlambda, v_bar1d, v_d, v_new1d, w
785 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: b, bt, qs, rinv
786 REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: btxlambda_new3d, btxlambda_old3d
787 TYPE(dct_type), POINTER :: dct_env
788 TYPE(dielectric_type), POINTER :: dielectric
789 TYPE(greens_fn_type), POINTER :: green
790 TYPE(ps_implicit_type), POINTER :: ps_implicit_env
791 TYPE(pw_grid_type), POINTER :: dct_pw_grid, pw_grid
792 TYPE(pw_pool_type), POINTER :: pw_pool, pw_pool_xpndd
793 TYPE(pw_r3d_rs_type) :: axvbar, density_xpndd, g, pxqainvxres, &
794 qainvxres, res_new, res_old, &
795 v_eps_xpndd, v_new_xpndd, v_old
796
797 CALL timeset(routinen, handle)
798
799 pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
800 pw_grid => pw_pool%pw_grid
801 dielectric => poisson_env%implicit_env%dielectric
802 green => poisson_env%green_fft
803 ps_implicit_env => poisson_env%implicit_env
804 dct_env => ps_implicit_env%dct_env
805
806 dct_pw_grid => poisson_env%dct_pw_grid
807 ngpts_local = dct_pw_grid%ngpts_local
808 ngpts = dct_pw_grid%ngpts
809 npts_local = dct_pw_grid%npts_local
810 bounds_local = dct_pw_grid%bounds_local
811 tol = poisson_env%parameters%ps_implicit_params%tol
812 omega = poisson_env%parameters%ps_implicit_params%omega
813 max_iter = poisson_env%parameters%ps_implicit_params%max_iter
814 use_zero_initial_guess = poisson_env%parameters%ps_implicit_params%zero_initial_guess
815 neumann_directions = poisson_env%parameters%ps_implicit_params%neumann_directions
816 times_called = ps_implicit_env%times_called
817
818 SELECT CASE (neumann_directions)
819 CASE (neumannxyz)
820 vol_scfac = 8.0_dp
821 CASE (neumannxy, neumannxz, neumannyz)
822 vol_scfac = 4.0_dp
823 CASE (neumannx, neumanny, neumannz)
824 vol_scfac = 2.0_dp
825 END SELECT
826
827 n_contacts = SIZE(ps_implicit_env%contacts)
828 n_tiles_tot = 0
829 DO j = 1, n_contacts
830 n_tiles_tot = n_tiles_tot + ps_implicit_env%contacts(j)%dirichlet_bc%n_tiles
831 END DO
832
833 IF (dct_pw_grid%para%blocked) THEN
834 data_size = product(npts_local)
835 ELSE IF (dct_pw_grid%para%ray_distribution) THEN
836 data_size = ngpts_local
837 ELSE ! parallel run with np = 1
838 data_size = product(npts_local)
839 END IF
840
841 CALL pw_pool_create(pw_pool_xpndd, pw_grid=dct_pw_grid)
842
843! check if this is the first scf iteration
844 IF (times_called == 0) CALL ps_implicit_initial_guess_create(ps_implicit_env, pw_pool_xpndd)
845
846 ALLOCATE (b(n_tiles_tot, data_size))
847 ALLOCATE (bt(data_size, n_tiles_tot))
848 ALLOCATE (qs(n_tiles_tot, n_tiles_tot))
849 ALLOCATE (rinv(n_tiles_tot + 1, n_tiles_tot + 1))
850
851 b(:, :) = ps_implicit_env%B
852 bt(:, :) = ps_implicit_env%Bt
853 qs(:, :) = ps_implicit_env%QS
854 rinv(:, :) = ps_implicit_env%Rinv
855 CALL get_voltage(poisson_env%parameters%dbc_params%time, ps_implicit_env%v_D, ps_implicit_env%osc_frac, &
856 ps_implicit_env%frequency, ps_implicit_env%phase, v_d)
857
858 lb1 = bounds_local(1, 1); ub1 = bounds_local(2, 1)
859 lb2 = bounds_local(1, 2); ub2 = bounds_local(2, 2)
860 lb3 = bounds_local(1, 3); ub3 = bounds_local(2, 3)
861
862 ALLOCATE (lambda0(n_tiles_tot), lambda_old(n_tiles_tot), lambda_new(n_tiles_tot))
863 ALLOCATE (btxlambda_old(data_size), btxlambda_new(data_size))
864 ALLOCATE (btxlambda_old3d(lb1:ub1, lb2:ub2, lb3:ub3), btxlambda_new3d(lb1:ub1, lb2:ub2, lb3:ub3))
865 ALLOCATE (qsxlambda(n_tiles_tot))
866 ALLOCATE (w(n_tiles_tot + 1))
867 ALLOCATE (lambda_newneta(n_tiles_tot + 1))
868 ALLOCATE (v_bar1d(data_size))
869 ALLOCATE (bxv_bar(n_tiles_tot))
870
871 ALLOCATE (v_new1d(data_size))
872 ALLOCATE (bxv_new(n_tiles_tot))
873
874 CALL pw_pool_xpndd%create_pw(g)
875 CALL pw_pool_xpndd%create_pw(v_old)
876 CALL pw_pool_xpndd%create_pw(res_old)
877 CALL pw_pool_xpndd%create_pw(res_new)
878 CALL pw_pool_xpndd%create_pw(qainvxres)
879 CALL pw_pool_xpndd%create_pw(pxqainvxres)
880 CALL pw_pool_xpndd%create_pw(axvbar)
881 CALL pw_pool_xpndd%create_pw(density_xpndd)
882 CALL pw_pool_xpndd%create_pw(v_new_xpndd)
883 CALL pw_pool_xpndd%create_pw(v_eps_xpndd)
884
885 IF (use_zero_initial_guess) THEN
886 CALL pw_zero(v_old)
887 lambda0 = 0.0_dp
888 ELSE
889 CALL pw_copy(ps_implicit_env%initial_guess, v_old)
890 lambda0(:) = ps_implicit_env%initial_lambda
891 END IF
892
893 CALL pw_expand(neumann_directions, &
894 dct_env%recv_msgs_bnds, dct_env%dests_expand, dct_env%srcs_expand, &
895 dct_env%flipg_stat, dct_env%bounds_shftd, density, density_xpndd)
896 CALL pw_expand(neumann_directions, &
897 dct_env%recv_msgs_bnds, dct_env%dests_expand, dct_env%srcs_expand, &
898 dct_env%flipg_stat, dct_env%bounds_shftd, v_new, v_new_xpndd)
899
900 g%array = fourpi*density_xpndd%array/dielectric%eps%array
901 g_avg = accurate_sum(g%array)/ngpts
902
903 lambda_old(:) = lambda0
904
905! res_old = g - \Delta(v_old) - P(v_old) - B^t * \lambda_old
906 CALL apply_poisson_operator_dct(pw_pool_xpndd, green, dielectric, v_old, res_old)
907 CALL pw_scale(res_old, -1.0_dp)
908 CALL pw_axpy(g, res_old)
909 IF (data_size /= 0) THEN
910 CALL dgemv('N', data_size, n_tiles_tot, 1.0_dp, bt, data_size, lambda_old, 1, 0.0_dp, btxlambda_old, 1)
911 END IF
912 CALL convert_1dto3d(ps_implicit_env%idx_1dto3d, btxlambda_old, btxlambda_old3d)
913 res_old%array = res_old%array - btxlambda_old3d
914
915! evaluate \Delta^-1(res_old)
916 CALL apply_inv_laplace_operator_dct(pw_pool_xpndd, green, res_old, qainvxres)
917
918 iter = 1
919 DO
920
921! v_new (v_bar) = v_old + \omega * QAinvxres_old
922 CALL pw_scale(qainvxres, omega)
923 CALL pw_copy(qainvxres, v_new_xpndd)
924 CALL pw_axpy(v_old, v_new_xpndd)
925
926! evaluate 1^t * (g - \Delta(\bar{v}) - P(\bar{v}))
927! = 1^t * (g - P(\bar{v}))
928 CALL apply_p_operator(pw_pool_xpndd, dielectric, v_new_xpndd, axvbar)
929 axvbar_avg = accurate_sum(axvbar%array)/ngpts
930 gminusaxvbar_avg = g_avg - axvbar_avg
931 CALL dct_pw_grid%para%group%sum(gminusaxvbar_avg)
932
933! evaluate Q_S * \lambda + v_D - B * \bar{v}
934 CALL dgemv('N', n_tiles_tot, n_tiles_tot, 1.0_dp, qs, n_tiles_tot, lambda_old, 1, 0.0_dp, qsxlambda, 1)
935 v_bar1d(ps_implicit_env%idx_1dto3d) = reshape(v_new_xpndd%array, [data_size])
936 CALL dgemv('N', n_tiles_tot, data_size, 1.0_dp, b, n_tiles_tot, v_bar1d, 1, 0.0_dp, bxv_bar, 1)
937 CALL dct_pw_grid%para%group%sum(bxv_bar)
938! solve R [\lambda; \eta] = [Q_S * \lambda + v_D - B * \bar{v}; 1^t * (g - \Delta(\bar{v}) - P(\bar{v}))]
939 w = 0.0_dp
940 w(:) = [qsxlambda + v_d - bxv_bar, gminusaxvbar_avg]
941 CALL dgemv('N', n_tiles_tot + 1, n_tiles_tot + 1, 1.0_dp, rinv, n_tiles_tot + 1, w, 1, 0.0_dp, lambda_newneta, 1)
942 lambda_new(:) = lambda_newneta(1:n_tiles_tot)
943 eta = lambda_newneta(n_tiles_tot + 1)
944
945! v_new = v_bar + 1 * \eta
946 v_new_xpndd%array = v_new_xpndd%array + eta/ngpts
947
948! evaluate B^t * \lambda_new
949 IF (data_size /= 0) THEN
950 CALL dgemv('N', data_size, n_tiles_tot, 1.0_dp, bt, data_size, lambda_new, 1, 0.0_dp, btxlambda_new, 1)
951 END IF
952 CALL convert_1dto3d(ps_implicit_env%idx_1dto3d, btxlambda_new, btxlambda_new3d)
953
954! res_new = res_old - \omega * ( \Delta(QAinvxres_old) + P(QAinvxres_old) ) - B^t * ( \lambda_new - \lambda_old )
955! = (1 - \omega) * res_old - \omega * P(QAinvxres_old) - B^t * ( \lambda_new - \lambda_old )
956 CALL pw_zero(res_new)
957 CALL apply_p_operator(pw_pool_xpndd, dielectric, qainvxres, pxqainvxres)
958 CALL pw_axpy(pxqainvxres, res_new, -1.0_dp)
959 CALL pw_axpy(res_old, res_new, 1.0_dp - omega)
960 res_new%array = res_new%array - btxlambda_new3d + btxlambda_old3d
961
962! compute the error
963 CALL ps_implicit_compute_error_dct(pw_pool_xpndd, green, res_new, v_old, v_new_xpndd, qainvxres, &
964 pres_error, nabs_error)
965! output
966 CALL ps_implicit_output(iter, pres_error, nabs_error, outp_unit)
967 IF (PRESENT(electric_enthalpy)) THEN
968 CALL ps_implicit_compute_ehartree(dielectric, density_xpndd, btxlambda_new3d, v_new_xpndd, &
969 ehartree, electric_enthalpy)
970 CALL ps_implicit_report_ehartree(ps_implicit_env, outp_unit, ehartree/vol_scfac)
971 ps_implicit_env%ehartree = ehartree/vol_scfac
972 ps_implicit_env%electric_enthalpy = electric_enthalpy/vol_scfac
973 ELSE
974 IF (outp_unit > 0) WRITE (outp_unit, '(A1,/)')
975 END IF
976
977! verbose output
978 IF (poisson_env%parameters%dbc_params%verbose_output) THEN
979 v_new1d(ps_implicit_env%idx_1dto3d) = reshape(v_new_xpndd%array, [data_size])
980 CALL dgemv('N', n_tiles_tot, data_size, 1.0_dp, b, n_tiles_tot, v_new1d, 1, 0.0_dp, bxv_new, 1)
981 CALL pw_grid%para%group%sum(bxv_new)
982 IF (outp_unit > 0) THEN
983 WRITE (outp_unit, '(T3,A)') "======== verbose "//repeat('=', 61)
984 WRITE (outp_unit, '(T20,A)') "Drgn tile vhartree lambda "
985 WRITE (outp_unit, '(T19,A)') repeat('-', 46)
986 nt_tot = 1
987 DO ng = 1, n_contacts
988 DO nt = 1, ps_implicit_env%contacts(ng)%dirichlet_bc%n_tiles
989 WRITE (outp_unit, '(T17,I6,5X,I6,3X,E13.4,E13.4)') ng, nt, bxv_new(nt_tot), lambda_new(nt_tot)
990 nt_tot = nt_tot + 1
991 END DO
992 END DO
993 WRITE (outp_unit, '(T3,A)') repeat('=', 78)
994 END IF
995 END IF
996
997! check the convergence
998 iter = iter + 1
999 reached_max_iter = iter > max_iter
1000 reached_tol = pres_error <= tol
1001 ps_implicit_env%times_called = ps_implicit_env%times_called + 1
1002 IF (pres_error > large_error) THEN
1003 cpabort("Poisson solver did not converge.")
1004 END IF
1005 IF (reached_max_iter .OR. reached_tol) EXIT
1006
1007! update
1008 CALL pw_copy(v_new_xpndd, v_old)
1009 lambda_old(:) = lambda_new
1010 CALL pw_copy(res_new, res_old)
1011 btxlambda_old3d(:, :, :) = btxlambda_new3d
1012
1013 END DO
1014 CALL ps_implicit_print_convergence_msg(iter, max_iter, outp_unit)
1015
1016 CALL pw_shrink(neumann_directions, dct_env%dests_shrink, dct_env%srcs_shrink, &
1017 dct_env%bounds_local_shftd, v_new_xpndd, v_new)
1018
1019 IF ((times_called /= 0) .AND. (.NOT. use_zero_initial_guess)) THEN
1020 CALL pw_copy(v_new_xpndd, ps_implicit_env%initial_guess)
1021 ps_implicit_env%initial_lambda(:) = lambda_new
1022 END IF
1023
1024 ps_implicit_env%cstr_charge%array = btxlambda_new3d
1025 IF (PRESENT(electric_enthalpy)) electric_enthalpy = ps_implicit_env%electric_enthalpy
1026! compute the extra contribution to the Hamiltonian due to the presence of dielectric
1027 CALL ps_implicit_compute_veps(pw_pool_xpndd, dielectric, v_new_xpndd, v_eps_xpndd)
1028 CALL pw_shrink(neumann_directions, dct_env%dests_shrink, dct_env%srcs_shrink, &
1029 dct_env%bounds_local_shftd, v_eps_xpndd, ps_implicit_env%v_eps)
1030
1031 CALL pw_pool_xpndd%give_back_pw(g)
1032 CALL pw_pool_xpndd%give_back_pw(v_old)
1033 CALL pw_pool_xpndd%give_back_pw(res_old)
1034 CALL pw_pool_xpndd%give_back_pw(res_new)
1035 CALL pw_pool_xpndd%give_back_pw(qainvxres)
1036 CALL pw_pool_xpndd%give_back_pw(pxqainvxres)
1037 CALL pw_pool_xpndd%give_back_pw(axvbar)
1038 CALL pw_pool_xpndd%give_back_pw(density_xpndd)
1039 CALL pw_pool_xpndd%give_back_pw(v_new_xpndd)
1040 CALL pw_pool_xpndd%give_back_pw(v_eps_xpndd)
1041 CALL pw_pool_release(pw_pool_xpndd)
1042
1043 CALL timestop(handle)
1044
1045 END SUBROUTINE implicit_poisson_solver_mixed
1046
1047! **************************************************************************************************
1048!> \brief allocates and zeroises initial guess for implicit (iterative) Poisson solver
1049!> \param ps_implicit_env the implicit env containing the initial guess
1050!> \param pw_pool pool of pw grid
1051!> \par History
1052!> 06.2014 created [Hossein Bani-Hashemian]
1053!> \author Mohammad Hossein Bani-Hashemian
1054! **************************************************************************************************
1055 SUBROUTINE ps_implicit_initial_guess_create(ps_implicit_env, pw_pool)
1056
1057 TYPE(ps_implicit_type), INTENT(INOUT), POINTER :: ps_implicit_env
1058 TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
1059
1060 CHARACTER(LEN=*), PARAMETER :: routinen = 'ps_implicit_initial_guess_create'
1061
1062 INTEGER :: handle, n_tiles_tot
1063
1064 CALL timeset(routinen, handle)
1065
1066 n_tiles_tot = SIZE(ps_implicit_env%v_D)
1067 NULLIFY (ps_implicit_env%initial_guess)
1068 ALLOCATE (ps_implicit_env%initial_guess)
1069 CALL pw_pool%create_pw(ps_implicit_env%initial_guess)
1070 CALL pw_zero(ps_implicit_env%initial_guess)
1071 ALLOCATE (ps_implicit_env%initial_lambda(n_tiles_tot))
1072 ps_implicit_env%initial_lambda = 0.0_dp
1073
1074 CALL timestop(handle)
1075
1076 END SUBROUTINE ps_implicit_initial_guess_create
1077
1078! **************************************************************************************************
1079!> \brief prepare blocks B, Bt, QS, R^-1, v_D
1080!> \param pw_pool_orig original pw grid
1081!> \param dct_pw_grid DCT (extended) grid
1082!> \param green green functions for FFT based inverse Laplacian
1083!> \param poisson_params paramaters of the poisson_env
1084!> \param ps_implicit_env the implicit_env that stores the blocks
1085!> \par History
1086!> 10.2014 created [Hossein Bani-Hashemian]
1087!> \author Mohammad Hossein Bani-Hashemian
1088! **************************************************************************************************
1089 SUBROUTINE ps_implicit_prepare_blocks(pw_pool_orig, dct_pw_grid, green, &
1090 poisson_params, ps_implicit_env)
1091
1092 TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool_orig
1093 TYPE(pw_grid_type), INTENT(IN), POINTER :: dct_pw_grid
1094 TYPE(greens_fn_type), INTENT(IN) :: green
1095 TYPE(pw_poisson_parameter_type), INTENT(IN) :: poisson_params
1096 TYPE(ps_implicit_type), INTENT(INOUT), POINTER :: ps_implicit_env
1097
1098 CHARACTER(LEN=*), PARAMETER :: routinen = 'ps_implicit_prepare_blocks'
1099
1100 INTEGER :: data_size, handle, i, indx1, indx2, info, j, k, l, lb1, lb2, lb3, n_contacts, &
1101 n_tiles, n_tiles_tot, neumann_directions, ngpts_local, ub1, ub2, ub3, unit_nr
1102 INTEGER(KIND=int_8) :: ngpts
1103 INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv
1104 INTEGER, DIMENSION(2, 3) :: bounds, bounds_local
1105 INTEGER, DIMENSION(3) :: npts, npts_local
1106 LOGICAL :: done_preparing
1107 REAL(dp) :: tile_volume, vol_scfac
1108 REAL(dp), ALLOCATABLE, DIMENSION(:) :: bxunit_vec, test_vec, work_arr
1109 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: qainvxbt, r
1110 TYPE(cp_logger_type), POINTER :: logger
1111 TYPE(dct_type), POINTER :: dct_env
1112 TYPE(pw_grid_type), POINTER :: pw_grid_orig
1113 TYPE(pw_pool_type), POINTER :: pw_pool_xpndd
1114 TYPE(pw_r3d_rs_type) :: pw_in, pw_out
1115
1116 CALL timeset(routinen, handle)
1117
1118 pw_grid_orig => pw_pool_orig%pw_grid
1119
1120 logger => cp_get_default_logger()
1121 IF (logger%para_env%is_source()) THEN
1122 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
1123 ELSE
1124 unit_nr = -1
1125 END IF
1126
1127 SELECT CASE (poisson_params%ps_implicit_params%boundary_condition)
1128 CASE (mixed_bc)
1129
1130 ngpts_local = dct_pw_grid%ngpts_local
1131 ngpts = dct_pw_grid%ngpts
1132 npts_local = dct_pw_grid%npts_local
1133 npts = dct_pw_grid%npts
1134 bounds_local = dct_pw_grid%bounds_local
1135 bounds = dct_pw_grid%bounds
1136 dct_env => ps_implicit_env%dct_env
1137
1138 neumann_directions = poisson_params%ps_implicit_params%neumann_directions
1139
1140 SELECT CASE (neumann_directions)
1141 CASE (neumannxyz)
1142 vol_scfac = 8.0_dp
1143 CASE (neumannxy, neumannxz, neumannyz)
1144 vol_scfac = 4.0_dp
1145 CASE (neumannx, neumanny, neumannz)
1146 vol_scfac = 2.0_dp
1147 END SELECT
1148
1149! evaluate indices for converting 3D arrays into 1D arrays
1150 lb1 = bounds_local(1, 1); ub1 = bounds_local(2, 1)
1151 lb2 = bounds_local(1, 2); ub2 = bounds_local(2, 2)
1152 lb3 = bounds_local(1, 3); ub3 = bounds_local(2, 3)
1153
1154 IF (dct_pw_grid%para%blocked) THEN
1155 data_size = product(npts_local)
1156 ELSE IF (dct_pw_grid%para%ray_distribution) THEN
1157 data_size = ngpts_local
1158 ELSE ! parallel run with np = 1
1159 data_size = product(npts_local)
1160 END IF
1161
1162 ALLOCATE (ps_implicit_env%idx_1dto3d(data_size))
1163 l = 1
1164 ! Suppress OpenMP (at least the Intel compiler has an issue here)
1165 ! An automatic OpenMP parallelization of this loop might be tricky
1166 ! because of the l incrementation
1167!$OMP PARALLEL IF(.FALSE.)
1168!$OMP DO
1169 DO k = lb3, ub3
1170 DO j = lb2, ub2
1171 DO i = lb1, ub1
1172 ps_implicit_env%idx_1dto3d(l) = (i - lb1 + 1) + &
1173 (j - lb2)*npts_local(1) + &
1174 (k - lb3)*npts_local(1)*npts_local(2)
1175 l = l + 1
1176 END DO
1177 END DO
1178 END DO
1179!$OMP END DO
1180!$OMP END PARALLEL
1181
1182 n_contacts = SIZE(ps_implicit_env%contacts)
1183 n_tiles_tot = 0
1184 DO j = 1, n_contacts
1185 n_tiles_tot = n_tiles_tot + ps_implicit_env%contacts(j)%dirichlet_bc%n_tiles
1186 END DO
1187
1188 ALLOCATE (ps_implicit_env%B(n_tiles_tot, data_size))
1189 ALLOCATE (ps_implicit_env%Bt(data_size, n_tiles_tot))
1190 ALLOCATE (ps_implicit_env%QS(n_tiles_tot, n_tiles_tot))
1191 ALLOCATE (ps_implicit_env%Rinv(n_tiles_tot + 1, n_tiles_tot + 1))
1192 ALLOCATE (ps_implicit_env%v_D(n_tiles_tot))
1193 ALLOCATE (ps_implicit_env%osc_frac(n_tiles_tot))
1194 ALLOCATE (ps_implicit_env%frequency(n_tiles_tot))
1195 ALLOCATE (ps_implicit_env%phase(n_tiles_tot))
1196
1197 ALLOCATE (qainvxbt(data_size, n_tiles_tot))
1198 ALLOCATE (bxunit_vec(n_tiles_tot))
1199 ALLOCATE (test_vec(n_tiles_tot))
1200 ALLOCATE (r(n_tiles_tot + 1, n_tiles_tot + 1))
1201 ALLOCATE (work_arr(n_tiles_tot + 1), ipiv(n_tiles_tot + 1)) ! LAPACK work and ipiv arrays
1202
1203! prepare pw_pool for evaluating inverse Laplacian of tile_pw's using DCT
1204 CALL pw_pool_create(pw_pool_xpndd, pw_grid=dct_pw_grid)
1205
1206! set up B, B^t, (\Delta^-1)*B^t
1207 indx1 = 1
1208 DO j = 1, n_contacts
1209 n_tiles = ps_implicit_env%contacts(j)%dirichlet_bc%n_tiles
1210 indx2 = indx1 + n_tiles - 1
1211 DO i = 1, n_tiles
1212
1213 CALL pw_pool_xpndd%create_pw(pw_in)
1214 CALL pw_expand(neumann_directions, &
1215 dct_env%recv_msgs_bnds, dct_env%dests_expand, dct_env%srcs_expand, &
1216 dct_env%flipg_stat, dct_env%bounds_shftd, &
1217 ps_implicit_env%contacts(j)%dirichlet_bc%tiles(i)%tile%tile_pw, pw_in)
1218
1219 tile_volume = ps_implicit_env%contacts(j)%dirichlet_bc%tiles(i)%tile%volume
1220 CALL pw_scale(pw_in, 1.0_dp/(vol_scfac*tile_volume)) ! normalize tile_pw
1221 ps_implicit_env%Bt(ps_implicit_env%idx_1dto3d, indx1 + i - 1) = reshape(pw_in%array, [data_size])
1222
1223 CALL pw_pool_xpndd%create_pw(pw_out)
1224 CALL apply_inv_laplace_operator_dct(pw_pool_xpndd, green, pw_in, pw_out)
1225 qainvxbt(ps_implicit_env%idx_1dto3d, indx1 + i - 1) = reshape(pw_out%array, [data_size])
1226 ! the electrostatic potential has opposite sign by internal convention
1227 ps_implicit_env%v_D(indx1 + i - 1) = -1.0_dp*ps_implicit_env%contacts(j)%dirichlet_bc%v_D
1228 ps_implicit_env%osc_frac(indx1 + i - 1) = ps_implicit_env%contacts(j)%dirichlet_bc%osc_frac
1229 ps_implicit_env%frequency(indx1 + i - 1) = ps_implicit_env%contacts(j)%dirichlet_bc%frequency
1230 ps_implicit_env%phase(indx1 + i - 1) = ps_implicit_env%contacts(j)%dirichlet_bc%phase
1231
1232 CALL pw_pool_xpndd%give_back_pw(pw_in)
1233 CALL pw_pool_xpndd%give_back_pw(pw_out)
1234 END DO
1235 indx1 = indx2 + 1
1236 END DO
1237 ps_implicit_env%B(:, :) = transpose(ps_implicit_env%Bt)
1238
1239! evaluate QS = - B*(\Delta^-1)*B^t
1240 IF (data_size /= 0) THEN
1241 CALL dgemm('N', 'N', n_tiles_tot, n_tiles_tot, data_size, &
1242 -1.0_dp, ps_implicit_env%B, n_tiles_tot, qainvxbt, &
1243 data_size, 0.0_dp, ps_implicit_env%QS, n_tiles_tot)
1244 END IF
1245 CALL pw_grid_orig%para%group%sum(ps_implicit_env%QS)
1246
1247! evaluate B*1
1248 bxunit_vec(:) = sum(ps_implicit_env%B, 2)/ngpts
1249 CALL pw_grid_orig%para%group%sum(bxunit_vec)
1250! set up R = [QS B*1; (B*1)^t 0]
1251 r = 0.0_dp
1252 r(1:n_tiles_tot, 1:n_tiles_tot) = ps_implicit_env%QS
1253 r(1:n_tiles_tot, n_tiles_tot + 1) = bxunit_vec
1254 r(n_tiles_tot + 1, 1:n_tiles_tot) = bxunit_vec
1255! evaluate R^(-1)
1256 ps_implicit_env%Rinv(:, :) = r
1257 CALL dgetrf(n_tiles_tot + 1, n_tiles_tot + 1, ps_implicit_env%Rinv, n_tiles_tot + 1, ipiv, info)
1258 IF (info /= 0) THEN
1259 CALL cp_abort(__location__, &
1260 "R is (nearly) singular! Either two Dirichlet constraints are identical or "// &
1261 "you need to reduce the number of tiles.")
1262 END IF
1263 CALL dgetri(n_tiles_tot + 1, ps_implicit_env%Rinv, n_tiles_tot + 1, ipiv, work_arr, n_tiles_tot + 1, info)
1264 IF (info /= 0) THEN
1265 cpabort("Inversion of R failed!")
1266 END IF
1267
1268 DEALLOCATE (qainvxbt, bxunit_vec, r, work_arr, ipiv)
1269 CALL pw_pool_release(pw_pool_xpndd)
1270
1271 done_preparing = .true.
1272 CALL pw_grid_orig%para%group%sum(done_preparing)
1273 IF ((unit_nr > 0) .AND. done_preparing) THEN
1274 WRITE (unit_nr, "(T3,A,/,T3,A,/,A)") "POISSON| ... Done. ", repeat('-', 78)
1275 END IF
1276
1277 CASE (mixed_periodic_bc)
1278
1279 ngpts_local = pw_grid_orig%ngpts_local
1280 ngpts = pw_grid_orig%ngpts
1281 npts_local = pw_grid_orig%npts_local
1282 npts = pw_grid_orig%npts
1283 bounds_local = pw_grid_orig%bounds_local
1284 bounds = pw_grid_orig%bounds
1285 dct_env => ps_implicit_env%dct_env
1286
1287! evaluate indices for converting 3D arrays into 1D arrays
1288 lb1 = bounds_local(1, 1); ub1 = bounds_local(2, 1)
1289 lb2 = bounds_local(1, 2); ub2 = bounds_local(2, 2)
1290 lb3 = bounds_local(1, 3); ub3 = bounds_local(2, 3)
1291
1292 IF (pw_grid_orig%para%blocked) THEN
1293 data_size = product(npts_local)
1294 ELSE IF (pw_grid_orig%para%ray_distribution) THEN
1295 data_size = ngpts_local
1296 ELSE ! parallel run with np = 1
1297 data_size = product(npts_local)
1298 END IF
1299
1300 ALLOCATE (ps_implicit_env%idx_1dto3d(data_size))
1301 l = 1
1302 ! Suppress OpenMP (at least the Intel compiler has an issue here)
1303 ! An automatic OpenMP parallelization of this loop might be tricky
1304 ! because of the l incrementation
1305!$OMP PARALLEL IF(.FALSE.)
1306!$OMP DO
1307 DO k = lb3, ub3
1308 DO j = lb2, ub2
1309 DO i = lb1, ub1
1310 ps_implicit_env%idx_1dto3d(l) = (i - lb1 + 1) + &
1311 (j - lb2)*npts_local(1) + &
1312 (k - lb3)*npts_local(1)*npts_local(2)
1313 l = l + 1
1314 END DO
1315 END DO
1316 END DO
1317!$OMP END DO
1318!$OMP END PARALLEL
1319
1320 n_contacts = SIZE(ps_implicit_env%contacts)
1321 n_tiles_tot = 0
1322 DO j = 1, n_contacts
1323 n_tiles_tot = n_tiles_tot + ps_implicit_env%contacts(j)%dirichlet_bc%n_tiles
1324 END DO
1325
1326 ALLOCATE (ps_implicit_env%B(n_tiles_tot, data_size))
1327 ALLOCATE (ps_implicit_env%Bt(data_size, n_tiles_tot))
1328 ALLOCATE (ps_implicit_env%QS(n_tiles_tot, n_tiles_tot))
1329 ALLOCATE (ps_implicit_env%Rinv(n_tiles_tot + 1, n_tiles_tot + 1))
1330 ALLOCATE (ps_implicit_env%v_D(n_tiles_tot))
1331 ALLOCATE (ps_implicit_env%osc_frac(n_tiles_tot))
1332 ALLOCATE (ps_implicit_env%frequency(n_tiles_tot))
1333 ALLOCATE (ps_implicit_env%phase(n_tiles_tot))
1334
1335 ALLOCATE (qainvxbt(data_size, n_tiles_tot))
1336 ALLOCATE (bxunit_vec(n_tiles_tot))
1337 ALLOCATE (test_vec(n_tiles_tot))
1338 ALLOCATE (r(n_tiles_tot + 1, n_tiles_tot + 1))
1339 ALLOCATE (work_arr(n_tiles_tot + 1), ipiv(n_tiles_tot + 1))
1340
1341! set up B, B^t, (\Delta^-1)*B^t
1342 indx1 = 1
1343 DO j = 1, n_contacts
1344 n_tiles = ps_implicit_env%contacts(j)%dirichlet_bc%n_tiles
1345 indx2 = indx1 + n_tiles - 1
1346 DO i = 1, n_tiles
1347 CALL pw_pool_orig%create_pw(pw_in)
1348 CALL pw_copy(ps_implicit_env%contacts(j)%dirichlet_bc%tiles(i)%tile%tile_pw, pw_in)
1349
1350 tile_volume = ps_implicit_env%contacts(j)%dirichlet_bc%tiles(i)%tile%volume
1351 CALL pw_scale(pw_in, 1.0_dp/tile_volume) ! normalize tile_pw
1352 ps_implicit_env%Bt(ps_implicit_env%idx_1dto3d, indx1 + i - 1) = reshape(pw_in%array, [data_size])
1353
1354 CALL pw_pool_orig%create_pw(pw_out)
1355 CALL apply_inv_laplace_operator_fft(pw_pool_orig, green, pw_in, pw_out)
1356 qainvxbt(ps_implicit_env%idx_1dto3d, indx1 + i - 1) = reshape(pw_out%array, [data_size])
1357 ! the electrostatic potential has opposite sign by internal convention
1358 ps_implicit_env%v_D(indx1 + i - 1) = -1.0_dp*ps_implicit_env%contacts(j)%dirichlet_bc%v_D
1359 ps_implicit_env%osc_frac(indx1 + i - 1) = ps_implicit_env%contacts(j)%dirichlet_bc%osc_frac
1360 ps_implicit_env%frequency(indx1 + i - 1) = ps_implicit_env%contacts(j)%dirichlet_bc%frequency
1361 ps_implicit_env%phase(indx1 + i - 1) = ps_implicit_env%contacts(j)%dirichlet_bc%phase
1362
1363 CALL pw_pool_orig%give_back_pw(pw_in)
1364 CALL pw_pool_orig%give_back_pw(pw_out)
1365 END DO
1366 indx1 = indx2 + 1
1367 END DO
1368 ps_implicit_env%B(:, :) = transpose(ps_implicit_env%Bt)
1369
1370! evaluate QS = - B*(\Delta^-1)*B^t
1371 IF (data_size /= 0) THEN
1372 CALL dgemm('N', 'N', n_tiles_tot, n_tiles_tot, data_size, &
1373 -1.0_dp, ps_implicit_env%B, n_tiles_tot, qainvxbt, &
1374 data_size, 0.0_dp, ps_implicit_env%QS, n_tiles_tot)
1375 END IF
1376 CALL pw_grid_orig%para%group%sum(ps_implicit_env%QS)
1377
1378! evaluate B*1
1379 bxunit_vec(:) = sum(ps_implicit_env%B, 2)/ngpts
1380 CALL pw_grid_orig%para%group%sum(bxunit_vec)
1381! set up R = [QS B*1; (B*1)^t 0]
1382 r = 0.0_dp
1383 r(1:n_tiles_tot, 1:n_tiles_tot) = ps_implicit_env%QS
1384 r(1:n_tiles_tot, n_tiles_tot + 1) = bxunit_vec
1385 r(n_tiles_tot + 1, 1:n_tiles_tot) = bxunit_vec
1386! evaluate R^(-1)
1387 ps_implicit_env%Rinv(:, :) = r
1388 CALL dgetrf(n_tiles_tot + 1, n_tiles_tot + 1, ps_implicit_env%Rinv, n_tiles_tot + 1, ipiv, info)
1389 IF (info /= 0) THEN
1390 CALL cp_abort(__location__, &
1391 "R is (nearly) singular! Either two Dirichlet constraints are identical or "// &
1392 "you need to reduce the number of tiles.")
1393 END IF
1394 CALL dgetri(n_tiles_tot + 1, ps_implicit_env%Rinv, n_tiles_tot + 1, ipiv, work_arr, n_tiles_tot + 1, info)
1395 IF (info /= 0) THEN
1396 cpabort("Inversion of R failed!")
1397 END IF
1398
1399 DEALLOCATE (qainvxbt, bxunit_vec, r, work_arr, ipiv)
1400
1401 done_preparing = .true.
1402 CALL pw_grid_orig%para%group%sum(done_preparing)
1403 IF ((unit_nr > 0) .AND. done_preparing) THEN
1404 WRITE (unit_nr, "(T3,A,/,T3,A,/,A)") "POISSON| ... Done. ", repeat('-', 78)
1405 END IF
1406
1407 CASE (periodic_bc, neumann_bc)
1408
1409 ALLOCATE (ps_implicit_env%idx_1dto3d(1))
1410 ALLOCATE (ps_implicit_env%B(1, 1))
1411 ALLOCATE (ps_implicit_env%Bt(1, 1))
1412 ALLOCATE (ps_implicit_env%QS(1, 1))
1413 ALLOCATE (ps_implicit_env%Rinv(1, 1))
1414 ALLOCATE (ps_implicit_env%v_D(1))
1415 ALLOCATE (ps_implicit_env%osc_frac(1))
1416 ALLOCATE (ps_implicit_env%frequency(1))
1417 ALLOCATE (ps_implicit_env%phase(1))
1418
1419 ps_implicit_env%idx_1dto3d = 1
1420 ps_implicit_env%B = 0.0_dp
1421 ps_implicit_env%Bt = 0.0_dp
1422 ps_implicit_env%QS = 0.0_dp
1423 ps_implicit_env%Rinv = 0.0_dp
1424 ps_implicit_env%v_D = 0.0_dp
1425
1426 CASE DEFAULT
1427 CALL cp_abort(__location__, &
1428 "Please specify the type of boundary conditions using the "// &
1429 "input file keyword BOUNDARY_CONDITIONS.")
1430 END SELECT
1431
1432 CALL timestop(handle)
1433
1434 END SUBROUTINE ps_implicit_prepare_blocks
1435
1436! **************************************************************************************************
1437!> \brief Evaluates the action of the operator P on a given matrix v, defined
1438!> as: P(v) := - \nabla_r(\ln(\eps)) \cdot \nabla_r(v)
1439!> \param pw_pool pool of pw grid
1440!> \param dielectric dielectric_type containing eps
1441!> \param v input matrix
1442!> \param Pxv action of the operator P on v
1443!> \par History
1444!> 07.2014 created [Hossein Bani-Hashemian]
1445!> \author Mohammad Hossein Bani-Hashemian
1446! **************************************************************************************************
1447 SUBROUTINE apply_p_operator(pw_pool, dielectric, v, Pxv)
1448
1449 TYPE(pw_pool_type), POINTER :: pw_pool
1450 TYPE(dielectric_type), INTENT(IN), POINTER :: dielectric
1451 TYPE(pw_r3d_rs_type), INTENT(IN) :: v
1452 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pxv
1453
1454 CHARACTER(LEN=*), PARAMETER :: routinen = 'apply_P_operator'
1455
1456 INTEGER :: handle, i
1457 TYPE(pw_r3d_rs_type), DIMENSION(3) :: dv
1458
1459 CALL timeset(routinen, handle)
1460
1461 DO i = 1, 3
1462 CALL pw_pool%create_pw(dv(i))
1463 END DO
1464
1465 CALL derive_fft(v, dv, pw_pool)
1466 associate(dln_eps => dielectric%dln_eps)
1467 pxv%array = -(dv(1)%array*dln_eps(1)%array + &
1468 dv(2)%array*dln_eps(2)%array + &
1469 dv(3)%array*dln_eps(3)%array)
1470 END associate
1471
1472 DO i = 1, 3
1473 CALL pw_pool%give_back_pw(dv(i))
1474 END DO
1475
1476 CALL timestop(handle)
1477
1478 END SUBROUTINE apply_p_operator
1479
1480! **************************************************************************************************
1481!> \brief Evaluates the action of the inverse of the Laplace operator on a given 3d matrix
1482!> \param pw_pool pool of pw grid
1483!> \param green green functions for FFT based inverse Laplacian
1484!> \param pw_in pw_in (density)
1485!> \param pw_out pw_out (potential)
1486!> \par History
1487!> 07.2014 created [Hossein Bani-Hashemian]
1488!> \author Mohammad Hossein Bani-Hashemian
1489! **************************************************************************************************
1490 SUBROUTINE apply_inv_laplace_operator_fft(pw_pool, green, pw_in, pw_out)
1491
1492 TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
1493 TYPE(greens_fn_type), INTENT(IN) :: green
1494 TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in
1495 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pw_out
1496
1497 CHARACTER(LEN=*), PARAMETER :: routinen = 'apply_inv_laplace_operator_fft'
1498
1499 INTEGER :: handle, ig, ng
1500 REAL(dp) :: prefactor
1501 TYPE(pw_c1d_gs_type) :: pw_in_gs
1502 TYPE(pw_grid_type), POINTER :: pw_grid
1503
1504 CALL timeset(routinen, handle)
1505
1506! here I divide by fourpi to cancel out the prefactor fourpi in influence_fn
1507 prefactor = 1.0_dp/fourpi
1508
1509 pw_grid => pw_pool%pw_grid
1510 ng = SIZE(pw_grid%gsq)
1511
1512 CALL pw_pool%create_pw(pw_in_gs)
1513
1514 CALL pw_transfer(pw_in, pw_in_gs)
1515 DO ig = 1, ng
1516 pw_in_gs%array(ig) = prefactor*pw_in_gs%array(ig)*green%influence_fn%array(ig)
1517 END DO
1518 CALL pw_transfer(pw_in_gs, pw_out)
1519
1520 CALL pw_pool%give_back_pw(pw_in_gs)
1521
1522 CALL timestop(handle)
1523
1524 END SUBROUTINE apply_inv_laplace_operator_fft
1525
1526! **************************************************************************************************
1527!> \brief Evaluates the action of the inverse of the Laplace operator on a given
1528!> 3d matrix using DCT-I
1529!> \param pw_pool pool of pw grid
1530!> \param green the greens_fn_type data holding a valid dct_influence_fn
1531!> \param pw_in pw_in (density)
1532!> \param pw_out pw_out (potential)
1533!> \par History
1534!> 07.2014 created [Hossein Bani-Hashemian]
1535!> 11.2015 revised [Hossein Bani-Hashemian]
1536!> \author Mohammad Hossein Bani-Hashemian
1537! **************************************************************************************************
1538 SUBROUTINE apply_inv_laplace_operator_dct(pw_pool, green, pw_in, pw_out)
1539
1540 TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
1541 TYPE(greens_fn_type), INTENT(IN) :: green
1542 TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in
1543 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pw_out
1544
1545 CHARACTER(LEN=*), PARAMETER :: routinen = 'apply_inv_laplace_operator_dct'
1546
1547 INTEGER :: handle, ig, ng
1548 REAL(dp) :: prefactor
1549 TYPE(pw_c1d_gs_type) :: pw_in_gs
1550 TYPE(pw_grid_type), POINTER :: pw_grid
1551
1552 CALL timeset(routinen, handle)
1553
1554! here I divide by fourpi to cancel out the prefactor fourpi in influence_fn
1555 prefactor = 1.0_dp/fourpi
1556
1557 pw_grid => pw_pool%pw_grid
1558 ng = SIZE(pw_grid%gsq)
1559
1560 CALL pw_pool%create_pw(pw_in_gs)
1561
1562 CALL pw_transfer(pw_in, pw_in_gs)
1563 DO ig = 1, ng
1564 pw_in_gs%array(ig) = prefactor*pw_in_gs%array(ig)*green%dct_influence_fn%array(ig)
1565 END DO
1566 CALL pw_transfer(pw_in_gs, pw_out)
1567
1568 CALL pw_pool%give_back_pw(pw_in_gs)
1569
1570 CALL timestop(handle)
1571
1572 END SUBROUTINE apply_inv_laplace_operator_dct
1573
1574! **************************************************************************************************
1575!> \brief Evaluates the action of the Laplace operator on a given 3d matrix
1576!> \param pw_pool pool of pw grid
1577!> \param green green functions for FFT based inverse Laplacian
1578!> \param pw_in pw_in (potential)
1579!> \param pw_out pw_out (density)
1580!> \par History
1581!> 07.2014 created [Hossein Bani-Hashemian]
1582!> \author Mohammad Hossein Bani-Hashemian
1583! **************************************************************************************************
1584 SUBROUTINE apply_laplace_operator_fft(pw_pool, green, pw_in, pw_out)
1585
1586 TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
1587 TYPE(greens_fn_type), INTENT(IN) :: green
1588 TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in
1589 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pw_out
1590
1591 CHARACTER(LEN=*), PARAMETER :: routinen = 'apply_laplace_operator_fft'
1592
1593 INTEGER :: g0_index, handle, ig, ng
1594 LOGICAL :: have_g0
1595 REAL(dp) :: prefactor
1596 TYPE(pw_c1d_gs_type) :: pw_in_gs
1597 TYPE(pw_grid_type), POINTER :: pw_grid
1598
1599 CALL timeset(routinen, handle)
1600
1601! here I multiply by fourpi to cancel out the prefactor fourpi in influence_fn
1602 prefactor = fourpi
1603
1604 pw_grid => pw_pool%pw_grid
1605 ng = SIZE(pw_in%pw_grid%gsq)
1606 have_g0 = green%influence_fn%pw_grid%have_g0
1607
1608 CALL pw_pool%create_pw(pw_in_gs)
1609
1610 CALL pw_transfer(pw_in, pw_in_gs)
1611
1612 IF (have_g0) THEN
1613 g0_index = green%influence_fn%pw_grid%first_gne0 - 1
1614 pw_in_gs%array(g0_index) = 0.0_dp
1615 END IF
1616 DO ig = green%influence_fn%pw_grid%first_gne0, ng
1617 pw_in_gs%array(ig) = prefactor*(pw_in_gs%array(ig)/green%influence_fn%array(ig))
1618 END DO
1619
1620 CALL pw_transfer(pw_in_gs, pw_out)
1621
1622 CALL pw_pool%give_back_pw(pw_in_gs)
1623
1624 CALL timestop(handle)
1625
1626 END SUBROUTINE apply_laplace_operator_fft
1627
1628! **************************************************************************************************
1629!> \brief Evaluates the action of the Laplace operator on a given 3d matrix using DCT-I
1630!> \param pw_pool pool of pw grid
1631!> \param green the greens_fn_type data holding a valid dct_influence_fn
1632!> \param pw_in pw_in (potential)
1633!> \param pw_out pw_out (density)
1634!> \par History
1635!> 07.2014 created [Hossein Bani-Hashemian]
1636!> 11.2015 revised [Hossein Bani-Hashemian]
1637!> \author Mohammad Hossein Bani-Hashemian
1638! **************************************************************************************************
1639 SUBROUTINE apply_laplace_operator_dct(pw_pool, green, pw_in, pw_out)
1640
1641 TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
1642 TYPE(greens_fn_type), INTENT(IN) :: green
1643 TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in
1644 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pw_out
1645
1646 CHARACTER(LEN=*), PARAMETER :: routinen = 'apply_laplace_operator_dct'
1647
1648 INTEGER :: g0_index, handle, ig, ng
1649 LOGICAL :: have_g0
1650 REAL(dp) :: prefactor
1651 TYPE(pw_c1d_gs_type) :: pw_in_gs
1652 TYPE(pw_grid_type), POINTER :: pw_grid
1653
1654 CALL timeset(routinen, handle)
1655
1656! here I multiply by fourpi to cancel out the prefactor fourpi in influence_fn
1657 prefactor = fourpi
1658
1659 pw_grid => pw_pool%pw_grid
1660 ng = SIZE(pw_in%pw_grid%gsq)
1661 have_g0 = green%dct_influence_fn%pw_grid%have_g0
1662
1663 CALL pw_pool%create_pw(pw_in_gs)
1664
1665 CALL pw_transfer(pw_in, pw_in_gs)
1666
1667 IF (have_g0) THEN
1668 g0_index = green%dct_influence_fn%pw_grid%first_gne0 - 1
1669 pw_in_gs%array(g0_index) = 0.0_dp
1670 END IF
1671 DO ig = green%dct_influence_fn%pw_grid%first_gne0, ng
1672 pw_in_gs%array(ig) = prefactor*(pw_in_gs%array(ig)/green%dct_influence_fn%array(ig))
1673 END DO
1674
1675 CALL pw_transfer(pw_in_gs, pw_out)
1676
1677 CALL pw_pool%give_back_pw(pw_in_gs)
1678
1679 CALL timestop(handle)
1680
1681 END SUBROUTINE apply_laplace_operator_dct
1682
1683! **************************************************************************************************
1684!> \brief Evaluates the action of the generalized Poisson operator on a given 3d matrix.
1685!> \param pw_pool pool of pw grid
1686!> \param green green functions for FFT based inverse Laplacian
1687!> \param dielectric dielectric environment
1688!> \param v potential
1689!> \param density density
1690!> \par History
1691!> 07.2014 created [Hossein Bani-Hashemian]
1692!> \author Mohammad Hossein Bani-Hashemian
1693! **************************************************************************************************
1694 SUBROUTINE apply_poisson_operator_fft(pw_pool, green, dielectric, v, density)
1695
1696 TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
1697 TYPE(greens_fn_type), INTENT(IN) :: green
1698 TYPE(dielectric_type), INTENT(IN), POINTER :: dielectric
1699 TYPE(pw_r3d_rs_type), INTENT(IN) :: v
1700 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: density
1701
1702 CHARACTER(LEN=*), PARAMETER :: routinen = 'apply_poisson_operator_fft'
1703
1704 INTEGER :: handle
1705 TYPE(pw_r3d_rs_type) :: pxv
1706
1707 CALL timeset(routinen, handle)
1708
1709 CALL pw_pool%create_pw(pxv)
1710
1711 CALL apply_p_operator(pw_pool, dielectric, v, pxv)
1712 CALL apply_laplace_operator_fft(pw_pool, green, v, density)
1713 CALL pw_axpy(pxv, density)
1714
1715 CALL pw_pool%give_back_pw(pxv)
1716
1717 CALL timestop(handle)
1718
1719 END SUBROUTINE apply_poisson_operator_fft
1720
1721! **************************************************************************************************
1722!> \brief Evaluates the action of the generalized Poisson operator on a given
1723!> 3d matrix using DCT-I.
1724!> \param pw_pool pool of pw grid
1725!> \param green the greens_fn_type data holding a valid dct_influence_fn
1726!> \param dielectric dielectric environment
1727!> \param v potential
1728!> \param density density
1729!> \par History
1730!> 07.2014 created [Hossein Bani-Hashemian]
1731!> 11.2015 revised [Hossein Bani-Hashemian]
1732!> \author Mohammad Hossein Bani-Hashemian
1733! **************************************************************************************************
1734 SUBROUTINE apply_poisson_operator_dct(pw_pool, green, dielectric, v, density)
1735
1736 TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
1737 TYPE(greens_fn_type), INTENT(IN) :: green
1738 TYPE(dielectric_type), INTENT(IN), POINTER :: dielectric
1739 TYPE(pw_r3d_rs_type), INTENT(IN) :: v
1740 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: density
1741
1742 CHARACTER(LEN=*), PARAMETER :: routinen = 'apply_poisson_operator_dct'
1743
1744 INTEGER :: handle
1745 TYPE(pw_r3d_rs_type) :: pxv
1746
1747 CALL timeset(routinen, handle)
1748
1749 CALL pw_pool%create_pw(pxv)
1750
1751 CALL apply_p_operator(pw_pool, dielectric, v, pxv)
1752 CALL apply_laplace_operator_dct(pw_pool, green, v, density)
1753 CALL pw_axpy(pxv, density)
1754
1755 CALL pw_pool%give_back_pw(pxv)
1756
1757 CALL timestop(handle)
1758
1759 END SUBROUTINE apply_poisson_operator_dct
1760
1761! **************************************************************************************************
1762!> \brief Computes the extra contribution (v_eps)
1763!> v_eps = - \frac{1}{8*\pi} * |\nabla_r(v)|^2 * \frac{d \eps}{d \rho}
1764!> to the functional derivative of the Hartree energy wrt the density, being
1765!> attributed to the dependency of the dielectric constant to the charge density.
1766!> [see V. M. Sanchez, M. Sued, and D. A. Scherlis, J. Chem. Phys. 131, 174108 (2009)]
1767!> \param pw_pool pool of the original plane-wave grid
1768!> \param dielectric dielectric environment
1769!> \param v Hartree potential
1770!> \param v_eps v_eps
1771!> \par History
1772!> 08.2014 created [Hossein Bani-Hashemian]
1773!> \author Mohammad Hossein Bani-Hashemian
1774! **************************************************************************************************
1775 SUBROUTINE ps_implicit_compute_veps(pw_pool, dielectric, v, v_eps)
1776
1777 TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
1778 TYPE(dielectric_type), INTENT(IN), POINTER :: dielectric
1779 TYPE(pw_r3d_rs_type), INTENT(IN) :: v
1780 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: v_eps
1781
1782 CHARACTER(LEN=*), PARAMETER :: routinen = 'ps_implicit_compute_veps'
1783
1784 INTEGER :: handle, i
1785 REAL(dp) :: eightpi
1786 TYPE(pw_r3d_rs_type) :: dv2
1787 TYPE(pw_r3d_rs_type), DIMENSION(3) :: dv
1788
1789 CALL timeset(routinen, handle)
1790
1791 eightpi = 2*fourpi
1792
1793 CALL pw_pool%create_pw(dv2)
1794 DO i = 1, 3
1795 CALL pw_pool%create_pw(dv(i))
1796 END DO
1797
1798 CALL derive_fft(v, dv, pw_pool)
1799
1800! evaluate |\nabla_r(v)|^2
1801 dv2%array = dv(1)%array**2 + dv(2)%array**2 + dv(3)%array**2
1802
1803 v_eps%array = -(1.0_dp/eightpi)*(dv2%array*dielectric%deps_drho%array)
1804
1805 CALL pw_pool%give_back_pw(dv2)
1806 DO i = 1, 3
1807 CALL pw_pool%give_back_pw(dv(i))
1808 END DO
1809
1810 CALL timestop(handle)
1811
1812 END SUBROUTINE ps_implicit_compute_veps
1813
1814! **************************************************************************************************
1815!> \brief Computes the Hartree energy
1816!> \param density electronic density
1817!> \param v Hartree potential
1818!> \param ehartree Hartree energy
1819!> \par History
1820!> 06.2015 created [Hossein Bani-Hashemian]
1821!> \author Mohammad Hossein Bani-Hashemian
1822! **************************************************************************************************
1823 SUBROUTINE compute_ehartree_periodic_bc(density, v, ehartree)
1824
1825 TYPE(pw_r3d_rs_type), INTENT(IN) :: density, v
1826 REAL(dp), INTENT(OUT) :: ehartree
1827
1828 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_ehartree_periodic_bc'
1829
1830 INTEGER :: handle
1831
1832 CALL timeset(routinen, handle)
1833
1834! E_H = \frac{1}{2} * \int \rho * v dr
1835 ehartree = 0.5_dp*pw_integral_ab(density, v)
1836
1837 CALL timestop(handle)
1838
1839 END SUBROUTINE compute_ehartree_periodic_bc
1840
1841! **************************************************************************************************
1842!> \brief Computes the Hartree energy
1843!> \param dielectric dielectric environment
1844!> \param density electronic density
1845!> \param Btxlambda B^t * \lambda (\lambda is the vector of Lagrange multipliers
1846!> and B^t is the transpose of the boundary operator
1847!> \param v Hartree potential
1848!> \param ehartree Hartree energy
1849!> \param electric_enthalpy electric enthalpy
1850!> \par History
1851!> 06.2015 created [Hossein Bani-Hashemian]
1852!> \author Mohammad Hossein Bani-Hashemian
1853! **************************************************************************************************
1854 SUBROUTINE compute_ehartree_mixed_bc(dielectric, density, Btxlambda, v, ehartree, electric_enthalpy)
1855
1856 TYPE(dielectric_type), INTENT(IN), POINTER :: dielectric
1857 TYPE(pw_r3d_rs_type), INTENT(IN) :: density
1858 REAL(dp), ALLOCATABLE, DIMENSION(:, :, :), &
1859 INTENT(IN) :: btxlambda
1860 TYPE(pw_r3d_rs_type), INTENT(IN) :: v
1861 REAL(dp), INTENT(OUT) :: ehartree, electric_enthalpy
1862
1863 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_ehartree_mixed_bc'
1864
1865 INTEGER :: handle
1866 REAL(dp) :: dvol, ehartree_rho, ehartree_rho_cstr
1867 TYPE(pw_grid_type), POINTER :: pw_grid
1868
1869 CALL timeset(routinen, handle)
1870
1871 pw_grid => v%pw_grid
1872
1873 dvol = pw_grid%dvol
1874
1875! E_H = \frac{1}{2} * \int \rho * v dr + \frac{1}{8 \pi} * \int Btxlambda * v dr
1876! the sign of the second term depends on the sign chosen for the Lagrange multiplier
1877! term in the variational form
1878 ehartree_rho = accurate_sum(density%array*v%array)
1879 ehartree_rho_cstr = accurate_sum(dielectric%eps%array*btxlambda*v%array/fourpi)
1880 ehartree_rho = 0.5_dp*ehartree_rho*dvol
1881 ehartree_rho_cstr = 0.5_dp*ehartree_rho_cstr*dvol
1882 CALL pw_grid%para%group%sum(ehartree_rho)
1883 CALL pw_grid%para%group%sum(ehartree_rho_cstr)
1884 electric_enthalpy = ehartree_rho + ehartree_rho_cstr
1885 ehartree = ehartree_rho - ehartree_rho_cstr
1886
1887 CALL timestop(handle)
1888
1889 END SUBROUTINE compute_ehartree_mixed_bc
1890
1891! **************************************************************************************************
1892!> \brief Computes the (normalized) preconditioned residual norm error and the
1893!> normalized absolute error
1894!> \param pw_pool pool of the original plane-wave grid
1895!> \param green greens functions for FFT based inverse Laplacian
1896!> \param res_new residual
1897!> \param v_old old v
1898!> \param v_new new v
1899!> \param QAinvxres_new Delta^-1(res_new)
1900!> \param pres_error preconditioned residual norm error
1901!> \param nabs_error normalized absolute error
1902!> \par History
1903!> 07.2014 created [Hossein Bani-Hashemian]
1904!> \author Mohammad Hossein Bani-Hashemian
1905! **************************************************************************************************
1906 SUBROUTINE ps_implicit_compute_error_fft(pw_pool, green, res_new, v_old, v_new, &
1907 QAinvxres_new, pres_error, nabs_error)
1908
1909 TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
1910 TYPE(greens_fn_type), INTENT(IN) :: green
1911 TYPE(pw_r3d_rs_type), INTENT(IN) :: res_new, v_old, v_new
1912 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: qainvxres_new
1913 REAL(dp), INTENT(OUT) :: pres_error, nabs_error
1914
1915 CHARACTER(LEN=*), PARAMETER :: routinen = 'ps_implicit_compute_error_fft'
1916
1917 INTEGER :: handle
1918 REAL(dp) :: vol
1919
1920 CALL timeset(routinen, handle)
1921
1922 vol = pw_pool%pw_grid%vol
1923
1924! evaluate \Delta^-1(res) = \Delta^-1 (g - \Delta(v_new) - P(v_new) + Bt \lambda)
1925 CALL apply_inv_laplace_operator_fft(pw_pool, green, res_new, qainvxres_new)
1926! (normalized) preconditioned residual norm error :
1927 pres_error = accurate_sum(qainvxres_new%array(:, :, :)**2)
1928 CALL pw_pool%pw_grid%para%group%sum(pres_error)
1929 pres_error = sqrt(pres_error)/vol
1930
1931! normalized absolute error :
1932! nabs_error := \frac{\| v_old - v_new \|}{volume}
1933 nabs_error = accurate_sum(abs(v_old%array - v_new%array)**2)
1934 CALL pw_pool%pw_grid%para%group%sum(nabs_error)
1935 nabs_error = sqrt(nabs_error)/vol
1936
1937 CALL timestop(handle)
1938
1939 END SUBROUTINE ps_implicit_compute_error_fft
1940
1941! **************************************************************************************************
1942!> \brief Computes the (normalized) preconditioned residual norm error and the
1943!> normalized absolute error
1944!> \param pw_pool pool of the original plane-wave grid
1945!> \param green the greens_fn_type data holding a valid dct_influence_fn
1946!> \param res_new residual
1947!> \param v_old old v
1948!> \param v_new new v
1949!> \param QAinvxres_new Delta^-1(res_new)
1950!> \param pres_error preconditioned residual norm error
1951!> \param nabs_error normalized absolute error
1952!> \par History
1953!> 07.2014 created [Hossein Bani-Hashemian]
1954!> 11.2015 revised [Hossein Bani-Hashemian]
1955!> \author Mohammad Hossein Bani-Hashemian
1956! **************************************************************************************************
1957 SUBROUTINE ps_implicit_compute_error_dct(pw_pool, green, res_new, v_old, v_new, &
1958 QAinvxres_new, pres_error, nabs_error)
1959
1960 TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
1961 TYPE(greens_fn_type), INTENT(IN) :: green
1962 TYPE(pw_r3d_rs_type), INTENT(IN) :: res_new, v_old, v_new
1963 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: qainvxres_new
1964 REAL(dp), INTENT(OUT) :: pres_error, nabs_error
1965
1966 CHARACTER(LEN=*), PARAMETER :: routinen = 'ps_implicit_compute_error_dct'
1967
1968 INTEGER :: handle
1969 REAL(dp) :: vol
1970
1971 CALL timeset(routinen, handle)
1972
1973 vol = pw_pool%pw_grid%vol
1974
1975! evaluate \Delta^-1(res) = \Delta^-1 (g - \Delta(v_new) - P(v_new) + Bt \lambda)
1976 CALL apply_inv_laplace_operator_dct(pw_pool, green, res_new, qainvxres_new)
1977! (normalized) preconditioned residual norm error :
1978 pres_error = accurate_sum(qainvxres_new%array(:, :, :)**2)
1979 CALL pw_pool%pw_grid%para%group%sum(pres_error)
1980 pres_error = sqrt(pres_error)/vol
1981
1982! normalized absolute error :
1983! nabs_error := \frac{\| v_old - v_new \|}{volume}
1984 nabs_error = accurate_sum(abs(v_old%array - v_new%array)**2)
1985 CALL pw_pool%pw_grid%para%group%sum(nabs_error)
1986 nabs_error = sqrt(nabs_error)/vol
1987
1988 CALL timestop(handle)
1989
1990 END SUBROUTINE ps_implicit_compute_error_dct
1991
1992! **************************************************************************************************
1993!> \brief output of the implicit (iterative) Poisson solver
1994!> \param iter current iteration
1995!> \param pres_error preconditioned residual norm error
1996!> \param nabs_error normalized absolute error
1997!> \param outp_unit output unit
1998!> \par History
1999!> 07.2014 created [Hossein Bani-Hashemian]
2000!> \author Mohammad Hossein Bani-Hashemian
2001! **************************************************************************************************
2002 SUBROUTINE ps_implicit_output(iter, pres_error, nabs_error, outp_unit)
2003
2004 INTEGER, INTENT(IN) :: iter
2005 REAL(dp), INTENT(IN) :: pres_error, nabs_error
2006 INTEGER, INTENT(OUT) :: outp_unit
2007
2008 CHARACTER(LEN=*), PARAMETER :: routinen = 'ps_implicit_output'
2009 INTEGER, PARAMETER :: low_print_level = 1
2010
2011 INTEGER :: handle
2012 TYPE(cp_logger_type), POINTER :: logger
2013
2014 CALL timeset(routinen, handle)
2015
2016 logger => cp_get_default_logger()
2017 IF (logger%para_env%is_source()) THEN
2018 outp_unit = cp_logger_get_default_unit_nr(logger, local=.true.)
2019 ELSE
2020 outp_unit = -1
2021 END IF
2022
2023 IF (logger%iter_info%print_level > low_print_level) THEN
2024 IF ((outp_unit > 0) .AND. (iter == 1)) THEN
2025 WRITE (outp_unit, '(T3,A)') &
2026 "POISSON| iter pres error nabs error E_hartree delta E"
2027 END IF
2028
2029 IF (outp_unit > 0) THEN
2030 WRITE (outp_unit, '(T3,A,I6,5X,E13.4,3X,E13.4)', advance='NO') &
2031 "POISSON| ", iter, pres_error, nabs_error
2032 END IF
2033 END IF
2034
2035 CALL timestop(handle)
2036
2037 END SUBROUTINE ps_implicit_output
2038
2039! **************************************************************************************************
2040!> \brief reports the Hartree energy in every iteration
2041!> \param ps_implicit_env the implicit poisson solver environment
2042!> \param outp_unit output unit
2043!> \param ehartree Hartree energy
2044!> \par History
2045!> 07.2014 created [Hossein Bani-Hashemian]
2046!> \author Mohammad Hossein Bani-Hashemian
2047! **************************************************************************************************
2048 SUBROUTINE ps_implicit_report_ehartree(ps_implicit_env, outp_unit, ehartree)
2049
2050 TYPE(ps_implicit_type) :: ps_implicit_env
2051 INTEGER, INTENT(IN) :: outp_unit
2052 REAL(dp), INTENT(IN) :: ehartree
2053
2054 CHARACTER(LEN=*), PARAMETER :: routinen = 'ps_implicit_report_ehartree'
2055 INTEGER, PARAMETER :: low_print_level = 1
2056
2057 INTEGER :: handle
2058 TYPE(cp_logger_type), POINTER :: logger
2059
2060 logger => cp_get_default_logger()
2061 CALL timeset(routinen, handle)
2062 IF (logger%iter_info%print_level > low_print_level) THEN
2063 IF (outp_unit > 0) WRITE (outp_unit, '(F19.10,E10.2)') &
2064 ehartree, ehartree - ps_implicit_env%ehartree
2065 END IF
2066 CALL timestop(handle)
2067
2068 END SUBROUTINE ps_implicit_report_ehartree
2069
2070! **************************************************************************************************
2071!> \brief reports the final number of iteration
2072!> \param iter the iteration number after exiting the main loop
2073!> \param max_iter maximum number of iterations
2074!> \param outp_unit output unit
2075!> \par History
2076!> 02.2016 created [Hossein Bani-Hashemian]
2077!> \author Mohammad Hossein Bani-Hashemian
2078! **************************************************************************************************
2079 SUBROUTINE ps_implicit_print_convergence_msg(iter, max_iter, outp_unit)
2080
2081 INTEGER, INTENT(IN) :: iter, max_iter, outp_unit
2082
2083 CHARACTER(LEN=*), PARAMETER :: routinen = 'ps_implicit_print_convergence_msg'
2084
2085 CHARACTER(LEN=12) :: msg
2086 INTEGER :: handle, last_iter
2087
2088 CALL timeset(routinen, handle)
2089
2090 last_iter = iter - 1
2091
2092 IF (outp_unit > 0) THEN
2093 IF (last_iter == max_iter) THEN
2094 WRITE (outp_unit, '(T3,A)') &
2095 "POISSON| No convergence achieved within the maximum number of iterations."
2096 END IF
2097 IF (last_iter < max_iter) THEN
2098 IF (last_iter == 1) THEN
2099 msg = " iteration."
2100 ELSE
2101 msg = " iterations."
2102 END IF
2103 WRITE (outp_unit, '(T3,A,I0,A)') &
2104 "POISSON| Poisson solver converged in ", last_iter, msg
2105 END IF
2106 END IF
2107 CALL timestop(handle)
2108
2109 END SUBROUTINE ps_implicit_print_convergence_msg
2110
2111! **************************************************************************************************
2112!> \brief converts a 1D array to a 3D array (contiguous layout)
2113!> \param idx_1dto3d mapping of indices
2114!> \param arr1d input 1D array
2115!> \param arr3d input 3D array
2116! **************************************************************************************************
2117 SUBROUTINE convert_1dto3d(idx_1dto3d, arr1d, arr3d)
2118
2119 INTEGER, DIMENSION(:), INTENT(IN) :: idx_1dto3d
2120 REAL(dp), DIMENSION(:), INTENT(IN) :: arr1d
2121 REAL(dp), ALLOCATABLE, DIMENSION(:, :, :), &
2122 INTENT(INOUT) :: arr3d
2123
2124 CHARACTER(LEN=*), PARAMETER :: routinen = 'convert_1dto3d'
2125
2126 INTEGER :: handle, i, j, k, l, lb1, lb2, lb3, &
2127 npts1, npts2, npts3, ub1, ub2, ub3
2128
2129 CALL timeset(routinen, handle)
2130
2131 lb1 = lbound(arr3d, 1); ub1 = ubound(arr3d, 1)
2132 lb2 = lbound(arr3d, 2); ub2 = ubound(arr3d, 2)
2133 lb3 = lbound(arr3d, 3); ub3 = ubound(arr3d, 3)
2134
2135 npts1 = ub1 - lb1 + 1
2136 npts2 = ub2 - lb2 + 1
2137 npts3 = ub3 - lb3 + 1
2138
2139 DO l = 1, SIZE(idx_1dto3d)
2140 k = ((idx_1dto3d(l) - 1)/(npts1*npts2)) + lb3
2141 j = ((idx_1dto3d(l) - 1) - (k - lb3)*npts1*npts2)/npts1 + lb2
2142 i = idx_1dto3d(l) - ((j - lb2)*npts1 + (k - lb3)*npts1*npts2) + lb1 - 1
2143 arr3d(i, j, k) = arr1d(l)
2144 END DO
2145
2146 CALL timestop(handle)
2147
2148 END SUBROUTINE convert_1dto3d
2149
2150! **************************************************************************************************
2151!> \brief Returns the voltage of a tile. In case an alternating field is used, the oltage is a function of time
2152!> \param time ...
2153!> \param v_D ...
2154!> \param osc_frac ...
2155!> \param frequency ...
2156!> \param phase ...
2157!> \param v_D_new ...
2158! **************************************************************************************************
2159 SUBROUTINE get_voltage(time, v_D, osc_frac, frequency, phase, v_D_new)
2160
2161 REAL(dp), INTENT(IN) :: time
2162 REAL(dp), DIMENSION(:), INTENT(IN) :: v_d, osc_frac, frequency, phase
2163 REAL(dp), ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: v_d_new
2164
2165 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_voltage'
2166
2167 INTEGER :: handle, i
2168
2169 CALL timeset(routinen, handle)
2170
2171 ALLOCATE (v_d_new(SIZE(v_d)))
2172
2173 DO i = 1, SIZE(v_d)
2174 v_d_new(i) = v_d(i)*(1 - osc_frac(i)) + &
2175 v_d(i)*osc_frac(i)*cos(2*pi*time*frequency(i) + phase(i))
2176 END DO
2177
2178 CALL timestop(handle)
2179
2180 END SUBROUTINE get_voltage
2181
2182END MODULE ps_implicit_methods
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:225
kahan_sum::accurate_sum
Definition kahan_sum.F:41
pw_methods::pw_axpy
Definition pw_methods.F:165
pw_methods::pw_copy
Definition pw_methods.F:146
pw_methods::pw_integral_ab
Definition pw_methods.F:222
pw_methods::pw_scale
Definition pw_methods.F:94
pw_methods::pw_transfer
Definition pw_methods.F:304
pw_methods::pw_zero
Definition pw_methods.F:83
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::banihashemian2016
integer, save, public banihashemian2016
Definition bibliography.F:36
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
cp_log_handling::cp_logger_get_default_unit_nr
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Definition cp_log_handling.F:576
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
dct
the type I Discrete Cosine Transform (DCT-I)
Definition dct.F:16
dct::neumannx
integer, parameter, public neumannx
Definition dct.F:63
dct::neumannxy
integer, parameter, public neumannxy
Definition dct.F:63
dct::dct_type_init
subroutine, public dct_type_init(pw_grid, neumann_directions, dct_env)
Initializes a dct_type.
Definition dct.F:83
dct::neumannxz
integer, parameter, public neumannxz
Definition dct.F:63
dct::pw_shrink
subroutine, public pw_shrink(neumann_directions, dests_shrink, srcs_shrink, bounds_local_shftd, pw_in, pw_shrinked)
shrinks an evenly symmetric pw_r3d_rs_type data to a pw_r3d_rs_type data that is 8 times smaller (the...
Definition dct.F:700
dct::neumannxyz
integer, parameter, public neumannxyz
Definition dct.F:63
dct::neumannz
integer, parameter, public neumannz
Definition dct.F:63
dct::neumannyz
integer, parameter, public neumannyz
Definition dct.F:63
dct::neumanny
integer, parameter, public neumanny
Definition dct.F:63
dct::pw_expand
subroutine, public pw_expand(neumann_directions, recv_msgs_bnds, dests_expand, srcs_expand, flipg_stat, bounds_shftd, pw_in, pw_expanded)
expands a pw_r3d_rs_type data to an evenly symmetric pw_r3d_rs_type data that is 8 times larger than ...
Definition dct.F:516
dielectric_methods
methods for evaluating the dielectric constant
Definition dielectric_methods.F:14
dielectric_methods::dielectric_create
subroutine, public dielectric_create(dielectric, pw_pool, dielectric_params)
allocates memory for a dielectric data type
Definition dielectric_methods.F:68
dielectric_methods::derive_fft
subroutine, public derive_fft(f, df, pw_pool)
computes the derivative of a function using FFT
Definition dielectric_methods.F:1202
dielectric_types
dielectric constant data type
Definition dielectric_types.F:14
dirichlet_bc_methods
subroutines for defining and creating Dirichlet type subdomains
Definition dirichlet_bc_methods.F:15
dirichlet_bc_methods::dirichlet_boundary_region_setup
subroutine, public dirichlet_boundary_region_setup(pw_pool, poisson_params, dbcs)
Sets up the Dirichlet boundary condition.
Definition dirichlet_bc_methods.F:65
dirichlet_bc_types
Dirichlet boundary condition data types.
Definition dirichlet_bc_types.F:15
dirichlet_bc_types::dbc_tile_release
subroutine, public dbc_tile_release(dbc, pw_pool)
releases tiles
Definition dirichlet_bc_types.F:195
kahan_sum
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Definition kahan_sum.F:29
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathconstants::fourpi
real(kind=dp), parameter, public fourpi
Definition mathconstants.F:113
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
pw_grid_types
Definition pw_grid_types.F:13
pw_methods
Definition pw_methods.F:25
pw_poisson_types
functions related to the poisson solver on regular grids
Definition pw_poisson_types.F:22
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_pool_release
subroutine, public pw_pool_release(pool)
releases the given pool (see cp2k/doc/ReferenceCounting.html)
Definition pw_pool_types.F:259
pw_pool_types::pw_pool_create
subroutine, public pw_pool_create(pool, pw_grid, max_cache)
creates a pool for pw
Definition pw_pool_types.F:173
pw_types
Definition pw_types.F:24
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
dct::dct_type
Definition dct.F:32
dielectric_types::dielectric_type
Definition dielectric_types.F:62
ps_implicit_types::ps_implicit_type
Definition ps_implicit_types.F:46
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_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