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