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