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 56 : 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 56 : CALL timeset(routineN, handle)
102 :
103 56 : CALL cite_reference(BaniHashemian2016)
104 :
105 56 : IF (.NOT. ASSOCIATED(ps_implicit_env)) THEN
106 2296 : ALLOCATE (ps_implicit_env)
107 :
108 56 : ps_implicit_env%do_dbc_cube = poisson_params%dbc_params%do_dbc_cube
109 56 : boundary_condition = poisson_params%ps_implicit_params%boundary_condition
110 56 : neumann_directions = poisson_params%ps_implicit_params%neumann_directions
111 :
112 : ! create dielectric
113 : NULLIFY (ps_implicit_env%dielectric)
114 34 : SELECT CASE (boundary_condition)
115 : CASE (PERIODIC_BC, MIXED_PERIODIC_BC)
116 34 : 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 78 : CALL pw_pool_release(pw_pool_xpndd)
121 : END SELECT
122 :
123 : ! initial guess
124 56 : NULLIFY (ps_implicit_env%initial_guess)
125 :
126 : ! v_eps
127 56 : NULLIFY (ps_implicit_env%v_eps)
128 56 : ALLOCATE (ps_implicit_env%v_eps)
129 56 : CALL pw_pool%create_pw(ps_implicit_env%v_eps)
130 56 : CALL pw_zero(ps_implicit_env%v_eps)
131 :
132 : ! constraint charge
133 56 : 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 72 : CALL pw_pool_release(pw_pool_xpndd)
145 : END SELECT
146 :
147 : ! initialize energies
148 56 : ps_implicit_env%ehartree = 0.0_dp
149 56 : ps_implicit_env%electric_enthalpy = 0.0_dp
150 : ! times called
151 56 : ps_implicit_env%times_called = 0
152 :
153 : ! dct env
154 56 : 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 56 : CALL dirichlet_boundary_region_setup(pw_pool, poisson_params, ps_implicit_env%contacts)
160 56 : 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 56 : IF ((boundary_condition == MIXED_PERIODIC_BC .OR. boundary_condition == MIXED_BC) .AND. &
163 : (.NOT. poisson_params%dbc_params%do_dbc_cube)) THEN
164 36 : n_contacts = SIZE(ps_implicit_env%contacts)
165 160 : DO j = 1, n_contacts
166 160 : 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 56 : CALL timestop(handle)
173 :
174 56 : 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 150 : 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 150 : CALL timeset(routineN, handle)
208 :
209 150 : pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool
210 150 : dielectric => poisson_env%implicit_env%dielectric
211 150 : green => poisson_env%green_fft
212 150 : ps_implicit_env => poisson_env%implicit_env
213 :
214 150 : tol = poisson_env%parameters%ps_implicit_params%tol
215 150 : omega = poisson_env%parameters%ps_implicit_params%omega
216 150 : max_iter = poisson_env%parameters%ps_implicit_params%max_iter
217 150 : use_zero_initial_guess = poisson_env%parameters%ps_implicit_params%zero_initial_guess
218 150 : times_called = ps_implicit_env%times_called
219 :
220 : ! check if this is the first scf iteration
221 150 : IF (times_called == 0) CALL ps_implicit_initial_guess_create(ps_implicit_env, pw_pool)
222 :
223 150 : CALL pw_pool%create_pw(g)
224 150 : CALL pw_pool%create_pw(v_old)
225 150 : CALL pw_pool%create_pw(res_old)
226 150 : CALL pw_pool%create_pw(res_new)
227 150 : CALL pw_pool%create_pw(QAinvxres)
228 150 : CALL pw_pool%create_pw(PxQAinvxres)
229 :
230 150 : IF (use_zero_initial_guess) THEN
231 0 : CALL pw_zero(v_old)
232 : ELSE
233 150 : CALL pw_copy(ps_implicit_env%initial_guess, v_old)
234 : END IF
235 :
236 22241170 : g%array = fourpi*density%array/dielectric%eps%array
237 :
238 : ! res_old = g - \Delta(v_old) - P(v_old)
239 150 : CALL apply_poisson_operator_fft(pw_pool, green, dielectric, v_old, res_old)
240 150 : CALL pw_scale(res_old, -1.0_dp)
241 150 : CALL pw_axpy(g, res_old)
242 :
243 : ! evaluate \Delta^-1(res_old)
244 150 : CALL apply_inv_laplace_operator_fft(pw_pool, green, res_old, QAinvxres)
245 :
246 150 : iter = 1
247 1472 : DO
248 :
249 : ! v_new = v_old + \omega * QAinvxres_old
250 736 : CALL pw_scale(QAinvxres, omega)
251 736 : CALL pw_copy(QAinvxres, v_new)
252 736 : 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 736 : CALL apply_P_operator(pw_pool, dielectric, QAinvxres, PxQAinvxres)
257 736 : CALL pw_copy(PxQAinvxres, res_new)
258 736 : CALL pw_scale(res_new, -1.0_dp)
259 736 : 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 736 : pres_error, nabs_error)
264 : ! output
265 736 : CALL ps_implicit_output(iter, pres_error, nabs_error, outp_unit)
266 736 : IF (PRESENT(ehartree)) THEN
267 736 : CALL ps_implicit_compute_ehartree(density, v_new, ehartree)
268 736 : CALL ps_implicit_report_ehartree(ps_implicit_env, outp_unit, ehartree)
269 736 : ps_implicit_env%ehartree = ehartree
270 : ELSE
271 0 : IF (outp_unit > 0) WRITE (outp_unit, '(A1,/)')
272 : END IF
273 :
274 736 : iter = iter + 1
275 736 : reached_max_iter = iter > max_iter
276 736 : reached_tol = pres_error <= tol
277 736 : IF (pres_error > large_error) THEN
278 0 : CPABORT("Poisson solver did not converge.")
279 : END IF
280 736 : ps_implicit_env%times_called = ps_implicit_env%times_called + 1
281 736 : 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 150 : CALL ps_implicit_print_convergence_msg(iter, max_iter, outp_unit)
289 :
290 150 : IF ((times_called /= 0) .AND. (.NOT. use_zero_initial_guess)) THEN
291 138 : CALL pw_copy(v_new, ps_implicit_env%initial_guess)
292 : END IF
293 :
294 150 : 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 150 : v_eps%pw_grid => ps_implicit_env%v_eps%pw_grid
299 150 : v_eps%array => ps_implicit_env%v_eps%array
300 150 : CALL ps_implicit_compute_veps(pw_pool, dielectric, v_new, v_eps)
301 : END BLOCK
302 :
303 150 : CALL pw_pool%give_back_pw(g)
304 150 : CALL pw_pool%give_back_pw(v_old)
305 150 : CALL pw_pool%give_back_pw(res_old)
306 150 : CALL pw_pool%give_back_pw(res_new)
307 150 : CALL pw_pool%give_back_pw(QAinvxres)
308 150 : CALL pw_pool%give_back_pw(PxQAinvxres)
309 :
310 150 : CALL timestop(handle)
311 :
312 150 : 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 56 : 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 56 : CALL timeset(routineN, handle)
1065 :
1066 56 : n_tiles_tot = SIZE(ps_implicit_env%v_D)
1067 56 : NULLIFY (ps_implicit_env%initial_guess)
1068 56 : ALLOCATE (ps_implicit_env%initial_guess)
1069 56 : CALL pw_pool%create_pw(ps_implicit_env%initial_guess)
1070 56 : CALL pw_zero(ps_implicit_env%initial_guess)
1071 168 : ALLOCATE (ps_implicit_env%initial_lambda(n_tiles_tot))
1072 298 : ps_implicit_env%initial_lambda = 0.0_dp
1073 :
1074 56 : CALL timestop(handle)
1075 :
1076 56 : 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 56 : 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 56 : 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 56 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: Bxunit_vec, test_vec, work_arr
1109 56 : 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 56 : CALL timeset(routineN, handle)
1117 :
1118 56 : pw_grid_orig => pw_pool_orig%pw_grid
1119 :
1120 56 : logger => cp_get_default_logger()
1121 56 : IF (logger%para_env%is_source()) THEN
1122 28 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
1123 : ELSE
1124 : unit_nr = -1
1125 : END IF
1126 :
1127 72 : 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 : ! Give up on OpenMP for now (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 1536 : DO k = lb3, ub3
1168 163376 : DO j = lb2, ub2
1169 7922620 : DO i = lb1, ub1
1170 : ps_implicit_env%idx_1dto3d(l) = (i - lb1 + 1) + &
1171 : (j - lb2)*npts_local(1) + &
1172 7759260 : (k - lb3)*npts_local(1)*npts_local(2)
1173 7921100 : l = l + 1
1174 : END DO
1175 : END DO
1176 : END DO
1177 :
1178 16 : n_contacts = SIZE(ps_implicit_env%contacts)
1179 16 : n_tiles_tot = 0
1180 56 : DO j = 1, n_contacts
1181 56 : n_tiles_tot = n_tiles_tot + ps_implicit_env%contacts(j)%dirichlet_bc%n_tiles
1182 : END DO
1183 :
1184 64 : ALLOCATE (ps_implicit_env%B(n_tiles_tot, data_size))
1185 48 : ALLOCATE (ps_implicit_env%Bt(data_size, n_tiles_tot))
1186 64 : ALLOCATE (ps_implicit_env%QS(n_tiles_tot, n_tiles_tot))
1187 64 : ALLOCATE (ps_implicit_env%Rinv(n_tiles_tot + 1, n_tiles_tot + 1))
1188 48 : ALLOCATE (ps_implicit_env%v_D(n_tiles_tot))
1189 32 : ALLOCATE (ps_implicit_env%osc_frac(n_tiles_tot))
1190 32 : ALLOCATE (ps_implicit_env%frequency(n_tiles_tot))
1191 32 : ALLOCATE (ps_implicit_env%phase(n_tiles_tot))
1192 :
1193 48 : ALLOCATE (QAinvxBt(data_size, n_tiles_tot))
1194 32 : ALLOCATE (Bxunit_vec(n_tiles_tot))
1195 32 : ALLOCATE (test_vec(n_tiles_tot))
1196 48 : ALLOCATE (R(n_tiles_tot + 1, n_tiles_tot + 1))
1197 80 : ALLOCATE (work_arr(n_tiles_tot + 1), ipiv(n_tiles_tot + 1)) ! LAPACK work and ipiv arrays
1198 :
1199 : ! prepare pw_pool for evaluating inverse Laplacian of tile_pw's using DCT
1200 16 : CALL pw_pool_create(pw_pool_xpndd, pw_grid=dct_pw_grid)
1201 :
1202 : ! set up B, B^t, (\Delta^-1)*B^t
1203 16 : indx1 = 1
1204 56 : DO j = 1, n_contacts
1205 40 : n_tiles = ps_implicit_env%contacts(j)%dirichlet_bc%n_tiles
1206 40 : indx2 = indx1 + n_tiles - 1
1207 90 : DO i = 1, n_tiles
1208 :
1209 50 : CALL pw_pool_xpndd%create_pw(pw_in)
1210 : CALL pw_expand(neumann_directions, &
1211 : dct_env%recv_msgs_bnds, dct_env%dests_expand, dct_env%srcs_expand, &
1212 : dct_env%flipg_stat, dct_env%bounds_shftd, &
1213 50 : ps_implicit_env%contacts(j)%dirichlet_bc%tiles(i)%tile%tile_pw, pw_in)
1214 :
1215 50 : tile_volume = ps_implicit_env%contacts(j)%dirichlet_bc%tiles(i)%tile%volume
1216 50 : CALL pw_scale(pw_in, 1.0_dp/(vol_scfac*tile_volume)) ! normalize tile_pw
1217 45985636 : ps_implicit_env%Bt(ps_implicit_env%idx_1dto3d, indx1 + i - 1) = RESHAPE(pw_in%array, [data_size])
1218 :
1219 50 : CALL pw_pool_xpndd%create_pw(pw_out)
1220 50 : CALL apply_inv_laplace_operator_dct(pw_pool_xpndd, green, pw_in, pw_out)
1221 45985636 : QAinvxBt(ps_implicit_env%idx_1dto3d, indx1 + i - 1) = RESHAPE(pw_out%array, [data_size])
1222 : ! the electrostatic potential has opposite sign by internal convention
1223 50 : ps_implicit_env%v_D(indx1 + i - 1) = -1.0_dp*ps_implicit_env%contacts(j)%dirichlet_bc%v_D
1224 50 : ps_implicit_env%osc_frac(indx1 + i - 1) = ps_implicit_env%contacts(j)%dirichlet_bc%osc_frac
1225 50 : ps_implicit_env%frequency(indx1 + i - 1) = ps_implicit_env%contacts(j)%dirichlet_bc%frequency
1226 50 : ps_implicit_env%phase(indx1 + i - 1) = ps_implicit_env%contacts(j)%dirichlet_bc%phase
1227 :
1228 50 : CALL pw_pool_xpndd%give_back_pw(pw_in)
1229 90 : CALL pw_pool_xpndd%give_back_pw(pw_out)
1230 : END DO
1231 56 : indx1 = indx2 + 1
1232 : END DO
1233 61504072 : ps_implicit_env%B(:, :) = TRANSPOSE(ps_implicit_env%Bt)
1234 :
1235 : ! evaluate QS = - B*(\Delta^-1)*B^t
1236 16 : IF (data_size /= 0) THEN
1237 : CALL DGEMM('N', 'N', n_tiles_tot, n_tiles_tot, data_size, &
1238 : -1.0_dp, ps_implicit_env%B, n_tiles_tot, QAinvxBt, &
1239 16 : data_size, 0.0_dp, ps_implicit_env%QS, n_tiles_tot)
1240 : END IF
1241 16 : CALL pw_grid_orig%para%group%sum(ps_implicit_env%QS)
1242 :
1243 : ! evaluate B*1
1244 22992834 : Bxunit_vec(:) = SUM(ps_implicit_env%B, 2)/ngpts
1245 16 : CALL pw_grid_orig%para%group%sum(Bxunit_vec)
1246 : ! set up R = [QS B*1; (B*1)^t 0]
1247 16 : R = 0.0_dp
1248 288 : R(1:n_tiles_tot, 1:n_tiles_tot) = ps_implicit_env%QS
1249 66 : R(1:n_tiles_tot, n_tiles_tot + 1) = Bxunit_vec
1250 66 : R(n_tiles_tot + 1, 1:n_tiles_tot) = Bxunit_vec
1251 : ! evaluate R^(-1)
1252 420 : ps_implicit_env%Rinv(:, :) = R
1253 16 : CALL DGETRF(n_tiles_tot + 1, n_tiles_tot + 1, ps_implicit_env%Rinv, n_tiles_tot + 1, ipiv, info)
1254 16 : IF (info /= 0) THEN
1255 : CALL cp_abort(__LOCATION__, &
1256 : "R is (nearly) singular! Either two Dirichlet constraints are identical or "// &
1257 0 : "you need to reduce the number of tiles.")
1258 : END IF
1259 16 : CALL DGETRI(n_tiles_tot + 1, ps_implicit_env%Rinv, n_tiles_tot + 1, ipiv, work_arr, n_tiles_tot + 1, info)
1260 16 : IF (info /= 0) THEN
1261 0 : CPABORT("Inversion of R failed!")
1262 : END IF
1263 :
1264 16 : DEALLOCATE (QAinvxBt, Bxunit_vec, R, work_arr, ipiv)
1265 16 : CALL pw_pool_release(pw_pool_xpndd)
1266 :
1267 16 : done_preparing = .TRUE.
1268 16 : CALL pw_grid_orig%para%group%sum(done_preparing)
1269 16 : IF ((unit_nr > 0) .AND. done_preparing) THEN
1270 8 : WRITE (unit_nr, "(T3,A,/,T3,A,/,A)") "POISSON| ... Done. ", REPEAT('-', 78)
1271 : END IF
1272 :
1273 : CASE (MIXED_PERIODIC_BC)
1274 :
1275 22 : ngpts_local = pw_grid_orig%ngpts_local
1276 22 : ngpts = pw_grid_orig%ngpts
1277 88 : npts_local = pw_grid_orig%npts_local
1278 88 : npts = pw_grid_orig%npts
1279 220 : bounds_local = pw_grid_orig%bounds_local
1280 220 : bounds = pw_grid_orig%bounds
1281 22 : dct_env => ps_implicit_env%dct_env
1282 :
1283 : ! evaluate indices for converting 3D arrays into 1D arrays
1284 22 : lb1 = bounds_local(1, 1); ub1 = bounds_local(2, 1)
1285 22 : lb2 = bounds_local(1, 2); ub2 = bounds_local(2, 2)
1286 22 : lb3 = bounds_local(1, 3); ub3 = bounds_local(2, 3)
1287 :
1288 22 : IF (pw_grid_orig%para%blocked) THEN
1289 0 : data_size = PRODUCT(npts_local)
1290 22 : ELSE IF (pw_grid_orig%para%ray_distribution) THEN
1291 22 : data_size = ngpts_local
1292 : ELSE ! parallel run with np = 1
1293 0 : data_size = PRODUCT(npts_local)
1294 : END IF
1295 :
1296 66 : ALLOCATE (ps_implicit_env%idx_1dto3d(data_size))
1297 22 : l = 1
1298 : ! Give up on OpenMP for now (at least the Intel compiler has an issue here)
1299 : ! An automatic OpenMP parallelization of this loop might be tricky
1300 : ! because of the l incrementation
1301 1444 : DO k = lb3, ub3
1302 95134 : DO j = lb2, ub2
1303 3231189 : DO i = lb1, ub1
1304 : ps_implicit_env%idx_1dto3d(l) = (i - lb1 + 1) + &
1305 : (j - lb2)*npts_local(1) + &
1306 3136077 : (k - lb3)*npts_local(1)*npts_local(2)
1307 3229767 : l = l + 1
1308 : END DO
1309 : END DO
1310 : END DO
1311 :
1312 22 : n_contacts = SIZE(ps_implicit_env%contacts)
1313 22 : n_tiles_tot = 0
1314 110 : DO j = 1, n_contacts
1315 110 : n_tiles_tot = n_tiles_tot + ps_implicit_env%contacts(j)%dirichlet_bc%n_tiles
1316 : END DO
1317 :
1318 88 : ALLOCATE (ps_implicit_env%B(n_tiles_tot, data_size))
1319 66 : ALLOCATE (ps_implicit_env%Bt(data_size, n_tiles_tot))
1320 88 : ALLOCATE (ps_implicit_env%QS(n_tiles_tot, n_tiles_tot))
1321 88 : ALLOCATE (ps_implicit_env%Rinv(n_tiles_tot + 1, n_tiles_tot + 1))
1322 66 : ALLOCATE (ps_implicit_env%v_D(n_tiles_tot))
1323 44 : ALLOCATE (ps_implicit_env%osc_frac(n_tiles_tot))
1324 44 : ALLOCATE (ps_implicit_env%frequency(n_tiles_tot))
1325 44 : ALLOCATE (ps_implicit_env%phase(n_tiles_tot))
1326 :
1327 66 : ALLOCATE (QAinvxBt(data_size, n_tiles_tot))
1328 44 : ALLOCATE (Bxunit_vec(n_tiles_tot))
1329 44 : ALLOCATE (test_vec(n_tiles_tot))
1330 66 : ALLOCATE (R(n_tiles_tot + 1, n_tiles_tot + 1))
1331 110 : ALLOCATE (work_arr(n_tiles_tot + 1), ipiv(n_tiles_tot + 1))
1332 :
1333 : ! set up B, B^t, (\Delta^-1)*B^t
1334 110 : indx1 = 1
1335 110 : DO j = 1, n_contacts
1336 88 : n_tiles = ps_implicit_env%contacts(j)%dirichlet_bc%n_tiles
1337 88 : indx2 = indx1 + n_tiles - 1
1338 262 : DO i = 1, n_tiles
1339 174 : CALL pw_pool_orig%create_pw(pw_in)
1340 174 : CALL pw_copy(ps_implicit_env%contacts(j)%dirichlet_bc%tiles(i)%tile%tile_pw, pw_in)
1341 :
1342 174 : tile_volume = ps_implicit_env%contacts(j)%dirichlet_bc%tiles(i)%tile%volume
1343 174 : CALL pw_scale(pw_in, 1.0_dp/tile_volume) ! normalize tile_pw
1344 40325064 : ps_implicit_env%Bt(ps_implicit_env%idx_1dto3d, indx1 + i - 1) = RESHAPE(pw_in%array, [data_size])
1345 :
1346 174 : CALL pw_pool_orig%create_pw(pw_out)
1347 174 : CALL apply_inv_laplace_operator_fft(pw_pool_orig, green, pw_in, pw_out)
1348 40325064 : QAinvxBt(ps_implicit_env%idx_1dto3d, indx1 + i - 1) = RESHAPE(pw_out%array, [data_size])
1349 : ! the electrostatic potential has opposite sign by internal convention
1350 174 : ps_implicit_env%v_D(indx1 + i - 1) = -1.0_dp*ps_implicit_env%contacts(j)%dirichlet_bc%v_D
1351 174 : ps_implicit_env%osc_frac(indx1 + i - 1) = ps_implicit_env%contacts(j)%dirichlet_bc%osc_frac
1352 174 : ps_implicit_env%frequency(indx1 + i - 1) = ps_implicit_env%contacts(j)%dirichlet_bc%frequency
1353 174 : ps_implicit_env%phase(indx1 + i - 1) = ps_implicit_env%contacts(j)%dirichlet_bc%phase
1354 :
1355 174 : CALL pw_pool_orig%give_back_pw(pw_in)
1356 262 : CALL pw_pool_orig%give_back_pw(pw_out)
1357 : END DO
1358 110 : indx1 = indx2 + 1
1359 : END DO
1360 46596892 : ps_implicit_env%B(:, :) = TRANSPOSE(ps_implicit_env%Bt)
1361 :
1362 : ! evaluate QS = - B*(\Delta^-1)*B^t
1363 22 : IF (data_size /= 0) THEN
1364 : CALL DGEMM('N', 'N', n_tiles_tot, n_tiles_tot, data_size, &
1365 : -1.0_dp, ps_implicit_env%B, n_tiles_tot, QAinvxBt, &
1366 22 : data_size, 0.0_dp, ps_implicit_env%QS, n_tiles_tot)
1367 : END IF
1368 22 : CALL pw_grid_orig%para%group%sum(ps_implicit_env%QS)
1369 :
1370 : ! evaluate B*1
1371 20162554 : Bxunit_vec(:) = SUM(ps_implicit_env%B, 2)/ngpts
1372 22 : CALL pw_grid_orig%para%group%sum(Bxunit_vec)
1373 : ! set up R = [QS B*1; (B*1)^t 0]
1374 22 : R = 0.0_dp
1375 3074 : R(1:n_tiles_tot, 1:n_tiles_tot) = ps_implicit_env%QS
1376 196 : R(1:n_tiles_tot, n_tiles_tot + 1) = Bxunit_vec
1377 196 : R(n_tiles_tot + 1, 1:n_tiles_tot) = Bxunit_vec
1378 : ! evaluate R^(-1)
1379 3466 : ps_implicit_env%Rinv(:, :) = R
1380 22 : CALL DGETRF(n_tiles_tot + 1, n_tiles_tot + 1, ps_implicit_env%Rinv, n_tiles_tot + 1, ipiv, info)
1381 22 : IF (info /= 0) THEN
1382 : CALL cp_abort(__LOCATION__, &
1383 : "R is (nearly) singular! Either two Dirichlet constraints are identical or "// &
1384 0 : "you need to reduce the number of tiles.")
1385 : END IF
1386 22 : CALL DGETRI(n_tiles_tot + 1, ps_implicit_env%Rinv, n_tiles_tot + 1, ipiv, work_arr, n_tiles_tot + 1, info)
1387 22 : IF (info /= 0) THEN
1388 0 : CPABORT("Inversion of R failed!")
1389 : END IF
1390 :
1391 22 : DEALLOCATE (QAinvxBt, Bxunit_vec, R, work_arr, ipiv)
1392 :
1393 22 : done_preparing = .TRUE.
1394 22 : CALL pw_grid_orig%para%group%sum(done_preparing)
1395 22 : IF ((unit_nr > 0) .AND. done_preparing) THEN
1396 11 : WRITE (unit_nr, "(T3,A,/,T3,A,/,A)") "POISSON| ... Done. ", REPEAT('-', 78)
1397 : END IF
1398 :
1399 : CASE (PERIODIC_BC, NEUMANN_BC)
1400 :
1401 18 : ALLOCATE (ps_implicit_env%idx_1dto3d(1))
1402 18 : ALLOCATE (ps_implicit_env%B(1, 1))
1403 18 : ALLOCATE (ps_implicit_env%Bt(1, 1))
1404 18 : ALLOCATE (ps_implicit_env%QS(1, 1))
1405 18 : ALLOCATE (ps_implicit_env%Rinv(1, 1))
1406 18 : ALLOCATE (ps_implicit_env%v_D(1))
1407 18 : ALLOCATE (ps_implicit_env%osc_frac(1))
1408 18 : ALLOCATE (ps_implicit_env%frequency(1))
1409 18 : ALLOCATE (ps_implicit_env%phase(1))
1410 :
1411 36 : ps_implicit_env%idx_1dto3d = 1
1412 54 : ps_implicit_env%B = 0.0_dp
1413 54 : ps_implicit_env%Bt = 0.0_dp
1414 54 : ps_implicit_env%QS = 0.0_dp
1415 54 : ps_implicit_env%Rinv = 0.0_dp
1416 36 : ps_implicit_env%v_D = 0.0_dp
1417 :
1418 : CASE DEFAULT
1419 : CALL cp_abort(__LOCATION__, &
1420 : "Please specify the type of boundary conditions using the "// &
1421 56 : "input file keyword BOUNDARY_CONDITIONS.")
1422 : END SELECT
1423 :
1424 56 : CALL timestop(handle)
1425 :
1426 112 : END SUBROUTINE ps_implicit_prepare_blocks
1427 :
1428 : ! **************************************************************************************************
1429 : !> \brief Evaluates the action of the operator P on a given matrix v, defined
1430 : !> as: P(v) := - \nabla_r(\ln(\eps)) \cdot \nabla_r(v)
1431 : !> \param pw_pool pool of pw grid
1432 : !> \param dielectric dielectric_type containing eps
1433 : !> \param v input matrix
1434 : !> \param Pxv action of the operator P on v
1435 : !> \par History
1436 : !> 07.2014 created [Hossein Bani-Hashemian]
1437 : !> \author Mohammad Hossein Bani-Hashemian
1438 : ! **************************************************************************************************
1439 4758 : SUBROUTINE apply_P_operator(pw_pool, dielectric, v, Pxv)
1440 :
1441 : TYPE(pw_pool_type), POINTER :: pw_pool
1442 : TYPE(dielectric_type), INTENT(IN), POINTER :: dielectric
1443 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v
1444 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: Pxv
1445 :
1446 : CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_P_operator'
1447 :
1448 : INTEGER :: handle, i
1449 19032 : TYPE(pw_r3d_rs_type), DIMENSION(3) :: dv
1450 :
1451 4758 : CALL timeset(routineN, handle)
1452 :
1453 19032 : DO i = 1, 3
1454 19032 : CALL pw_pool%create_pw(dv(i))
1455 : END DO
1456 :
1457 4758 : CALL derive_fft(v, dv, pw_pool)
1458 : ASSOCIATE (dln_eps => dielectric%dln_eps)
1459 : Pxv%array = -(dv(1)%array*dln_eps(1)%array + &
1460 : dv(2)%array*dln_eps(2)%array + &
1461 850866506 : dv(3)%array*dln_eps(3)%array)
1462 : END ASSOCIATE
1463 :
1464 19032 : DO i = 1, 3
1465 19032 : CALL pw_pool%give_back_pw(dv(i))
1466 : END DO
1467 :
1468 4758 : CALL timestop(handle)
1469 :
1470 4758 : END SUBROUTINE apply_P_operator
1471 :
1472 : ! **************************************************************************************************
1473 : !> \brief Evaluates the action of the inverse of the Laplace operator on a given 3d matrix
1474 : !> \param pw_pool pool of pw grid
1475 : !> \param green green functions for FFT based inverse Laplacian
1476 : !> \param pw_in pw_in (density)
1477 : !> \param pw_out pw_out (potential)
1478 : !> \par History
1479 : !> 07.2014 created [Hossein Bani-Hashemian]
1480 : !> \author Mohammad Hossein Bani-Hashemian
1481 : ! **************************************************************************************************
1482 2770 : SUBROUTINE apply_inv_laplace_operator_fft(pw_pool, green, pw_in, pw_out)
1483 :
1484 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
1485 : TYPE(greens_fn_type), INTENT(IN) :: green
1486 : TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in
1487 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pw_out
1488 :
1489 : CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_inv_laplace_operator_fft'
1490 :
1491 : INTEGER :: handle, ig, ng
1492 : REAL(dp) :: prefactor
1493 : TYPE(pw_c1d_gs_type) :: pw_in_gs
1494 : TYPE(pw_grid_type), POINTER :: pw_grid
1495 :
1496 2770 : CALL timeset(routineN, handle)
1497 :
1498 : ! here I divide by fourpi to cancel out the prefactor fourpi in influence_fn
1499 2770 : prefactor = 1.0_dp/fourpi
1500 :
1501 2770 : pw_grid => pw_pool%pw_grid
1502 2770 : ng = SIZE(pw_grid%gsq)
1503 :
1504 2770 : CALL pw_pool%create_pw(pw_in_gs)
1505 :
1506 2770 : CALL pw_transfer(pw_in, pw_in_gs)
1507 385725972 : DO ig = 1, ng
1508 385725972 : pw_in_gs%array(ig) = prefactor*pw_in_gs%array(ig)*green%influence_fn%array(ig)
1509 : END DO
1510 2770 : CALL pw_transfer(pw_in_gs, pw_out)
1511 :
1512 2770 : CALL pw_pool%give_back_pw(pw_in_gs)
1513 :
1514 2770 : CALL timestop(handle)
1515 :
1516 2770 : END SUBROUTINE apply_inv_laplace_operator_fft
1517 :
1518 : ! **************************************************************************************************
1519 : !> \brief Evaluates the action of the inverse of the Laplace operator on a given
1520 : !> 3d matrix using DCT-I
1521 : !> \param pw_pool pool of pw grid
1522 : !> \param green the greens_fn_type data holding a valid dct_influence_fn
1523 : !> \param pw_in pw_in (density)
1524 : !> \param pw_out pw_out (potential)
1525 : !> \par History
1526 : !> 07.2014 created [Hossein Bani-Hashemian]
1527 : !> 11.2015 revised [Hossein Bani-Hashemian]
1528 : !> \author Mohammad Hossein Bani-Hashemian
1529 : ! **************************************************************************************************
1530 474 : SUBROUTINE apply_inv_laplace_operator_dct(pw_pool, green, pw_in, pw_out)
1531 :
1532 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
1533 : TYPE(greens_fn_type), INTENT(IN) :: green
1534 : TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in
1535 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pw_out
1536 :
1537 : CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_inv_laplace_operator_dct'
1538 :
1539 : INTEGER :: handle, ig, ng
1540 : REAL(dp) :: prefactor
1541 : TYPE(pw_c1d_gs_type) :: pw_in_gs
1542 : TYPE(pw_grid_type), POINTER :: pw_grid
1543 :
1544 474 : CALL timeset(routineN, handle)
1545 :
1546 : ! here I divide by fourpi to cancel out the prefactor fourpi in influence_fn
1547 474 : prefactor = 1.0_dp/fourpi
1548 :
1549 474 : pw_grid => pw_pool%pw_grid
1550 474 : ng = SIZE(pw_grid%gsq)
1551 :
1552 474 : CALL pw_pool%create_pw(pw_in_gs)
1553 :
1554 474 : CALL pw_transfer(pw_in, pw_in_gs)
1555 205618218 : DO ig = 1, ng
1556 205618218 : pw_in_gs%array(ig) = prefactor*pw_in_gs%array(ig)*green%dct_influence_fn%array(ig)
1557 : END DO
1558 474 : CALL pw_transfer(pw_in_gs, pw_out)
1559 :
1560 474 : CALL pw_pool%give_back_pw(pw_in_gs)
1561 :
1562 474 : CALL timestop(handle)
1563 :
1564 474 : END SUBROUTINE apply_inv_laplace_operator_dct
1565 :
1566 : ! **************************************************************************************************
1567 : !> \brief Evaluates the action of the Laplace operator on a given 3d matrix
1568 : !> \param pw_pool pool of pw grid
1569 : !> \param green green functions for FFT based inverse Laplacian
1570 : !> \param pw_in pw_in (potential)
1571 : !> \param pw_out pw_out (density)
1572 : !> \par History
1573 : !> 07.2014 created [Hossein Bani-Hashemian]
1574 : !> \author Mohammad Hossein Bani-Hashemian
1575 : ! **************************************************************************************************
1576 342 : SUBROUTINE apply_laplace_operator_fft(pw_pool, green, pw_in, pw_out)
1577 :
1578 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
1579 : TYPE(greens_fn_type), INTENT(IN) :: green
1580 : TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in
1581 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pw_out
1582 :
1583 : CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_laplace_operator_fft'
1584 :
1585 : INTEGER :: g0_index, handle, ig, ng
1586 : LOGICAL :: have_g0
1587 : REAL(dp) :: prefactor
1588 : TYPE(pw_c1d_gs_type) :: pw_in_gs
1589 : TYPE(pw_grid_type), POINTER :: pw_grid
1590 :
1591 342 : CALL timeset(routineN, handle)
1592 :
1593 : ! here I multiply by fourpi to cancel out the prefactor fourpi in influence_fn
1594 342 : prefactor = fourpi
1595 :
1596 342 : pw_grid => pw_pool%pw_grid
1597 342 : ng = SIZE(pw_in%pw_grid%gsq)
1598 342 : have_g0 = green%influence_fn%pw_grid%have_g0
1599 :
1600 342 : CALL pw_pool%create_pw(pw_in_gs)
1601 :
1602 342 : CALL pw_transfer(pw_in, pw_in_gs)
1603 :
1604 342 : IF (have_g0) THEN
1605 171 : g0_index = green%influence_fn%pw_grid%first_gne0 - 1
1606 171 : pw_in_gs%array(g0_index) = 0.0_dp
1607 : END IF
1608 48044571 : DO ig = green%influence_fn%pw_grid%first_gne0, ng
1609 48044571 : pw_in_gs%array(ig) = prefactor*(pw_in_gs%array(ig)/green%influence_fn%array(ig))
1610 : END DO
1611 :
1612 342 : CALL pw_transfer(pw_in_gs, pw_out)
1613 :
1614 342 : CALL pw_pool%give_back_pw(pw_in_gs)
1615 :
1616 342 : CALL timestop(handle)
1617 :
1618 342 : END SUBROUTINE apply_laplace_operator_fft
1619 :
1620 : ! **************************************************************************************************
1621 : !> \brief Evaluates the action of the Laplace operator on a given 3d matrix using DCT-I
1622 : !> \param pw_pool pool of pw grid
1623 : !> \param green the greens_fn_type data holding a valid dct_influence_fn
1624 : !> \param pw_in pw_in (potential)
1625 : !> \param pw_out pw_out (density)
1626 : !> \par History
1627 : !> 07.2014 created [Hossein Bani-Hashemian]
1628 : !> 11.2015 revised [Hossein Bani-Hashemian]
1629 : !> \author Mohammad Hossein Bani-Hashemian
1630 : ! **************************************************************************************************
1631 156 : SUBROUTINE apply_laplace_operator_dct(pw_pool, green, pw_in, pw_out)
1632 :
1633 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
1634 : TYPE(greens_fn_type), INTENT(IN) :: green
1635 : TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in
1636 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pw_out
1637 :
1638 : CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_laplace_operator_dct'
1639 :
1640 : INTEGER :: g0_index, handle, ig, ng
1641 : LOGICAL :: have_g0
1642 : REAL(dp) :: prefactor
1643 : TYPE(pw_c1d_gs_type) :: pw_in_gs
1644 : TYPE(pw_grid_type), POINTER :: pw_grid
1645 :
1646 156 : CALL timeset(routineN, handle)
1647 :
1648 : ! here I multiply by fourpi to cancel out the prefactor fourpi in influence_fn
1649 156 : prefactor = fourpi
1650 :
1651 156 : pw_grid => pw_pool%pw_grid
1652 156 : ng = SIZE(pw_in%pw_grid%gsq)
1653 156 : have_g0 = green%dct_influence_fn%pw_grid%have_g0
1654 :
1655 156 : CALL pw_pool%create_pw(pw_in_gs)
1656 :
1657 156 : CALL pw_transfer(pw_in, pw_in_gs)
1658 :
1659 156 : IF (have_g0) THEN
1660 78 : g0_index = green%dct_influence_fn%pw_grid%first_gne0 - 1
1661 78 : pw_in_gs%array(g0_index) = 0.0_dp
1662 : END IF
1663 65185854 : DO ig = green%dct_influence_fn%pw_grid%first_gne0, ng
1664 65185854 : pw_in_gs%array(ig) = prefactor*(pw_in_gs%array(ig)/green%dct_influence_fn%array(ig))
1665 : END DO
1666 :
1667 156 : CALL pw_transfer(pw_in_gs, pw_out)
1668 :
1669 156 : CALL pw_pool%give_back_pw(pw_in_gs)
1670 :
1671 156 : CALL timestop(handle)
1672 :
1673 156 : END SUBROUTINE apply_laplace_operator_dct
1674 :
1675 : ! **************************************************************************************************
1676 : !> \brief Evaluates the action of the generalized Poisson operator on a given 3d matrix.
1677 : !> \param pw_pool pool of pw grid
1678 : !> \param green green functions for FFT based inverse Laplacian
1679 : !> \param dielectric dielectric environment
1680 : !> \param v potential
1681 : !> \param density density
1682 : !> \par History
1683 : !> 07.2014 created [Hossein Bani-Hashemian]
1684 : !> \author Mohammad Hossein Bani-Hashemian
1685 : ! **************************************************************************************************
1686 342 : SUBROUTINE apply_poisson_operator_fft(pw_pool, green, dielectric, v, density)
1687 :
1688 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
1689 : TYPE(greens_fn_type), INTENT(IN) :: green
1690 : TYPE(dielectric_type), INTENT(IN), POINTER :: dielectric
1691 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v
1692 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: density
1693 :
1694 : CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_poisson_operator_fft'
1695 :
1696 : INTEGER :: handle
1697 : TYPE(pw_r3d_rs_type) :: Pxv
1698 :
1699 342 : CALL timeset(routineN, handle)
1700 :
1701 342 : CALL pw_pool%create_pw(Pxv)
1702 :
1703 342 : CALL apply_P_operator(pw_pool, dielectric, v, Pxv)
1704 342 : CALL apply_laplace_operator_fft(pw_pool, green, v, density)
1705 342 : CALL pw_axpy(Pxv, density)
1706 :
1707 342 : CALL pw_pool%give_back_pw(Pxv)
1708 :
1709 342 : CALL timestop(handle)
1710 :
1711 342 : END SUBROUTINE apply_poisson_operator_fft
1712 :
1713 : ! **************************************************************************************************
1714 : !> \brief Evaluates the action of the generalized Poisson operator on a given
1715 : !> 3d matrix using DCT-I.
1716 : !> \param pw_pool pool of pw grid
1717 : !> \param green the greens_fn_type data holding a valid dct_influence_fn
1718 : !> \param dielectric dielectric environment
1719 : !> \param v potential
1720 : !> \param density density
1721 : !> \par History
1722 : !> 07.2014 created [Hossein Bani-Hashemian]
1723 : !> 11.2015 revised [Hossein Bani-Hashemian]
1724 : !> \author Mohammad Hossein Bani-Hashemian
1725 : ! **************************************************************************************************
1726 156 : SUBROUTINE apply_poisson_operator_dct(pw_pool, green, dielectric, v, density)
1727 :
1728 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
1729 : TYPE(greens_fn_type), INTENT(IN) :: green
1730 : TYPE(dielectric_type), INTENT(IN), POINTER :: dielectric
1731 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v
1732 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: density
1733 :
1734 : CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_poisson_operator_dct'
1735 :
1736 : INTEGER :: handle
1737 : TYPE(pw_r3d_rs_type) :: Pxv
1738 :
1739 156 : CALL timeset(routineN, handle)
1740 :
1741 156 : CALL pw_pool%create_pw(Pxv)
1742 :
1743 156 : CALL apply_P_operator(pw_pool, dielectric, v, Pxv)
1744 156 : CALL apply_laplace_operator_dct(pw_pool, green, v, density)
1745 156 : CALL pw_axpy(Pxv, density)
1746 :
1747 156 : CALL pw_pool%give_back_pw(Pxv)
1748 :
1749 156 : CALL timestop(handle)
1750 :
1751 156 : END SUBROUTINE apply_poisson_operator_dct
1752 :
1753 : ! **************************************************************************************************
1754 : !> \brief Computes the extra contribution (v_eps)
1755 : !> v_eps = - \frac{1}{8*\pi} * |\nabla_r(v)|^2 * \frac{d \eps}{d \rho}
1756 : !> to the functional derivative of the Hartree energy wrt the density, being
1757 : !> attributed to the dependency of the dielectric constant to the charge density.
1758 : !> [see V. M. Sanchez, M. Sued, and D. A. Scherlis, J. Chem. Phys. 131, 174108 (2009)]
1759 : !> \param pw_pool pool of the original plane-wave grid
1760 : !> \param dielectric dielectric environment
1761 : !> \param v Hartree potential
1762 : !> \param v_eps v_eps
1763 : !> \par History
1764 : !> 08.2014 created [Hossein Bani-Hashemian]
1765 : !> \author Mohammad Hossein Bani-Hashemian
1766 : ! **************************************************************************************************
1767 498 : SUBROUTINE ps_implicit_compute_veps(pw_pool, dielectric, v, v_eps)
1768 :
1769 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
1770 : TYPE(dielectric_type), INTENT(IN), POINTER :: dielectric
1771 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v
1772 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: v_eps
1773 :
1774 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ps_implicit_compute_veps'
1775 :
1776 : INTEGER :: handle, i
1777 : REAL(dp) :: eightpi
1778 : TYPE(pw_r3d_rs_type) :: dv2
1779 1992 : TYPE(pw_r3d_rs_type), DIMENSION(3) :: dv
1780 :
1781 498 : CALL timeset(routineN, handle)
1782 :
1783 498 : eightpi = 2*fourpi
1784 :
1785 498 : CALL pw_pool%create_pw(dv2)
1786 1992 : DO i = 1, 3
1787 1992 : CALL pw_pool%create_pw(dv(i))
1788 : END DO
1789 :
1790 498 : CALL derive_fft(v, dv, pw_pool)
1791 :
1792 : ! evaluate |\nabla_r(v)|^2
1793 115979154 : dv2%array = dv(1)%array**2 + dv(2)%array**2 + dv(3)%array**2
1794 :
1795 115979154 : v_eps%array = -(1.0_dp/eightpi)*(dv2%array*dielectric%deps_drho%array)
1796 :
1797 498 : CALL pw_pool%give_back_pw(dv2)
1798 1992 : DO i = 1, 3
1799 1992 : CALL pw_pool%give_back_pw(dv(i))
1800 : END DO
1801 :
1802 498 : CALL timestop(handle)
1803 :
1804 498 : END SUBROUTINE ps_implicit_compute_veps
1805 :
1806 : ! **************************************************************************************************
1807 : !> \brief Computes the Hartree energy
1808 : !> \param density electronic density
1809 : !> \param v Hartree potential
1810 : !> \param ehartree Hartree energy
1811 : !> \par History
1812 : !> 06.2015 created [Hossein Bani-Hashemian]
1813 : !> \author Mohammad Hossein Bani-Hashemian
1814 : ! **************************************************************************************************
1815 784 : SUBROUTINE compute_ehartree_periodic_bc(density, v, ehartree)
1816 :
1817 : TYPE(pw_r3d_rs_type), INTENT(IN) :: density, v
1818 : REAL(dp), INTENT(OUT) :: ehartree
1819 :
1820 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_ehartree_periodic_bc'
1821 :
1822 : INTEGER :: handle
1823 :
1824 784 : CALL timeset(routineN, handle)
1825 :
1826 : ! E_H = \frac{1}{2} * \int \rho * v dr
1827 784 : ehartree = 0.5_dp*pw_integral_ab(density, v)
1828 :
1829 784 : CALL timestop(handle)
1830 :
1831 784 : END SUBROUTINE compute_ehartree_periodic_bc
1832 :
1833 : ! **************************************************************************************************
1834 : !> \brief Computes the Hartree energy
1835 : !> \param dielectric dielectric environment
1836 : !> \param density electronic density
1837 : !> \param Btxlambda B^t * \lambda (\lambda is the vector of Lagrange multipliers
1838 : !> and B^t is the transpose of the boundary operator
1839 : !> \param v Hartree potential
1840 : !> \param ehartree Hartree energy
1841 : !> \param electric_enthalpy electric enthalpy
1842 : !> \par History
1843 : !> 06.2015 created [Hossein Bani-Hashemian]
1844 : !> \author Mohammad Hossein Bani-Hashemian
1845 : ! **************************************************************************************************
1846 1738 : SUBROUTINE compute_ehartree_mixed_bc(dielectric, density, Btxlambda, v, ehartree, electric_enthalpy)
1847 :
1848 : TYPE(dielectric_type), INTENT(IN), POINTER :: dielectric
1849 : TYPE(pw_r3d_rs_type), INTENT(IN) :: density
1850 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :), &
1851 : INTENT(IN) :: Btxlambda
1852 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v
1853 : REAL(dp), INTENT(OUT) :: ehartree, electric_enthalpy
1854 :
1855 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_ehartree_mixed_bc'
1856 :
1857 : INTEGER :: handle
1858 : REAL(dp) :: dvol, ehartree_rho, ehartree_rho_cstr
1859 : TYPE(pw_grid_type), POINTER :: pw_grid
1860 :
1861 1738 : CALL timeset(routineN, handle)
1862 :
1863 1738 : pw_grid => v%pw_grid
1864 :
1865 1738 : dvol = pw_grid%dvol
1866 :
1867 : ! E_H = \frac{1}{2} * \int \rho * v dr + \frac{1}{8 \pi} * \int Btxlambda * v dr
1868 : ! the sign of the second term depends on the sign chosen for the Lagrange multiplier
1869 : ! term in the variational form
1870 287894894 : ehartree_rho = accurate_sum(density%array*v%array)
1871 287894894 : ehartree_rho_cstr = accurate_sum(dielectric%eps%array*Btxlambda*v%array/fourpi)
1872 1738 : ehartree_rho = 0.5_dp*ehartree_rho*dvol
1873 1738 : ehartree_rho_cstr = 0.5_dp*ehartree_rho_cstr*dvol
1874 1738 : CALL pw_grid%para%group%sum(ehartree_rho)
1875 1738 : CALL pw_grid%para%group%sum(ehartree_rho_cstr)
1876 1738 : electric_enthalpy = ehartree_rho + ehartree_rho_cstr
1877 1738 : ehartree = ehartree_rho - ehartree_rho_cstr
1878 :
1879 1738 : CALL timestop(handle)
1880 :
1881 1738 : END SUBROUTINE compute_ehartree_mixed_bc
1882 :
1883 : ! **************************************************************************************************
1884 : !> \brief Computes the (normalized) preconditioned residual norm error and the
1885 : !> normalized absolute error
1886 : !> \param pw_pool pool of the original plane-wave grid
1887 : !> \param green greens functions for FFT based inverse Laplacian
1888 : !> \param res_new residual
1889 : !> \param v_old old v
1890 : !> \param v_new new v
1891 : !> \param QAinvxres_new Delta^-1(res_new)
1892 : !> \param pres_error preconditioned residual norm error
1893 : !> \param nabs_error normalized absolute error
1894 : !> \par History
1895 : !> 07.2014 created [Hossein Bani-Hashemian]
1896 : !> \author Mohammad Hossein Bani-Hashemian
1897 : ! **************************************************************************************************
1898 2254 : SUBROUTINE ps_implicit_compute_error_fft(pw_pool, green, res_new, v_old, v_new, &
1899 : QAinvxres_new, pres_error, nabs_error)
1900 :
1901 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
1902 : TYPE(greens_fn_type), INTENT(IN) :: green
1903 : TYPE(pw_r3d_rs_type), INTENT(IN) :: res_new, v_old, v_new
1904 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: QAinvxres_new
1905 : REAL(dp), INTENT(OUT) :: pres_error, nabs_error
1906 :
1907 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ps_implicit_compute_error_fft'
1908 :
1909 : INTEGER :: handle
1910 : REAL(dp) :: vol
1911 :
1912 2254 : CALL timeset(routineN, handle)
1913 :
1914 2254 : vol = pw_pool%pw_grid%vol
1915 :
1916 : ! evaluate \Delta^-1(res) = \Delta^-1 (g - \Delta(v_new) - P(v_new) + Bt \lambda)
1917 2254 : CALL apply_inv_laplace_operator_fft(pw_pool, green, res_new, QAinvxres_new)
1918 : ! (normalized) preconditioned residual norm error :
1919 327225790 : pres_error = accurate_sum(QAinvxres_new%array(:, :, :)**2)
1920 2254 : CALL pw_pool%pw_grid%para%group%sum(pres_error)
1921 2254 : pres_error = SQRT(pres_error)/vol
1922 :
1923 : ! normalized absolute error :
1924 : ! nabs_error := \frac{\| v_old - v_new \|}{volume}
1925 327225790 : nabs_error = accurate_sum(ABS(v_old%array - v_new%array)**2)
1926 2254 : CALL pw_pool%pw_grid%para%group%sum(nabs_error)
1927 2254 : nabs_error = SQRT(nabs_error)/vol
1928 :
1929 2254 : CALL timestop(handle)
1930 :
1931 2254 : END SUBROUTINE ps_implicit_compute_error_fft
1932 :
1933 : ! **************************************************************************************************
1934 : !> \brief Computes the (normalized) preconditioned residual norm error and the
1935 : !> normalized absolute error
1936 : !> \param pw_pool pool of the original plane-wave grid
1937 : !> \param green the greens_fn_type data holding a valid dct_influence_fn
1938 : !> \param res_new residual
1939 : !> \param v_old old v
1940 : !> \param v_new new v
1941 : !> \param QAinvxres_new Delta^-1(res_new)
1942 : !> \param pres_error preconditioned residual norm error
1943 : !> \param nabs_error normalized absolute error
1944 : !> \par History
1945 : !> 07.2014 created [Hossein Bani-Hashemian]
1946 : !> 11.2015 revised [Hossein Bani-Hashemian]
1947 : !> \author Mohammad Hossein Bani-Hashemian
1948 : ! **************************************************************************************************
1949 268 : SUBROUTINE ps_implicit_compute_error_dct(pw_pool, green, res_new, v_old, v_new, &
1950 : QAinvxres_new, pres_error, nabs_error)
1951 :
1952 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
1953 : TYPE(greens_fn_type), INTENT(IN) :: green
1954 : TYPE(pw_r3d_rs_type), INTENT(IN) :: res_new, v_old, v_new
1955 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: QAinvxres_new
1956 : REAL(dp), INTENT(OUT) :: pres_error, nabs_error
1957 :
1958 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ps_implicit_compute_error_dct'
1959 :
1960 : INTEGER :: handle
1961 : REAL(dp) :: vol
1962 :
1963 268 : CALL timeset(routineN, handle)
1964 :
1965 268 : vol = pw_pool%pw_grid%vol
1966 :
1967 : ! evaluate \Delta^-1(res) = \Delta^-1 (g - \Delta(v_new) - P(v_new) + Bt \lambda)
1968 268 : CALL apply_inv_laplace_operator_dct(pw_pool, green, res_new, QAinvxres_new)
1969 : ! (normalized) preconditioned residual norm error :
1970 119766668 : pres_error = accurate_sum(QAinvxres_new%array(:, :, :)**2)
1971 268 : CALL pw_pool%pw_grid%para%group%sum(pres_error)
1972 268 : pres_error = SQRT(pres_error)/vol
1973 :
1974 : ! normalized absolute error :
1975 : ! nabs_error := \frac{\| v_old - v_new \|}{volume}
1976 119766668 : nabs_error = accurate_sum(ABS(v_old%array - v_new%array)**2)
1977 268 : CALL pw_pool%pw_grid%para%group%sum(nabs_error)
1978 268 : nabs_error = SQRT(nabs_error)/vol
1979 :
1980 268 : CALL timestop(handle)
1981 :
1982 268 : END SUBROUTINE ps_implicit_compute_error_dct
1983 :
1984 : ! **************************************************************************************************
1985 : !> \brief output of the implicit (iterative) Poisson solver
1986 : !> \param iter current iteration
1987 : !> \param pres_error preconditioned residual norm error
1988 : !> \param nabs_error normalized absolute error
1989 : !> \param outp_unit output unit
1990 : !> \par History
1991 : !> 07.2014 created [Hossein Bani-Hashemian]
1992 : !> \author Mohammad Hossein Bani-Hashemian
1993 : ! **************************************************************************************************
1994 2522 : SUBROUTINE ps_implicit_output(iter, pres_error, nabs_error, outp_unit)
1995 :
1996 : INTEGER, INTENT(IN) :: iter
1997 : REAL(dp), INTENT(IN) :: pres_error, nabs_error
1998 : INTEGER, INTENT(OUT) :: outp_unit
1999 :
2000 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ps_implicit_output'
2001 : INTEGER, PARAMETER :: low_print_level = 1
2002 :
2003 : INTEGER :: handle
2004 : TYPE(cp_logger_type), POINTER :: logger
2005 :
2006 2522 : CALL timeset(routineN, handle)
2007 :
2008 2522 : logger => cp_get_default_logger()
2009 2522 : IF (logger%para_env%is_source()) THEN
2010 1261 : outp_unit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
2011 : ELSE
2012 1261 : outp_unit = -1
2013 : END IF
2014 :
2015 2522 : IF (logger%iter_info%print_level > low_print_level) THEN
2016 2476 : IF ((outp_unit > 0) .AND. (iter == 1)) THEN
2017 : WRITE (outp_unit, '(T3,A)') &
2018 226 : "POISSON| iter pres error nabs error E_hartree delta E"
2019 : END IF
2020 :
2021 2476 : IF (outp_unit > 0) THEN
2022 : WRITE (outp_unit, '(T3,A,I6,5X,E13.4,3X,E13.4)', ADVANCE='NO') &
2023 1238 : "POISSON| ", iter, pres_error, nabs_error
2024 : END IF
2025 : END IF
2026 :
2027 2522 : CALL timestop(handle)
2028 :
2029 2522 : END SUBROUTINE ps_implicit_output
2030 :
2031 : ! **************************************************************************************************
2032 : !> \brief reports the Hartree energy in every iteration
2033 : !> \param ps_implicit_env the implicit poisson solver environment
2034 : !> \param outp_unit output unit
2035 : !> \param ehartree Hartree energy
2036 : !> \par History
2037 : !> 07.2014 created [Hossein Bani-Hashemian]
2038 : !> \author Mohammad Hossein Bani-Hashemian
2039 : ! **************************************************************************************************
2040 5044 : SUBROUTINE ps_implicit_report_ehartree(ps_implicit_env, outp_unit, ehartree)
2041 :
2042 : TYPE(ps_implicit_type) :: ps_implicit_env
2043 : INTEGER, INTENT(IN) :: outp_unit
2044 : REAL(dp), INTENT(IN) :: ehartree
2045 :
2046 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ps_implicit_report_ehartree'
2047 : INTEGER, PARAMETER :: low_print_level = 1
2048 :
2049 : INTEGER :: handle
2050 : TYPE(cp_logger_type), POINTER :: logger
2051 :
2052 2522 : logger => cp_get_default_logger()
2053 2522 : CALL timeset(routineN, handle)
2054 2522 : IF (logger%iter_info%print_level > low_print_level) THEN
2055 2476 : IF (outp_unit > 0) WRITE (outp_unit, '(F19.10,E10.2)') &
2056 1238 : ehartree, ehartree - ps_implicit_env%ehartree
2057 : END IF
2058 2522 : CALL timestop(handle)
2059 :
2060 2522 : END SUBROUTINE ps_implicit_report_ehartree
2061 :
2062 : ! **************************************************************************************************
2063 : !> \brief reports the final number of iteration
2064 : !> \param iter the iteration number after exiting the main loop
2065 : !> \param max_iter maximum number of iterations
2066 : !> \param outp_unit output unit
2067 : !> \par History
2068 : !> 02.2016 created [Hossein Bani-Hashemian]
2069 : !> \author Mohammad Hossein Bani-Hashemian
2070 : ! **************************************************************************************************
2071 498 : SUBROUTINE ps_implicit_print_convergence_msg(iter, max_iter, outp_unit)
2072 :
2073 : INTEGER, INTENT(IN) :: iter, max_iter, outp_unit
2074 :
2075 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ps_implicit_print_convergence_msg'
2076 :
2077 : CHARACTER(LEN=12) :: msg
2078 : INTEGER :: handle, last_iter
2079 :
2080 498 : CALL timeset(routineN, handle)
2081 :
2082 498 : last_iter = iter - 1
2083 :
2084 498 : IF (outp_unit > 0) THEN
2085 249 : IF (last_iter == max_iter) THEN
2086 : WRITE (outp_unit, '(T3,A)') &
2087 0 : "POISSON| No convergence achieved within the maximum number of iterations."
2088 : END IF
2089 249 : IF (last_iter < max_iter) THEN
2090 249 : IF (last_iter == 1) THEN
2091 117 : msg = " iteration."
2092 : ELSE
2093 132 : msg = " iterations."
2094 : END IF
2095 : WRITE (outp_unit, '(T3,A,I0,A)') &
2096 249 : "POISSON| Poisson solver converged in ", last_iter, msg
2097 : END IF
2098 : END IF
2099 498 : CALL timestop(handle)
2100 :
2101 498 : END SUBROUTINE ps_implicit_print_convergence_msg
2102 :
2103 : ! **************************************************************************************************
2104 : !> \brief converts a 1D array to a 3D array (contiguous layout)
2105 : !> \param idx_1dto3d mapping of indices
2106 : !> \param arr1d input 1D array
2107 : !> \param arr3d input 3D array
2108 : ! **************************************************************************************************
2109 2062 : SUBROUTINE convert_1dto3d(idx_1dto3d, arr1d, arr3d)
2110 :
2111 : INTEGER, DIMENSION(:), INTENT(IN) :: idx_1dto3d
2112 : REAL(dp), DIMENSION(:), INTENT(IN) :: arr1d
2113 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :), &
2114 : INTENT(INOUT) :: arr3d
2115 :
2116 : CHARACTER(LEN=*), PARAMETER :: routineN = 'convert_1dto3d'
2117 :
2118 : INTEGER :: handle, i, j, k, l, lb1, lb2, lb3, &
2119 : npts1, npts2, npts3, ub1, ub2, ub3
2120 :
2121 2062 : CALL timeset(routineN, handle)
2122 :
2123 2062 : lb1 = LBOUND(arr3d, 1); ub1 = UBOUND(arr3d, 1)
2124 2062 : lb2 = LBOUND(arr3d, 2); ub2 = UBOUND(arr3d, 2)
2125 2062 : lb3 = LBOUND(arr3d, 3); ub3 = UBOUND(arr3d, 3)
2126 :
2127 2062 : npts1 = ub1 - lb1 + 1
2128 2062 : npts2 = ub2 - lb2 + 1
2129 2062 : npts3 = ub3 - lb3 + 1
2130 :
2131 363318274 : DO l = 1, SIZE(idx_1dto3d)
2132 363316212 : k = ((idx_1dto3d(l) - 1)/(npts1*npts2)) + lb3
2133 363316212 : j = ((idx_1dto3d(l) - 1) - (k - lb3)*npts1*npts2)/npts1 + lb2
2134 363316212 : i = idx_1dto3d(l) - ((j - lb2)*npts1 + (k - lb3)*npts1*npts2) + lb1 - 1
2135 363318274 : arr3d(i, j, k) = arr1d(l)
2136 : END DO
2137 :
2138 2062 : CALL timestop(handle)
2139 :
2140 2062 : END SUBROUTINE convert_1dto3d
2141 :
2142 : ! **************************************************************************************************
2143 : !> \brief Returns the voltage of a tile. In case an alternating field is used, the oltage is a function of time
2144 : !> \param time ...
2145 : !> \param v_D ...
2146 : !> \param osc_frac ...
2147 : !> \param frequency ...
2148 : !> \param phase ...
2149 : !> \param v_D_new ...
2150 : ! **************************************************************************************************
2151 324 : SUBROUTINE get_voltage(time, v_D, osc_frac, frequency, phase, v_D_new)
2152 :
2153 : REAL(dp), INTENT(IN) :: time
2154 : REAL(dp), DIMENSION(:), INTENT(IN) :: v_D, osc_frac, frequency, phase
2155 : REAL(dp), ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: v_D_new
2156 :
2157 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_voltage'
2158 :
2159 : INTEGER :: handle, i
2160 :
2161 324 : CALL timeset(routineN, handle)
2162 :
2163 972 : ALLOCATE (v_D_new(SIZE(v_D)))
2164 :
2165 2256 : DO i = 1, SIZE(v_D)
2166 : v_D_new(i) = v_D(i)*(1 - osc_frac(i)) + &
2167 2256 : v_D(i)*osc_frac(i)*COS(2*pi*time*frequency(i) + phase(i))
2168 : END DO
2169 :
2170 324 : CALL timestop(handle)
2171 :
2172 324 : END SUBROUTINE get_voltage
2173 :
2174 : END MODULE ps_implicit_methods
|
