(git:b31c023)
Loading...
Searching...
No Matches
dgs.F
Go to the documentation of this file.
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! **************************************************************************************************
13MODULE 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
81CONTAINS
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 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 CALL dg_find_cutoff(b_cell_hmat, npts_s, cutoff_radius, bounds_s, bounds_b, cutoff)
115
116 ecut = 0.5_dp*cutoff*cutoff
117 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 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 rs_dims=rs_dims, iounit=iounit, fft_usage=fft_usage)
123 END IF
124
125 CALL dg_find_basis(grid_b%npts, b_cell_hmat, unit_cell_hmat)
126
127 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 cutoff=ecut, iounit=iounit, fft_usage=fft_usage)
131
132 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 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 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 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 cell_hmat(3, 3)*cell_hmat(3, 3))
154 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 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 b_cell_deth = abs(det_3x3(b_cell_hmat))
179 IF (b_cell_deth < 1.0e-10_dp) THEN
180 CALL cp_abort(__location__, &
181 "An invalid set of cell vectors was specified. "// &
182 "The determinant det(h) is too small")
183 END IF
184 b_cell_h_inv = inv_3x3(b_cell_hmat)
185
186 CALL fft_radix_operations(npts_s(1), nout(1), &
187 operation=fft_radix_next_odd)
188 CALL fft_radix_operations(npts_s(1), nout(2), &
189 operation=fft_radix_next_odd)
190 CALL fft_radix_operations(npts_s(1), nout(3), &
191 operation=fft_radix_next_odd)
192
193 cell_lengths = get_cell_lengths(b_cell_hmat)
194
195 CALL dg_get_spacing(nout, cutoff_radius, dr)
196 CALL dg_find_radix(dr, cell_lengths, npts_b)
197
198! In-line code to set grid_b % npts = npts_s if necessary
199 IF (nout(1) > npts_b(1)) THEN
200 npts_b(1) = nout(1)
201 dr(1) = cell_lengths(1)/real(nout(1), kind=dp)
202 END IF
203 IF (nout(2) > npts_b(2)) THEN
204 npts_b(2) = nout(2)
205 dr(2) = cell_lengths(2)/real(nout(2), kind=dp)
206 END IF
207 IF (nout(3) > npts_b(3)) THEN
208 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 bounds_b(1, :) = -npts_b/2
214 bounds_b(2, :) = +(npts_b - 1)/2
215! small grid bounds
216 bounds_s(1, :) = -nout(:)/2
217 bounds_s(2, :) = (+nout(:) - 1)/2
218
219 cutoff = pw_find_cutoff(npts_b, b_cell_h_inv)
220
221 END SUBROUTINE dg_find_cutoff
222
223! **************************************************************************************************
224!> \brief ...
225!> \param npts ...
226!> \param cutoff_radius ...
227!> \param dr ...
228! **************************************************************************************************
229 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 dr(:) = cutoff_radius/(real(npts(:), kind=dp)/2.0_dp)
236
237 END SUBROUTINE dg_get_spacing
238
239! **************************************************************************************************
240!> \brief ...
241!> \param b_cell_hmat ...
242!> \param grid_b ...
243!> \param grid_s ...
244! **************************************************************************************************
245 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 CALL dg_find_basis(grid_b%npts, b_cell_hmat, unit_cell_hmat)
252 CALL dg_set_cell(grid_s%npts, unit_cell_hmat, s_cell_hmat)
253 CALL pw_grid_change(s_cell_hmat, grid_s)
254 END SUBROUTINE dg_grid_change
255
256! **************************************************************************************************
257!> \brief ...
258!> \param dr ...
259!> \param cell_lengths ...
260!> \param npts ...
261! **************************************************************************************************
262 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 nin(:) = nint(cell_lengths(:)/dr(:))
271 CALL fft_radix_operations(nin(1), npts(1), &
272 operation=fft_radix_closest)
273 CALL fft_radix_operations(nin(2), npts(2), &
274 operation=fft_radix_closest)
275 CALL fft_radix_operations(nin(3), npts(3), &
276 operation=fft_radix_closest)
277 dr(:) = cell_lengths(:)/real(npts(:), kind=dp)
278
279 END SUBROUTINE dg_find_radix
280
281! **************************************************************************************************
282!> \brief ...
283!> \param npts ...
284!> \param cell_hmat ...
285!> \param unit_cell_hmat ...
286! **************************************************************************************************
287 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 DO i = 1, 3
295 unit_cell_hmat(:, i) = cell_hmat(:, i)/real(npts(:), kind=dp)
296 END DO
297
298 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 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 cell_hmat(:, 1) = unit_cell_hmat(:, 1)*npts(1)
316 cell_hmat(:, 2) = unit_cell_hmat(:, 2)*npts(2)
317 cell_hmat(:, 3) = unit_cell_hmat(:, 3)*npts(3)
318
319 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 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 CALL dg_get_delta(cell_hmat, r, npts_s, npts_b, centre, delta)
346
347 CALL structure_factor_evaluate(delta, lb, ex, ey, ez)
348
349 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 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 cell_deth = abs(det_3x3(cell_hmat))
372 IF (cell_deth < 1.0e-10_dp) THEN
373 CALL cp_abort(__location__, &
374 "An invalid set of cell vectors was specified. "// &
375 "The determinant det(h) is too small")
376 END IF
377
378 cell_h_inv = inv_3x3(cell_hmat)
379
380! compute the scaled coordinate of atomi
381 s = matmul(cell_h_inv, r)
382 s = s - nint(s)
383
384! find the continuous ``grid'' point (on big grid)
385 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 centre(:) = nint(grid_i(:))
389
390! find the distance vector
391 delta(:) = (grid_i(:) - centre(:))/real(npts_s(:), kind=dp)
392
393 centre(:) = centre(:) + npts_b(:)/2
394 centre(:) = modulo(centre(:), npts_b(:))
395 centre(:) = centre(:) - npts_b(:)/2
396
397 END SUBROUTINE dg_get_delta
398
399! **************************************************************************************************
400!> \brief ...
401!> \param rs ...
402!> \param rhos ...
403!> \param center ...
404! **************************************************************************************************
405 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 folded = .false.
416
417 DO i = rhos%pw_grid%bounds(1, 1), rhos%pw_grid%bounds(2, 1)
418 ia = i - rhos%pw_grid%bounds(1, 1) + 1
419 ii = center(1) + i - rs%lb_local(1)
420 IF (ii < 0) THEN
421 rs%px(ia) = ii + rs%npts_local(1) + 1
422 folded = .true.
423 ELSEIF (ii >= rs%npts_local(1)) THEN
424 rs%px(ia) = ii - rs%npts_local(1) + 1
425 folded = .true.
426 ELSE
427 rs%px(ia) = ii + 1
428 END IF
429 END DO
430 DO i = rhos%pw_grid%bounds(1, 2), rhos%pw_grid%bounds(2, 2)
431 ia = i - rhos%pw_grid%bounds(1, 2) + 1
432 ii = center(2) + i - rs%lb_local(2)
433 IF (ii < 0) THEN
434 rs%py(ia) = ii + rs%npts_local(2) + 1
435 folded = .true.
436 ELSEIF (ii >= rs%npts_local(2)) THEN
437 rs%py(ia) = ii - rs%npts_local(2) + 1
438 folded = .true.
439 ELSE
440 rs%py(ia) = ii + 1
441 END IF
442 END DO
443 DO i = rhos%pw_grid%bounds(1, 3), rhos%pw_grid%bounds(2, 3)
444 ia = i - rhos%pw_grid%bounds(1, 3) + 1
445 ii = center(3) + i - rs%lb_local(3)
446 IF (ii < 0) THEN
447 rs%pz(ia) = ii + rs%npts_local(3) + 1
448 folded = .true.
449 ELSEIF (ii >= rs%npts_local(3)) THEN
450 rs%pz(ia) = ii - rs%npts_local(3) + 1
451 folded = .true.
452 ELSE
453 rs%pz(ia) = ii + 1
454 END IF
455 END DO
456
457 IF (folded) THEN
458 CALL dg_add_patch(rs%r, rhos%array, rhos%pw_grid%npts, &
459 rs%px, rs%py, rs%pz)
460 ELSE
461 nc(1) = rs%px(1) - 1
462 nc(2) = rs%py(1) - 1
463 nc(3) = rs%pz(1) - 1
464 CALL dg_add_patch(rs%r, rhos%array, rhos%pw_grid%npts, nc)
465 END IF
466
467 END SUBROUTINE dg_sum_patch_coef
468
469! **************************************************************************************************
470!> \brief ...
471!> \param rs ...
472!> \param rhos ...
473!> \param center ...
474! **************************************************************************************************
475 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 ns(1) = SIZE(rhos, 1)
486 ns(2) = SIZE(rhos, 2)
487 ns(3) = SIZE(rhos, 3)
488 lb = -(ns - 1)/2
489 ub = lb + ns - 1
490 folded = .false.
491
492 DO i = lb(1), ub(1)
493 ia = i - lb(1) + 1
494 ii = center(1) + i - rs%lb_local(1)
495 IF (ii < 0) THEN
496 rs%px(ia) = ii + rs%npts_local(1) + 1
497 folded = .true.
498 ELSEIF (ii >= rs%npts_local(1)) THEN
499 rs%px(ia) = ii - rs%npts_local(1) + 1
500 folded = .true.
501 ELSE
502 rs%px(ia) = ii + 1
503 END IF
504 END DO
505 DO i = lb(2), ub(2)
506 ia = i - lb(2) + 1
507 ii = center(2) + i - rs%lb_local(2)
508 IF (ii < 0) THEN
509 rs%py(ia) = ii + rs%npts_local(2) + 1
510 folded = .true.
511 ELSEIF (ii >= rs%npts_local(2)) THEN
512 rs%py(ia) = ii - rs%npts_local(2) + 1
513 folded = .true.
514 ELSE
515 rs%py(ia) = ii + 1
516 END IF
517 END DO
518 DO i = lb(3), ub(3)
519 ia = i - lb(3) + 1
520 ii = center(3) + i - rs%lb_local(3)
521 IF (ii < 0) THEN
522 rs%pz(ia) = ii + rs%npts_local(3) + 1
523 folded = .true.
524 ELSEIF (ii >= rs%npts_local(3)) THEN
525 rs%pz(ia) = ii - rs%npts_local(3) + 1
526 folded = .true.
527 ELSE
528 rs%pz(ia) = ii + 1
529 END IF
530 END DO
531
532 IF (folded) THEN
533 CALL dg_add_patch(rs%r, rhos, ns, rs%px, rs%py, rs%pz)
534 ELSE
535 nc(1) = rs%px(1) - 1
536 nc(2) = rs%py(1) - 1
537 nc(3) = rs%pz(1) - 1
538 CALL dg_add_patch(rs%r, rhos, ns, nc)
539 END IF
540
541 END SUBROUTINE dg_sum_patch_arr
542
543! **************************************************************************************************
544!> \brief ...
545!> \param drpot ...
546!> \param rhos ...
547!> \param center ...
548!> \param force ...
549! **************************************************************************************************
550 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 ns(1) = SIZE(rhos, 1)
563 ns(2) = SIZE(rhos, 2)
564 ns(3) = SIZE(rhos, 3)
565 lb = -(ns - 1)/2
566 ub = lb + ns - 1
567 folded = .false.
568
569 DO i = lb(1), ub(1)
570 ia = i - lb(1) + 1
571 ii = center(1) + i - drpot(1)%lb_local(1)
572 IF (ii < 0) THEN
573 drpot(1)%px(ia) = ii + drpot(1)%npts_local(1) + 1
574 folded = .true.
575 ELSEIF (ii >= drpot(1)%npts_local(1)) THEN
576 drpot(1)%px(ia) = ii - drpot(1)%npts_local(1) + 1
577 folded = .true.
578 ELSE
579 drpot(1)%px(ia) = ii + 1
580 END IF
581 END DO
582 DO i = lb(2), ub(2)
583 ia = i - lb(2) + 1
584 ii = center(2) + i - drpot(1)%lb_local(2)
585 IF (ii < 0) THEN
586 drpot(1)%py(ia) = ii + drpot(1)%npts_local(2) + 1
587 folded = .true.
588 ELSEIF (ii >= drpot(1)%npts_local(2)) THEN
589 drpot(1)%py(ia) = ii - drpot(1)%npts_local(2) + 1
590 folded = .true.
591 ELSE
592 drpot(1)%py(ia) = ii + 1
593 END IF
594 END DO
595 DO i = lb(3), ub(3)
596 ia = i - lb(3) + 1
597 ii = center(3) + i - drpot(1)%lb_local(3)
598 IF (ii < 0) THEN
599 drpot(1)%pz(ia) = ii + drpot(1)%npts_local(3) + 1
600 folded = .true.
601 ELSEIF (ii >= drpot(1)%npts_local(3)) THEN
602 drpot(1)%pz(ia) = ii - drpot(1)%npts_local(3) + 1
603 folded = .true.
604 ELSE
605 drpot(1)%pz(ia) = ii + 1
606 END IF
607 END DO
608
609 IF (folded) THEN
610 CALL dg_int_patch_3d(drpot(1)%r, drpot(2)%r, &
611 drpot(3)%r, rhos, force, ns, &
612 drpot(1)%px, drpot(1)%py, drpot(1)%pz)
613 ELSE
614 nc(1) = drpot(1)%px(1) - 1
615 nc(2) = drpot(1)%py(1) - 1
616 nc(3) = drpot(1)%pz(1) - 1
617 CALL dg_int_patch_3d(drpot(1)%r, drpot(2)%r, &
618 drpot(3)%r, rhos, force, ns, nc)
619 END IF
620
621 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 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 ns(1) = SIZE(rhos, 1)
642 ns(2) = SIZE(rhos, 2)
643 ns(3) = SIZE(rhos, 3)
644 lb = -(ns - 1)/2
645 ub = lb + ns - 1
646 folded = .false.
647
648 DO i = lb(1), ub(1)
649 ia = i - lb(1) + 1
650 ii = center(1) + i - drpot%lb_local(1)
651 IF (ii < 0) THEN
652 drpot%px(ia) = ii + drpot%npts_local(1) + 1
653 folded = .true.
654 ELSEIF (ii >= drpot%desc%npts(1)) THEN
655 drpot%px(ia) = ii - drpot%npts_local(1) + 1
656 folded = .true.
657 ELSE
658 drpot%px(ia) = ii + 1
659 END IF
660 END DO
661 DO i = lb(2), ub(2)
662 ia = i - lb(2) + 1
663 ii = center(2) + i - drpot%lb_local(2)
664 IF (ii < 0) THEN
665 drpot%py(ia) = ii + drpot%npts_local(2) + 1
666 folded = .true.
667 ELSEIF (ii >= drpot%desc%npts(2)) THEN
668 drpot%py(ia) = ii - drpot%npts_local(2) + 1
669 folded = .true.
670 ELSE
671 drpot%py(ia) = ii + 1
672 END IF
673 END DO
674 DO i = lb(3), ub(3)
675 ia = i - lb(3) + 1
676 ii = center(3) + i - drpot%lb_local(3)
677 IF (ii < 0) THEN
678 drpot%pz(ia) = ii + drpot%npts_local(3) + 1
679 folded = .true.
680 ELSEIF (ii >= drpot%desc%npts(3)) THEN
681 drpot%pz(ia) = ii - drpot%npts_local(3) + 1
682 folded = .true.
683 ELSE
684 drpot%pz(ia) = ii + 1
685 END IF
686 END DO
687
688 IF (folded) THEN
689 CALL dg_int_patch_1d(drpot%r, rhos, force, ns, &
690 drpot%px, drpot%py, drpot%pz)
691 ELSE
692 nc(1) = drpot%px(1) - 1
693 nc(2) = drpot%py(1) - 1
694 nc(3) = drpot%pz(1) - 1
695 CALL dg_int_patch_1d(drpot%r, rhos, force, ns, nc)
696 END IF
697
698 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 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 folded = .false.
720
721 DO i = rhos%pw_grid%bounds(1, 1), rhos%pw_grid%bounds(2, 1)
722 ia = i - rhos%pw_grid%bounds(1, 1) + 1
723 ii = center(1) + i - drpot(1)%lb_local(1)
724 IF (ii < 0) THEN
725 drpot(1)%px(ia) = ii + drpot(1)%desc%npts(1) + 1
726 folded = .true.
727 ELSEIF (ii >= drpot(1)%desc%npts(1)) THEN
728 drpot(1)%px(ia) = ii - drpot(1)%desc%npts(1) + 1
729 folded = .true.
730 ELSE
731 drpot(1)%px(ia) = ii + 1
732 END IF
733 END DO
734 DO i = rhos%pw_grid%bounds(1, 2), rhos%pw_grid%bounds(2, 2)
735 ia = i - rhos%pw_grid%bounds(1, 2) + 1
736 ii = center(2) + i - drpot(1)%lb_local(2)
737 IF (ii < 0) THEN
738 drpot(1)%py(ia) = ii + drpot(1)%desc%npts(2) + 1
739 folded = .true.
740 ELSEIF (ii >= drpot(1)%desc%npts(2)) THEN
741 drpot(1)%py(ia) = ii - drpot(1)%desc%npts(2) + 1
742 folded = .true.
743 ELSE
744 drpot(1)%py(ia) = ii + 1
745 END IF
746 END DO
747 DO i = rhos%pw_grid%bounds(1, 3), rhos%pw_grid%bounds(2, 3)
748 ia = i - rhos%pw_grid%bounds(1, 3) + 1
749 ii = center(3) + i - drpot(1)%lb_local(3)
750 IF (ii < 0) THEN
751 drpot(1)%pz(ia) = ii + drpot(1)%desc%npts(3) + 1
752 folded = .true.
753 ELSEIF (ii >= drpot(1)%desc%npts(3)) THEN
754 drpot(1)%pz(ia) = ii - drpot(1)%desc%npts(3) + 1
755 folded = .true.
756 ELSE
757 drpot(1)%pz(ia) = ii + 1
758 END IF
759 END DO
760
761 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 drpot(1)%px, drpot(1)%py, drpot(1)%pz)
765 ELSE
766 nc(1) = drpot(1)%px(1) - 1
767 nc(2) = drpot(1)%py(1) - 1
768 nc(3) = drpot(1)%pz(1) - 1
769 CALL dg_int_patch_3d(drpot(1)%r, drpot(2)%r, &
770 drpot(3)%r, rhos%array, force, rhos%pw_grid%npts, nc)
771 END IF
772
773 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 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 folded = .false.
794
795 DO i = rhos%pw_grid%bounds(1, 1), rhos%pw_grid%bounds(2, 1)
796 ia = i - rhos%pw_grid%bounds(1, 1) + 1
797 ii = center(1) + i - drpot%lb_local(1)
798 IF (ii < 0) THEN
799 drpot%px(ia) = ii + drpot%desc%npts(1) + 1
800 folded = .true.
801 ELSEIF (ii >= drpot%desc%npts(1)) THEN
802 drpot%px(ia) = ii - drpot%desc%npts(1) + 1
803 folded = .true.
804 ELSE
805 drpot%px(ia) = ii + 1
806 END IF
807 END DO
808 DO i = rhos%pw_grid%bounds(1, 2), rhos%pw_grid%bounds(2, 2)
809 ia = i - rhos%pw_grid%bounds(1, 2) + 1
810 ii = center(2) + i - drpot%lb_local(2)
811 IF (ii < 0) THEN
812 drpot%py(ia) = ii + drpot%desc%npts(2) + 1
813 folded = .true.
814 ELSEIF (ii >= drpot%desc%npts(2)) THEN
815 drpot%py(ia) = ii - drpot%desc%npts(2) + 1
816 folded = .true.
817 ELSE
818 drpot%py(ia) = ii + 1
819 END IF
820 END DO
821 DO i = rhos%pw_grid%bounds(1, 3), rhos%pw_grid%bounds(2, 3)
822 ia = i - rhos%pw_grid%bounds(1, 3) + 1
823 ii = center(3) + i - drpot%lb_local(3)
824 IF (ii < 0) THEN
825 drpot%pz(ia) = ii + drpot%desc%npts(3) + 1
826 folded = .true.
827 ELSEIF (ii >= drpot%desc%npts(3)) THEN
828 drpot%pz(ia) = ii - drpot%desc%npts(3) + 1
829 folded = .true.
830 ELSE
831 drpot%pz(ia) = ii + 1
832 END IF
833 END DO
834
835 IF (folded) THEN
836 CALL dg_int_patch_1d(drpot%r, rhos%array, force, &
837 rhos%pw_grid%npts, drpot%px, drpot%py, drpot%pz)
838 ELSE
839 nc(1) = drpot%px(1) - 1
840 nc(2) = drpot%py(1) - 1
841 nc(3) = drpot%pz(1) - 1
842 CALL dg_int_patch_1d(drpot%r, rhos%array, force, rhos%pw_grid%npts, nc)
843 END IF
844
845 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 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 nd = rhos1%pw_grid%npts
869
870 ALLOCATE (zs(nd(1)*nd(2)))
871 zs = 0.0_dp
872 CALL cd%create(rho0%pw_grid)
873 CALL pw_zero(cd)
874
875 za = z_zero
876 zb = cmplx(charge1, 0.0_dp, kind=dp)
877 CALL rankup(nd, za, cd%array, zb, ex1, ey1, ez1, zs)
878 CALL pw_multiply_with(cd, rho0)
879 CALL fft3d(bwfft, nd, cd%array)
880 CALL pw_copy(cd, rhos1)
881
882 DEALLOCATE (zs)
883 CALL cd%release()
884
885 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 SUBROUTINE dg_get_patch_2(rho0, rhos1, rhos2, charge1, charge2, &
902 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 nd = rhos1%pw_grid%npts
915
916 ALLOCATE (zs(nd(1)*nd(2)))
917 zs = 0.0_dp
918 CALL cd%create(rhos1%pw_grid)
919 CALL pw_zero(cd)
920
921 za = z_zero
922 zb = cmplx(charge2, 0.0_dp, kind=dp)
923 CALL rankup(nd, za, cd%array, zb, ex2, ey2, ez2, zs)
924 za = gaussi
925 zb = cmplx(charge1, 0.0_dp, kind=dp)
926 CALL rankup(nd, za, cd%array, zb, ex1, ey1, ez1, zs)
927 CALL pw_multiply_with(cd, rho0)
928 CALL fft3d(bwfft, nd, cd%array)
929 CALL copy_cri(cd%array, rhos1%array, rhos2%array)
930
931 DEALLOCATE (zs)
932 CALL cd%release()
933
934 END SUBROUTINE dg_get_patch_2
935
936! **************************************************************************************************
937!> \brief ...
938!> \param rb ...
939!> \param rs ...
940!> \param ns ...
941!> \param nc ...
942! **************************************************************************************************
943 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 DO k = 1, ns(3)
952 kk = nc(3) + k
953 DO j = 1, ns(2)
954 jj = nc(2) + j
955 DO i = 1, ns(1)
956 ii = nc(1) + i
957 rb(ii, jj, kk) = rb(ii, jj, kk) + rs(i, j, k)
958 END DO
959 END DO
960 END DO
961
962 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 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 DO k = 1, ns(3)
983 kk = pz(k)
984 DO j = 1, ns(2)
985 jj = py(j)
986 DO i = 1, ns(1)
987 ii = px(i)
988 rb(ii, jj, kk) = rb(ii, jj, kk) + rs(i, j, k)
989 END DO
990 END DO
991 END DO
992
993 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 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 f = 0.0_dp
1015 DO k = 1, ns(3)
1016 kk = nc(3) + k
1017 DO j = 1, ns(2)
1018 jj = nc(2) + j
1019 DO i = 1, ns(1)
1020 ii = nc(1) + i
1021 s = rs(i, j, k)
1022 f(1) = f(1) + s*rbx(ii, jj, kk)
1023 f(2) = f(2) + s*rby(ii, jj, kk)
1024 f(3) = f(3) + s*rbz(ii, jj, kk)
1025 END DO
1026 END DO
1027 END DO
1028
1029 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 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 f = 0.0_dp
1049 DO k = 1, ns(3)
1050 kk = nc(3) + k
1051 DO j = 1, ns(2)
1052 jj = nc(2) + j
1053 DO i = 1, ns(1)
1054 ii = nc(1) + i
1055 s = rs(i, j, k)
1056 f = f + s*rb(ii, jj, kk)
1057 END DO
1058 END DO
1059 END DO
1060
1061 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 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 f = 0.0_dp
1086 DO k = 1, ns(3)
1087 kk = pz(k)
1088 DO j = 1, ns(2)
1089 jj = py(j)
1090 DO i = 1, ns(1)
1091 ii = px(i)
1092 s = rs(i, j, k)
1093 f(1) = f(1) + s*rbx(ii, jj, kk)
1094 f(2) = f(2) + s*rby(ii, jj, kk)
1095 f(3) = f(3) + s*rbz(ii, jj, kk)
1096 END DO
1097 END DO
1098 END DO
1099
1100 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 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 f = 0.0_dp
1123 DO k = 1, ns(3)
1124 kk = pz(k)
1125 DO j = 1, ns(2)
1126 jj = py(j)
1127 DO i = 1, ns(1)
1128 ii = px(i)
1129 s = rs(i, j, k)
1130 f = f + s*rb(ii, jj, kk)
1131 END DO
1132 END DO
1133 END DO
1134
1135 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 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 n2 = n(1)*n(2)
1164 n3 = n2*n(3)
1165 scr(1:n2) = z_zero
1166 CALL zgeru(n(1), n(2), zb, ex, 1, ey, 1, scr, n(1))
1167 CALL zscal(n3, za, cmat, 1)
1168 CALL zgeru(n2, n(3), z_one, scr, 1, ez, 1, cmat, n2)
1169
1170 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 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!$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 END SUBROUTINE copy_cri
1193
1194END MODULE dgs
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:120
dgs::dg_get_patch
Definition dgs.F:65
dgs::dg_sum_patch_force_1d
Definition dgs.F:61
dgs::dg_sum_patch_force_3d
Definition dgs.F:57
dgs::dg_sum_patch
Definition dgs.F:53
fft_tools::fft3d
Definition fft_tools.F:166
mathlib::det_3x3
Definition mathlib.F:66
pw_grids::pw_grid_create
Definition pw_grids.F:81
pw_methods::pw_copy
Definition pw_methods.F:146
pw_methods::pw_multiply_with
Definition pw_methods.F:203
pw_methods::pw_zero
Definition pw_methods.F:83
dgs
Definition dgs.F:13
dgs::dg_pme_grid_setup
subroutine, public dg_pme_grid_setup(b_cell_hmat, npts_s, cutoff_radius, grid_s, grid_b, mp_comm, grid_ref, rs_dims, iounit, fft_usage)
...
Definition dgs.F:98
dgs::dg_grid_change
subroutine, public dg_grid_change(b_cell_hmat, grid_b, grid_s)
...
Definition dgs.F:246
dgs::dg_get_strucfac
subroutine, public dg_get_strucfac(cell_hmat, r, npts_s, npts_b, centre, lb, ex, ey, ez)
...
Definition dgs.F:334
fft_tools
Definition fft_tools.F:29
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:146
fft_tools::fft_radix_closest
integer, parameter, public fft_radix_closest
Definition fft_tools.F:147
fft_tools::fft_radix_next_odd
integer, parameter, public fft_radix_next_odd
Definition fft_tools.F:149
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::gaussi
complex(kind=dp), parameter, public gaussi
Definition mathconstants.F:142
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
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
message_passing::mp_comm_self
type(mp_comm_type), parameter, public mp_comm_self
Definition message_passing.F:812
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:26
pw_grids
This module defines the grid data type and some basic operations on it.
Definition pw_grids.F:36
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:2098
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
message_passing::mp_comm_type
Definition message_passing.F:189
pw_grid_types::pw_grid_type
Definition pw_grid_types.F:52
pw_types::pw_c3d_rs_type
Definition pw_types.F:88
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
realspace_grid_types::realspace_grid_type
Definition realspace_grid_types.F:136