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