(git:34ef472)
dgs.F
Go to the documentation of this file.
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  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  CALL dg_find_cutoff(b_cell_hmat, npts_s, cutoff_radius, grid_s, grid_b, cutoff)
109 
110  ecut = 0.5_dp*cutoff*cutoff
111  bo = grid_b%bounds
112  IF (PRESENT(grid_ref)) THEN
113  CALL pw_grid_setup(b_cell_hmat, grid_b, bounds=bo, cutoff=ecut, spherical=.true., &
114  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  rs_dims=rs_dims, iounit=iounit, fft_usage=fft_usage)
118  END IF
119 
120  CALL dg_find_basis(grid_b%npts, b_cell_hmat, unit_cell_hmat)
121 
122  CALL dg_set_cell(grid_s%npts, unit_cell_hmat, s_cell_hmat)
123 
124  bo = grid_s%bounds
125  CALL pw_grid_setup(s_cell_hmat, grid_s, bounds=bo, cutoff=ecut, iounit=iounit, fft_usage=fft_usage)
126 
127  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  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  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  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  cell_hmat(3, 3)*cell_hmat(3, 3))
149  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  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  b_cell_deth = abs(det_3x3(b_cell_hmat))
174  IF (b_cell_deth < 1.0e-10_dp) THEN
175  CALL cp_abort(__location__, &
176  "An invalid set of cell vectors was specified. "// &
177  "The determinant det(h) is too small")
178  END IF
179  b_cell_h_inv = inv_3x3(b_cell_hmat)
180 
181  CALL fft_radix_operations(npts_s(1), nout(1), &
182  operation=fft_radix_next_odd)
183  CALL fft_radix_operations(npts_s(1), nout(2), &
184  operation=fft_radix_next_odd)
185  CALL fft_radix_operations(npts_s(1), nout(3), &
186  operation=fft_radix_next_odd)
187 
188  cell_lengths = get_cell_lengths(b_cell_hmat)
189 
190  CALL dg_get_spacing(nout, cutoff_radius, dr)
191  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  IF (nout(1) > grid_b%npts(1)) THEN
195  grid_b%npts(1) = nout(1)
196  dr(1) = cell_lengths(1)/real(nout(1), kind=dp)
197  END IF
198  IF (nout(2) > grid_b%npts(2)) THEN
199  grid_b%npts(2) = nout(2)
200  dr(2) = cell_lengths(2)/real(nout(2), kind=dp)
201  END IF
202  IF (nout(3) > grid_b%npts(3)) THEN
203  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  grid_b%bounds(1, :) = -grid_b%npts/2
209  grid_b%bounds(2, :) = +(grid_b%npts - 1)/2
210  grid_b%grid_span = halfspace
211 ! small grid bounds
212  grid_s%bounds(1, :) = -nout(:)/2
213  grid_s%bounds(2, :) = (+nout(:) - 1)/2
214  grid_s%grid_span = halfspace
215  grid_s%npts = nout
216 
217  cutoff = pw_find_cutoff(grid_b%npts, b_cell_h_inv)
218 
219  END SUBROUTINE dg_find_cutoff
220 
221 ! **************************************************************************************************
222 !> \brief ...
223 !> \param npts ...
224 !> \param cutoff_radius ...
225 !> \param dr ...
226 ! **************************************************************************************************
227  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  dr(:) = cutoff_radius/(real(npts(:), kind=dp)/2.0_dp)
234 
235  END SUBROUTINE dg_get_spacing
236 
237 ! **************************************************************************************************
238 !> \brief ...
239 !> \param b_cell_hmat ...
240 !> \param grid_b ...
241 !> \param grid_s ...
242 ! **************************************************************************************************
243  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  CALL dg_find_basis(grid_b%npts, b_cell_hmat, unit_cell_hmat)
250  CALL dg_set_cell(grid_s%npts, unit_cell_hmat, s_cell_hmat)
251  CALL pw_grid_change(s_cell_hmat, grid_s)
252  END SUBROUTINE dg_grid_change
253 
254 ! **************************************************************************************************
255 !> \brief ...
256 !> \param dr ...
257 !> \param cell_lengths ...
258 !> \param npts ...
259 ! **************************************************************************************************
260  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  nin(:) = nint(cell_lengths(:)/dr(:))
269  CALL fft_radix_operations(nin(1), npts(1), &
270  operation=fft_radix_closest)
271  CALL fft_radix_operations(nin(2), npts(2), &
272  operation=fft_radix_closest)
273  CALL fft_radix_operations(nin(3), npts(3), &
274  operation=fft_radix_closest)
275  dr(:) = cell_lengths(:)/real(npts(:), kind=dp)
276 
277  END SUBROUTINE dg_find_radix
278 
279 ! **************************************************************************************************
280 !> \brief ...
281 !> \param npts ...
282 !> \param cell_hmat ...
283 !> \param unit_cell_hmat ...
284 ! **************************************************************************************************
285  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  DO i = 1, 3
293  unit_cell_hmat(:, i) = cell_hmat(:, i)/real(npts(:), kind=dp)
294  END DO
295 
296  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  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  cell_hmat(:, 1) = unit_cell_hmat(:, 1)*npts(1)
314  cell_hmat(:, 2) = unit_cell_hmat(:, 2)*npts(2)
315  cell_hmat(:, 3) = unit_cell_hmat(:, 3)*npts(3)
316 
317  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  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  CALL dg_get_delta(cell_hmat, r, npts_s, npts_b, centre, delta)
344 
345  CALL structure_factor_evaluate(delta, lb, ex, ey, ez)
346 
347  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  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  cell_deth = abs(det_3x3(cell_hmat))
370  IF (cell_deth < 1.0e-10_dp) THEN
371  CALL cp_abort(__location__, &
372  "An invalid set of cell vectors was specified. "// &
373  "The determinant det(h) is too small")
374  END IF
375 
376  cell_h_inv = inv_3x3(cell_hmat)
377 
378 ! compute the scaled coordinate of atomi
379  s = matmul(cell_h_inv, r)
380  s = s - nint(s)
381 
382 ! find the continuous ``grid'' point (on big grid)
383  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  centre(:) = nint(grid_i(:))
387 
388 ! find the distance vector
389  delta(:) = (grid_i(:) - centre(:))/real(npts_s(:), kind=dp)
390 
391  centre(:) = centre(:) + npts_b(:)/2
392  centre(:) = modulo(centre(:), npts_b(:))
393  centre(:) = centre(:) - npts_b(:)/2
394 
395  END SUBROUTINE dg_get_delta
396 
397 ! **************************************************************************************************
398 !> \brief ...
399 !> \param rs ...
400 !> \param rhos ...
401 !> \param center ...
402 ! **************************************************************************************************
403  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  folded = .false.
414 
415  DO i = rhos%pw_grid%bounds(1, 1), rhos%pw_grid%bounds(2, 1)
416  ia = i - rhos%pw_grid%bounds(1, 1) + 1
417  ii = center(1) + i - rs%lb_local(1)
418  IF (ii < 0) THEN
419  rs%px(ia) = ii + rs%npts_local(1) + 1
420  folded = .true.
421  ELSEIF (ii >= rs%npts_local(1)) THEN
422  rs%px(ia) = ii - rs%npts_local(1) + 1
423  folded = .true.
424  ELSE
425  rs%px(ia) = ii + 1
426  END IF
427  END DO
428  DO i = rhos%pw_grid%bounds(1, 2), rhos%pw_grid%bounds(2, 2)
429  ia = i - rhos%pw_grid%bounds(1, 2) + 1
430  ii = center(2) + i - rs%lb_local(2)
431  IF (ii < 0) THEN
432  rs%py(ia) = ii + rs%npts_local(2) + 1
433  folded = .true.
434  ELSEIF (ii >= rs%npts_local(2)) THEN
435  rs%py(ia) = ii - rs%npts_local(2) + 1
436  folded = .true.
437  ELSE
438  rs%py(ia) = ii + 1
439  END IF
440  END DO
441  DO i = rhos%pw_grid%bounds(1, 3), rhos%pw_grid%bounds(2, 3)
442  ia = i - rhos%pw_grid%bounds(1, 3) + 1
443  ii = center(3) + i - rs%lb_local(3)
444  IF (ii < 0) THEN
445  rs%pz(ia) = ii + rs%npts_local(3) + 1
446  folded = .true.
447  ELSEIF (ii >= rs%npts_local(3)) THEN
448  rs%pz(ia) = ii - rs%npts_local(3) + 1
449  folded = .true.
450  ELSE
451  rs%pz(ia) = ii + 1
452  END IF
453  END DO
454 
455  IF (folded) THEN
456  CALL dg_add_patch(rs%r, rhos%array, rhos%pw_grid%npts, &
457  rs%px, rs%py, rs%pz)
458  ELSE
459  nc(1) = rs%px(1) - 1
460  nc(2) = rs%py(1) - 1
461  nc(3) = rs%pz(1) - 1
462  CALL dg_add_patch(rs%r, rhos%array, rhos%pw_grid%npts, nc)
463  END IF
464 
465  END SUBROUTINE dg_sum_patch_coef
466 
467 ! **************************************************************************************************
468 !> \brief ...
469 !> \param rs ...
470 !> \param rhos ...
471 !> \param center ...
472 ! **************************************************************************************************
473  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  ns(1) = SIZE(rhos, 1)
484  ns(2) = SIZE(rhos, 2)
485  ns(3) = SIZE(rhos, 3)
486  lb = -(ns - 1)/2
487  ub = lb + ns - 1
488  folded = .false.
489 
490  DO i = lb(1), ub(1)
491  ia = i - lb(1) + 1
492  ii = center(1) + i - rs%lb_local(1)
493  IF (ii < 0) THEN
494  rs%px(ia) = ii + rs%npts_local(1) + 1
495  folded = .true.
496  ELSEIF (ii >= rs%npts_local(1)) THEN
497  rs%px(ia) = ii - rs%npts_local(1) + 1
498  folded = .true.
499  ELSE
500  rs%px(ia) = ii + 1
501  END IF
502  END DO
503  DO i = lb(2), ub(2)
504  ia = i - lb(2) + 1
505  ii = center(2) + i - rs%lb_local(2)
506  IF (ii < 0) THEN
507  rs%py(ia) = ii + rs%npts_local(2) + 1
508  folded = .true.
509  ELSEIF (ii >= rs%npts_local(2)) THEN
510  rs%py(ia) = ii - rs%npts_local(2) + 1
511  folded = .true.
512  ELSE
513  rs%py(ia) = ii + 1
514  END IF
515  END DO
516  DO i = lb(3), ub(3)
517  ia = i - lb(3) + 1
518  ii = center(3) + i - rs%lb_local(3)
519  IF (ii < 0) THEN
520  rs%pz(ia) = ii + rs%npts_local(3) + 1
521  folded = .true.
522  ELSEIF (ii >= rs%npts_local(3)) THEN
523  rs%pz(ia) = ii - rs%npts_local(3) + 1
524  folded = .true.
525  ELSE
526  rs%pz(ia) = ii + 1
527  END IF
528  END DO
529 
530  IF (folded) THEN
531  CALL dg_add_patch(rs%r, rhos, ns, rs%px, rs%py, rs%pz)
532  ELSE
533  nc(1) = rs%px(1) - 1
534  nc(2) = rs%py(1) - 1
535  nc(3) = rs%pz(1) - 1
536  CALL dg_add_patch(rs%r, rhos, ns, nc)
537  END IF
538 
539  END SUBROUTINE dg_sum_patch_arr
540 
541 ! **************************************************************************************************
542 !> \brief ...
543 !> \param drpot ...
544 !> \param rhos ...
545 !> \param center ...
546 !> \param force ...
547 ! **************************************************************************************************
548  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  ns(1) = SIZE(rhos, 1)
561  ns(2) = SIZE(rhos, 2)
562  ns(3) = SIZE(rhos, 3)
563  lb = -(ns - 1)/2
564  ub = lb + ns - 1
565  folded = .false.
566 
567  DO i = lb(1), ub(1)
568  ia = i - lb(1) + 1
569  ii = center(1) + i - drpot(1)%lb_local(1)
570  IF (ii < 0) THEN
571  drpot(1)%px(ia) = ii + drpot(1)%npts_local(1) + 1
572  folded = .true.
573  ELSEIF (ii >= drpot(1)%npts_local(1)) THEN
574  drpot(1)%px(ia) = ii - drpot(1)%npts_local(1) + 1
575  folded = .true.
576  ELSE
577  drpot(1)%px(ia) = ii + 1
578  END IF
579  END DO
580  DO i = lb(2), ub(2)
581  ia = i - lb(2) + 1
582  ii = center(2) + i - drpot(1)%lb_local(2)
583  IF (ii < 0) THEN
584  drpot(1)%py(ia) = ii + drpot(1)%npts_local(2) + 1
585  folded = .true.
586  ELSEIF (ii >= drpot(1)%npts_local(2)) THEN
587  drpot(1)%py(ia) = ii - drpot(1)%npts_local(2) + 1
588  folded = .true.
589  ELSE
590  drpot(1)%py(ia) = ii + 1
591  END IF
592  END DO
593  DO i = lb(3), ub(3)
594  ia = i - lb(3) + 1
595  ii = center(3) + i - drpot(1)%lb_local(3)
596  IF (ii < 0) THEN
597  drpot(1)%pz(ia) = ii + drpot(1)%npts_local(3) + 1
598  folded = .true.
599  ELSEIF (ii >= drpot(1)%npts_local(3)) THEN
600  drpot(1)%pz(ia) = ii - drpot(1)%npts_local(3) + 1
601  folded = .true.
602  ELSE
603  drpot(1)%pz(ia) = ii + 1
604  END IF
605  END DO
606 
607  IF (folded) THEN
608  CALL dg_int_patch_3d(drpot(1)%r, drpot(2)%r, &
609  drpot(3)%r, rhos, force, ns, &
610  drpot(1)%px, drpot(1)%py, drpot(1)%pz)
611  ELSE
612  nc(1) = drpot(1)%px(1) - 1
613  nc(2) = drpot(1)%py(1) - 1
614  nc(3) = drpot(1)%pz(1) - 1
615  CALL dg_int_patch_3d(drpot(1)%r, drpot(2)%r, &
616  drpot(3)%r, rhos, force, ns, nc)
617  END IF
618 
619  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  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  ns(1) = SIZE(rhos, 1)
640  ns(2) = SIZE(rhos, 2)
641  ns(3) = SIZE(rhos, 3)
642  lb = -(ns - 1)/2
643  ub = lb + ns - 1
644  folded = .false.
645 
646  DO i = lb(1), ub(1)
647  ia = i - lb(1) + 1
648  ii = center(1) + i - drpot%lb_local(1)
649  IF (ii < 0) THEN
650  drpot%px(ia) = ii + drpot%npts_local(1) + 1
651  folded = .true.
652  ELSEIF (ii >= drpot%desc%npts(1)) THEN
653  drpot%px(ia) = ii - drpot%npts_local(1) + 1
654  folded = .true.
655  ELSE
656  drpot%px(ia) = ii + 1
657  END IF
658  END DO
659  DO i = lb(2), ub(2)
660  ia = i - lb(2) + 1
661  ii = center(2) + i - drpot%lb_local(2)
662  IF (ii < 0) THEN
663  drpot%py(ia) = ii + drpot%npts_local(2) + 1
664  folded = .true.
665  ELSEIF (ii >= drpot%desc%npts(2)) THEN
666  drpot%py(ia) = ii - drpot%npts_local(2) + 1
667  folded = .true.
668  ELSE
669  drpot%py(ia) = ii + 1
670  END IF
671  END DO
672  DO i = lb(3), ub(3)
673  ia = i - lb(3) + 1
674  ii = center(3) + i - drpot%lb_local(3)
675  IF (ii < 0) THEN
676  drpot%pz(ia) = ii + drpot%npts_local(3) + 1
677  folded = .true.
678  ELSEIF (ii >= drpot%desc%npts(3)) THEN
679  drpot%pz(ia) = ii - drpot%npts_local(3) + 1
680  folded = .true.
681  ELSE
682  drpot%pz(ia) = ii + 1
683  END IF
684  END DO
685 
686  IF (folded) THEN
687  CALL dg_int_patch_1d(drpot%r, rhos, force, ns, &
688  drpot%px, drpot%py, drpot%pz)
689  ELSE
690  nc(1) = drpot%px(1) - 1
691  nc(2) = drpot%py(1) - 1
692  nc(3) = drpot%pz(1) - 1
693  CALL dg_int_patch_1d(drpot%r, rhos, force, ns, nc)
694  END IF
695 
696  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  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  folded = .false.
718 
719  DO i = rhos%pw_grid%bounds(1, 1), rhos%pw_grid%bounds(2, 1)
720  ia = i - rhos%pw_grid%bounds(1, 1) + 1
721  ii = center(1) + i - drpot(1)%lb_local(1)
722  IF (ii < 0) THEN
723  drpot(1)%px(ia) = ii + drpot(1)%desc%npts(1) + 1
724  folded = .true.
725  ELSEIF (ii >= drpot(1)%desc%npts(1)) THEN
726  drpot(1)%px(ia) = ii - drpot(1)%desc%npts(1) + 1
727  folded = .true.
728  ELSE
729  drpot(1)%px(ia) = ii + 1
730  END IF
731  END DO
732  DO i = rhos%pw_grid%bounds(1, 2), rhos%pw_grid%bounds(2, 2)
733  ia = i - rhos%pw_grid%bounds(1, 2) + 1
734  ii = center(2) + i - drpot(1)%lb_local(2)
735  IF (ii < 0) THEN
736  drpot(1)%py(ia) = ii + drpot(1)%desc%npts(2) + 1
737  folded = .true.
738  ELSEIF (ii >= drpot(1)%desc%npts(2)) THEN
739  drpot(1)%py(ia) = ii - drpot(1)%desc%npts(2) + 1
740  folded = .true.
741  ELSE
742  drpot(1)%py(ia) = ii + 1
743  END IF
744  END DO
745  DO i = rhos%pw_grid%bounds(1, 3), rhos%pw_grid%bounds(2, 3)
746  ia = i - rhos%pw_grid%bounds(1, 3) + 1
747  ii = center(3) + i - drpot(1)%lb_local(3)
748  IF (ii < 0) THEN
749  drpot(1)%pz(ia) = ii + drpot(1)%desc%npts(3) + 1
750  folded = .true.
751  ELSEIF (ii >= drpot(1)%desc%npts(3)) THEN
752  drpot(1)%pz(ia) = ii - drpot(1)%desc%npts(3) + 1
753  folded = .true.
754  ELSE
755  drpot(1)%pz(ia) = ii + 1
756  END IF
757  END DO
758 
759  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  drpot(1)%px, drpot(1)%py, drpot(1)%pz)
763  ELSE
764  nc(1) = drpot(1)%px(1) - 1
765  nc(2) = drpot(1)%py(1) - 1
766  nc(3) = drpot(1)%pz(1) - 1
767  CALL dg_int_patch_3d(drpot(1)%r, drpot(2)%r, &
768  drpot(3)%r, rhos%array, force, rhos%pw_grid%npts, nc)
769  END IF
770 
771  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  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  folded = .false.
792 
793  DO i = rhos%pw_grid%bounds(1, 1), rhos%pw_grid%bounds(2, 1)
794  ia = i - rhos%pw_grid%bounds(1, 1) + 1
795  ii = center(1) + i - drpot%lb_local(1)
796  IF (ii < 0) THEN
797  drpot%px(ia) = ii + drpot%desc%npts(1) + 1
798  folded = .true.
799  ELSEIF (ii >= drpot%desc%npts(1)) THEN
800  drpot%px(ia) = ii - drpot%desc%npts(1) + 1
801  folded = .true.
802  ELSE
803  drpot%px(ia) = ii + 1
804  END IF
805  END DO
806  DO i = rhos%pw_grid%bounds(1, 2), rhos%pw_grid%bounds(2, 2)
807  ia = i - rhos%pw_grid%bounds(1, 2) + 1
808  ii = center(2) + i - drpot%lb_local(2)
809  IF (ii < 0) THEN
810  drpot%py(ia) = ii + drpot%desc%npts(2) + 1
811  folded = .true.
812  ELSEIF (ii >= drpot%desc%npts(2)) THEN
813  drpot%py(ia) = ii - drpot%desc%npts(2) + 1
814  folded = .true.
815  ELSE
816  drpot%py(ia) = ii + 1
817  END IF
818  END DO
819  DO i = rhos%pw_grid%bounds(1, 3), rhos%pw_grid%bounds(2, 3)
820  ia = i - rhos%pw_grid%bounds(1, 3) + 1
821  ii = center(3) + i - drpot%lb_local(3)
822  IF (ii < 0) THEN
823  drpot%pz(ia) = ii + drpot%desc%npts(3) + 1
824  folded = .true.
825  ELSEIF (ii >= drpot%desc%npts(3)) THEN
826  drpot%pz(ia) = ii - drpot%desc%npts(3) + 1
827  folded = .true.
828  ELSE
829  drpot%pz(ia) = ii + 1
830  END IF
831  END DO
832 
833  IF (folded) THEN
834  CALL dg_int_patch_1d(drpot%r, rhos%array, force, &
835  rhos%pw_grid%npts, drpot%px, drpot%py, drpot%pz)
836  ELSE
837  nc(1) = drpot%px(1) - 1
838  nc(2) = drpot%py(1) - 1
839  nc(3) = drpot%pz(1) - 1
840  CALL dg_int_patch_1d(drpot%r, rhos%array, force, rhos%pw_grid%npts, nc)
841  END IF
842 
843  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  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  nd = rhos1%pw_grid%npts
867 
868  ALLOCATE (zs(nd(1)*nd(2)))
869  zs = 0.0_dp
870  CALL cd%create(rho0%pw_grid)
871  CALL pw_zero(cd)
872 
873  za = cmplx(0.0_dp, 0.0_dp, kind=dp)
874  zb = cmplx(charge1, 0.0_dp, kind=dp)
875  CALL rankup(nd, za, cd%array, zb, ex1, ey1, ez1, zs)
876  CALL pw_multiply_with(cd, rho0)
877  CALL fft3d(bwfft, nd, cd%array)
878  CALL pw_copy(cd, rhos1)
879 
880  DEALLOCATE (zs)
881  CALL cd%release()
882 
883  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  SUBROUTINE dg_get_patch_2(rho0, rhos1, rhos2, charge1, charge2, &
900  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  nd = rhos1%pw_grid%npts
913 
914  ALLOCATE (zs(nd(1)*nd(2)))
915  zs = 0.0_dp
916  CALL cd%create(rhos1%pw_grid)
917  CALL pw_zero(cd)
918 
919  za = cmplx(0.0_dp, 0.0_dp, kind=dp)
920  zb = cmplx(charge2, 0.0_dp, kind=dp)
921  CALL rankup(nd, za, cd%array, zb, ex2, ey2, ez2, zs)
922  za = cmplx(0.0_dp, 1.0_dp, kind=dp)
923  zb = cmplx(charge1, 0.0_dp, kind=dp)
924  CALL rankup(nd, za, cd%array, zb, ex1, ey1, ez1, zs)
925  CALL pw_multiply_with(cd, rho0)
926  CALL fft3d(bwfft, nd, cd%array)
927  CALL copy_cri(cd%array, rhos1%array, rhos2%array)
928 
929  DEALLOCATE (zs)
930  CALL cd%release()
931 
932  END SUBROUTINE dg_get_patch_2
933 
934 ! **************************************************************************************************
935 !> \brief ...
936 !> \param rb ...
937 !> \param rs ...
938 !> \param ns ...
939 !> \param nc ...
940 ! **************************************************************************************************
941  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  DO k = 1, ns(3)
950  kk = nc(3) + k
951  DO j = 1, ns(2)
952  jj = nc(2) + j
953  DO i = 1, ns(1)
954  ii = nc(1) + i
955  rb(ii, jj, kk) = rb(ii, jj, kk) + rs(i, j, k)
956  END DO
957  END DO
958  END DO
959 
960  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  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  DO k = 1, ns(3)
981  kk = pz(k)
982  DO j = 1, ns(2)
983  jj = py(j)
984  DO i = 1, ns(1)
985  ii = px(i)
986  rb(ii, jj, kk) = rb(ii, jj, kk) + rs(i, j, k)
987  END DO
988  END DO
989  END DO
990 
991  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  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  f = 0.0_dp
1013  DO k = 1, ns(3)
1014  kk = nc(3) + k
1015  DO j = 1, ns(2)
1016  jj = nc(2) + j
1017  DO i = 1, ns(1)
1018  ii = nc(1) + i
1019  s = rs(i, j, k)
1020  f(1) = f(1) + s*rbx(ii, jj, kk)
1021  f(2) = f(2) + s*rby(ii, jj, kk)
1022  f(3) = f(3) + s*rbz(ii, jj, kk)
1023  END DO
1024  END DO
1025  END DO
1026 
1027  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  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  f = 0.0_dp
1047  DO k = 1, ns(3)
1048  kk = nc(3) + k
1049  DO j = 1, ns(2)
1050  jj = nc(2) + j
1051  DO i = 1, ns(1)
1052  ii = nc(1) + i
1053  s = rs(i, j, k)
1054  f = f + s*rb(ii, jj, kk)
1055  END DO
1056  END DO
1057  END DO
1058 
1059  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  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  f = 0.0_dp
1084  DO k = 1, ns(3)
1085  kk = pz(k)
1086  DO j = 1, ns(2)
1087  jj = py(j)
1088  DO i = 1, ns(1)
1089  ii = px(i)
1090  s = rs(i, j, k)
1091  f(1) = f(1) + s*rbx(ii, jj, kk)
1092  f(2) = f(2) + s*rby(ii, jj, kk)
1093  f(3) = f(3) + s*rbz(ii, jj, kk)
1094  END DO
1095  END DO
1096  END DO
1097 
1098  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  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  f = 0.0_dp
1121  DO k = 1, ns(3)
1122  kk = pz(k)
1123  DO j = 1, ns(2)
1124  jj = py(j)
1125  DO i = 1, ns(1)
1126  ii = px(i)
1127  s = rs(i, j, k)
1128  f = f + s*rb(ii, jj, kk)
1129  END DO
1130  END DO
1131  END DO
1132 
1133  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  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  n2 = n(1)*n(2)
1162  n3 = n2*n(3)
1163  scr(1:n2) = z_zero
1164  CALL zgeru(n(1), n(2), zb, ex, 1, ey, 1, scr, n(1))
1165  CALL zscal(n3, za, cmat, 1)
1166  CALL zgeru(n2, n(3), z_one, scr, 1, ez, 1, cmat, n2)
1167 
1168  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  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 !$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  END SUBROUTINE copy_cri
1191 
1192 END MODULE dgs
det_3x3
real(kind=dp) function det_3x3(a)
...
Definition: dumpdcd.F:1091
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition: grid_common.h:117
dgs
Definition: dgs.F:13
dgs::dg_grid_change
subroutine, public dg_grid_change(b_cell_hmat, grid_b, grid_s)
...
Definition: dgs.F:244
dgs::dg_pme_grid_setup
subroutine, public dg_pme_grid_setup(b_cell_hmat, npts_s, cutoff_radius, grid_s, grid_b, grid_ref, rs_dims, iounit, fft_usage)
...
Definition: dgs.F:94
dgs::dg_get_strucfac
subroutine, public dg_get_strucfac(cell_hmat, r, npts_s, npts_b, centre, lb, ex, ey, ez)
...
Definition: dgs.F:332
fft_tools
Definition: fft_tools.F:28
fft_tools::fft_radix_operations
subroutine, public fft_radix_operations(radix_in, radix_out, operation)
Determine the allowed lengths of FFT's '''.
Definition: fft_tools.F:239
fft_tools::bwfft
integer, parameter, public bwfft
Definition: fft_tools.F:145
fft_tools::fft_radix_closest
integer, parameter, public fft_radix_closest
Definition: fft_tools.F:146
fft_tools::fft_radix_next_odd
integer, parameter, public fft_radix_next_odd
Definition: fft_tools.F:148
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::z_one
complex(kind=dp), parameter, public z_one
Definition: mathconstants.F:143
mathconstants::z_zero
complex(kind=dp), parameter, public z_zero
Definition: mathconstants.F:143
mathlib
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
mathlib::inv_3x3
pure real(kind=dp) function, dimension(3, 3), public inv_3x3(a)
Returns the inverse of the 3 x 3 matrix a.
Definition: mathlib.F:516
pw_grid_info
This module returns additional info on PW grids.
Definition: pw_grid_info.F:14
pw_grid_info::pw_find_cutoff
real(kind=dp) function, public pw_find_cutoff(npts, h_inv)
Given a grid and a box, calculate the corresponding cutoff *** This routine calculates the cutoff in ...
Definition: pw_grid_info.F:272
pw_grid_types
Definition: pw_grid_types.F:13
pw_grid_types::halfspace
integer, parameter, public halfspace
Definition: pw_grid_types.F:28
pw_grids
This module defines the grid data type and some basic operations on it.
Definition: pw_grids.F:36
pw_grids::pw_grid_setup
subroutine, public pw_grid_setup(cell_hmat, pw_grid, grid_span, cutoff, bounds, bounds_local, npts, spherical, odd, fft_usage, ncommensurate, icommensurate, blocked, ref_grid, rs_dims, iounit)
sets up a pw_grid
Definition: pw_grids.F:286
pw_grids::pw_grid_change
subroutine, public pw_grid_change(cell_hmat, pw_grid)
Recalculate the g-vectors after a change of the box.
Definition: pw_grids.F:2068
pw_methods
Definition: pw_methods.F:25
pw_types
Definition: pw_types.F:24
realspace_grid_types
Definition: realspace_grid_types.F:18
structure_factors
Definition: structure_factors.F:12
structure_factors::structure_factor_evaluate
subroutine, public structure_factor_evaluate(delta, lb, ex, ey, ez)
...
Definition: structure_factors.F:151