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 : !> \par History
10 : !> JGH (15-Mar-2001) : Update small grid when cell changes
11 : !> with dg_grid_change
12 : ! **************************************************************************************************
13 : MODULE dgs
14 :
15 : USE fft_tools, ONLY: BWFFT,&
16 : FFT_RADIX_CLOSEST,&
17 : FFT_RADIX_NEXT_ODD,&
18 : fft3d,&
19 : fft_radix_operations
20 : USE kinds, ONLY: dp
21 : USE mathconstants, ONLY: z_one,&
22 : z_zero
23 : USE mathlib, ONLY: det_3x3,&
24 : inv_3x3
25 : USE message_passing, ONLY: mp_comm_self,&
26 : mp_comm_type
27 : USE pw_grid_info, ONLY: pw_find_cutoff
28 : USE pw_grid_types, ONLY: HALFSPACE,&
29 : pw_grid_type
30 : USE pw_grids, ONLY: pw_grid_change,&
31 : pw_grid_create
32 : USE pw_methods, ONLY: pw_copy,&
33 : pw_multiply_with,&
34 : pw_zero
35 : USE pw_types, ONLY: pw_c3d_rs_type,&
36 : pw_r3d_rs_type
37 : USE realspace_grid_types, ONLY: realspace_grid_type
38 : USE structure_factors, ONLY: structure_factor_evaluate
39 : #include "../base/base_uses.f90"
40 :
41 : IMPLICIT NONE
42 :
43 : PRIVATE
44 :
45 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dgs'
46 :
47 : PUBLIC :: dg_get_patch
48 : PUBLIC :: dg_pme_grid_setup, &
49 : dg_sum_patch, dg_sum_patch_force_3d, dg_sum_patch_force_1d, &
50 : dg_get_strucfac, dg_grid_change
51 :
52 : INTERFACE dg_sum_patch
53 : MODULE PROCEDURE dg_sum_patch_coef, dg_sum_patch_arr
54 : END INTERFACE
55 :
56 : INTERFACE dg_sum_patch_force_3d
57 : MODULE PROCEDURE dg_sum_patch_force_coef_3d, dg_sum_patch_force_arr_3d
58 : END INTERFACE
59 :
60 : INTERFACE dg_sum_patch_force_1d
61 : MODULE PROCEDURE dg_sum_patch_force_coef_1d, dg_sum_patch_force_arr_1d
62 : END INTERFACE
63 :
64 : INTERFACE dg_get_patch
65 : MODULE PROCEDURE dg_get_patch_1, dg_get_patch_2
66 : END INTERFACE
67 :
68 : INTERFACE dg_add_patch
69 : MODULE PROCEDURE dg_add_patch_simple, dg_add_patch_folded
70 : END INTERFACE
71 :
72 : INTERFACE dg_int_patch_3d
73 : MODULE PROCEDURE dg_int_patch_simple_3d, dg_int_patch_folded_3d
74 : END INTERFACE
75 :
76 : INTERFACE dg_int_patch_1d
77 : MODULE PROCEDURE dg_int_patch_simple_1d, dg_int_patch_folded_1d
78 : END INTERFACE
79 :
80 : CONTAINS
81 :
82 : ! **************************************************************************************************
83 : !> \brief ...
84 : !> \param b_cell_hmat ...
85 : !> \param npts_s ...
86 : !> \param cutoff_radius ...
87 : !> \param grid_s ...
88 : !> \param grid_b ...
89 : !> \param mp_comm ...
90 : !> \param grid_ref ...
91 : !> \param rs_dims ...
92 : !> \param iounit ...
93 : !> \param fft_usage ...
94 : ! **************************************************************************************************
95 86 : SUBROUTINE dg_pme_grid_setup(b_cell_hmat, npts_s, cutoff_radius, grid_s, grid_b, &
96 : mp_comm, grid_ref, rs_dims, iounit, fft_usage)
97 :
98 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: b_cell_hmat
99 : INTEGER, DIMENSION(:), INTENT(IN) :: npts_s
100 : REAL(KIND=dp), INTENT(IN) :: cutoff_radius
101 : TYPE(pw_grid_type), POINTER :: grid_s, grid_b
102 :
103 : CLASS(mp_comm_type), INTENT(IN) :: mp_comm
104 : TYPE(pw_grid_type), INTENT(IN), OPTIONAL :: grid_ref
105 : INTEGER, DIMENSION(2), INTENT(in), OPTIONAL :: rs_dims
106 : INTEGER, INTENT(IN), OPTIONAL :: iounit
107 : LOGICAL, INTENT(IN), OPTIONAL :: fft_usage
108 :
109 : INTEGER, DIMENSION(2, 3) :: bounds_b, bounds_s
110 : REAL(KIND=dp) :: cutoff, ecut
111 : REAL(KIND=dp), DIMENSION(3, 3) :: s_cell_hmat, unit_cell_hmat
112 :
113 86 : CALL dg_find_cutoff(b_cell_hmat, npts_s, cutoff_radius, bounds_s, bounds_b, cutoff)
114 :
115 86 : ecut = 0.5_dp*cutoff*cutoff
116 86 : IF (PRESENT(grid_ref)) THEN
117 : CALL pw_grid_create(grid_b, mp_comm, b_cell_hmat, bounds=bounds_b, grid_span=HALFSPACE, cutoff=ecut, spherical=.TRUE., &
118 0 : ref_grid=grid_ref, rs_dims=rs_dims, iounit=iounit, fft_usage=fft_usage)
119 : ELSE
120 : CALL pw_grid_create(grid_b, mp_comm, b_cell_hmat, bounds=bounds_b, grid_span=HALFSPACE, cutoff=ecut, spherical=.TRUE., &
121 86 : rs_dims=rs_dims, iounit=iounit, fft_usage=fft_usage)
122 : END IF
123 :
124 86 : CALL dg_find_basis(grid_b%npts, b_cell_hmat, unit_cell_hmat)
125 :
126 344 : CALL dg_set_cell(bounds_s(2, :) - bounds_s(1, :) + 1, unit_cell_hmat, s_cell_hmat)
127 :
128 : CALL pw_grid_create(grid_s, mp_comm_self, s_cell_hmat, bounds=bounds_s, grid_span=HALFSPACE, &
129 86 : cutoff=ecut, iounit=iounit, fft_usage=fft_usage)
130 :
131 86 : END SUBROUTINE dg_pme_grid_setup
132 :
133 : ! **************************************************************************************************
134 : !> \brief Calculate the lengths of the cell vectors a, b, and c
135 : !> \param cell_hmat ...
136 : !> \return ...
137 : !> \par History 01.2014 copied from cell_types::get_cell()
138 : !> \author Ole Schuett
139 : ! **************************************************************************************************
140 86 : PURE FUNCTION get_cell_lengths(cell_hmat) RESULT(abc)
141 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: cell_hmat
142 : REAL(KIND=dp), DIMENSION(3) :: abc
143 :
144 : abc(1) = SQRT(cell_hmat(1, 1)*cell_hmat(1, 1) + &
145 : cell_hmat(2, 1)*cell_hmat(2, 1) + &
146 86 : cell_hmat(3, 1)*cell_hmat(3, 1))
147 : abc(2) = SQRT(cell_hmat(1, 2)*cell_hmat(1, 2) + &
148 : cell_hmat(2, 2)*cell_hmat(2, 2) + &
149 86 : cell_hmat(3, 2)*cell_hmat(3, 2))
150 : abc(3) = SQRT(cell_hmat(1, 3)*cell_hmat(1, 3) + &
151 : cell_hmat(2, 3)*cell_hmat(2, 3) + &
152 86 : cell_hmat(3, 3)*cell_hmat(3, 3))
153 86 : END FUNCTION get_cell_lengths
154 :
155 : ! **************************************************************************************************
156 : !> \brief ...
157 : !> \param b_cell_hmat ...
158 : !> \param npts_s ...
159 : !> \param cutoff_radius ...
160 : !> \param bounds_s ...
161 : !> \param bounds_b ...
162 : !> \param cutoff ...
163 : ! **************************************************************************************************
164 86 : SUBROUTINE dg_find_cutoff(b_cell_hmat, npts_s, cutoff_radius, &
165 : bounds_s, bounds_b, cutoff)
166 :
167 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: b_cell_hmat
168 : INTEGER, DIMENSION(:), INTENT(IN) :: npts_s
169 : REAL(KIND=dp), INTENT(IN) :: cutoff_radius
170 : INTEGER, DIMENSION(2, 3), INTENT(OUT) :: bounds_s, bounds_b
171 : REAL(KIND=dp), INTENT(OUT) :: cutoff
172 :
173 : INTEGER, DIMENSION(3) :: nout, npts_b
174 : REAL(KIND=dp) :: b_cell_deth, cell_lengths(3), dr(3)
175 : REAL(KIND=dp), DIMENSION(3, 3) :: b_cell_h_inv
176 :
177 86 : b_cell_deth = ABS(det_3x3(b_cell_hmat))
178 86 : IF (b_cell_deth < 1.0E-10_dp) THEN
179 : CALL cp_abort(__LOCATION__, &
180 : "An invalid set of cell vectors was specified. "// &
181 0 : "The determinant det(h) is too small")
182 : END IF
183 86 : b_cell_h_inv = inv_3x3(b_cell_hmat)
184 :
185 : CALL fft_radix_operations(npts_s(1), nout(1), &
186 86 : operation=FFT_RADIX_NEXT_ODD)
187 : CALL fft_radix_operations(npts_s(1), nout(2), &
188 86 : operation=FFT_RADIX_NEXT_ODD)
189 : CALL fft_radix_operations(npts_s(1), nout(3), &
190 86 : operation=FFT_RADIX_NEXT_ODD)
191 :
192 86 : cell_lengths = get_cell_lengths(b_cell_hmat)
193 :
194 86 : CALL dg_get_spacing(nout, cutoff_radius, dr)
195 86 : CALL dg_find_radix(dr, cell_lengths, npts_b)
196 :
197 : ! In-line code to set grid_b % npts = npts_s if necessary
198 86 : IF (nout(1) > npts_b(1)) THEN
199 4 : npts_b(1) = nout(1)
200 : dr(1) = cell_lengths(1)/REAL(nout(1), KIND=dp)
201 : END IF
202 86 : IF (nout(2) > npts_b(2)) THEN
203 4 : npts_b(2) = nout(2)
204 : dr(2) = cell_lengths(2)/REAL(nout(2), KIND=dp)
205 : END IF
206 86 : IF (nout(3) > npts_b(3)) THEN
207 4 : npts_b(3) = nout(3)
208 : dr(3) = cell_lengths(3)/REAL(nout(3), KIND=dp)
209 : END IF
210 :
211 : ! big grid bounds
212 344 : bounds_b(1, :) = -npts_b/2
213 344 : bounds_b(2, :) = +(npts_b - 1)/2
214 : ! small grid bounds
215 344 : bounds_s(1, :) = -nout(:)/2
216 344 : bounds_s(2, :) = (+nout(:) - 1)/2
217 :
218 86 : cutoff = pw_find_cutoff(npts_b, b_cell_h_inv)
219 :
220 86 : END SUBROUTINE dg_find_cutoff
221 :
222 : ! **************************************************************************************************
223 : !> \brief ...
224 : !> \param npts ...
225 : !> \param cutoff_radius ...
226 : !> \param dr ...
227 : ! **************************************************************************************************
228 86 : SUBROUTINE dg_get_spacing(npts, cutoff_radius, dr)
229 :
230 : INTEGER, DIMENSION(:), INTENT(IN) :: npts
231 : REAL(KIND=dp), INTENT(IN) :: cutoff_radius
232 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: dr
233 :
234 344 : dr(:) = cutoff_radius/(REAL(npts(:), KIND=dp)/2.0_dp)
235 :
236 86 : END SUBROUTINE dg_get_spacing
237 :
238 : ! **************************************************************************************************
239 : !> \brief ...
240 : !> \param b_cell_hmat ...
241 : !> \param grid_b ...
242 : !> \param grid_s ...
243 : ! **************************************************************************************************
244 86 : SUBROUTINE dg_grid_change(b_cell_hmat, grid_b, grid_s)
245 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: b_cell_hmat
246 : TYPE(pw_grid_type), POINTER :: grid_b, grid_s
247 :
248 : REAL(KIND=dp), DIMENSION(3, 3) :: s_cell_hmat, unit_cell_hmat
249 :
250 86 : CALL dg_find_basis(grid_b%npts, b_cell_hmat, unit_cell_hmat)
251 86 : CALL dg_set_cell(grid_s%npts, unit_cell_hmat, s_cell_hmat)
252 86 : CALL pw_grid_change(s_cell_hmat, grid_s)
253 86 : END SUBROUTINE dg_grid_change
254 :
255 : ! **************************************************************************************************
256 : !> \brief ...
257 : !> \param dr ...
258 : !> \param cell_lengths ...
259 : !> \param npts ...
260 : ! **************************************************************************************************
261 86 : SUBROUTINE dg_find_radix(dr, cell_lengths, npts)
262 :
263 : REAL(KIND=dp), INTENT(INOUT) :: dr(3)
264 : REAL(KIND=dp), INTENT(IN) :: cell_lengths(3)
265 : INTEGER, DIMENSION(:), INTENT(OUT) :: npts
266 :
267 : INTEGER, DIMENSION(3) :: nin
268 :
269 344 : nin(:) = NINT(cell_lengths(:)/dr(:))
270 : CALL fft_radix_operations(nin(1), npts(1), &
271 86 : operation=FFT_RADIX_CLOSEST)
272 : CALL fft_radix_operations(nin(2), npts(2), &
273 86 : operation=FFT_RADIX_CLOSEST)
274 : CALL fft_radix_operations(nin(3), npts(3), &
275 86 : operation=FFT_RADIX_CLOSEST)
276 344 : dr(:) = cell_lengths(:)/REAL(npts(:), KIND=dp)
277 :
278 86 : END SUBROUTINE dg_find_radix
279 :
280 : ! **************************************************************************************************
281 : !> \brief ...
282 : !> \param npts ...
283 : !> \param cell_hmat ...
284 : !> \param unit_cell_hmat ...
285 : ! **************************************************************************************************
286 172 : PURE SUBROUTINE dg_find_basis(npts, cell_hmat, unit_cell_hmat)
287 : INTEGER, DIMENSION(:), INTENT(IN) :: npts
288 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: cell_hmat
289 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: unit_cell_hmat
290 :
291 : INTEGER :: i
292 :
293 688 : DO i = 1, 3
294 2236 : unit_cell_hmat(:, i) = cell_hmat(:, i)/REAL(npts(:), KIND=dp)
295 : END DO
296 :
297 172 : END SUBROUTINE dg_find_basis
298 :
299 : !! Calculation of the basis on the mesh 'box'
300 :
301 : ! **************************************************************************************************
302 : !> \brief ...
303 : !> \param npts ...
304 : !> \param unit_cell_hmat ...
305 : !> \param cell_hmat ...
306 : ! **************************************************************************************************
307 172 : PURE SUBROUTINE dg_set_cell(npts, unit_cell_hmat, cell_hmat)
308 : INTEGER, DIMENSION(:), INTENT(IN) :: npts
309 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: unit_cell_hmat
310 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: cell_hmat
311 :
312 : ! computing the unit vector along a, b, c and scaling it to length dr:
313 :
314 688 : cell_hmat(:, 1) = unit_cell_hmat(:, 1)*npts(1)
315 688 : cell_hmat(:, 2) = unit_cell_hmat(:, 2)*npts(2)
316 688 : cell_hmat(:, 3) = unit_cell_hmat(:, 3)*npts(3)
317 :
318 172 : END SUBROUTINE dg_set_cell
319 :
320 : ! **************************************************************************************************
321 : !> \brief ...
322 : !> \param cell_hmat ...
323 : !> \param r ...
324 : !> \param npts_s ...
325 : !> \param npts_b ...
326 : !> \param centre ...
327 : !> \param lb ...
328 : !> \param ex ...
329 : !> \param ey ...
330 : !> \param ez ...
331 : ! **************************************************************************************************
332 24942 : SUBROUTINE dg_get_strucfac(cell_hmat, r, npts_s, npts_b, centre, lb, ex, ey, ez)
333 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: cell_hmat
334 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: r
335 : INTEGER, DIMENSION(:), INTENT(IN) :: npts_s, npts_b
336 : INTEGER, INTENT(OUT) :: centre(3)
337 : INTEGER, INTENT(IN) :: lb(3)
338 : COMPLEX(KIND=dp), DIMENSION(lb(1):), INTENT(OUT) :: ex
339 : COMPLEX(KIND=dp), DIMENSION(lb(2):), INTENT(OUT) :: ey
340 : COMPLEX(KIND=dp), DIMENSION(lb(3):), INTENT(OUT) :: ez
341 :
342 : REAL(KIND=dp) :: delta(3)
343 :
344 24942 : CALL dg_get_delta(cell_hmat, r, npts_s, npts_b, centre, delta)
345 :
346 24942 : CALL structure_factor_evaluate(delta, lb, ex, ey, ez)
347 :
348 24942 : END SUBROUTINE dg_get_strucfac
349 :
350 : ! **************************************************************************************************
351 : !> \brief ...
352 : !> \param cell_hmat ...
353 : !> \param r ...
354 : !> \param npts_s ...
355 : !> \param npts_b ...
356 : !> \param centre ...
357 : !> \param delta ...
358 : ! **************************************************************************************************
359 24942 : SUBROUTINE dg_get_delta(cell_hmat, r, npts_s, npts_b, centre, delta)
360 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: cell_hmat
361 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: r
362 : INTEGER, DIMENSION(:), INTENT(IN) :: npts_s, npts_b
363 : INTEGER, DIMENSION(:), INTENT(OUT) :: centre
364 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: delta
365 :
366 : REAL(KIND=dp) :: cell_deth
367 : REAL(KIND=dp), DIMENSION(3) :: grid_i, s
368 : REAL(KIND=dp), DIMENSION(3, 3) :: cell_h_inv
369 :
370 24942 : cell_deth = ABS(det_3x3(cell_hmat))
371 24942 : IF (cell_deth < 1.0E-10_dp) THEN
372 : CALL cp_abort(__LOCATION__, &
373 : "An invalid set of cell vectors was specified. "// &
374 0 : "The determinant det(h) is too small")
375 : END IF
376 :
377 24942 : cell_h_inv = inv_3x3(cell_hmat)
378 :
379 : ! compute the scaled coordinate of atomi
380 324246 : s = MATMUL(cell_h_inv, r)
381 99768 : s = s - NINT(s)
382 :
383 : ! find the continuous ``grid'' point (on big grid)
384 99768 : grid_i(1:3) = REAL(npts_b(1:3), KIND=dp)*s(1:3)
385 :
386 : ! find the closest grid point (on big grid)
387 99768 : centre(:) = NINT(grid_i(:))
388 :
389 : ! find the distance vector
390 99768 : delta(:) = (grid_i(:) - centre(:))/REAL(npts_s(:), KIND=dp)
391 :
392 99768 : centre(:) = centre(:) + npts_b(:)/2
393 99768 : centre(:) = MODULO(centre(:), npts_b(:))
394 99768 : centre(:) = centre(:) - npts_b(:)/2
395 :
396 24942 : END SUBROUTINE dg_get_delta
397 :
398 : ! **************************************************************************************************
399 : !> \brief ...
400 : !> \param rs ...
401 : !> \param rhos ...
402 : !> \param center ...
403 : ! **************************************************************************************************
404 24639 : PURE SUBROUTINE dg_sum_patch_coef(rs, rhos, center)
405 :
406 : TYPE(realspace_grid_type), INTENT(INOUT) :: rs
407 : TYPE(pw_r3d_rs_type), INTENT(IN) :: rhos
408 : INTEGER, DIMENSION(3), INTENT(IN) :: center
409 :
410 : INTEGER :: i, ia, ii
411 : INTEGER, DIMENSION(3) :: nc
412 : LOGICAL :: folded
413 :
414 24639 : folded = .FALSE.
415 :
416 616464 : DO i = rhos%pw_grid%bounds(1, 1), rhos%pw_grid%bounds(2, 1)
417 591825 : ia = i - rhos%pw_grid%bounds(1, 1) + 1
418 591825 : ii = center(1) + i - rs%lb_local(1)
419 616464 : IF (ii < 0) THEN
420 24659 : rs%px(ia) = ii + rs%npts_local(1) + 1
421 24659 : folded = .TRUE.
422 567166 : ELSEIF (ii >= rs%npts_local(1)) THEN
423 19957 : rs%px(ia) = ii - rs%npts_local(1) + 1
424 19957 : folded = .TRUE.
425 : ELSE
426 547209 : rs%px(ia) = ii + 1
427 : END IF
428 : END DO
429 616464 : DO i = rhos%pw_grid%bounds(1, 2), rhos%pw_grid%bounds(2, 2)
430 591825 : ia = i - rhos%pw_grid%bounds(1, 2) + 1
431 591825 : ii = center(2) + i - rs%lb_local(2)
432 616464 : IF (ii < 0) THEN
433 6199 : rs%py(ia) = ii + rs%npts_local(2) + 1
434 6199 : folded = .TRUE.
435 585626 : ELSEIF (ii >= rs%npts_local(2)) THEN
436 7161 : rs%py(ia) = ii - rs%npts_local(2) + 1
437 7161 : folded = .TRUE.
438 : ELSE
439 578465 : rs%py(ia) = ii + 1
440 : END IF
441 : END DO
442 616464 : DO i = rhos%pw_grid%bounds(1, 3), rhos%pw_grid%bounds(2, 3)
443 591825 : ia = i - rhos%pw_grid%bounds(1, 3) + 1
444 591825 : ii = center(3) + i - rs%lb_local(3)
445 616464 : IF (ii < 0) THEN
446 25463 : rs%pz(ia) = ii + rs%npts_local(3) + 1
447 25463 : folded = .TRUE.
448 566362 : ELSEIF (ii >= rs%npts_local(3)) THEN
449 4783 : rs%pz(ia) = ii - rs%npts_local(3) + 1
450 4783 : folded = .TRUE.
451 : ELSE
452 561579 : rs%pz(ia) = ii + 1
453 : END IF
454 : END DO
455 :
456 24639 : IF (folded) THEN
457 : CALL dg_add_patch(rs%r, rhos%array, rhos%pw_grid%npts, &
458 13854 : rs%px, rs%py, rs%pz)
459 : ELSE
460 10785 : nc(1) = rs%px(1) - 1
461 10785 : nc(2) = rs%py(1) - 1
462 10785 : nc(3) = rs%pz(1) - 1
463 10785 : CALL dg_add_patch(rs%r, rhos%array, rhos%pw_grid%npts, nc)
464 : END IF
465 :
466 24639 : END SUBROUTINE dg_sum_patch_coef
467 :
468 : ! **************************************************************************************************
469 : !> \brief ...
470 : !> \param rs ...
471 : !> \param rhos ...
472 : !> \param center ...
473 : ! **************************************************************************************************
474 4177545 : PURE SUBROUTINE dg_sum_patch_arr(rs, rhos, center)
475 :
476 : TYPE(realspace_grid_type), INTENT(INOUT) :: rs
477 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: rhos
478 : INTEGER, DIMENSION(3), INTENT(IN) :: center
479 :
480 : INTEGER :: i, ia, ii
481 : INTEGER, DIMENSION(3) :: lb, nc, ns, ub
482 : LOGICAL :: folded
483 :
484 4177545 : ns(1) = SIZE(rhos, 1)
485 4177545 : ns(2) = SIZE(rhos, 2)
486 4177545 : ns(3) = SIZE(rhos, 3)
487 16710180 : lb = -(ns - 1)/2
488 16710180 : ub = lb + ns - 1
489 4177545 : folded = .FALSE.
490 :
491 26176507 : DO i = lb(1), ub(1)
492 21998962 : ia = i - lb(1) + 1
493 21998962 : ii = center(1) + i - rs%lb_local(1)
494 26176507 : IF (ii < 0) THEN
495 348169 : rs%px(ia) = ii + rs%npts_local(1) + 1
496 348169 : folded = .TRUE.
497 21650793 : ELSEIF (ii >= rs%npts_local(1)) THEN
498 757735 : rs%px(ia) = ii - rs%npts_local(1) + 1
499 757735 : folded = .TRUE.
500 : ELSE
501 20893058 : rs%px(ia) = ii + 1
502 : END IF
503 : END DO
504 26176507 : DO i = lb(2), ub(2)
505 21998962 : ia = i - lb(2) + 1
506 21998962 : ii = center(2) + i - rs%lb_local(2)
507 26176507 : IF (ii < 0) THEN
508 281355 : rs%py(ia) = ii + rs%npts_local(2) + 1
509 281355 : folded = .TRUE.
510 21717607 : ELSEIF (ii >= rs%npts_local(2)) THEN
511 634464 : rs%py(ia) = ii - rs%npts_local(2) + 1
512 634464 : folded = .TRUE.
513 : ELSE
514 21083143 : rs%py(ia) = ii + 1
515 : END IF
516 : END DO
517 26176507 : DO i = lb(3), ub(3)
518 21998962 : ia = i - lb(3) + 1
519 21998962 : ii = center(3) + i - rs%lb_local(3)
520 26176507 : IF (ii < 0) THEN
521 283067 : rs%pz(ia) = ii + rs%npts_local(3) + 1
522 283067 : folded = .TRUE.
523 21715895 : ELSEIF (ii >= rs%npts_local(3)) THEN
524 771820 : rs%pz(ia) = ii - rs%npts_local(3) + 1
525 771820 : folded = .TRUE.
526 : ELSE
527 20944075 : rs%pz(ia) = ii + 1
528 : END IF
529 : END DO
530 :
531 4177545 : IF (folded) THEN
532 1518211 : CALL dg_add_patch(rs%r, rhos, ns, rs%px, rs%py, rs%pz)
533 : ELSE
534 2659334 : nc(1) = rs%px(1) - 1
535 2659334 : nc(2) = rs%py(1) - 1
536 2659334 : nc(3) = rs%pz(1) - 1
537 2659334 : CALL dg_add_patch(rs%r, rhos, ns, nc)
538 : END IF
539 :
540 4177545 : END SUBROUTINE dg_sum_patch_arr
541 :
542 : ! **************************************************************************************************
543 : !> \brief ...
544 : !> \param drpot ...
545 : !> \param rhos ...
546 : !> \param center ...
547 : !> \param force ...
548 : ! **************************************************************************************************
549 4002407 : PURE SUBROUTINE dg_sum_patch_force_arr_3d(drpot, rhos, center, force)
550 :
551 : TYPE(realspace_grid_type), DIMENSION(:), &
552 : INTENT(INOUT) :: drpot
553 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: rhos
554 : INTEGER, DIMENSION(3), INTENT(IN) :: center
555 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: force
556 :
557 : INTEGER :: i, ia, ii
558 : INTEGER, DIMENSION(3) :: lb, nc, ns, ub
559 : LOGICAL :: folded
560 :
561 4002407 : ns(1) = SIZE(rhos, 1)
562 4002407 : ns(2) = SIZE(rhos, 2)
563 4002407 : ns(3) = SIZE(rhos, 3)
564 16009628 : lb = -(ns - 1)/2
565 16009628 : ub = lb + ns - 1
566 4002407 : folded = .FALSE.
567 :
568 25059085 : DO i = lb(1), ub(1)
569 21056678 : ia = i - lb(1) + 1
570 21056678 : ii = center(1) + i - drpot(1)%lb_local(1)
571 25059085 : IF (ii < 0) THEN
572 337618 : drpot(1)%px(ia) = ii + drpot(1)%npts_local(1) + 1
573 337618 : folded = .TRUE.
574 20719060 : ELSEIF (ii >= drpot(1)%npts_local(1)) THEN
575 719176 : drpot(1)%px(ia) = ii - drpot(1)%npts_local(1) + 1
576 719176 : folded = .TRUE.
577 : ELSE
578 19999884 : drpot(1)%px(ia) = ii + 1
579 : END IF
580 : END DO
581 25059085 : DO i = lb(2), ub(2)
582 21056678 : ia = i - lb(2) + 1
583 21056678 : ii = center(2) + i - drpot(1)%lb_local(2)
584 25059085 : IF (ii < 0) THEN
585 268679 : drpot(1)%py(ia) = ii + drpot(1)%npts_local(2) + 1
586 268679 : folded = .TRUE.
587 20787999 : ELSEIF (ii >= drpot(1)%npts_local(2)) THEN
588 601437 : drpot(1)%py(ia) = ii - drpot(1)%npts_local(2) + 1
589 601437 : folded = .TRUE.
590 : ELSE
591 20186562 : drpot(1)%py(ia) = ii + 1
592 : END IF
593 : END DO
594 25059085 : DO i = lb(3), ub(3)
595 21056678 : ia = i - lb(3) + 1
596 21056678 : ii = center(3) + i - drpot(1)%lb_local(3)
597 25059085 : IF (ii < 0) THEN
598 275584 : drpot(1)%pz(ia) = ii + drpot(1)%npts_local(3) + 1
599 275584 : folded = .TRUE.
600 20781094 : ELSEIF (ii >= drpot(1)%npts_local(3)) THEN
601 740585 : drpot(1)%pz(ia) = ii - drpot(1)%npts_local(3) + 1
602 740585 : folded = .TRUE.
603 : ELSE
604 20040509 : drpot(1)%pz(ia) = ii + 1
605 : END IF
606 : END DO
607 :
608 4002407 : IF (folded) THEN
609 : CALL dg_int_patch_3d(drpot(1)%r, drpot(2)%r, &
610 : drpot(3)%r, rhos, force, ns, &
611 1464331 : drpot(1)%px, drpot(1)%py, drpot(1)%pz)
612 : ELSE
613 2538076 : nc(1) = drpot(1)%px(1) - 1
614 2538076 : nc(2) = drpot(1)%py(1) - 1
615 2538076 : nc(3) = drpot(1)%pz(1) - 1
616 : CALL dg_int_patch_3d(drpot(1)%r, drpot(2)%r, &
617 2538076 : drpot(3)%r, rhos, force, ns, nc)
618 : END IF
619 :
620 4002407 : END SUBROUTINE dg_sum_patch_force_arr_3d
621 :
622 : ! **************************************************************************************************
623 : !> \brief ...
624 : !> \param drpot ...
625 : !> \param rhos ...
626 : !> \param center ...
627 : !> \param force ...
628 : ! **************************************************************************************************
629 186989 : PURE SUBROUTINE dg_sum_patch_force_arr_1d(drpot, rhos, center, force)
630 :
631 : TYPE(realspace_grid_type), INTENT(INOUT) :: drpot
632 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: rhos
633 : INTEGER, DIMENSION(3), INTENT(IN) :: center
634 : REAL(KIND=dp), INTENT(OUT) :: force
635 :
636 : INTEGER :: i, ia, ii
637 : INTEGER, DIMENSION(3) :: lb, nc, ns, ub
638 : LOGICAL :: folded
639 :
640 186989 : ns(1) = SIZE(rhos, 1)
641 186989 : ns(2) = SIZE(rhos, 2)
642 186989 : ns(3) = SIZE(rhos, 3)
643 747956 : lb = -(ns - 1)/2
644 747956 : ub = lb + ns - 1
645 186989 : folded = .FALSE.
646 :
647 1200237 : DO i = lb(1), ub(1)
648 1013248 : ia = i - lb(1) + 1
649 1013248 : ii = center(1) + i - drpot%lb_local(1)
650 1200237 : IF (ii < 0) THEN
651 12013 : drpot%px(ia) = ii + drpot%npts_local(1) + 1
652 12013 : folded = .TRUE.
653 1001235 : ELSEIF (ii >= drpot%desc%npts(1)) THEN
654 41275 : drpot%px(ia) = ii - drpot%npts_local(1) + 1
655 41275 : folded = .TRUE.
656 : ELSE
657 959960 : drpot%px(ia) = ii + 1
658 : END IF
659 : END DO
660 1200237 : DO i = lb(2), ub(2)
661 1013248 : ia = i - lb(2) + 1
662 1013248 : ii = center(2) + i - drpot%lb_local(2)
663 1200237 : IF (ii < 0) THEN
664 14546 : drpot%py(ia) = ii + drpot%npts_local(2) + 1
665 14546 : folded = .TRUE.
666 998702 : ELSEIF (ii >= drpot%desc%npts(2)) THEN
667 35658 : drpot%py(ia) = ii - drpot%npts_local(2) + 1
668 35658 : folded = .TRUE.
669 : ELSE
670 963044 : drpot%py(ia) = ii + 1
671 : END IF
672 : END DO
673 1200237 : DO i = lb(3), ub(3)
674 1013248 : ia = i - lb(3) + 1
675 1013248 : ii = center(3) + i - drpot%lb_local(3)
676 1200237 : IF (ii < 0) THEN
677 8476 : drpot%pz(ia) = ii + drpot%npts_local(3) + 1
678 8476 : folded = .TRUE.
679 1004772 : ELSEIF (ii >= drpot%desc%npts(3)) THEN
680 33896 : drpot%pz(ia) = ii - drpot%npts_local(3) + 1
681 33896 : folded = .TRUE.
682 : ELSE
683 970876 : drpot%pz(ia) = ii + 1
684 : END IF
685 : END DO
686 :
687 186989 : IF (folded) THEN
688 : CALL dg_int_patch_1d(drpot%r, rhos, force, ns, &
689 59309 : drpot%px, drpot%py, drpot%pz)
690 : ELSE
691 127680 : nc(1) = drpot%px(1) - 1
692 127680 : nc(2) = drpot%py(1) - 1
693 127680 : nc(3) = drpot%pz(1) - 1
694 127680 : CALL dg_int_patch_1d(drpot%r, rhos, force, ns, nc)
695 : END IF
696 :
697 186989 : END SUBROUTINE dg_sum_patch_force_arr_1d
698 :
699 : ! **************************************************************************************************
700 : !> \brief ...
701 : !> \param drpot ...
702 : !> \param rhos ...
703 : !> \param center ...
704 : !> \param force ...
705 : ! **************************************************************************************************
706 24639 : PURE SUBROUTINE dg_sum_patch_force_coef_3d(drpot, rhos, center, force)
707 :
708 : TYPE(realspace_grid_type), DIMENSION(:), &
709 : INTENT(INOUT) :: drpot
710 : TYPE(pw_r3d_rs_type), INTENT(IN) :: rhos
711 : INTEGER, DIMENSION(3), INTENT(IN) :: center
712 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: force
713 :
714 : INTEGER :: i, ia, ii
715 : INTEGER, DIMENSION(3) :: nc
716 : LOGICAL :: folded
717 :
718 24639 : folded = .FALSE.
719 :
720 616464 : DO i = rhos%pw_grid%bounds(1, 1), rhos%pw_grid%bounds(2, 1)
721 591825 : ia = i - rhos%pw_grid%bounds(1, 1) + 1
722 591825 : ii = center(1) + i - drpot(1)%lb_local(1)
723 616464 : IF (ii < 0) THEN
724 24659 : drpot(1)%px(ia) = ii + drpot(1)%desc%npts(1) + 1
725 24659 : folded = .TRUE.
726 567166 : ELSEIF (ii >= drpot(1)%desc%npts(1)) THEN
727 19957 : drpot(1)%px(ia) = ii - drpot(1)%desc%npts(1) + 1
728 19957 : folded = .TRUE.
729 : ELSE
730 547209 : drpot(1)%px(ia) = ii + 1
731 : END IF
732 : END DO
733 616464 : DO i = rhos%pw_grid%bounds(1, 2), rhos%pw_grid%bounds(2, 2)
734 591825 : ia = i - rhos%pw_grid%bounds(1, 2) + 1
735 591825 : ii = center(2) + i - drpot(1)%lb_local(2)
736 616464 : IF (ii < 0) THEN
737 6199 : drpot(1)%py(ia) = ii + drpot(1)%desc%npts(2) + 1
738 6199 : folded = .TRUE.
739 585626 : ELSEIF (ii >= drpot(1)%desc%npts(2)) THEN
740 7161 : drpot(1)%py(ia) = ii - drpot(1)%desc%npts(2) + 1
741 7161 : folded = .TRUE.
742 : ELSE
743 578465 : drpot(1)%py(ia) = ii + 1
744 : END IF
745 : END DO
746 616464 : DO i = rhos%pw_grid%bounds(1, 3), rhos%pw_grid%bounds(2, 3)
747 591825 : ia = i - rhos%pw_grid%bounds(1, 3) + 1
748 591825 : ii = center(3) + i - drpot(1)%lb_local(3)
749 616464 : IF (ii < 0) THEN
750 25463 : drpot(1)%pz(ia) = ii + drpot(1)%desc%npts(3) + 1
751 25463 : folded = .TRUE.
752 566362 : ELSEIF (ii >= drpot(1)%desc%npts(3)) THEN
753 4783 : drpot(1)%pz(ia) = ii - drpot(1)%desc%npts(3) + 1
754 4783 : folded = .TRUE.
755 : ELSE
756 561579 : drpot(1)%pz(ia) = ii + 1
757 : END IF
758 : END DO
759 :
760 24639 : IF (folded) THEN
761 : CALL dg_int_patch_3d(drpot(1)%r, drpot(2)%r, &
762 : drpot(3)%r, rhos%array, force, rhos%pw_grid%npts, &
763 13854 : drpot(1)%px, drpot(1)%py, drpot(1)%pz)
764 : ELSE
765 10785 : nc(1) = drpot(1)%px(1) - 1
766 10785 : nc(2) = drpot(1)%py(1) - 1
767 10785 : nc(3) = drpot(1)%pz(1) - 1
768 : CALL dg_int_patch_3d(drpot(1)%r, drpot(2)%r, &
769 10785 : drpot(3)%r, rhos%array, force, rhos%pw_grid%npts, nc)
770 : END IF
771 :
772 24639 : END SUBROUTINE dg_sum_patch_force_coef_3d
773 :
774 : ! **************************************************************************************************
775 : !> \brief ...
776 : !> \param drpot ...
777 : !> \param rhos ...
778 : !> \param center ...
779 : !> \param force ...
780 : ! **************************************************************************************************
781 303 : PURE SUBROUTINE dg_sum_patch_force_coef_1d(drpot, rhos, center, force)
782 :
783 : TYPE(realspace_grid_type), INTENT(INOUT) :: drpot
784 : TYPE(pw_r3d_rs_type), INTENT(IN) :: rhos
785 : INTEGER, DIMENSION(3), INTENT(IN) :: center
786 : REAL(KIND=dp), INTENT(OUT) :: force
787 :
788 : INTEGER :: i, ia, ii
789 : INTEGER, DIMENSION(3) :: nc
790 : LOGICAL :: folded
791 :
792 303 : folded = .FALSE.
793 :
794 4848 : DO i = rhos%pw_grid%bounds(1, 1), rhos%pw_grid%bounds(2, 1)
795 4545 : ia = i - rhos%pw_grid%bounds(1, 1) + 1
796 4545 : ii = center(1) + i - drpot%lb_local(1)
797 4848 : IF (ii < 0) THEN
798 0 : drpot%px(ia) = ii + drpot%desc%npts(1) + 1
799 0 : folded = .TRUE.
800 4545 : ELSEIF (ii >= drpot%desc%npts(1)) THEN
801 0 : drpot%px(ia) = ii - drpot%desc%npts(1) + 1
802 0 : folded = .TRUE.
803 : ELSE
804 4545 : drpot%px(ia) = ii + 1
805 : END IF
806 : END DO
807 4848 : DO i = rhos%pw_grid%bounds(1, 2), rhos%pw_grid%bounds(2, 2)
808 4545 : ia = i - rhos%pw_grid%bounds(1, 2) + 1
809 4545 : ii = center(2) + i - drpot%lb_local(2)
810 4848 : IF (ii < 0) THEN
811 0 : drpot%py(ia) = ii + drpot%desc%npts(2) + 1
812 0 : folded = .TRUE.
813 4545 : ELSEIF (ii >= drpot%desc%npts(2)) THEN
814 0 : drpot%py(ia) = ii - drpot%desc%npts(2) + 1
815 0 : folded = .TRUE.
816 : ELSE
817 4545 : drpot%py(ia) = ii + 1
818 : END IF
819 : END DO
820 4848 : DO i = rhos%pw_grid%bounds(1, 3), rhos%pw_grid%bounds(2, 3)
821 4545 : ia = i - rhos%pw_grid%bounds(1, 3) + 1
822 4545 : ii = center(3) + i - drpot%lb_local(3)
823 4848 : IF (ii < 0) THEN
824 0 : drpot%pz(ia) = ii + drpot%desc%npts(3) + 1
825 0 : folded = .TRUE.
826 4545 : ELSEIF (ii >= drpot%desc%npts(3)) THEN
827 0 : drpot%pz(ia) = ii - drpot%desc%npts(3) + 1
828 0 : folded = .TRUE.
829 : ELSE
830 4545 : drpot%pz(ia) = ii + 1
831 : END IF
832 : END DO
833 :
834 303 : IF (folded) THEN
835 : CALL dg_int_patch_1d(drpot%r, rhos%array, force, &
836 0 : rhos%pw_grid%npts, drpot%px, drpot%py, drpot%pz)
837 : ELSE
838 303 : nc(1) = drpot%px(1) - 1
839 303 : nc(2) = drpot%py(1) - 1
840 303 : nc(3) = drpot%pz(1) - 1
841 303 : CALL dg_int_patch_1d(drpot%r, rhos%array, force, rhos%pw_grid%npts, nc)
842 : END IF
843 :
844 303 : END SUBROUTINE dg_sum_patch_force_coef_1d
845 :
846 : ! **************************************************************************************************
847 : !> \brief ...
848 : !> \param rho0 ...
849 : !> \param rhos1 ...
850 : !> \param charge1 ...
851 : !> \param ex1 ...
852 : !> \param ey1 ...
853 : !> \param ez1 ...
854 : ! **************************************************************************************************
855 2351 : SUBROUTINE dg_get_patch_1(rho0, rhos1, charge1, ex1, ey1, ez1)
856 :
857 : TYPE(pw_r3d_rs_type), INTENT(IN) :: rho0
858 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: rhos1
859 : REAL(KIND=dp), INTENT(IN) :: charge1
860 : COMPLEX(KIND=dp), DIMENSION(:), INTENT(IN) :: ex1, ey1, ez1
861 :
862 : COMPLEX(KIND=dp) :: za, zb
863 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: zs
864 : INTEGER :: nd(3)
865 : TYPE(pw_c3d_rs_type) :: cd
866 :
867 9404 : nd = rhos1%pw_grid%npts
868 :
869 7053 : ALLOCATE (zs(nd(1)*nd(2)))
870 1350526 : zs = 0.0_dp
871 2351 : CALL cd%create(rho0%pw_grid)
872 2351 : CALL pw_zero(cd)
873 :
874 2351 : za = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
875 2351 : zb = CMPLX(charge1, 0.0_dp, KIND=dp)
876 2351 : CALL rankup(nd, za, cd%array, zb, ex1, ey1, ez1, zs)
877 2351 : CALL pw_multiply_with(cd, rho0)
878 2351 : CALL fft3d(BWFFT, nd, cd%array)
879 2351 : CALL pw_copy(cd, rhos1)
880 :
881 2351 : DEALLOCATE (zs)
882 2351 : CALL cd%release()
883 :
884 2351 : END SUBROUTINE dg_get_patch_1
885 :
886 : ! **************************************************************************************************
887 : !> \brief ...
888 : !> \param rho0 ...
889 : !> \param rhos1 ...
890 : !> \param rhos2 ...
891 : !> \param charge1 ...
892 : !> \param charge2 ...
893 : !> \param ex1 ...
894 : !> \param ey1 ...
895 : !> \param ez1 ...
896 : !> \param ex2 ...
897 : !> \param ey2 ...
898 : !> \param ez2 ...
899 : ! **************************************************************************************************
900 23615 : SUBROUTINE dg_get_patch_2(rho0, rhos1, rhos2, charge1, charge2, &
901 23615 : ex1, ey1, ez1, ex2, ey2, ez2)
902 :
903 : TYPE(pw_r3d_rs_type), INTENT(IN) :: rho0
904 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: rhos1, rhos2
905 : REAL(KIND=dp), INTENT(IN) :: charge1, charge2
906 : COMPLEX(KIND=dp), DIMENSION(:), INTENT(IN) :: ex1, ey1, ez1, ex2, ey2, ez2
907 :
908 : COMPLEX(KIND=dp) :: za, zb
909 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: zs
910 : INTEGER :: nd(3)
911 : TYPE(pw_c3d_rs_type) :: cd
912 :
913 94460 : nd = rhos1%pw_grid%npts
914 :
915 70845 : ALLOCATE (zs(nd(1)*nd(2)))
916 13816990 : zs = 0.0_dp
917 23615 : CALL cd%create(rhos1%pw_grid)
918 23615 : CALL pw_zero(cd)
919 :
920 23615 : za = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
921 23615 : zb = CMPLX(charge2, 0.0_dp, KIND=dp)
922 23615 : CALL rankup(nd, za, cd%array, zb, ex2, ey2, ez2, zs)
923 23615 : za = CMPLX(0.0_dp, 1.0_dp, KIND=dp)
924 23615 : zb = CMPLX(charge1, 0.0_dp, KIND=dp)
925 23615 : CALL rankup(nd, za, cd%array, zb, ex1, ey1, ez1, zs)
926 23615 : CALL pw_multiply_with(cd, rho0)
927 23615 : CALL fft3d(BWFFT, nd, cd%array)
928 23615 : CALL copy_cri(cd%array, rhos1%array, rhos2%array)
929 :
930 23615 : DEALLOCATE (zs)
931 23615 : CALL cd%release()
932 :
933 23615 : END SUBROUTINE dg_get_patch_2
934 :
935 : ! **************************************************************************************************
936 : !> \brief ...
937 : !> \param rb ...
938 : !> \param rs ...
939 : !> \param ns ...
940 : !> \param nc ...
941 : ! **************************************************************************************************
942 2670119 : PURE SUBROUTINE dg_add_patch_simple(rb, rs, ns, nc)
943 :
944 : REAL(KIND=dp), INTENT(INOUT) :: rb(:, :, :)
945 : REAL(KIND=dp), INTENT(IN) :: rs(:, :, :)
946 : INTEGER, INTENT(IN) :: ns(3), nc(3)
947 :
948 : INTEGER :: i, ii, j, jj, k, kk
949 :
950 16707401 : DO k = 1, ns(3)
951 14037282 : kk = nc(3) + k
952 97231857 : DO j = 1, ns(2)
953 80524456 : jj = nc(2) + j
954 668172880 : DO i = 1, ns(1)
955 573611142 : ii = nc(1) + i
956 654135598 : rb(ii, jj, kk) = rb(ii, jj, kk) + rs(i, j, k)
957 : END DO
958 : END DO
959 : END DO
960 :
961 2670119 : END SUBROUTINE dg_add_patch_simple
962 :
963 : ! **************************************************************************************************
964 : !> \brief ...
965 : !> \param rb ...
966 : !> \param rs ...
967 : !> \param ns ...
968 : !> \param px ...
969 : !> \param py ...
970 : !> \param pz ...
971 : ! **************************************************************************************************
972 1532065 : PURE SUBROUTINE dg_add_patch_folded(rb, rs, ns, px, py, pz)
973 :
974 : REAL(KIND=dp), INTENT(INOUT) :: rb(:, :, :)
975 : REAL(KIND=dp), INTENT(IN) :: rs(:, :, :)
976 : INTEGER, INTENT(IN) :: ns(:)
977 : INTEGER, DIMENSION(:), INTENT(IN) :: px, py, pz
978 :
979 : INTEGER :: i, ii, j, jj, k, kk
980 :
981 10085570 : DO k = 1, ns(3)
982 8553505 : kk = pz(k)
983 63764667 : DO j = 1, ns(2)
984 53679097 : jj = py(j)
985 514198887 : DO i = 1, ns(1)
986 451966285 : ii = px(i)
987 505645382 : rb(ii, jj, kk) = rb(ii, jj, kk) + rs(i, j, k)
988 : END DO
989 : END DO
990 : END DO
991 :
992 1532065 : END SUBROUTINE dg_add_patch_folded
993 :
994 : ! **************************************************************************************************
995 : !> \brief ...
996 : !> \param rbx ...
997 : !> \param rby ...
998 : !> \param rbz ...
999 : !> \param rs ...
1000 : !> \param f ...
1001 : !> \param ns ...
1002 : !> \param nc ...
1003 : ! **************************************************************************************************
1004 2548861 : PURE SUBROUTINE dg_int_patch_simple_3d(rbx, rby, rbz, rs, f, ns, nc)
1005 :
1006 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: rbx, rby, rbz, rs
1007 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: f
1008 : INTEGER, INTENT(IN) :: ns(3), nc(3)
1009 :
1010 : INTEGER :: i, ii, j, jj, k, kk
1011 : REAL(KIND=dp) :: s
1012 :
1013 2548861 : f = 0.0_dp
1014 15936809 : DO k = 1, ns(3)
1015 13387948 : kk = nc(3) + k
1016 92862551 : DO j = 1, ns(2)
1017 76925742 : jj = nc(2) + j
1018 643325726 : DO i = 1, ns(1)
1019 553012036 : ii = nc(1) + i
1020 553012036 : s = rs(i, j, k)
1021 553012036 : f(1) = f(1) + s*rbx(ii, jj, kk)
1022 553012036 : f(2) = f(2) + s*rby(ii, jj, kk)
1023 629937778 : f(3) = f(3) + s*rbz(ii, jj, kk)
1024 : END DO
1025 : END DO
1026 : END DO
1027 :
1028 2548861 : END SUBROUTINE dg_int_patch_simple_3d
1029 :
1030 : ! **************************************************************************************************
1031 : !> \brief ...
1032 : !> \param rb ...
1033 : !> \param rs ...
1034 : !> \param f ...
1035 : !> \param ns ...
1036 : !> \param nc ...
1037 : ! **************************************************************************************************
1038 127983 : PURE SUBROUTINE dg_int_patch_simple_1d(rb, rs, f, ns, nc)
1039 :
1040 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: rb, rs
1041 : REAL(KIND=dp), INTENT(OUT) :: f
1042 : INTEGER, INTENT(IN) :: ns(3), nc(3)
1043 :
1044 : INTEGER :: i, ii, j, jj, k, kk
1045 : REAL(KIND=dp) :: s
1046 :
1047 127983 : f = 0.0_dp
1048 820195 : DO k = 1, ns(3)
1049 692212 : kk = nc(3) + k
1050 4722897 : DO j = 1, ns(2)
1051 3902702 : jj = nc(2) + j
1052 27712246 : DO i = 1, ns(1)
1053 23117332 : ii = nc(1) + i
1054 23117332 : s = rs(i, j, k)
1055 27020034 : f = f + s*rb(ii, jj, kk)
1056 : END DO
1057 : END DO
1058 : END DO
1059 :
1060 127983 : END SUBROUTINE dg_int_patch_simple_1d
1061 :
1062 : ! **************************************************************************************************
1063 : !> \brief ...
1064 : !> \param rbx ...
1065 : !> \param rby ...
1066 : !> \param rbz ...
1067 : !> \param rs ...
1068 : !> \param f ...
1069 : !> \param ns ...
1070 : !> \param px ...
1071 : !> \param py ...
1072 : !> \param pz ...
1073 : ! **************************************************************************************************
1074 1478185 : PURE SUBROUTINE dg_int_patch_folded_3d(rbx, rby, rbz, rs, f, ns, px, py, pz)
1075 :
1076 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: rbx, rby, rbz, rs
1077 : REAL(KIND=dp), DIMENSION(3), INTENT(INOUT) :: f
1078 : INTEGER, INTENT(IN) :: ns(3)
1079 : INTEGER, DIMENSION(:), INTENT(IN) :: px, py, pz
1080 :
1081 : INTEGER :: i, ii, j, jj, k, kk
1082 : REAL(KIND=dp) :: s
1083 :
1084 1478185 : f = 0.0_dp
1085 9738740 : DO k = 1, ns(3)
1086 8260555 : kk = pz(k)
1087 61766559 : DO j = 1, ns(2)
1088 52027819 : jj = py(j)
1089 502615893 : DO i = 1, ns(1)
1090 442327519 : ii = px(i)
1091 442327519 : s = rs(i, j, k)
1092 442327519 : f(1) = f(1) + s*rbx(ii, jj, kk)
1093 442327519 : f(2) = f(2) + s*rby(ii, jj, kk)
1094 494355338 : f(3) = f(3) + s*rbz(ii, jj, kk)
1095 : END DO
1096 : END DO
1097 : END DO
1098 :
1099 1478185 : END SUBROUTINE dg_int_patch_folded_3d
1100 :
1101 : ! **************************************************************************************************
1102 : !> \brief ...
1103 : !> \param rb ...
1104 : !> \param rs ...
1105 : !> \param f ...
1106 : !> \param ns ...
1107 : !> \param px ...
1108 : !> \param py ...
1109 : !> \param pz ...
1110 : ! **************************************************************************************************
1111 59309 : PURE SUBROUTINE dg_int_patch_folded_1d(rb, rs, f, ns, px, py, pz)
1112 :
1113 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: rb, rs
1114 : REAL(KIND=dp), INTENT(INOUT) :: f
1115 : INTEGER, INTENT(IN) :: ns(3)
1116 : INTEGER, DIMENSION(:), INTENT(IN) :: px, py, pz
1117 :
1118 : INTEGER :: i, ii, j, jj, k, kk
1119 : REAL(KIND=dp) :: s
1120 :
1121 59309 : f = 0.0_dp
1122 384890 : DO k = 1, ns(3)
1123 325581 : kk = pz(k)
1124 2236745 : DO j = 1, ns(2)
1125 1851855 : jj = py(j)
1126 13078707 : DO i = 1, ns(1)
1127 10901271 : ii = px(i)
1128 10901271 : s = rs(i, j, k)
1129 12753126 : f = f + s*rb(ii, jj, kk)
1130 : END DO
1131 : END DO
1132 : END DO
1133 :
1134 59309 : END SUBROUTINE dg_int_patch_folded_1d
1135 :
1136 : ! **************************************************************************************************
1137 : !> \brief ...
1138 : !> \param n ...
1139 : !> \param za ...
1140 : !> \param cmat ...
1141 : !> \param zb ...
1142 : !> \param ex ...
1143 : !> \param ey ...
1144 : !> \param ez ...
1145 : !> \param scr ...
1146 : ! **************************************************************************************************
1147 49581 : SUBROUTINE rankup(n, za, cmat, zb, ex, ey, ez, scr)
1148 : !
1149 : ! cmat(i,j,k) <- za * cmat(i,j,k) + ex(i) * ey(j) * ez(k)
1150 : !
1151 :
1152 : INTEGER, DIMENSION(3), INTENT(IN) :: n
1153 : COMPLEX(KIND=dp), INTENT(IN) :: za
1154 : COMPLEX(KIND=dp), DIMENSION(:, :, :), &
1155 : INTENT(INOUT) :: cmat
1156 : COMPLEX(KIND=dp), INTENT(IN) :: zb
1157 : COMPLEX(KIND=dp), DIMENSION(:), INTENT(IN) :: ex, ey, ez
1158 : COMPLEX(KIND=dp), DIMENSION(:), INTENT(INOUT) :: scr
1159 :
1160 : INTEGER :: n2, n3
1161 :
1162 49581 : n2 = n(1)*n(2)
1163 49581 : n3 = n2*n(3)
1164 28984506 : scr(1:n2) = z_zero
1165 49581 : CALL zgeru(n(1), n(2), zb, ex, 1, ey, 1, scr, n(1))
1166 49581 : CALL zscal(n3, za, cmat, 1)
1167 49581 : CALL zgeru(n2, n(3), z_one, scr, 1, ez, 1, cmat, n2)
1168 :
1169 49581 : END SUBROUTINE rankup
1170 :
1171 : ! **************************************************************************************************
1172 : !> \brief Copy a the real and imag. parts of a complex 3D array into two real arrays
1173 : !> \param z the complex array
1174 : !> \param r1 the real array for the real part
1175 : !> \param r2 the real array for the imaginary part
1176 : ! **************************************************************************************************
1177 23615 : SUBROUTINE copy_cri(z, r1, r2)
1178 : !
1179 : ! r1 = real ( z )
1180 : ! r2 = imag ( z )
1181 : !
1182 :
1183 : COMPLEX(KIND=dp), INTENT(IN) :: z(:, :, :)
1184 : REAL(KIND=dp), INTENT(INOUT) :: r1(:, :, :), r2(:, :, :)
1185 :
1186 23615 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE), SHARED(r1,r2,z)
1187 : r1(:, :, :) = REAL(z(:, :, :), KIND=dp)
1188 : r2(:, :, :) = AIMAG(z(:, :, :))
1189 : !$OMP END PARALLEL WORKSHARE
1190 :
1191 23615 : END SUBROUTINE copy_cri
1192 :
1193 : END MODULE dgs
|
