(git:e7e05ae)
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 ! **************************************************************************************************
16 MODULE 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 
74 CONTAINS
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 
2172 END 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
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