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