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