(git:374b731)
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-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! **************************************************************************************************
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 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
78CONTAINS
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
1192END 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:117
dgs::dg_get_patch
Definition dgs.F:62
dgs::dg_sum_patch_force_1d
Definition dgs.F:58
dgs::dg_sum_patch_force_3d
Definition dgs.F:54
dgs::dg_sum_patch
Definition dgs.F:50
fft_tools::fft3d
Definition fft_tools.F:165
mathlib::det_3x3
Definition mathlib.F:65
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_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
pw_grid_types::pw_grid_type
Definition pw_grid_types.F:62
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:137