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