61#include "../base/base_uses.f90"
78 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'pw_methods'
79 LOGICAL,
PARAMETER,
PRIVATE :: debug_this_module = .false.
84 MODULE PROCEDURE pw_zero_r1d_rs
85 MODULE PROCEDURE pw_zero_r3d_rs
86 MODULE PROCEDURE pw_zero_c1d_rs
87 MODULE PROCEDURE pw_zero_c3d_rs
88 MODULE PROCEDURE pw_zero_r1d_gs
89 MODULE PROCEDURE pw_zero_r3d_gs
90 MODULE PROCEDURE pw_zero_c1d_gs
91 MODULE PROCEDURE pw_zero_c3d_gs
95 MODULE PROCEDURE pw_scale_r1d_rs
96 MODULE PROCEDURE pw_scale_r3d_rs
97 MODULE PROCEDURE pw_scale_c1d_rs
98 MODULE PROCEDURE pw_scale_c3d_rs
99 MODULE PROCEDURE pw_scale_r1d_gs
100 MODULE PROCEDURE pw_scale_r3d_gs
101 MODULE PROCEDURE pw_scale_c1d_gs
102 MODULE PROCEDURE pw_scale_c3d_gs
106 MODULE PROCEDURE pw_write_r1d_rs
107 MODULE PROCEDURE pw_write_r3d_rs
108 MODULE PROCEDURE pw_write_c1d_rs
109 MODULE PROCEDURE pw_write_c3d_rs
110 MODULE PROCEDURE pw_write_r1d_gs
111 MODULE PROCEDURE pw_write_r3d_gs
112 MODULE PROCEDURE pw_write_c1d_gs
113 MODULE PROCEDURE pw_write_c3d_gs
117 MODULE PROCEDURE pw_integrate_function_r1d_rs
118 MODULE PROCEDURE pw_integrate_function_r3d_rs
119 MODULE PROCEDURE pw_integrate_function_c1d_rs
120 MODULE PROCEDURE pw_integrate_function_c3d_rs
121 MODULE PROCEDURE pw_integrate_function_r1d_gs
122 MODULE PROCEDURE pw_integrate_function_r3d_gs
123 MODULE PROCEDURE pw_integrate_function_c1d_gs
124 MODULE PROCEDURE pw_integrate_function_c3d_gs
128 MODULE PROCEDURE pw_set_value_r1d_rs
129 MODULE PROCEDURE pw_zero_r1d_rs
130 MODULE PROCEDURE pw_set_value_r3d_rs
131 MODULE PROCEDURE pw_zero_r3d_rs
132 MODULE PROCEDURE pw_set_value_c1d_rs
133 MODULE PROCEDURE pw_zero_c1d_rs
134 MODULE PROCEDURE pw_set_value_c3d_rs
135 MODULE PROCEDURE pw_zero_c3d_rs
136 MODULE PROCEDURE pw_set_value_r1d_gs
137 MODULE PROCEDURE pw_zero_r1d_gs
138 MODULE PROCEDURE pw_set_value_r3d_gs
139 MODULE PROCEDURE pw_zero_r3d_gs
140 MODULE PROCEDURE pw_set_value_c1d_gs
141 MODULE PROCEDURE pw_zero_c1d_gs
142 MODULE PROCEDURE pw_set_value_c3d_gs
143 MODULE PROCEDURE pw_zero_c3d_gs
147 MODULE PROCEDURE pw_copy_r1d_r1d_rs
148 MODULE PROCEDURE pw_copy_r1d_c1d_rs
149 MODULE PROCEDURE pw_copy_r3d_r3d_rs
150 MODULE PROCEDURE pw_copy_r3d_c3d_rs
151 MODULE PROCEDURE pw_copy_c1d_r1d_rs
152 MODULE PROCEDURE pw_copy_c1d_c1d_rs
153 MODULE PROCEDURE pw_copy_c3d_r3d_rs
154 MODULE PROCEDURE pw_copy_c3d_c3d_rs
155 MODULE PROCEDURE pw_copy_r1d_r1d_gs
156 MODULE PROCEDURE pw_copy_r1d_c1d_gs
157 MODULE PROCEDURE pw_copy_r3d_r3d_gs
158 MODULE PROCEDURE pw_copy_r3d_c3d_gs
159 MODULE PROCEDURE pw_copy_c1d_r1d_gs
160 MODULE PROCEDURE pw_copy_c1d_c1d_gs
161 MODULE PROCEDURE pw_copy_c3d_r3d_gs
162 MODULE PROCEDURE pw_copy_c3d_c3d_gs
166 MODULE PROCEDURE pw_axpy_r1d_r1d_rs
167 MODULE PROCEDURE pw_axpy_r1d_c1d_rs
168 MODULE PROCEDURE pw_axpy_r3d_r3d_rs
169 MODULE PROCEDURE pw_axpy_r3d_c3d_rs
170 MODULE PROCEDURE pw_axpy_c1d_r1d_rs
171 MODULE PROCEDURE pw_axpy_c1d_c1d_rs
172 MODULE PROCEDURE pw_axpy_c3d_r3d_rs
173 MODULE PROCEDURE pw_axpy_c3d_c3d_rs
174 MODULE PROCEDURE pw_axpy_r1d_r1d_gs
175 MODULE PROCEDURE pw_axpy_r1d_c1d_gs
176 MODULE PROCEDURE pw_axpy_r3d_r3d_gs
177 MODULE PROCEDURE pw_axpy_r3d_c3d_gs
178 MODULE PROCEDURE pw_axpy_c1d_r1d_gs
179 MODULE PROCEDURE pw_axpy_c1d_c1d_gs
180 MODULE PROCEDURE pw_axpy_c3d_r3d_gs
181 MODULE PROCEDURE pw_axpy_c3d_c3d_gs
185 MODULE PROCEDURE pw_multiply_r1d_r1d_rs
186 MODULE PROCEDURE pw_multiply_r1d_c1d_rs
187 MODULE PROCEDURE pw_multiply_r3d_r3d_rs
188 MODULE PROCEDURE pw_multiply_r3d_c3d_rs
189 MODULE PROCEDURE pw_multiply_c1d_r1d_rs
190 MODULE PROCEDURE pw_multiply_c1d_c1d_rs
191 MODULE PROCEDURE pw_multiply_c3d_r3d_rs
192 MODULE PROCEDURE pw_multiply_c3d_c3d_rs
193 MODULE PROCEDURE pw_multiply_r1d_r1d_gs
194 MODULE PROCEDURE pw_multiply_r1d_c1d_gs
195 MODULE PROCEDURE pw_multiply_r3d_r3d_gs
196 MODULE PROCEDURE pw_multiply_r3d_c3d_gs
197 MODULE PROCEDURE pw_multiply_c1d_r1d_gs
198 MODULE PROCEDURE pw_multiply_c1d_c1d_gs
199 MODULE PROCEDURE pw_multiply_c3d_r3d_gs
200 MODULE PROCEDURE pw_multiply_c3d_c3d_gs
204 MODULE PROCEDURE pw_multiply_with_r1d_r1d_rs
205 MODULE PROCEDURE pw_multiply_with_r1d_c1d_rs
206 MODULE PROCEDURE pw_multiply_with_r3d_r3d_rs
207 MODULE PROCEDURE pw_multiply_with_r3d_c3d_rs
208 MODULE PROCEDURE pw_multiply_with_c1d_r1d_rs
209 MODULE PROCEDURE pw_multiply_with_c1d_c1d_rs
210 MODULE PROCEDURE pw_multiply_with_c3d_r3d_rs
211 MODULE PROCEDURE pw_multiply_with_c3d_c3d_rs
212 MODULE PROCEDURE pw_multiply_with_r1d_r1d_gs
213 MODULE PROCEDURE pw_multiply_with_r1d_c1d_gs
214 MODULE PROCEDURE pw_multiply_with_r3d_r3d_gs
215 MODULE PROCEDURE pw_multiply_with_r3d_c3d_gs
216 MODULE PROCEDURE pw_multiply_with_c1d_r1d_gs
217 MODULE PROCEDURE pw_multiply_with_c1d_c1d_gs
218 MODULE PROCEDURE pw_multiply_with_c3d_r3d_gs
219 MODULE PROCEDURE pw_multiply_with_c3d_c3d_gs
223 MODULE PROCEDURE pw_integral_ab_r1d_r1d_rs
224 MODULE PROCEDURE pw_integral_ab_r1d_c1d_rs
225 MODULE PROCEDURE pw_integral_ab_r3d_r3d_rs
226 MODULE PROCEDURE pw_integral_ab_r3d_c3d_rs
227 MODULE PROCEDURE pw_integral_ab_c1d_r1d_rs
228 MODULE PROCEDURE pw_integral_ab_c1d_c1d_rs
229 MODULE PROCEDURE pw_integral_ab_c3d_r3d_rs
230 MODULE PROCEDURE pw_integral_ab_c3d_c3d_rs
231 MODULE PROCEDURE pw_integral_ab_r1d_r1d_gs
232 MODULE PROCEDURE pw_integral_ab_r1d_c1d_gs
233 MODULE PROCEDURE pw_integral_ab_r3d_r3d_gs
234 MODULE PROCEDURE pw_integral_ab_r3d_c3d_gs
235 MODULE PROCEDURE pw_integral_ab_c1d_r1d_gs
236 MODULE PROCEDURE pw_integral_ab_c1d_c1d_gs
237 MODULE PROCEDURE pw_integral_ab_c3d_r3d_gs
238 MODULE PROCEDURE pw_integral_ab_c3d_c3d_gs
242 MODULE PROCEDURE pw_integral_a2b_r1d_r1d
243 MODULE PROCEDURE pw_integral_a2b_r1d_c1d
244 MODULE PROCEDURE pw_integral_a2b_c1d_r1d
245 MODULE PROCEDURE pw_integral_a2b_c1d_c1d
249 MODULE PROCEDURE pw_gather_p_r1d
250 MODULE PROCEDURE pw_gather_p_c1d
251 MODULE PROCEDURE pw_gather_s_r1d_r3d
252 MODULE PROCEDURE pw_gather_s_r1d_c3d
253 MODULE PROCEDURE pw_gather_s_c1d_r3d
254 MODULE PROCEDURE pw_gather_s_c1d_c3d
258 MODULE PROCEDURE pw_scatter_p_r1d
259 MODULE PROCEDURE pw_scatter_p_c1d
260 MODULE PROCEDURE pw_scatter_s_r1d_r3d
261 MODULE PROCEDURE pw_scatter_s_r1d_c3d
262 MODULE PROCEDURE pw_scatter_s_c1d_r3d
263 MODULE PROCEDURE pw_scatter_s_c1d_c3d
267 MODULE PROCEDURE pw_copy_to_array_r1d_r1d_rs
268 MODULE PROCEDURE pw_copy_to_array_r1d_c1d_rs
269 MODULE PROCEDURE pw_copy_to_array_r3d_r3d_rs
270 MODULE PROCEDURE pw_copy_to_array_r3d_c3d_rs
271 MODULE PROCEDURE pw_copy_to_array_c1d_r1d_rs
272 MODULE PROCEDURE pw_copy_to_array_c1d_c1d_rs
273 MODULE PROCEDURE pw_copy_to_array_c3d_r3d_rs
274 MODULE PROCEDURE pw_copy_to_array_c3d_c3d_rs
275 MODULE PROCEDURE pw_copy_to_array_r1d_r1d_gs
276 MODULE PROCEDURE pw_copy_to_array_r1d_c1d_gs
277 MODULE PROCEDURE pw_copy_to_array_r3d_r3d_gs
278 MODULE PROCEDURE pw_copy_to_array_r3d_c3d_gs
279 MODULE PROCEDURE pw_copy_to_array_c1d_r1d_gs
280 MODULE PROCEDURE pw_copy_to_array_c1d_c1d_gs
281 MODULE PROCEDURE pw_copy_to_array_c3d_r3d_gs
282 MODULE PROCEDURE pw_copy_to_array_c3d_c3d_gs
286 MODULE PROCEDURE pw_copy_from_array_r1d_r1d_rs
287 MODULE PROCEDURE pw_copy_from_array_r1d_c1d_rs
288 MODULE PROCEDURE pw_copy_from_array_r3d_r3d_rs
289 MODULE PROCEDURE pw_copy_from_array_r3d_c3d_rs
290 MODULE PROCEDURE pw_copy_from_array_c1d_r1d_rs
291 MODULE PROCEDURE pw_copy_from_array_c1d_c1d_rs
292 MODULE PROCEDURE pw_copy_from_array_c3d_r3d_rs
293 MODULE PROCEDURE pw_copy_from_array_c3d_c3d_rs
294 MODULE PROCEDURE pw_copy_from_array_r1d_r1d_gs
295 MODULE PROCEDURE pw_copy_from_array_r1d_c1d_gs
296 MODULE PROCEDURE pw_copy_from_array_r3d_r3d_gs
297 MODULE PROCEDURE pw_copy_from_array_r3d_c3d_gs
298 MODULE PROCEDURE pw_copy_from_array_c1d_r1d_gs
299 MODULE PROCEDURE pw_copy_from_array_c1d_c1d_gs
300 MODULE PROCEDURE pw_copy_from_array_c3d_r3d_gs
301 MODULE PROCEDURE pw_copy_from_array_c3d_c3d_gs
305 MODULE PROCEDURE pw_copy_r1d_r1d_rs
306 MODULE PROCEDURE pw_copy_r1d_r1d_gs
307 MODULE PROCEDURE pw_gather_s_r1d_r3d_2
308 MODULE PROCEDURE pw_scatter_s_r1d_r3d_2
309 MODULE PROCEDURE pw_copy_r1d_c1d_rs
310 MODULE PROCEDURE pw_copy_r1d_c1d_gs
311 MODULE PROCEDURE pw_gather_s_r1d_c3d_2
312 MODULE PROCEDURE pw_scatter_s_r1d_c3d_2
313 MODULE PROCEDURE pw_copy_r3d_r3d_rs
314 MODULE PROCEDURE pw_copy_r3d_r3d_gs
315 MODULE PROCEDURE fft_wrap_pw1pw2_r3d_c1d_rs_gs
316 MODULE PROCEDURE pw_copy_r3d_c3d_rs
317 MODULE PROCEDURE pw_copy_r3d_c3d_gs
318 MODULE PROCEDURE fft_wrap_pw1pw2_r3d_c3d_rs_gs
319 MODULE PROCEDURE pw_copy_c1d_r1d_rs
320 MODULE PROCEDURE pw_copy_c1d_r1d_gs
321 MODULE PROCEDURE pw_gather_s_c1d_r3d_2
322 MODULE PROCEDURE pw_scatter_s_c1d_r3d_2
323 MODULE PROCEDURE fft_wrap_pw1pw2_c1d_r3d_gs_rs
324 MODULE PROCEDURE pw_copy_c1d_c1d_rs
325 MODULE PROCEDURE pw_copy_c1d_c1d_gs
326 MODULE PROCEDURE pw_gather_s_c1d_c3d_2
327 MODULE PROCEDURE pw_scatter_s_c1d_c3d_2
328 MODULE PROCEDURE fft_wrap_pw1pw2_c1d_c3d_gs_rs
329 MODULE PROCEDURE pw_copy_c3d_r3d_rs
330 MODULE PROCEDURE pw_copy_c3d_r3d_gs
331 MODULE PROCEDURE fft_wrap_pw1pw2_c3d_r3d_gs_rs
332 MODULE PROCEDURE fft_wrap_pw1pw2_c3d_c1d_rs_gs
333 MODULE PROCEDURE pw_copy_c3d_c3d_rs
334 MODULE PROCEDURE pw_copy_c3d_c3d_gs
335 MODULE PROCEDURE fft_wrap_pw1pw2_c3d_c3d_rs_gs
336 MODULE PROCEDURE fft_wrap_pw1pw2_c3d_c3d_gs_rs
347 SUBROUTINE pw_zero_r1d_rs (pw)
351 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_zero'
355 CALL timeset(routinen, handle)
357#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
365 CALL timestop(handle)
367 END SUBROUTINE pw_zero_r1d_rs
376 SUBROUTINE pw_scale_r1d_rs (pw, a)
379 REAL(KIND=
dp),
INTENT(IN) :: a
381 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_scale'
385 CALL timeset(routinen, handle)
388 pw%array = a*pw%array
391 CALL timestop(handle)
393 END SUBROUTINE pw_scale_r1d_rs
405 SUBROUTINE pw_write_r1d_rs (pw, unit_nr)
408 INTEGER,
INTENT(in) :: unit_nr
412 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
414 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=r1d"
415 IF (
ASSOCIATED(pw%array))
THEN
416 WRITE (unit=unit_nr, fmt=
"(' array=<r(',i8,':',i8,')>')") &
417 lbound(pw%array, 1), ubound(pw%array, 1)
419 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
422 END SUBROUTINE pw_write_r1d_rs
431 FUNCTION pw_integrate_function_r1d_rs (fun, isign, oprt)
RESULT(total_fun)
434 INTEGER,
INTENT(IN),
OPTIONAL :: isign
435 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
436 REAL(kind=
dp) :: total_fun
442 IF (
PRESENT(oprt))
THEN
447 cpabort(
"Unknown operator")
455 total_fun = fun%pw_grid%dvol*
accurate_sum(abs(fun%array))
461 CALL fun%pw_grid%para%group%sum(total_fun)
464 IF (
PRESENT(isign))
THEN
465 total_fun = total_fun*sign(1._dp, real(isign,
dp))
468 END FUNCTION pw_integrate_function_r1d_rs
475 SUBROUTINE pw_set_value_r1d_rs (pw, value)
477 REAL(KIND=
dp),
INTENT(IN) ::
value
479 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_set_value'
483 CALL timeset(routinen, handle)
489 CALL timestop(handle)
491 END SUBROUTINE pw_set_value_r1d_rs
499 SUBROUTINE pw_zero_r1d_gs (pw)
503 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_zero'
507 CALL timeset(routinen, handle)
509#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
517 CALL timestop(handle)
519 END SUBROUTINE pw_zero_r1d_gs
528 SUBROUTINE pw_scale_r1d_gs (pw, a)
531 REAL(KIND=
dp),
INTENT(IN) :: a
533 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_scale'
537 CALL timeset(routinen, handle)
540 pw%array = a*pw%array
543 CALL timestop(handle)
545 END SUBROUTINE pw_scale_r1d_gs
557 SUBROUTINE pw_write_r1d_gs (pw, unit_nr)
560 INTEGER,
INTENT(in) :: unit_nr
564 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
566 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=r1d"
567 IF (
ASSOCIATED(pw%array))
THEN
568 WRITE (unit=unit_nr, fmt=
"(' array=<r(',i8,':',i8,')>')") &
569 lbound(pw%array, 1), ubound(pw%array, 1)
571 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
574 END SUBROUTINE pw_write_r1d_gs
583 FUNCTION pw_integrate_function_r1d_gs (fun, isign, oprt)
RESULT(total_fun)
586 INTEGER,
INTENT(IN),
OPTIONAL :: isign
587 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
588 REAL(kind=
dp) :: total_fun
594 IF (
PRESENT(oprt))
THEN
599 cpabort(
"Unknown operator")
606 cpabort(
"Operator ABS not implemented")
607 IF (fun%pw_grid%have_g0) total_fun = fun%pw_grid%vol* fun%array(1)
610 CALL fun%pw_grid%para%group%sum(total_fun)
613 IF (
PRESENT(isign))
THEN
614 total_fun = total_fun*sign(1._dp, real(isign,
dp))
617 END FUNCTION pw_integrate_function_r1d_gs
624 SUBROUTINE pw_set_value_r1d_gs (pw, value)
626 REAL(KIND=
dp),
INTENT(IN) ::
value
628 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_set_value'
632 CALL timeset(routinen, handle)
638 CALL timestop(handle)
640 END SUBROUTINE pw_set_value_r1d_gs
648 SUBROUTINE pw_gather_p_r1d (pw, c, scale)
651 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
INTENT(IN) :: c
652 REAL(KIND=
dp),
INTENT(IN),
OPTIONAL :: scale
654 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_gather_p'
656 INTEGER :: gpt, handle, l, m, mn, n
658 CALL timeset(routinen, handle)
661 cpabort(
"This grid type is not distributed")
664 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
665 ngpts =>
SIZE(pw%pw_grid%gsq), ghat => pw%pw_grid%g_hat, yzq => pw%pw_grid%para%yzq)
667 IF (
PRESENT(scale))
THEN
672 l = mapl(ghat(1, gpt)) + 1
673 m = mapm(ghat(2, gpt)) + 1
674 n = mapn(ghat(3, gpt)) + 1
676 pw%array(gpt) = scale* real(c(l, mn), kind=
dp)
684 l = mapl(ghat(1, gpt)) + 1
685 m = mapm(ghat(2, gpt)) + 1
686 n = mapn(ghat(3, gpt)) + 1
688 pw%array(gpt) = real(c(l, mn), kind=
dp)
695 CALL timestop(handle)
697 END SUBROUTINE pw_gather_p_r1d
705 SUBROUTINE pw_scatter_p_r1d (pw, c, scale)
706 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw
707 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
INTENT(INOUT) :: c
708 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
710 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_scatter_p'
712 INTEGER :: gpt, handle, l, m, mn, n
714 CALL timeset(routinen, handle)
716 IF (pw%pw_grid%para%mode /= pw_mode_distributed)
THEN
717 cpabort(
"This grid type is not distributed")
720 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
721 ghat => pw%pw_grid%g_hat, yzq => pw%pw_grid%para%yzq, ngpts =>
SIZE(pw%pw_grid%gsq))
723 IF (.NOT.
PRESENT(scale)) c = z_zero
725 IF (
PRESENT(scale))
THEN
730 l = mapl(ghat(1, gpt)) + 1
731 m = mapm(ghat(2, gpt)) + 1
732 n = mapn(ghat(3, gpt)) + 1
734 c(l, mn) = cmplx(scale*pw%array(gpt), 0.0_dp, kind=dp)
742 l = mapl(ghat(1, gpt)) + 1
743 m = mapm(ghat(2, gpt)) + 1
744 n = mapn(ghat(3, gpt)) + 1
746 c(l, mn) = cmplx(pw%array(gpt), 0.0_dp, kind=dp)
753 IF (pw%pw_grid%grid_span == halfspace)
THEN
755 associate(mapm => pw%pw_grid%mapm%neg, mapn => pw%pw_grid%mapn%neg, mapl => pw%pw_grid%mapl%neg, &
756 ghat => pw%pw_grid%g_hat, ngpts =>
SIZE(pw%pw_grid%gsq), yzq => pw%pw_grid%para%yzq)
758 IF (
PRESENT(scale))
THEN
763 l = mapl(ghat(1, gpt)) + 1
764 m = mapm(ghat(2, gpt)) + 1
765 n = mapn(ghat(3, gpt)) + 1
767 c(l, mn) = scale*( cmplx(pw%array(gpt), 0.0_dp, kind=dp))
775 l = mapl(ghat(1, gpt)) + 1
776 m = mapm(ghat(2, gpt)) + 1
777 n = mapn(ghat(3, gpt)) + 1
779 c(l, mn) = ( cmplx(pw%array(gpt), 0.0_dp, kind=dp))
786 CALL timestop(handle)
788 END SUBROUTINE pw_scatter_p_r1d
796 SUBROUTINE pw_zero_r3d_rs (pw)
798 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw
800 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_zero'
804 CALL timeset(routinen, handle)
806#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
814 CALL timestop(handle)
816 END SUBROUTINE pw_zero_r3d_rs
825 SUBROUTINE pw_scale_r3d_rs (pw, a)
827 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw
828 REAL(KIND=dp),
INTENT(IN) :: a
830 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_scale'
834 CALL timeset(routinen, handle)
837 pw%array = a*pw%array
840 CALL timestop(handle)
842 END SUBROUTINE pw_scale_r3d_rs
854 SUBROUTINE pw_write_r3d_rs (pw, unit_nr)
856 TYPE(pw_r3d_rs_type),
INTENT(in) :: pw
857 INTEGER,
INTENT(in) :: unit_nr
861 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
863 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=r3d"
864 IF (
ASSOCIATED(pw%array))
THEN
865 WRITE (unit=unit_nr, fmt=
"(' array=<r(',i8,':',i8,',',i8,':',i8,',',i8,':',i8,')>')") &
866 lbound(pw%array, 1), ubound(pw%array, 1), lbound(pw%array, 2), ubound(pw%array, 2), &
867 lbound(pw%array, 3), ubound(pw%array, 3)
869 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
872 END SUBROUTINE pw_write_r3d_rs
881 FUNCTION pw_integrate_function_r3d_rs (fun, isign, oprt)
RESULT(total_fun)
883 TYPE(pw_r3d_rs_type),
INTENT(IN) :: fun
884 INTEGER,
INTENT(IN),
OPTIONAL :: isign
885 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
886 REAL(kind=dp) :: total_fun
892 IF (
PRESENT(oprt))
THEN
897 cpabort(
"Unknown operator")
905 total_fun = fun%pw_grid%dvol*accurate_sum(abs(fun%array))
907 total_fun = fun%pw_grid%dvol*accurate_sum( fun%array)
910 IF (fun%pw_grid%para%mode /= pw_mode_local)
THEN
911 CALL fun%pw_grid%para%group%sum(total_fun)
914 IF (
PRESENT(isign))
THEN
915 total_fun = total_fun*sign(1._dp, real(isign, dp))
918 END FUNCTION pw_integrate_function_r3d_rs
925 SUBROUTINE pw_set_value_r3d_rs (pw, value)
926 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw
927 REAL(KIND=dp),
INTENT(IN) ::
value
929 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_set_value'
933 CALL timeset(routinen, handle)
939 CALL timestop(handle)
941 END SUBROUTINE pw_set_value_r3d_rs
949 SUBROUTINE pw_zero_r3d_gs (pw)
951 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw
953 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_zero'
957 CALL timeset(routinen, handle)
959#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
967 CALL timestop(handle)
969 END SUBROUTINE pw_zero_r3d_gs
978 SUBROUTINE pw_scale_r3d_gs (pw, a)
980 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw
981 REAL(KIND=dp),
INTENT(IN) :: a
983 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_scale'
987 CALL timeset(routinen, handle)
990 pw%array = a*pw%array
993 CALL timestop(handle)
995 END SUBROUTINE pw_scale_r3d_gs
1007 SUBROUTINE pw_write_r3d_gs (pw, unit_nr)
1009 TYPE(pw_r3d_gs_type),
INTENT(in) :: pw
1010 INTEGER,
INTENT(in) :: unit_nr
1014 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
1016 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=r3d"
1017 IF (
ASSOCIATED(pw%array))
THEN
1018 WRITE (unit=unit_nr, fmt=
"(' array=<r(',i8,':',i8,',',i8,':',i8,',',i8,':',i8,')>')") &
1019 lbound(pw%array, 1), ubound(pw%array, 1), lbound(pw%array, 2), ubound(pw%array, 2), &
1020 lbound(pw%array, 3), ubound(pw%array, 3)
1022 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
1025 END SUBROUTINE pw_write_r3d_gs
1034 FUNCTION pw_integrate_function_r3d_gs (fun, isign, oprt)
RESULT(total_fun)
1036 TYPE(pw_r3d_gs_type),
INTENT(IN) :: fun
1037 INTEGER,
INTENT(IN),
OPTIONAL :: isign
1038 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
1039 REAL(kind=dp) :: total_fun
1045 IF (
PRESENT(oprt))
THEN
1050 cpabort(
"Unknown operator")
1057 cpabort(
"Operator ABS not implemented")
1058 cpabort(
"Reciprocal space integration for 3D grids not implemented!")
1060 IF (fun%pw_grid%para%mode /= pw_mode_local)
THEN
1061 CALL fun%pw_grid%para%group%sum(total_fun)
1064 IF (
PRESENT(isign))
THEN
1065 total_fun = total_fun*sign(1._dp, real(isign, dp))
1068 END FUNCTION pw_integrate_function_r3d_gs
1075 SUBROUTINE pw_set_value_r3d_gs (pw, value)
1076 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw
1077 REAL(KIND=dp),
INTENT(IN) ::
value
1079 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_set_value'
1083 CALL timeset(routinen, handle)
1089 CALL timestop(handle)
1091 END SUBROUTINE pw_set_value_r3d_gs
1100 SUBROUTINE pw_zero_c1d_rs (pw)
1102 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw
1104 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_zero'
1108 CALL timeset(routinen, handle)
1110#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
1112 pw%array = cmplx(0.0_dp, 0.0_dp, kind=dp)
1115 pw%array = cmplx(0.0_dp, 0.0_dp, kind=dp)
1118 CALL timestop(handle)
1120 END SUBROUTINE pw_zero_c1d_rs
1129 SUBROUTINE pw_scale_c1d_rs (pw, a)
1131 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw
1132 REAL(KIND=dp),
INTENT(IN) :: a
1134 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_scale'
1138 CALL timeset(routinen, handle)
1141 pw%array = a*pw%array
1144 CALL timestop(handle)
1146 END SUBROUTINE pw_scale_c1d_rs
1158 SUBROUTINE pw_write_c1d_rs (pw, unit_nr)
1160 TYPE(pw_c1d_rs_type),
INTENT(in) :: pw
1161 INTEGER,
INTENT(in) :: unit_nr
1165 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
1167 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=c1d"
1168 IF (
ASSOCIATED(pw%array))
THEN
1169 WRITE (unit=unit_nr, fmt=
"(' array=<c(',i8,':',i8,')>')") &
1170 lbound(pw%array, 1), ubound(pw%array, 1)
1172 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
1175 END SUBROUTINE pw_write_c1d_rs
1184 FUNCTION pw_integrate_function_c1d_rs (fun, isign, oprt)
RESULT(total_fun)
1186 TYPE(pw_c1d_rs_type),
INTENT(IN) :: fun
1187 INTEGER,
INTENT(IN),
OPTIONAL :: isign
1188 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
1189 REAL(kind=dp) :: total_fun
1195 IF (
PRESENT(oprt))
THEN
1200 cpabort(
"Unknown operator")
1208 total_fun = fun%pw_grid%dvol*accurate_sum(abs(fun%array))
1210 total_fun = fun%pw_grid%dvol*accurate_sum( real(fun%array, kind=dp))
1213 IF (fun%pw_grid%para%mode /= pw_mode_local)
THEN
1214 CALL fun%pw_grid%para%group%sum(total_fun)
1217 IF (
PRESENT(isign))
THEN
1218 total_fun = total_fun*sign(1._dp, real(isign, dp))
1221 END FUNCTION pw_integrate_function_c1d_rs
1228 SUBROUTINE pw_set_value_c1d_rs (pw, value)
1229 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw
1230 REAL(KIND=dp),
INTENT(IN) ::
value
1232 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_set_value'
1236 CALL timeset(routinen, handle)
1239 pw%array = cmplx(
value, 0.0_dp, kind=dp)
1242 CALL timestop(handle)
1244 END SUBROUTINE pw_set_value_c1d_rs
1252 SUBROUTINE pw_zero_c1d_gs (pw)
1254 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
1256 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_zero'
1260 CALL timeset(routinen, handle)
1262#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
1264 pw%array = cmplx(0.0_dp, 0.0_dp, kind=dp)
1267 pw%array = cmplx(0.0_dp, 0.0_dp, kind=dp)
1270 CALL timestop(handle)
1272 END SUBROUTINE pw_zero_c1d_gs
1281 SUBROUTINE pw_scale_c1d_gs (pw, a)
1283 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
1284 REAL(KIND=dp),
INTENT(IN) :: a
1286 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_scale'
1290 CALL timeset(routinen, handle)
1293 pw%array = a*pw%array
1296 CALL timestop(handle)
1298 END SUBROUTINE pw_scale_c1d_gs
1310 SUBROUTINE pw_write_c1d_gs (pw, unit_nr)
1312 TYPE(pw_c1d_gs_type),
INTENT(in) :: pw
1313 INTEGER,
INTENT(in) :: unit_nr
1317 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
1319 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=c1d"
1320 IF (
ASSOCIATED(pw%array))
THEN
1321 WRITE (unit=unit_nr, fmt=
"(' array=<c(',i8,':',i8,')>')") &
1322 lbound(pw%array, 1), ubound(pw%array, 1)
1324 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
1327 END SUBROUTINE pw_write_c1d_gs
1336 FUNCTION pw_integrate_function_c1d_gs (fun, isign, oprt)
RESULT(total_fun)
1338 TYPE(pw_c1d_gs_type),
INTENT(IN) :: fun
1339 INTEGER,
INTENT(IN),
OPTIONAL :: isign
1340 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
1341 REAL(kind=dp) :: total_fun
1347 IF (
PRESENT(oprt))
THEN
1352 cpabort(
"Unknown operator")
1359 cpabort(
"Operator ABS not implemented")
1360 IF (fun%pw_grid%have_g0) total_fun = fun%pw_grid%vol* real(fun%array(1), kind=dp)
1362 IF (fun%pw_grid%para%mode /= pw_mode_local)
THEN
1363 CALL fun%pw_grid%para%group%sum(total_fun)
1366 IF (
PRESENT(isign))
THEN
1367 total_fun = total_fun*sign(1._dp, real(isign, dp))
1370 END FUNCTION pw_integrate_function_c1d_gs
1377 SUBROUTINE pw_set_value_c1d_gs (pw, value)
1378 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
1379 REAL(KIND=dp),
INTENT(IN) ::
value
1381 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_set_value'
1385 CALL timeset(routinen, handle)
1388 pw%array = cmplx(
value, 0.0_dp, kind=dp)
1391 CALL timestop(handle)
1393 END SUBROUTINE pw_set_value_c1d_gs
1401 SUBROUTINE pw_gather_p_c1d (pw, c, scale)
1403 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
1404 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
INTENT(IN) :: c
1405 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
1407 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_gather_p'
1409 INTEGER :: gpt, handle, l, m, mn, n
1411 CALL timeset(routinen, handle)
1413 IF (pw%pw_grid%para%mode /= pw_mode_distributed)
THEN
1414 cpabort(
"This grid type is not distributed")
1417 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
1418 ngpts =>
SIZE(pw%pw_grid%gsq), ghat => pw%pw_grid%g_hat, yzq => pw%pw_grid%para%yzq)
1420 IF (
PRESENT(scale))
THEN
1425 l = mapl(ghat(1, gpt)) + 1
1426 m = mapm(ghat(2, gpt)) + 1
1427 n = mapn(ghat(3, gpt)) + 1
1429 pw%array(gpt) = scale* c(l, mn)
1437 l = mapl(ghat(1, gpt)) + 1
1438 m = mapm(ghat(2, gpt)) + 1
1439 n = mapn(ghat(3, gpt)) + 1
1441 pw%array(gpt) = c(l, mn)
1448 CALL timestop(handle)
1450 END SUBROUTINE pw_gather_p_c1d
1458 SUBROUTINE pw_scatter_p_c1d (pw, c, scale)
1459 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
1460 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
INTENT(INOUT) :: c
1461 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
1463 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_scatter_p'
1465 INTEGER :: gpt, handle, l, m, mn, n
1467 CALL timeset(routinen, handle)
1469 IF (pw%pw_grid%para%mode /= pw_mode_distributed)
THEN
1470 cpabort(
"This grid type is not distributed")
1473 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
1474 ghat => pw%pw_grid%g_hat, yzq => pw%pw_grid%para%yzq, ngpts =>
SIZE(pw%pw_grid%gsq))
1476 IF (.NOT.
PRESENT(scale)) c = z_zero
1478 IF (
PRESENT(scale))
THEN
1483 l = mapl(ghat(1, gpt)) + 1
1484 m = mapm(ghat(2, gpt)) + 1
1485 n = mapn(ghat(3, gpt)) + 1
1487 c(l, mn) = scale*pw%array(gpt)
1495 l = mapl(ghat(1, gpt)) + 1
1496 m = mapm(ghat(2, gpt)) + 1
1497 n = mapn(ghat(3, gpt)) + 1
1499 c(l, mn) = pw%array(gpt)
1506 IF (pw%pw_grid%grid_span == halfspace)
THEN
1508 associate(mapm => pw%pw_grid%mapm%neg, mapn => pw%pw_grid%mapn%neg, mapl => pw%pw_grid%mapl%neg, &
1509 ghat => pw%pw_grid%g_hat, ngpts =>
SIZE(pw%pw_grid%gsq), yzq => pw%pw_grid%para%yzq)
1511 IF (
PRESENT(scale))
THEN
1516 l = mapl(ghat(1, gpt)) + 1
1517 m = mapm(ghat(2, gpt)) + 1
1518 n = mapn(ghat(3, gpt)) + 1
1520 c(l, mn) = scale*conjg( pw%array(gpt))
1528 l = mapl(ghat(1, gpt)) + 1
1529 m = mapm(ghat(2, gpt)) + 1
1530 n = mapn(ghat(3, gpt)) + 1
1532 c(l, mn) = conjg( pw%array(gpt))
1539 CALL timestop(handle)
1541 END SUBROUTINE pw_scatter_p_c1d
1549 SUBROUTINE pw_zero_c3d_rs (pw)
1551 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw
1553 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_zero'
1557 CALL timeset(routinen, handle)
1559#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
1561 pw%array = cmplx(0.0_dp, 0.0_dp, kind=dp)
1564 pw%array = cmplx(0.0_dp, 0.0_dp, kind=dp)
1567 CALL timestop(handle)
1569 END SUBROUTINE pw_zero_c3d_rs
1578 SUBROUTINE pw_scale_c3d_rs (pw, a)
1580 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw
1581 REAL(KIND=dp),
INTENT(IN) :: a
1583 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_scale'
1587 CALL timeset(routinen, handle)
1590 pw%array = a*pw%array
1593 CALL timestop(handle)
1595 END SUBROUTINE pw_scale_c3d_rs
1607 SUBROUTINE pw_write_c3d_rs (pw, unit_nr)
1609 TYPE(pw_c3d_rs_type),
INTENT(in) :: pw
1610 INTEGER,
INTENT(in) :: unit_nr
1614 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
1616 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=c3d"
1617 IF (
ASSOCIATED(pw%array))
THEN
1618 WRITE (unit=unit_nr, fmt=
"(' array=<c(',i8,':',i8,',',i8,':',i8,',',i8,':',i8,')>')") &
1619 lbound(pw%array, 1), ubound(pw%array, 1), lbound(pw%array, 2), ubound(pw%array, 2), &
1620 lbound(pw%array, 3), ubound(pw%array, 3)
1622 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
1625 END SUBROUTINE pw_write_c3d_rs
1634 FUNCTION pw_integrate_function_c3d_rs (fun, isign, oprt)
RESULT(total_fun)
1636 TYPE(pw_c3d_rs_type),
INTENT(IN) :: fun
1637 INTEGER,
INTENT(IN),
OPTIONAL :: isign
1638 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
1639 REAL(kind=dp) :: total_fun
1645 IF (
PRESENT(oprt))
THEN
1650 cpabort(
"Unknown operator")
1658 total_fun = fun%pw_grid%dvol*accurate_sum(abs(fun%array))
1660 total_fun = fun%pw_grid%dvol*accurate_sum( real(fun%array, kind=dp))
1663 IF (fun%pw_grid%para%mode /= pw_mode_local)
THEN
1664 CALL fun%pw_grid%para%group%sum(total_fun)
1667 IF (
PRESENT(isign))
THEN
1668 total_fun = total_fun*sign(1._dp, real(isign, dp))
1671 END FUNCTION pw_integrate_function_c3d_rs
1678 SUBROUTINE pw_set_value_c3d_rs (pw, value)
1679 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw
1680 REAL(KIND=dp),
INTENT(IN) ::
value
1682 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_set_value'
1686 CALL timeset(routinen, handle)
1689 pw%array = cmplx(
value, 0.0_dp, kind=dp)
1692 CALL timestop(handle)
1694 END SUBROUTINE pw_set_value_c3d_rs
1702 SUBROUTINE pw_zero_c3d_gs (pw)
1704 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw
1706 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_zero'
1710 CALL timeset(routinen, handle)
1712#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
1714 pw%array = cmplx(0.0_dp, 0.0_dp, kind=dp)
1717 pw%array = cmplx(0.0_dp, 0.0_dp, kind=dp)
1720 CALL timestop(handle)
1722 END SUBROUTINE pw_zero_c3d_gs
1731 SUBROUTINE pw_scale_c3d_gs (pw, a)
1733 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw
1734 REAL(KIND=dp),
INTENT(IN) :: a
1736 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_scale'
1740 CALL timeset(routinen, handle)
1743 pw%array = a*pw%array
1746 CALL timestop(handle)
1748 END SUBROUTINE pw_scale_c3d_gs
1760 SUBROUTINE pw_write_c3d_gs (pw, unit_nr)
1762 TYPE(pw_c3d_gs_type),
INTENT(in) :: pw
1763 INTEGER,
INTENT(in) :: unit_nr
1767 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
1769 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=c3d"
1770 IF (
ASSOCIATED(pw%array))
THEN
1771 WRITE (unit=unit_nr, fmt=
"(' array=<c(',i8,':',i8,',',i8,':',i8,',',i8,':',i8,')>')") &
1772 lbound(pw%array, 1), ubound(pw%array, 1), lbound(pw%array, 2), ubound(pw%array, 2), &
1773 lbound(pw%array, 3), ubound(pw%array, 3)
1775 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
1778 END SUBROUTINE pw_write_c3d_gs
1787 FUNCTION pw_integrate_function_c3d_gs (fun, isign, oprt)
RESULT(total_fun)
1789 TYPE(pw_c3d_gs_type),
INTENT(IN) :: fun
1790 INTEGER,
INTENT(IN),
OPTIONAL :: isign
1791 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
1792 REAL(kind=dp) :: total_fun
1798 IF (
PRESENT(oprt))
THEN
1803 cpabort(
"Unknown operator")
1810 cpabort(
"Operator ABS not implemented")
1811 cpabort(
"Reciprocal space integration for 3D grids not implemented!")
1813 IF (fun%pw_grid%para%mode /= pw_mode_local)
THEN
1814 CALL fun%pw_grid%para%group%sum(total_fun)
1817 IF (
PRESENT(isign))
THEN
1818 total_fun = total_fun*sign(1._dp, real(isign, dp))
1821 END FUNCTION pw_integrate_function_c3d_gs
1828 SUBROUTINE pw_set_value_c3d_gs (pw, value)
1829 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw
1830 REAL(KIND=dp),
INTENT(IN) ::
value
1832 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_set_value'
1836 CALL timeset(routinen, handle)
1839 pw%array = cmplx(
value, 0.0_dp, kind=dp)
1842 CALL timestop(handle)
1844 END SUBROUTINE pw_set_value_c3d_gs
1860 SUBROUTINE pw_copy_r1d_r1d_rs (pw1, pw2)
1862 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
1863 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw2
1865 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
1868 INTEGER :: i, j, ng, ng1, ng2, ns
1870 CALL timeset(routinen, handle)
1872 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
1873 cpabort(
"Both grids must be either spherical or non-spherical!")
1874 IF (pw1%pw_grid%spherical) &
1875 cpabort(
"Spherical grids only exist in reciprocal space!")
1877 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
1878 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
1879 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
1880 ng1 =
SIZE(pw1%array)
1881 ng2 =
SIZE(pw2%array)
1884 pw2%array(1:ng) = pw1%array(1:ng)
1888 pw2%array(ng + 1:ng2) = 0.0_dp
1892 cpabort(
"Copies between spherical grids require compatible grids!")
1895 ng1 =
SIZE(pw1%array)
1896 ng2 =
SIZE(pw2%array)
1897 ns = 2*max(ng1, ng2)
1899 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
1900 IF (ng1 >= ng2)
THEN
1905 j = pw2%pw_grid%gidx(i)
1906 pw2%array(i) = pw1%array(j)
1915 j = pw2%pw_grid%gidx(i)
1916 pw2%array(j) = pw1%array(i)
1920 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
1921 IF (ng1 >= ng2)
THEN
1924 j = pw1%pw_grid%gidx(i)
1925 pw2%array(i) = pw1%array(j)
1932 j = pw1%pw_grid%gidx(i)
1933 pw2%array(j) = pw1%array(i)
1938 cpabort(
"Copy not implemented!")
1945 pw2%array = pw1%array
1949 CALL timestop(handle)
1951 END SUBROUTINE pw_copy_r1d_r1d_rs
1958 SUBROUTINE pw_copy_to_array_r1d_r1d_rs (pw, array)
1959 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw
1960 REAL(KIND=dp),
DIMENSION(:),
INTENT(INOUT) :: array
1962 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
1966 CALL timeset(routinen, handle)
1969 array(:) = pw%array(:)
1972 CALL timestop(handle)
1973 END SUBROUTINE pw_copy_to_array_r1d_r1d_rs
1980 SUBROUTINE pw_copy_from_array_r1d_r1d_rs (pw, array)
1981 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw
1982 REAL(KIND=dp),
DIMENSION(:),
INTENT(IN) :: array
1984 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
1988 CALL timeset(routinen, handle)
1994 CALL timestop(handle)
1995 END SUBROUTINE pw_copy_from_array_r1d_r1d_rs
2013 SUBROUTINE pw_axpy_r1d_r1d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
2015 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
2016 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw2
2017 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
2018 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
2020 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
2023 LOGICAL :: my_allow_noncompatible_grids
2024 REAL(KIND=dp) :: my_alpha, my_beta
2025 INTEGER :: i, j, ng, ng1, ng2
2027 CALL timeset(routinen, handle)
2030 IF (
PRESENT(alpha)) my_alpha = alpha
2033 IF (
PRESENT(beta)) my_beta = beta
2035 my_allow_noncompatible_grids = .false.
2036 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
2038 IF (my_beta /= 1.0_dp)
THEN
2039 IF (my_beta == 0.0_dp)
THEN
2043 pw2%array = pw2%array*my_beta
2048 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
2050 IF (my_alpha == 1.0_dp)
THEN
2052 pw2%array = pw2%array + pw1%array
2054 ELSE IF (my_alpha /= 0.0_dp)
THEN
2056 pw2%array = pw2%array + my_alpha* pw1%array
2060 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
2062 ng1 =
SIZE(pw1%array)
2063 ng2 =
SIZE(pw2%array)
2066 IF (pw1%pw_grid%spherical)
THEN
2067 IF (my_alpha == 1.0_dp)
THEN
2070 pw2%array(i) = pw2%array(i) + pw1%array(i)
2073 ELSE IF (my_alpha /= 0.0_dp)
THEN
2076 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(i)
2080 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
2081 IF (ng1 >= ng2)
THEN
2082 IF (my_alpha == 1.0_dp)
THEN
2085 j = pw2%pw_grid%gidx(i)
2086 pw2%array(i) = pw2%array(i) + pw1%array(j)
2089 ELSE IF (my_alpha /= 0.0_dp)
THEN
2092 j = pw2%pw_grid%gidx(i)
2093 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
2098 IF (my_alpha == 1.0_dp)
THEN
2101 j = pw2%pw_grid%gidx(i)
2102 pw2%array(j) = pw2%array(j) + pw1%array(i)
2105 ELSE IF (my_alpha /= 0.0_dp)
THEN
2108 j = pw2%pw_grid%gidx(i)
2109 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
2114 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
2115 IF (ng1 >= ng2)
THEN
2116 IF (my_alpha == 1.0_dp)
THEN
2119 j = pw1%pw_grid%gidx(i)
2120 pw2%array(i) = pw2%array(i) + pw1%array(j)
2123 ELSE IF (my_alpha /= 0.0_dp)
THEN
2126 j = pw1%pw_grid%gidx(i)
2127 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
2132 IF (my_alpha == 1.0_dp)
THEN
2135 j = pw1%pw_grid%gidx(i)
2136 pw2%array(j) = pw2%array(j) + pw1%array(i)
2139 ELSE IF (my_alpha /= 0.0_dp)
THEN
2142 j = pw1%pw_grid%gidx(i)
2143 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
2149 cpabort(
"Grids not compatible")
2154 cpabort(
"Grids not compatible")
2158 CALL timestop(handle)
2160 END SUBROUTINE pw_axpy_r1d_r1d_rs
2171 SUBROUTINE pw_multiply_r1d_r1d_rs (pw_out, pw1, pw2, alpha)
2173 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw_out
2174 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
2175 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw2
2176 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
2178 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
2181 REAL(KIND=dp) :: my_alpha
2183 CALL timeset(routinen, handle)
2186 IF (
PRESENT(alpha)) my_alpha = alpha
2188 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
2189 cpabort(
"pw_multiply not implemented for non-identical grids!")
2191#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
2192 IF (my_alpha == 1.0_dp)
THEN
2194 pw_out%array = pw_out%array + pw1%array* pw2%array
2198 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
2202 IF (my_alpha == 1.0_dp)
THEN
2203 pw_out%array = pw_out%array + pw1%array* pw2%array
2205 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
2209 CALL timestop(handle)
2211 END SUBROUTINE pw_multiply_r1d_r1d_rs
2218 SUBROUTINE pw_multiply_with_r1d_r1d_rs (pw1, pw2)
2219 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw1
2220 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw2
2222 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
2226 CALL timeset(routinen, handle)
2228 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
2229 cpabort(
"Incompatible grids!")
2232 pw1%array = pw1%array* pw2%array
2235 CALL timestop(handle)
2237 END SUBROUTINE pw_multiply_with_r1d_r1d_rs
2252 FUNCTION pw_integral_ab_r1d_r1d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
2254 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
2255 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw2
2256 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
2257 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
2258 REAL(kind=dp) :: integral_value
2260 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_r1d_r1d_rs'
2262 INTEGER :: handle, loc_sumtype
2263 LOGICAL :: my_just_sum, my_local_only
2265 CALL timeset(routinen, handle)
2268 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
2270 my_local_only = .false.
2271 IF (
PRESENT(local_only)) my_local_only = local_only
2273 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
2274 cpabort(
"Grids incompatible")
2277 my_just_sum = .false.
2278 IF (
PRESENT(just_sum)) my_just_sum = just_sum
2285 integral_value = dot_product(pw1%array, pw2%array)
2290 integral_value = accurate_dot_product(pw1%array, pw2%array)
2294 IF (.NOT. my_just_sum)
THEN
2295 integral_value = integral_value*pw1%pw_grid%dvol
2298 IF (pw1%pw_grid%grid_span == halfspace)
THEN
2299 integral_value = 2.0_dp*integral_value
2300 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
2301 pw1%array(1)*pw2%array(1)
2304 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
2305 CALL pw1%pw_grid%para%group%sum(integral_value)
2307 CALL timestop(handle)
2309 END FUNCTION pw_integral_ab_r1d_r1d_rs
2323 SUBROUTINE pw_copy_r1d_r1d_gs (pw1, pw2)
2325 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
2326 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw2
2328 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
2331 INTEGER :: i, j, ng, ng1, ng2, ns
2333 CALL timeset(routinen, handle)
2335 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
2336 cpabort(
"Both grids must be either spherical or non-spherical!")
2338 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
2339 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
2340 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
2341 ng1 =
SIZE(pw1%array)
2342 ng2 =
SIZE(pw2%array)
2345 pw2%array(1:ng) = pw1%array(1:ng)
2349 pw2%array(ng + 1:ng2) = 0.0_dp
2353 cpabort(
"Copies between spherical grids require compatible grids!")
2356 ng1 =
SIZE(pw1%array)
2357 ng2 =
SIZE(pw2%array)
2358 ns = 2*max(ng1, ng2)
2360 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
2361 IF (ng1 >= ng2)
THEN
2366 j = pw2%pw_grid%gidx(i)
2367 pw2%array(i) = pw1%array(j)
2376 j = pw2%pw_grid%gidx(i)
2377 pw2%array(j) = pw1%array(i)
2381 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
2382 IF (ng1 >= ng2)
THEN
2385 j = pw1%pw_grid%gidx(i)
2386 pw2%array(i) = pw1%array(j)
2393 j = pw1%pw_grid%gidx(i)
2394 pw2%array(j) = pw1%array(i)
2399 cpabort(
"Copy not implemented!")
2406 pw2%array = pw1%array
2410 CALL timestop(handle)
2412 END SUBROUTINE pw_copy_r1d_r1d_gs
2419 SUBROUTINE pw_copy_to_array_r1d_r1d_gs (pw, array)
2420 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw
2421 REAL(KIND=dp),
DIMENSION(:),
INTENT(INOUT) :: array
2423 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
2427 CALL timeset(routinen, handle)
2430 array(:) = pw%array(:)
2433 CALL timestop(handle)
2434 END SUBROUTINE pw_copy_to_array_r1d_r1d_gs
2441 SUBROUTINE pw_copy_from_array_r1d_r1d_gs (pw, array)
2442 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw
2443 REAL(KIND=dp),
DIMENSION(:),
INTENT(IN) :: array
2445 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
2449 CALL timeset(routinen, handle)
2455 CALL timestop(handle)
2456 END SUBROUTINE pw_copy_from_array_r1d_r1d_gs
2474 SUBROUTINE pw_axpy_r1d_r1d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
2476 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
2477 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw2
2478 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
2479 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
2481 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
2484 LOGICAL :: my_allow_noncompatible_grids
2485 REAL(KIND=dp) :: my_alpha, my_beta
2486 INTEGER :: i, j, ng, ng1, ng2
2488 CALL timeset(routinen, handle)
2491 IF (
PRESENT(alpha)) my_alpha = alpha
2494 IF (
PRESENT(beta)) my_beta = beta
2496 my_allow_noncompatible_grids = .false.
2497 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
2499 IF (my_beta /= 1.0_dp)
THEN
2500 IF (my_beta == 0.0_dp)
THEN
2504 pw2%array = pw2%array*my_beta
2509 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
2511 IF (my_alpha == 1.0_dp)
THEN
2513 pw2%array = pw2%array + pw1%array
2515 ELSE IF (my_alpha /= 0.0_dp)
THEN
2517 pw2%array = pw2%array + my_alpha* pw1%array
2521 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
2523 ng1 =
SIZE(pw1%array)
2524 ng2 =
SIZE(pw2%array)
2527 IF (pw1%pw_grid%spherical)
THEN
2528 IF (my_alpha == 1.0_dp)
THEN
2531 pw2%array(i) = pw2%array(i) + pw1%array(i)
2534 ELSE IF (my_alpha /= 0.0_dp)
THEN
2537 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(i)
2541 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
2542 IF (ng1 >= ng2)
THEN
2543 IF (my_alpha == 1.0_dp)
THEN
2546 j = pw2%pw_grid%gidx(i)
2547 pw2%array(i) = pw2%array(i) + pw1%array(j)
2550 ELSE IF (my_alpha /= 0.0_dp)
THEN
2553 j = pw2%pw_grid%gidx(i)
2554 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
2559 IF (my_alpha == 1.0_dp)
THEN
2562 j = pw2%pw_grid%gidx(i)
2563 pw2%array(j) = pw2%array(j) + pw1%array(i)
2566 ELSE IF (my_alpha /= 0.0_dp)
THEN
2569 j = pw2%pw_grid%gidx(i)
2570 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
2575 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
2576 IF (ng1 >= ng2)
THEN
2577 IF (my_alpha == 1.0_dp)
THEN
2580 j = pw1%pw_grid%gidx(i)
2581 pw2%array(i) = pw2%array(i) + pw1%array(j)
2584 ELSE IF (my_alpha /= 0.0_dp)
THEN
2587 j = pw1%pw_grid%gidx(i)
2588 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
2593 IF (my_alpha == 1.0_dp)
THEN
2596 j = pw1%pw_grid%gidx(i)
2597 pw2%array(j) = pw2%array(j) + pw1%array(i)
2600 ELSE IF (my_alpha /= 0.0_dp)
THEN
2603 j = pw1%pw_grid%gidx(i)
2604 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
2610 cpabort(
"Grids not compatible")
2615 cpabort(
"Grids not compatible")
2619 CALL timestop(handle)
2621 END SUBROUTINE pw_axpy_r1d_r1d_gs
2632 SUBROUTINE pw_multiply_r1d_r1d_gs (pw_out, pw1, pw2, alpha)
2634 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw_out
2635 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
2636 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
2637 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
2639 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
2642 REAL(KIND=dp) :: my_alpha
2644 CALL timeset(routinen, handle)
2647 IF (
PRESENT(alpha)) my_alpha = alpha
2649 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
2650 cpabort(
"pw_multiply not implemented for non-identical grids!")
2652#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
2653 IF (my_alpha == 1.0_dp)
THEN
2655 pw_out%array = pw_out%array + pw1%array* pw2%array
2659 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
2663 IF (my_alpha == 1.0_dp)
THEN
2664 pw_out%array = pw_out%array + pw1%array* pw2%array
2666 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
2670 CALL timestop(handle)
2672 END SUBROUTINE pw_multiply_r1d_r1d_gs
2679 SUBROUTINE pw_multiply_with_r1d_r1d_gs (pw1, pw2)
2680 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw1
2681 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
2683 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
2687 CALL timeset(routinen, handle)
2689 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
2690 cpabort(
"Incompatible grids!")
2693 pw1%array = pw1%array* pw2%array
2696 CALL timestop(handle)
2698 END SUBROUTINE pw_multiply_with_r1d_r1d_gs
2713 FUNCTION pw_integral_ab_r1d_r1d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
2715 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
2716 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
2717 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
2718 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
2719 REAL(kind=dp) :: integral_value
2721 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_r1d_r1d_gs'
2723 INTEGER :: handle, loc_sumtype
2724 LOGICAL :: my_just_sum, my_local_only
2726 CALL timeset(routinen, handle)
2729 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
2731 my_local_only = .false.
2732 IF (
PRESENT(local_only)) my_local_only = local_only
2734 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
2735 cpabort(
"Grids incompatible")
2738 my_just_sum = .false.
2739 IF (
PRESENT(just_sum)) my_just_sum = just_sum
2746 integral_value = dot_product(pw1%array, pw2%array)
2751 integral_value = accurate_dot_product(pw1%array, pw2%array)
2755 IF (.NOT. my_just_sum)
THEN
2756 integral_value = integral_value*pw1%pw_grid%vol
2759 IF (pw1%pw_grid%grid_span == halfspace)
THEN
2760 integral_value = 2.0_dp*integral_value
2761 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
2762 pw1%array(1)*pw2%array(1)
2765 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
2766 CALL pw1%pw_grid%para%group%sum(integral_value)
2768 CALL timestop(handle)
2770 END FUNCTION pw_integral_ab_r1d_r1d_gs
2778 FUNCTION pw_integral_a2b_r1d_r1d (pw1, pw2)
RESULT(integral_value)
2780 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
2781 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
2782 REAL(kind=dp) :: integral_value
2784 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_a2b'
2788 CALL timeset(routinen, handle)
2790 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
2791 cpabort(
"Grids incompatible")
2794 integral_value = accurate_sum(pw1%array*pw2%array*pw1%pw_grid%gsq)
2795 IF (pw1%pw_grid%grid_span == halfspace)
THEN
2796 integral_value = 2.0_dp*integral_value
2799 integral_value = integral_value*pw1%pw_grid%vol
2801 IF (pw1%pw_grid%para%mode == pw_mode_distributed) &
2802 CALL pw1%pw_grid%para%group%sum(integral_value)
2803 CALL timestop(handle)
2805 END FUNCTION pw_integral_a2b_r1d_r1d
2819 SUBROUTINE pw_copy_r1d_c1d_rs (pw1, pw2)
2821 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
2822 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw2
2824 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
2827 INTEGER :: i, j, ng, ng1, ng2, ns
2829 CALL timeset(routinen, handle)
2831 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
2832 cpabort(
"Both grids must be either spherical or non-spherical!")
2833 IF (pw1%pw_grid%spherical) &
2834 cpabort(
"Spherical grids only exist in reciprocal space!")
2836 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
2837 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
2838 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
2839 ng1 =
SIZE(pw1%array)
2840 ng2 =
SIZE(pw2%array)
2843 pw2%array(1:ng) = cmplx(pw1%array(1:ng), 0.0_dp, kind=dp)
2847 pw2%array(ng + 1:ng2) = cmplx(0.0_dp, 0.0_dp, kind=dp)
2851 cpabort(
"Copies between spherical grids require compatible grids!")
2854 ng1 =
SIZE(pw1%array)
2855 ng2 =
SIZE(pw2%array)
2856 ns = 2*max(ng1, ng2)
2858 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
2859 IF (ng1 >= ng2)
THEN
2864 j = pw2%pw_grid%gidx(i)
2865 pw2%array(i) = cmplx(pw1%array(j), 0.0_dp, kind=dp)
2874 j = pw2%pw_grid%gidx(i)
2875 pw2%array(j) = cmplx(pw1%array(i), 0.0_dp, kind=dp)
2879 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
2880 IF (ng1 >= ng2)
THEN
2883 j = pw1%pw_grid%gidx(i)
2884 pw2%array(i) = cmplx(pw1%array(j), 0.0_dp, kind=dp)
2891 j = pw1%pw_grid%gidx(i)
2892 pw2%array(j) = cmplx(pw1%array(i), 0.0_dp, kind=dp)
2897 cpabort(
"Copy not implemented!")
2904 pw2%array = cmplx(pw1%array, 0.0_dp, kind=dp)
2908 CALL timestop(handle)
2910 END SUBROUTINE pw_copy_r1d_c1d_rs
2917 SUBROUTINE pw_copy_to_array_r1d_c1d_rs (pw, array)
2918 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw
2919 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(INOUT) :: array
2921 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
2925 CALL timeset(routinen, handle)
2928 array(:) = cmplx(pw%array(:), 0.0_dp, kind=dp)
2931 CALL timestop(handle)
2932 END SUBROUTINE pw_copy_to_array_r1d_c1d_rs
2939 SUBROUTINE pw_copy_from_array_r1d_c1d_rs (pw, array)
2940 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw
2941 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(IN) :: array
2943 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
2947 CALL timeset(routinen, handle)
2950 pw%array = real(array, kind=dp)
2953 CALL timestop(handle)
2954 END SUBROUTINE pw_copy_from_array_r1d_c1d_rs
2972 SUBROUTINE pw_axpy_r1d_c1d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
2974 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
2975 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw2
2976 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
2977 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
2979 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
2982 LOGICAL :: my_allow_noncompatible_grids
2983 REAL(KIND=dp) :: my_alpha, my_beta
2984 INTEGER :: i, j, ng, ng1, ng2
2986 CALL timeset(routinen, handle)
2989 IF (
PRESENT(alpha)) my_alpha = alpha
2992 IF (
PRESENT(beta)) my_beta = beta
2994 my_allow_noncompatible_grids = .false.
2995 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
2997 IF (my_beta /= 1.0_dp)
THEN
2998 IF (my_beta == 0.0_dp)
THEN
3002 pw2%array = pw2%array*my_beta
3007 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
3009 IF (my_alpha == 1.0_dp)
THEN
3011 pw2%array = pw2%array + cmplx(pw1%array, 0.0_dp, kind=dp)
3013 ELSE IF (my_alpha /= 0.0_dp)
THEN
3015 pw2%array = pw2%array + my_alpha* cmplx(pw1%array, 0.0_dp, kind=dp)
3019 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
3021 ng1 =
SIZE(pw1%array)
3022 ng2 =
SIZE(pw2%array)
3025 IF (pw1%pw_grid%spherical)
THEN
3026 IF (my_alpha == 1.0_dp)
THEN
3029 pw2%array(i) = pw2%array(i) + cmplx(pw1%array(i), 0.0_dp, kind=dp)
3032 ELSE IF (my_alpha /= 0.0_dp)
THEN
3035 pw2%array(i) = pw2%array(i) + my_alpha* cmplx(pw1%array(i), 0.0_dp, kind=dp)
3039 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
3040 IF (ng1 >= ng2)
THEN
3041 IF (my_alpha == 1.0_dp)
THEN
3044 j = pw2%pw_grid%gidx(i)
3045 pw2%array(i) = pw2%array(i) + cmplx(pw1%array(j), 0.0_dp, kind=dp)
3048 ELSE IF (my_alpha /= 0.0_dp)
THEN
3051 j = pw2%pw_grid%gidx(i)
3052 pw2%array(i) = pw2%array(i) + my_alpha* cmplx(pw1%array(j), 0.0_dp, kind=dp)
3057 IF (my_alpha == 1.0_dp)
THEN
3060 j = pw2%pw_grid%gidx(i)
3061 pw2%array(j) = pw2%array(j) + cmplx(pw1%array(i), 0.0_dp, kind=dp)
3064 ELSE IF (my_alpha /= 0.0_dp)
THEN
3067 j = pw2%pw_grid%gidx(i)
3068 pw2%array(j) = pw2%array(j) + my_alpha* cmplx(pw1%array(i), 0.0_dp, kind=dp)
3073 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
3074 IF (ng1 >= ng2)
THEN
3075 IF (my_alpha == 1.0_dp)
THEN
3078 j = pw1%pw_grid%gidx(i)
3079 pw2%array(i) = pw2%array(i) + cmplx(pw1%array(j), 0.0_dp, kind=dp)
3082 ELSE IF (my_alpha /= 0.0_dp)
THEN
3085 j = pw1%pw_grid%gidx(i)
3086 pw2%array(i) = pw2%array(i) + my_alpha* cmplx(pw1%array(j), 0.0_dp, kind=dp)
3091 IF (my_alpha == 1.0_dp)
THEN
3094 j = pw1%pw_grid%gidx(i)
3095 pw2%array(j) = pw2%array(j) + cmplx(pw1%array(i), 0.0_dp, kind=dp)
3098 ELSE IF (my_alpha /= 0.0_dp)
THEN
3101 j = pw1%pw_grid%gidx(i)
3102 pw2%array(j) = pw2%array(j) + my_alpha* cmplx(pw1%array(i), 0.0_dp, kind=dp)
3108 cpabort(
"Grids not compatible")
3113 cpabort(
"Grids not compatible")
3117 CALL timestop(handle)
3119 END SUBROUTINE pw_axpy_r1d_c1d_rs
3130 SUBROUTINE pw_multiply_r1d_c1d_rs (pw_out, pw1, pw2, alpha)
3132 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw_out
3133 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
3134 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw2
3135 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
3137 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
3140 REAL(KIND=dp) :: my_alpha
3142 CALL timeset(routinen, handle)
3145 IF (
PRESENT(alpha)) my_alpha = alpha
3147 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
3148 cpabort(
"pw_multiply not implemented for non-identical grids!")
3150#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
3151 IF (my_alpha == 1.0_dp)
THEN
3153 pw_out%array = pw_out%array + pw1%array* real(pw2%array, kind=dp)
3157 pw_out%array = pw_out%array + my_alpha*pw1%array* real(pw2%array, kind=dp)
3161 IF (my_alpha == 1.0_dp)
THEN
3162 pw_out%array = pw_out%array + pw1%array* real(pw2%array, kind=dp)
3164 pw_out%array = pw_out%array + my_alpha*pw1%array* real(pw2%array, kind=dp)
3168 CALL timestop(handle)
3170 END SUBROUTINE pw_multiply_r1d_c1d_rs
3177 SUBROUTINE pw_multiply_with_r1d_c1d_rs (pw1, pw2)
3178 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw1
3179 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw2
3181 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
3185 CALL timeset(routinen, handle)
3187 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
3188 cpabort(
"Incompatible grids!")
3191 pw1%array = pw1%array* real(pw2%array, kind=dp)
3194 CALL timestop(handle)
3196 END SUBROUTINE pw_multiply_with_r1d_c1d_rs
3211 FUNCTION pw_integral_ab_r1d_c1d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
3213 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
3214 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw2
3215 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
3216 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
3217 REAL(kind=dp) :: integral_value
3219 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_r1d_c1d_rs'
3221 INTEGER :: handle, loc_sumtype
3222 LOGICAL :: my_just_sum, my_local_only
3224 CALL timeset(routinen, handle)
3227 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
3229 my_local_only = .false.
3230 IF (
PRESENT(local_only)) my_local_only = local_only
3232 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
3233 cpabort(
"Grids incompatible")
3236 my_just_sum = .false.
3237 IF (
PRESENT(just_sum)) my_just_sum = just_sum
3244 integral_value = sum(pw1%array*real(pw2%array, kind=dp))
3249 integral_value = accurate_sum(pw1%array*real(pw2%array, kind=dp))
3253 IF (.NOT. my_just_sum)
THEN
3254 integral_value = integral_value*pw1%pw_grid%dvol
3257 IF (pw1%pw_grid%grid_span == halfspace)
THEN
3258 integral_value = 2.0_dp*integral_value
3259 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
3260 pw1%array(1)*real(pw2%array(1), kind=dp)
3263 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
3264 CALL pw1%pw_grid%para%group%sum(integral_value)
3266 CALL timestop(handle)
3268 END FUNCTION pw_integral_ab_r1d_c1d_rs
3282 SUBROUTINE pw_copy_r1d_c1d_gs (pw1, pw2)
3284 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
3285 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
3287 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
3290 INTEGER :: i, j, ng, ng1, ng2, ns
3292 CALL timeset(routinen, handle)
3294 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
3295 cpabort(
"Both grids must be either spherical or non-spherical!")
3297 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
3298 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
3299 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
3300 ng1 =
SIZE(pw1%array)
3301 ng2 =
SIZE(pw2%array)
3304 pw2%array(1:ng) = cmplx(pw1%array(1:ng), 0.0_dp, kind=dp)
3308 pw2%array(ng + 1:ng2) = cmplx(0.0_dp, 0.0_dp, kind=dp)
3312 cpabort(
"Copies between spherical grids require compatible grids!")
3315 ng1 =
SIZE(pw1%array)
3316 ng2 =
SIZE(pw2%array)
3317 ns = 2*max(ng1, ng2)
3319 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
3320 IF (ng1 >= ng2)
THEN
3325 j = pw2%pw_grid%gidx(i)
3326 pw2%array(i) = cmplx(pw1%array(j), 0.0_dp, kind=dp)
3335 j = pw2%pw_grid%gidx(i)
3336 pw2%array(j) = cmplx(pw1%array(i), 0.0_dp, kind=dp)
3340 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
3341 IF (ng1 >= ng2)
THEN
3344 j = pw1%pw_grid%gidx(i)
3345 pw2%array(i) = cmplx(pw1%array(j), 0.0_dp, kind=dp)
3352 j = pw1%pw_grid%gidx(i)
3353 pw2%array(j) = cmplx(pw1%array(i), 0.0_dp, kind=dp)
3358 cpabort(
"Copy not implemented!")
3365 pw2%array = cmplx(pw1%array, 0.0_dp, kind=dp)
3369 CALL timestop(handle)
3371 END SUBROUTINE pw_copy_r1d_c1d_gs
3378 SUBROUTINE pw_copy_to_array_r1d_c1d_gs (pw, array)
3379 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw
3380 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(INOUT) :: array
3382 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
3386 CALL timeset(routinen, handle)
3389 array(:) = cmplx(pw%array(:), 0.0_dp, kind=dp)
3392 CALL timestop(handle)
3393 END SUBROUTINE pw_copy_to_array_r1d_c1d_gs
3400 SUBROUTINE pw_copy_from_array_r1d_c1d_gs (pw, array)
3401 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw
3402 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(IN) :: array
3404 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
3408 CALL timeset(routinen, handle)
3411 pw%array = real(array, kind=dp)
3414 CALL timestop(handle)
3415 END SUBROUTINE pw_copy_from_array_r1d_c1d_gs
3433 SUBROUTINE pw_axpy_r1d_c1d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
3435 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
3436 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
3437 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
3438 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
3440 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
3443 LOGICAL :: my_allow_noncompatible_grids
3444 REAL(KIND=dp) :: my_alpha, my_beta
3445 INTEGER :: i, j, ng, ng1, ng2
3447 CALL timeset(routinen, handle)
3450 IF (
PRESENT(alpha)) my_alpha = alpha
3453 IF (
PRESENT(beta)) my_beta = beta
3455 my_allow_noncompatible_grids = .false.
3456 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
3458 IF (my_beta /= 1.0_dp)
THEN
3459 IF (my_beta == 0.0_dp)
THEN
3463 pw2%array = pw2%array*my_beta
3468 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
3470 IF (my_alpha == 1.0_dp)
THEN
3472 pw2%array = pw2%array + cmplx(pw1%array, 0.0_dp, kind=dp)
3474 ELSE IF (my_alpha /= 0.0_dp)
THEN
3476 pw2%array = pw2%array + my_alpha* cmplx(pw1%array, 0.0_dp, kind=dp)
3480 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
3482 ng1 =
SIZE(pw1%array)
3483 ng2 =
SIZE(pw2%array)
3486 IF (pw1%pw_grid%spherical)
THEN
3487 IF (my_alpha == 1.0_dp)
THEN
3490 pw2%array(i) = pw2%array(i) + cmplx(pw1%array(i), 0.0_dp, kind=dp)
3493 ELSE IF (my_alpha /= 0.0_dp)
THEN
3496 pw2%array(i) = pw2%array(i) + my_alpha* cmplx(pw1%array(i), 0.0_dp, kind=dp)
3500 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
3501 IF (ng1 >= ng2)
THEN
3502 IF (my_alpha == 1.0_dp)
THEN
3505 j = pw2%pw_grid%gidx(i)
3506 pw2%array(i) = pw2%array(i) + cmplx(pw1%array(j), 0.0_dp, kind=dp)
3509 ELSE IF (my_alpha /= 0.0_dp)
THEN
3512 j = pw2%pw_grid%gidx(i)
3513 pw2%array(i) = pw2%array(i) + my_alpha* cmplx(pw1%array(j), 0.0_dp, kind=dp)
3518 IF (my_alpha == 1.0_dp)
THEN
3521 j = pw2%pw_grid%gidx(i)
3522 pw2%array(j) = pw2%array(j) + cmplx(pw1%array(i), 0.0_dp, kind=dp)
3525 ELSE IF (my_alpha /= 0.0_dp)
THEN
3528 j = pw2%pw_grid%gidx(i)
3529 pw2%array(j) = pw2%array(j) + my_alpha* cmplx(pw1%array(i), 0.0_dp, kind=dp)
3534 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
3535 IF (ng1 >= ng2)
THEN
3536 IF (my_alpha == 1.0_dp)
THEN
3539 j = pw1%pw_grid%gidx(i)
3540 pw2%array(i) = pw2%array(i) + cmplx(pw1%array(j), 0.0_dp, kind=dp)
3543 ELSE IF (my_alpha /= 0.0_dp)
THEN
3546 j = pw1%pw_grid%gidx(i)
3547 pw2%array(i) = pw2%array(i) + my_alpha* cmplx(pw1%array(j), 0.0_dp, kind=dp)
3552 IF (my_alpha == 1.0_dp)
THEN
3555 j = pw1%pw_grid%gidx(i)
3556 pw2%array(j) = pw2%array(j) + cmplx(pw1%array(i), 0.0_dp, kind=dp)
3559 ELSE IF (my_alpha /= 0.0_dp)
THEN
3562 j = pw1%pw_grid%gidx(i)
3563 pw2%array(j) = pw2%array(j) + my_alpha* cmplx(pw1%array(i), 0.0_dp, kind=dp)
3569 cpabort(
"Grids not compatible")
3574 cpabort(
"Grids not compatible")
3578 CALL timestop(handle)
3580 END SUBROUTINE pw_axpy_r1d_c1d_gs
3591 SUBROUTINE pw_multiply_r1d_c1d_gs (pw_out, pw1, pw2, alpha)
3593 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw_out
3594 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
3595 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
3596 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
3598 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
3601 REAL(KIND=dp) :: my_alpha
3603 CALL timeset(routinen, handle)
3606 IF (
PRESENT(alpha)) my_alpha = alpha
3608 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
3609 cpabort(
"pw_multiply not implemented for non-identical grids!")
3611#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
3612 IF (my_alpha == 1.0_dp)
THEN
3614 pw_out%array = pw_out%array + pw1%array* real(pw2%array, kind=dp)
3618 pw_out%array = pw_out%array + my_alpha*pw1%array* real(pw2%array, kind=dp)
3622 IF (my_alpha == 1.0_dp)
THEN
3623 pw_out%array = pw_out%array + pw1%array* real(pw2%array, kind=dp)
3625 pw_out%array = pw_out%array + my_alpha*pw1%array* real(pw2%array, kind=dp)
3629 CALL timestop(handle)
3631 END SUBROUTINE pw_multiply_r1d_c1d_gs
3638 SUBROUTINE pw_multiply_with_r1d_c1d_gs (pw1, pw2)
3639 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw1
3640 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
3642 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
3646 CALL timeset(routinen, handle)
3648 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
3649 cpabort(
"Incompatible grids!")
3652 pw1%array = pw1%array* real(pw2%array, kind=dp)
3655 CALL timestop(handle)
3657 END SUBROUTINE pw_multiply_with_r1d_c1d_gs
3672 FUNCTION pw_integral_ab_r1d_c1d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
3674 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
3675 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
3676 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
3677 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
3678 REAL(kind=dp) :: integral_value
3680 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_r1d_c1d_gs'
3682 INTEGER :: handle, loc_sumtype
3683 LOGICAL :: my_just_sum, my_local_only
3685 CALL timeset(routinen, handle)
3688 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
3690 my_local_only = .false.
3691 IF (
PRESENT(local_only)) my_local_only = local_only
3693 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
3694 cpabort(
"Grids incompatible")
3697 my_just_sum = .false.
3698 IF (
PRESENT(just_sum)) my_just_sum = just_sum
3705 integral_value = sum(pw1%array*real(pw2%array, kind=dp))
3710 integral_value = accurate_sum(pw1%array*real(pw2%array, kind=dp))
3714 IF (.NOT. my_just_sum)
THEN
3715 integral_value = integral_value*pw1%pw_grid%vol
3718 IF (pw1%pw_grid%grid_span == halfspace)
THEN
3719 integral_value = 2.0_dp*integral_value
3720 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
3721 pw1%array(1)*real(pw2%array(1), kind=dp)
3724 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
3725 CALL pw1%pw_grid%para%group%sum(integral_value)
3727 CALL timestop(handle)
3729 END FUNCTION pw_integral_ab_r1d_c1d_gs
3737 FUNCTION pw_integral_a2b_r1d_c1d (pw1, pw2)
RESULT(integral_value)
3739 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
3740 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
3741 REAL(kind=dp) :: integral_value
3743 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_a2b'
3747 CALL timeset(routinen, handle)
3749 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
3750 cpabort(
"Grids incompatible")
3753 integral_value = accurate_sum(pw1%array*real(pw2%array, kind=dp)*pw1%pw_grid%gsq)
3754 IF (pw1%pw_grid%grid_span == halfspace)
THEN
3755 integral_value = 2.0_dp*integral_value
3758 integral_value = integral_value*pw1%pw_grid%vol
3760 IF (pw1%pw_grid%para%mode == pw_mode_distributed) &
3761 CALL pw1%pw_grid%para%group%sum(integral_value)
3762 CALL timestop(handle)
3764 END FUNCTION pw_integral_a2b_r1d_c1d
3778 SUBROUTINE pw_copy_r3d_r3d_rs (pw1, pw2)
3780 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
3781 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw2
3783 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
3787 CALL timeset(routinen, handle)
3789 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
3790 cpabort(
"Both grids must be either spherical or non-spherical!")
3791 IF (pw1%pw_grid%spherical) &
3792 cpabort(
"Spherical grids only exist in reciprocal space!")
3794 IF (any(shape(pw2%array) /= shape(pw1%array))) &
3795 cpabort(
"3D grids must be compatible!")
3796 IF (pw1%pw_grid%spherical) &
3797 cpabort(
"3D grids must not be spherical!")
3799 pw2%array(:, :, :) = pw1%array(:, :, :)
3802 CALL timestop(handle)
3804 END SUBROUTINE pw_copy_r3d_r3d_rs
3811 SUBROUTINE pw_copy_to_array_r3d_r3d_rs (pw, array)
3812 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw
3813 REAL(KIND=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
3815 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
3819 CALL timeset(routinen, handle)
3822 array(:, :, :) = pw%array(:, :, :)
3825 CALL timestop(handle)
3826 END SUBROUTINE pw_copy_to_array_r3d_r3d_rs
3833 SUBROUTINE pw_copy_from_array_r3d_r3d_rs (pw, array)
3834 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw
3835 REAL(KIND=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
3837 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
3841 CALL timeset(routinen, handle)
3847 CALL timestop(handle)
3848 END SUBROUTINE pw_copy_from_array_r3d_r3d_rs
3866 SUBROUTINE pw_axpy_r3d_r3d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
3868 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
3869 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw2
3870 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
3871 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
3873 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
3876 LOGICAL :: my_allow_noncompatible_grids
3877 REAL(KIND=dp) :: my_alpha, my_beta
3879 CALL timeset(routinen, handle)
3882 IF (
PRESENT(alpha)) my_alpha = alpha
3885 IF (
PRESENT(beta)) my_beta = beta
3887 my_allow_noncompatible_grids = .false.
3888 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
3890 IF (my_beta /= 1.0_dp)
THEN
3891 IF (my_beta == 0.0_dp)
THEN
3895 pw2%array = pw2%array*my_beta
3900 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
3901 IF (my_alpha == 1.0_dp)
THEN
3903 pw2%array = pw2%array + pw1%array
3905 ELSE IF (my_alpha /= 0.0_dp)
THEN
3907 pw2%array = pw2%array + my_alpha* pw1%array
3911 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
3913 IF (any(shape(pw1%array) /= shape(pw2%array))) &
3914 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
3916 IF (my_alpha == 1.0_dp)
THEN
3918 pw2%array = pw2%array + pw1%array
3922 pw2%array = pw2%array + my_alpha* pw1%array
3928 cpabort(
"Grids not compatible")
3932 CALL timestop(handle)
3934 END SUBROUTINE pw_axpy_r3d_r3d_rs
3945 SUBROUTINE pw_multiply_r3d_r3d_rs (pw_out, pw1, pw2, alpha)
3947 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw_out
3948 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
3949 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw2
3950 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
3952 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
3955 REAL(KIND=dp) :: my_alpha
3957 CALL timeset(routinen, handle)
3960 IF (
PRESENT(alpha)) my_alpha = alpha
3962 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
3963 cpabort(
"pw_multiply not implemented for non-identical grids!")
3965#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
3966 IF (my_alpha == 1.0_dp)
THEN
3968 pw_out%array = pw_out%array + pw1%array* pw2%array
3972 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
3976 IF (my_alpha == 1.0_dp)
THEN
3977 pw_out%array = pw_out%array + pw1%array* pw2%array
3979 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
3983 CALL timestop(handle)
3985 END SUBROUTINE pw_multiply_r3d_r3d_rs
3992 SUBROUTINE pw_multiply_with_r3d_r3d_rs (pw1, pw2)
3993 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw1
3994 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw2
3996 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
4000 CALL timeset(routinen, handle)
4002 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
4003 cpabort(
"Incompatible grids!")
4006 pw1%array = pw1%array* pw2%array
4009 CALL timestop(handle)
4011 END SUBROUTINE pw_multiply_with_r3d_r3d_rs
4026 FUNCTION pw_integral_ab_r3d_r3d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
4028 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
4029 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw2
4030 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
4031 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
4032 REAL(kind=dp) :: integral_value
4034 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_r3d_r3d_rs'
4036 INTEGER :: handle, loc_sumtype
4037 LOGICAL :: my_just_sum, my_local_only
4039 CALL timeset(routinen, handle)
4042 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
4044 my_local_only = .false.
4045 IF (
PRESENT(local_only)) my_local_only = local_only
4047 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
4048 cpabort(
"Grids incompatible")
4051 my_just_sum = .false.
4052 IF (
PRESENT(just_sum)) my_just_sum = just_sum
4059 integral_value = sum(pw1%array*pw2%array)
4064 integral_value = accurate_dot_product(pw1%array, pw2%array)
4068 IF (.NOT. my_just_sum)
THEN
4069 integral_value = integral_value*pw1%pw_grid%dvol
4073 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
4074 CALL pw1%pw_grid%para%group%sum(integral_value)
4076 CALL timestop(handle)
4078 END FUNCTION pw_integral_ab_r3d_r3d_rs
4092 SUBROUTINE pw_copy_r3d_r3d_gs (pw1, pw2)
4094 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4095 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw2
4097 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
4101 CALL timeset(routinen, handle)
4103 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
4104 cpabort(
"Both grids must be either spherical or non-spherical!")
4106 IF (any(shape(pw2%array) /= shape(pw1%array))) &
4107 cpabort(
"3D grids must be compatible!")
4108 IF (pw1%pw_grid%spherical) &
4109 cpabort(
"3D grids must not be spherical!")
4111 pw2%array(:, :, :) = pw1%array(:, :, :)
4114 CALL timestop(handle)
4116 END SUBROUTINE pw_copy_r3d_r3d_gs
4123 SUBROUTINE pw_copy_to_array_r3d_r3d_gs (pw, array)
4124 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw
4125 REAL(KIND=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
4127 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
4131 CALL timeset(routinen, handle)
4134 array(:, :, :) = pw%array(:, :, :)
4137 CALL timestop(handle)
4138 END SUBROUTINE pw_copy_to_array_r3d_r3d_gs
4145 SUBROUTINE pw_copy_from_array_r3d_r3d_gs (pw, array)
4146 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw
4147 REAL(KIND=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
4149 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
4153 CALL timeset(routinen, handle)
4159 CALL timestop(handle)
4160 END SUBROUTINE pw_copy_from_array_r3d_r3d_gs
4178 SUBROUTINE pw_axpy_r3d_r3d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
4180 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4181 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw2
4182 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
4183 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
4185 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
4188 LOGICAL :: my_allow_noncompatible_grids
4189 REAL(KIND=dp) :: my_alpha, my_beta
4191 CALL timeset(routinen, handle)
4194 IF (
PRESENT(alpha)) my_alpha = alpha
4197 IF (
PRESENT(beta)) my_beta = beta
4199 my_allow_noncompatible_grids = .false.
4200 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
4202 IF (my_beta /= 1.0_dp)
THEN
4203 IF (my_beta == 0.0_dp)
THEN
4207 pw2%array = pw2%array*my_beta
4212 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
4213 IF (my_alpha == 1.0_dp)
THEN
4215 pw2%array = pw2%array + pw1%array
4217 ELSE IF (my_alpha /= 0.0_dp)
THEN
4219 pw2%array = pw2%array + my_alpha* pw1%array
4223 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
4225 IF (any(shape(pw1%array) /= shape(pw2%array))) &
4226 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
4228 IF (my_alpha == 1.0_dp)
THEN
4230 pw2%array = pw2%array + pw1%array
4234 pw2%array = pw2%array + my_alpha* pw1%array
4240 cpabort(
"Grids not compatible")
4244 CALL timestop(handle)
4246 END SUBROUTINE pw_axpy_r3d_r3d_gs
4257 SUBROUTINE pw_multiply_r3d_r3d_gs (pw_out, pw1, pw2, alpha)
4259 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw_out
4260 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4261 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw2
4262 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
4264 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
4267 REAL(KIND=dp) :: my_alpha
4269 CALL timeset(routinen, handle)
4272 IF (
PRESENT(alpha)) my_alpha = alpha
4274 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
4275 cpabort(
"pw_multiply not implemented for non-identical grids!")
4277#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
4278 IF (my_alpha == 1.0_dp)
THEN
4280 pw_out%array = pw_out%array + pw1%array* pw2%array
4284 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
4288 IF (my_alpha == 1.0_dp)
THEN
4289 pw_out%array = pw_out%array + pw1%array* pw2%array
4291 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
4295 CALL timestop(handle)
4297 END SUBROUTINE pw_multiply_r3d_r3d_gs
4304 SUBROUTINE pw_multiply_with_r3d_r3d_gs (pw1, pw2)
4305 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw1
4306 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw2
4308 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
4312 CALL timeset(routinen, handle)
4314 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
4315 cpabort(
"Incompatible grids!")
4318 pw1%array = pw1%array* pw2%array
4321 CALL timestop(handle)
4323 END SUBROUTINE pw_multiply_with_r3d_r3d_gs
4338 FUNCTION pw_integral_ab_r3d_r3d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
4340 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4341 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw2
4342 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
4343 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
4344 REAL(kind=dp) :: integral_value
4346 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_r3d_r3d_gs'
4348 INTEGER :: handle, loc_sumtype
4349 LOGICAL :: my_just_sum, my_local_only
4351 CALL timeset(routinen, handle)
4354 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
4356 my_local_only = .false.
4357 IF (
PRESENT(local_only)) my_local_only = local_only
4359 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
4360 cpabort(
"Grids incompatible")
4363 my_just_sum = .false.
4364 IF (
PRESENT(just_sum)) my_just_sum = just_sum
4371 integral_value = sum(pw1%array*pw2%array)
4376 integral_value = accurate_dot_product(pw1%array, pw2%array)
4380 IF (.NOT. my_just_sum)
THEN
4381 integral_value = integral_value*pw1%pw_grid%vol
4385 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
4386 CALL pw1%pw_grid%para%group%sum(integral_value)
4388 CALL timestop(handle)
4390 END FUNCTION pw_integral_ab_r3d_r3d_gs
4405 SUBROUTINE pw_copy_r3d_c3d_rs (pw1, pw2)
4407 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
4408 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw2
4410 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
4414 CALL timeset(routinen, handle)
4416 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
4417 cpabort(
"Both grids must be either spherical or non-spherical!")
4418 IF (pw1%pw_grid%spherical) &
4419 cpabort(
"Spherical grids only exist in reciprocal space!")
4421 IF (any(shape(pw2%array) /= shape(pw1%array))) &
4422 cpabort(
"3D grids must be compatible!")
4423 IF (pw1%pw_grid%spherical) &
4424 cpabort(
"3D grids must not be spherical!")
4426 pw2%array(:, :, :) = cmplx(pw1%array(:, :, :), 0.0_dp, kind=dp)
4429 CALL timestop(handle)
4431 END SUBROUTINE pw_copy_r3d_c3d_rs
4438 SUBROUTINE pw_copy_to_array_r3d_c3d_rs (pw, array)
4439 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw
4440 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
4442 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
4446 CALL timeset(routinen, handle)
4449 array(:, :, :) = cmplx(pw%array(:, :, :), 0.0_dp, kind=dp)
4452 CALL timestop(handle)
4453 END SUBROUTINE pw_copy_to_array_r3d_c3d_rs
4460 SUBROUTINE pw_copy_from_array_r3d_c3d_rs (pw, array)
4461 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw
4462 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
4464 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
4468 CALL timeset(routinen, handle)
4471 pw%array = real(array, kind=dp)
4474 CALL timestop(handle)
4475 END SUBROUTINE pw_copy_from_array_r3d_c3d_rs
4493 SUBROUTINE pw_axpy_r3d_c3d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
4495 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
4496 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw2
4497 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
4498 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
4500 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
4503 LOGICAL :: my_allow_noncompatible_grids
4504 REAL(KIND=dp) :: my_alpha, my_beta
4506 CALL timeset(routinen, handle)
4509 IF (
PRESENT(alpha)) my_alpha = alpha
4512 IF (
PRESENT(beta)) my_beta = beta
4514 my_allow_noncompatible_grids = .false.
4515 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
4517 IF (my_beta /= 1.0_dp)
THEN
4518 IF (my_beta == 0.0_dp)
THEN
4522 pw2%array = pw2%array*my_beta
4527 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
4528 IF (my_alpha == 1.0_dp)
THEN
4530 pw2%array = pw2%array + cmplx(pw1%array, 0.0_dp, kind=dp)
4532 ELSE IF (my_alpha /= 0.0_dp)
THEN
4534 pw2%array = pw2%array + my_alpha* cmplx(pw1%array, 0.0_dp, kind=dp)
4538 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
4540 IF (any(shape(pw1%array) /= shape(pw2%array))) &
4541 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
4543 IF (my_alpha == 1.0_dp)
THEN
4545 pw2%array = pw2%array + cmplx(pw1%array, 0.0_dp, kind=dp)
4549 pw2%array = pw2%array + my_alpha* cmplx(pw1%array, 0.0_dp, kind=dp)
4555 cpabort(
"Grids not compatible")
4559 CALL timestop(handle)
4561 END SUBROUTINE pw_axpy_r3d_c3d_rs
4572 SUBROUTINE pw_multiply_r3d_c3d_rs (pw_out, pw1, pw2, alpha)
4574 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw_out
4575 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
4576 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw2
4577 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
4579 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
4582 REAL(KIND=dp) :: my_alpha
4584 CALL timeset(routinen, handle)
4587 IF (
PRESENT(alpha)) my_alpha = alpha
4589 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
4590 cpabort(
"pw_multiply not implemented for non-identical grids!")
4592#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
4593 IF (my_alpha == 1.0_dp)
THEN
4595 pw_out%array = pw_out%array + pw1%array* real(pw2%array, kind=dp)
4599 pw_out%array = pw_out%array + my_alpha*pw1%array* real(pw2%array, kind=dp)
4603 IF (my_alpha == 1.0_dp)
THEN
4604 pw_out%array = pw_out%array + pw1%array* real(pw2%array, kind=dp)
4606 pw_out%array = pw_out%array + my_alpha*pw1%array* real(pw2%array, kind=dp)
4610 CALL timestop(handle)
4612 END SUBROUTINE pw_multiply_r3d_c3d_rs
4619 SUBROUTINE pw_multiply_with_r3d_c3d_rs (pw1, pw2)
4620 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw1
4621 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw2
4623 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
4627 CALL timeset(routinen, handle)
4629 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
4630 cpabort(
"Incompatible grids!")
4633 pw1%array = pw1%array* real(pw2%array, kind=dp)
4636 CALL timestop(handle)
4638 END SUBROUTINE pw_multiply_with_r3d_c3d_rs
4653 FUNCTION pw_integral_ab_r3d_c3d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
4655 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
4656 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw2
4657 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
4658 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
4659 REAL(kind=dp) :: integral_value
4661 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_r3d_c3d_rs'
4663 INTEGER :: handle, loc_sumtype
4664 LOGICAL :: my_just_sum, my_local_only
4666 CALL timeset(routinen, handle)
4669 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
4671 my_local_only = .false.
4672 IF (
PRESENT(local_only)) my_local_only = local_only
4674 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
4675 cpabort(
"Grids incompatible")
4678 my_just_sum = .false.
4679 IF (
PRESENT(just_sum)) my_just_sum = just_sum
4686 integral_value = sum(pw1%array*real(pw2%array, kind=dp))
4691 integral_value = accurate_sum(pw1%array*real(pw2%array, kind=dp))
4695 IF (.NOT. my_just_sum)
THEN
4696 integral_value = integral_value*pw1%pw_grid%dvol
4700 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
4701 CALL pw1%pw_grid%para%group%sum(integral_value)
4703 CALL timestop(handle)
4705 END FUNCTION pw_integral_ab_r3d_c3d_rs
4719 SUBROUTINE pw_copy_r3d_c3d_gs (pw1, pw2)
4721 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4722 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
4724 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
4728 CALL timeset(routinen, handle)
4730 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
4731 cpabort(
"Both grids must be either spherical or non-spherical!")
4733 IF (any(shape(pw2%array) /= shape(pw1%array))) &
4734 cpabort(
"3D grids must be compatible!")
4735 IF (pw1%pw_grid%spherical) &
4736 cpabort(
"3D grids must not be spherical!")
4738 pw2%array(:, :, :) = cmplx(pw1%array(:, :, :), 0.0_dp, kind=dp)
4741 CALL timestop(handle)
4743 END SUBROUTINE pw_copy_r3d_c3d_gs
4750 SUBROUTINE pw_copy_to_array_r3d_c3d_gs (pw, array)
4751 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw
4752 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
4754 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
4758 CALL timeset(routinen, handle)
4761 array(:, :, :) = cmplx(pw%array(:, :, :), 0.0_dp, kind=dp)
4764 CALL timestop(handle)
4765 END SUBROUTINE pw_copy_to_array_r3d_c3d_gs
4772 SUBROUTINE pw_copy_from_array_r3d_c3d_gs (pw, array)
4773 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw
4774 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
4776 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
4780 CALL timeset(routinen, handle)
4783 pw%array = real(array, kind=dp)
4786 CALL timestop(handle)
4787 END SUBROUTINE pw_copy_from_array_r3d_c3d_gs
4805 SUBROUTINE pw_axpy_r3d_c3d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
4807 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4808 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
4809 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
4810 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
4812 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
4815 LOGICAL :: my_allow_noncompatible_grids
4816 REAL(KIND=dp) :: my_alpha, my_beta
4818 CALL timeset(routinen, handle)
4821 IF (
PRESENT(alpha)) my_alpha = alpha
4824 IF (
PRESENT(beta)) my_beta = beta
4826 my_allow_noncompatible_grids = .false.
4827 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
4829 IF (my_beta /= 1.0_dp)
THEN
4830 IF (my_beta == 0.0_dp)
THEN
4834 pw2%array = pw2%array*my_beta
4839 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
4840 IF (my_alpha == 1.0_dp)
THEN
4842 pw2%array = pw2%array + cmplx(pw1%array, 0.0_dp, kind=dp)
4844 ELSE IF (my_alpha /= 0.0_dp)
THEN
4846 pw2%array = pw2%array + my_alpha* cmplx(pw1%array, 0.0_dp, kind=dp)
4850 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
4852 IF (any(shape(pw1%array) /= shape(pw2%array))) &
4853 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
4855 IF (my_alpha == 1.0_dp)
THEN
4857 pw2%array = pw2%array + cmplx(pw1%array, 0.0_dp, kind=dp)
4861 pw2%array = pw2%array + my_alpha* cmplx(pw1%array, 0.0_dp, kind=dp)
4867 cpabort(
"Grids not compatible")
4871 CALL timestop(handle)
4873 END SUBROUTINE pw_axpy_r3d_c3d_gs
4884 SUBROUTINE pw_multiply_r3d_c3d_gs (pw_out, pw1, pw2, alpha)
4886 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw_out
4887 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4888 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw2
4889 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
4891 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
4894 REAL(KIND=dp) :: my_alpha
4896 CALL timeset(routinen, handle)
4899 IF (
PRESENT(alpha)) my_alpha = alpha
4901 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
4902 cpabort(
"pw_multiply not implemented for non-identical grids!")
4904#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
4905 IF (my_alpha == 1.0_dp)
THEN
4907 pw_out%array = pw_out%array + pw1%array* real(pw2%array, kind=dp)
4911 pw_out%array = pw_out%array + my_alpha*pw1%array* real(pw2%array, kind=dp)
4915 IF (my_alpha == 1.0_dp)
THEN
4916 pw_out%array = pw_out%array + pw1%array* real(pw2%array, kind=dp)
4918 pw_out%array = pw_out%array + my_alpha*pw1%array* real(pw2%array, kind=dp)
4922 CALL timestop(handle)
4924 END SUBROUTINE pw_multiply_r3d_c3d_gs
4931 SUBROUTINE pw_multiply_with_r3d_c3d_gs (pw1, pw2)
4932 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw1
4933 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw2
4935 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
4939 CALL timeset(routinen, handle)
4941 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
4942 cpabort(
"Incompatible grids!")
4945 pw1%array = pw1%array* real(pw2%array, kind=dp)
4948 CALL timestop(handle)
4950 END SUBROUTINE pw_multiply_with_r3d_c3d_gs
4965 FUNCTION pw_integral_ab_r3d_c3d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
4967 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4968 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw2
4969 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
4970 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
4971 REAL(kind=dp) :: integral_value
4973 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_r3d_c3d_gs'
4975 INTEGER :: handle, loc_sumtype
4976 LOGICAL :: my_just_sum, my_local_only
4978 CALL timeset(routinen, handle)
4981 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
4983 my_local_only = .false.
4984 IF (
PRESENT(local_only)) my_local_only = local_only
4986 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
4987 cpabort(
"Grids incompatible")
4990 my_just_sum = .false.
4991 IF (
PRESENT(just_sum)) my_just_sum = just_sum
4998 integral_value = sum(pw1%array*real(pw2%array, kind=dp))
5003 integral_value = accurate_sum(pw1%array*real(pw2%array, kind=dp))
5007 IF (.NOT. my_just_sum)
THEN
5008 integral_value = integral_value*pw1%pw_grid%vol
5012 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
5013 CALL pw1%pw_grid%para%group%sum(integral_value)
5015 CALL timestop(handle)
5017 END FUNCTION pw_integral_ab_r3d_c3d_gs
5032 SUBROUTINE pw_copy_c1d_r1d_rs (pw1, pw2)
5034 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
5035 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw2
5037 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
5040 INTEGER :: i, j, ng, ng1, ng2, ns
5042 CALL timeset(routinen, handle)
5044 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
5045 cpabort(
"Both grids must be either spherical or non-spherical!")
5046 IF (pw1%pw_grid%spherical) &
5047 cpabort(
"Spherical grids only exist in reciprocal space!")
5049 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
5050 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
5051 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
5052 ng1 =
SIZE(pw1%array)
5053 ng2 =
SIZE(pw2%array)
5056 pw2%array(1:ng) = real(pw1%array(1:ng), kind=dp)
5060 pw2%array(ng + 1:ng2) = 0.0_dp
5064 cpabort(
"Copies between spherical grids require compatible grids!")
5067 ng1 =
SIZE(pw1%array)
5068 ng2 =
SIZE(pw2%array)
5069 ns = 2*max(ng1, ng2)
5071 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
5072 IF (ng1 >= ng2)
THEN
5077 j = pw2%pw_grid%gidx(i)
5078 pw2%array(i) = real(pw1%array(j), kind=dp)
5087 j = pw2%pw_grid%gidx(i)
5088 pw2%array(j) = real(pw1%array(i), kind=dp)
5092 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
5093 IF (ng1 >= ng2)
THEN
5096 j = pw1%pw_grid%gidx(i)
5097 pw2%array(i) = real(pw1%array(j), kind=dp)
5104 j = pw1%pw_grid%gidx(i)
5105 pw2%array(j) = real(pw1%array(i), kind=dp)
5110 cpabort(
"Copy not implemented!")
5117 pw2%array = real(pw1%array, kind=dp)
5121 CALL timestop(handle)
5123 END SUBROUTINE pw_copy_c1d_r1d_rs
5130 SUBROUTINE pw_copy_to_array_c1d_r1d_rs (pw, array)
5131 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw
5132 REAL(KIND=dp),
DIMENSION(:),
INTENT(INOUT) :: array
5134 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
5138 CALL timeset(routinen, handle)
5141 array(:) = real(pw%array(:), kind=dp)
5144 CALL timestop(handle)
5145 END SUBROUTINE pw_copy_to_array_c1d_r1d_rs
5152 SUBROUTINE pw_copy_from_array_c1d_r1d_rs (pw, array)
5153 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw
5154 REAL(KIND=dp),
DIMENSION(:),
INTENT(IN) :: array
5156 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
5160 CALL timeset(routinen, handle)
5163 pw%array = cmplx(array, 0.0_dp, kind=dp)
5166 CALL timestop(handle)
5167 END SUBROUTINE pw_copy_from_array_c1d_r1d_rs
5185 SUBROUTINE pw_axpy_c1d_r1d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
5187 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
5188 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw2
5189 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
5190 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
5192 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
5195 LOGICAL :: my_allow_noncompatible_grids
5196 REAL(KIND=dp) :: my_alpha, my_beta
5197 INTEGER :: i, j, ng, ng1, ng2
5199 CALL timeset(routinen, handle)
5202 IF (
PRESENT(alpha)) my_alpha = alpha
5205 IF (
PRESENT(beta)) my_beta = beta
5207 my_allow_noncompatible_grids = .false.
5208 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
5210 IF (my_beta /= 1.0_dp)
THEN
5211 IF (my_beta == 0.0_dp)
THEN
5215 pw2%array = pw2%array*my_beta
5220 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
5222 IF (my_alpha == 1.0_dp)
THEN
5224 pw2%array = pw2%array + real(pw1%array, kind=dp)
5226 ELSE IF (my_alpha /= 0.0_dp)
THEN
5228 pw2%array = pw2%array + my_alpha* real(pw1%array, kind=dp)
5232 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
5234 ng1 =
SIZE(pw1%array)
5235 ng2 =
SIZE(pw2%array)
5238 IF (pw1%pw_grid%spherical)
THEN
5239 IF (my_alpha == 1.0_dp)
THEN
5242 pw2%array(i) = pw2%array(i) + real(pw1%array(i), kind=dp)
5245 ELSE IF (my_alpha /= 0.0_dp)
THEN
5248 pw2%array(i) = pw2%array(i) + my_alpha* real(pw1%array(i), kind=dp)
5252 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
5253 IF (ng1 >= ng2)
THEN
5254 IF (my_alpha == 1.0_dp)
THEN
5257 j = pw2%pw_grid%gidx(i)
5258 pw2%array(i) = pw2%array(i) + real(pw1%array(j), kind=dp)
5261 ELSE IF (my_alpha /= 0.0_dp)
THEN
5264 j = pw2%pw_grid%gidx(i)
5265 pw2%array(i) = pw2%array(i) + my_alpha* real(pw1%array(j), kind=dp)
5270 IF (my_alpha == 1.0_dp)
THEN
5273 j = pw2%pw_grid%gidx(i)
5274 pw2%array(j) = pw2%array(j) + real(pw1%array(i), kind=dp)
5277 ELSE IF (my_alpha /= 0.0_dp)
THEN
5280 j = pw2%pw_grid%gidx(i)
5281 pw2%array(j) = pw2%array(j) + my_alpha* real(pw1%array(i), kind=dp)
5286 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
5287 IF (ng1 >= ng2)
THEN
5288 IF (my_alpha == 1.0_dp)
THEN
5291 j = pw1%pw_grid%gidx(i)
5292 pw2%array(i) = pw2%array(i) + real(pw1%array(j), kind=dp)
5295 ELSE IF (my_alpha /= 0.0_dp)
THEN
5298 j = pw1%pw_grid%gidx(i)
5299 pw2%array(i) = pw2%array(i) + my_alpha* real(pw1%array(j), kind=dp)
5304 IF (my_alpha == 1.0_dp)
THEN
5307 j = pw1%pw_grid%gidx(i)
5308 pw2%array(j) = pw2%array(j) + real(pw1%array(i), kind=dp)
5311 ELSE IF (my_alpha /= 0.0_dp)
THEN
5314 j = pw1%pw_grid%gidx(i)
5315 pw2%array(j) = pw2%array(j) + my_alpha* real(pw1%array(i), kind=dp)
5321 cpabort(
"Grids not compatible")
5326 cpabort(
"Grids not compatible")
5330 CALL timestop(handle)
5332 END SUBROUTINE pw_axpy_c1d_r1d_rs
5343 SUBROUTINE pw_multiply_c1d_r1d_rs (pw_out, pw1, pw2, alpha)
5345 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw_out
5346 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
5347 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw2
5348 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
5350 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
5353 REAL(KIND=dp) :: my_alpha
5355 CALL timeset(routinen, handle)
5358 IF (
PRESENT(alpha)) my_alpha = alpha
5360 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
5361 cpabort(
"pw_multiply not implemented for non-identical grids!")
5363#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
5364 IF (my_alpha == 1.0_dp)
THEN
5366 pw_out%array = pw_out%array + pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5370 pw_out%array = pw_out%array + my_alpha*pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5374 IF (my_alpha == 1.0_dp)
THEN
5375 pw_out%array = pw_out%array + pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5377 pw_out%array = pw_out%array + my_alpha*pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5381 CALL timestop(handle)
5383 END SUBROUTINE pw_multiply_c1d_r1d_rs
5390 SUBROUTINE pw_multiply_with_c1d_r1d_rs (pw1, pw2)
5391 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw1
5392 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw2
5394 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
5398 CALL timeset(routinen, handle)
5400 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
5401 cpabort(
"Incompatible grids!")
5404 pw1%array = pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5407 CALL timestop(handle)
5409 END SUBROUTINE pw_multiply_with_c1d_r1d_rs
5424 FUNCTION pw_integral_ab_c1d_r1d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
5426 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
5427 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw2
5428 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
5429 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
5430 REAL(kind=dp) :: integral_value
5432 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_c1d_r1d_rs'
5434 INTEGER :: handle, loc_sumtype
5435 LOGICAL :: my_just_sum, my_local_only
5437 CALL timeset(routinen, handle)
5440 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
5442 my_local_only = .false.
5443 IF (
PRESENT(local_only)) my_local_only = local_only
5445 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
5446 cpabort(
"Grids incompatible")
5449 my_just_sum = .false.
5450 IF (
PRESENT(just_sum)) my_just_sum = just_sum
5457 integral_value = sum(real(pw1%array, kind=dp)*pw2%array)
5462 integral_value = accurate_sum(real(pw1%array, kind=dp)*pw2%array)
5466 IF (.NOT. my_just_sum)
THEN
5467 integral_value = integral_value*pw1%pw_grid%dvol
5470 IF (pw1%pw_grid%grid_span == halfspace)
THEN
5471 integral_value = 2.0_dp*integral_value
5472 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
5473 REAL(conjg(pw1%array(1))*pw2%array(1), kind=dp)
5476 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
5477 CALL pw1%pw_grid%para%group%sum(integral_value)
5479 CALL timestop(handle)
5481 END FUNCTION pw_integral_ab_c1d_r1d_rs
5495 SUBROUTINE pw_copy_c1d_r1d_gs (pw1, pw2)
5497 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
5498 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw2
5500 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
5503 INTEGER :: i, j, ng, ng1, ng2, ns
5505 CALL timeset(routinen, handle)
5507 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
5508 cpabort(
"Both grids must be either spherical or non-spherical!")
5510 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
5511 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
5512 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
5513 ng1 =
SIZE(pw1%array)
5514 ng2 =
SIZE(pw2%array)
5517 pw2%array(1:ng) = real(pw1%array(1:ng), kind=dp)
5521 pw2%array(ng + 1:ng2) = 0.0_dp
5525 cpabort(
"Copies between spherical grids require compatible grids!")
5528 ng1 =
SIZE(pw1%array)
5529 ng2 =
SIZE(pw2%array)
5530 ns = 2*max(ng1, ng2)
5532 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
5533 IF (ng1 >= ng2)
THEN
5538 j = pw2%pw_grid%gidx(i)
5539 pw2%array(i) = real(pw1%array(j), kind=dp)
5548 j = pw2%pw_grid%gidx(i)
5549 pw2%array(j) = real(pw1%array(i), kind=dp)
5553 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
5554 IF (ng1 >= ng2)
THEN
5557 j = pw1%pw_grid%gidx(i)
5558 pw2%array(i) = real(pw1%array(j), kind=dp)
5565 j = pw1%pw_grid%gidx(i)
5566 pw2%array(j) = real(pw1%array(i), kind=dp)
5571 cpabort(
"Copy not implemented!")
5578 pw2%array = real(pw1%array, kind=dp)
5582 CALL timestop(handle)
5584 END SUBROUTINE pw_copy_c1d_r1d_gs
5591 SUBROUTINE pw_copy_to_array_c1d_r1d_gs (pw, array)
5592 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
5593 REAL(KIND=dp),
DIMENSION(:),
INTENT(INOUT) :: array
5595 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
5599 CALL timeset(routinen, handle)
5602 array(:) = real(pw%array(:), kind=dp)
5605 CALL timestop(handle)
5606 END SUBROUTINE pw_copy_to_array_c1d_r1d_gs
5613 SUBROUTINE pw_copy_from_array_c1d_r1d_gs (pw, array)
5614 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
5615 REAL(KIND=dp),
DIMENSION(:),
INTENT(IN) :: array
5617 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
5621 CALL timeset(routinen, handle)
5624 pw%array = cmplx(array, 0.0_dp, kind=dp)
5627 CALL timestop(handle)
5628 END SUBROUTINE pw_copy_from_array_c1d_r1d_gs
5646 SUBROUTINE pw_axpy_c1d_r1d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
5648 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
5649 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw2
5650 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
5651 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
5653 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
5656 LOGICAL :: my_allow_noncompatible_grids
5657 REAL(KIND=dp) :: my_alpha, my_beta
5658 INTEGER :: i, j, ng, ng1, ng2
5660 CALL timeset(routinen, handle)
5663 IF (
PRESENT(alpha)) my_alpha = alpha
5666 IF (
PRESENT(beta)) my_beta = beta
5668 my_allow_noncompatible_grids = .false.
5669 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
5671 IF (my_beta /= 1.0_dp)
THEN
5672 IF (my_beta == 0.0_dp)
THEN
5676 pw2%array = pw2%array*my_beta
5681 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
5683 IF (my_alpha == 1.0_dp)
THEN
5685 pw2%array = pw2%array + real(pw1%array, kind=dp)
5687 ELSE IF (my_alpha /= 0.0_dp)
THEN
5689 pw2%array = pw2%array + my_alpha* real(pw1%array, kind=dp)
5693 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
5695 ng1 =
SIZE(pw1%array)
5696 ng2 =
SIZE(pw2%array)
5699 IF (pw1%pw_grid%spherical)
THEN
5700 IF (my_alpha == 1.0_dp)
THEN
5703 pw2%array(i) = pw2%array(i) + real(pw1%array(i), kind=dp)
5706 ELSE IF (my_alpha /= 0.0_dp)
THEN
5709 pw2%array(i) = pw2%array(i) + my_alpha* real(pw1%array(i), kind=dp)
5713 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
5714 IF (ng1 >= ng2)
THEN
5715 IF (my_alpha == 1.0_dp)
THEN
5718 j = pw2%pw_grid%gidx(i)
5719 pw2%array(i) = pw2%array(i) + real(pw1%array(j), kind=dp)
5722 ELSE IF (my_alpha /= 0.0_dp)
THEN
5725 j = pw2%pw_grid%gidx(i)
5726 pw2%array(i) = pw2%array(i) + my_alpha* real(pw1%array(j), kind=dp)
5731 IF (my_alpha == 1.0_dp)
THEN
5734 j = pw2%pw_grid%gidx(i)
5735 pw2%array(j) = pw2%array(j) + real(pw1%array(i), kind=dp)
5738 ELSE IF (my_alpha /= 0.0_dp)
THEN
5741 j = pw2%pw_grid%gidx(i)
5742 pw2%array(j) = pw2%array(j) + my_alpha* real(pw1%array(i), kind=dp)
5747 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
5748 IF (ng1 >= ng2)
THEN
5749 IF (my_alpha == 1.0_dp)
THEN
5752 j = pw1%pw_grid%gidx(i)
5753 pw2%array(i) = pw2%array(i) + real(pw1%array(j), kind=dp)
5756 ELSE IF (my_alpha /= 0.0_dp)
THEN
5759 j = pw1%pw_grid%gidx(i)
5760 pw2%array(i) = pw2%array(i) + my_alpha* real(pw1%array(j), kind=dp)
5765 IF (my_alpha == 1.0_dp)
THEN
5768 j = pw1%pw_grid%gidx(i)
5769 pw2%array(j) = pw2%array(j) + real(pw1%array(i), kind=dp)
5772 ELSE IF (my_alpha /= 0.0_dp)
THEN
5775 j = pw1%pw_grid%gidx(i)
5776 pw2%array(j) = pw2%array(j) + my_alpha* real(pw1%array(i), kind=dp)
5782 cpabort(
"Grids not compatible")
5787 cpabort(
"Grids not compatible")
5791 CALL timestop(handle)
5793 END SUBROUTINE pw_axpy_c1d_r1d_gs
5804 SUBROUTINE pw_multiply_c1d_r1d_gs (pw_out, pw1, pw2, alpha)
5806 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw_out
5807 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
5808 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
5809 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
5811 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
5814 REAL(KIND=dp) :: my_alpha
5816 CALL timeset(routinen, handle)
5819 IF (
PRESENT(alpha)) my_alpha = alpha
5821 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
5822 cpabort(
"pw_multiply not implemented for non-identical grids!")
5824#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
5825 IF (my_alpha == 1.0_dp)
THEN
5827 pw_out%array = pw_out%array + pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5831 pw_out%array = pw_out%array + my_alpha*pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5835 IF (my_alpha == 1.0_dp)
THEN
5836 pw_out%array = pw_out%array + pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5838 pw_out%array = pw_out%array + my_alpha*pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5842 CALL timestop(handle)
5844 END SUBROUTINE pw_multiply_c1d_r1d_gs
5851 SUBROUTINE pw_multiply_with_c1d_r1d_gs (pw1, pw2)
5852 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw1
5853 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
5855 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
5859 CALL timeset(routinen, handle)
5861 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
5862 cpabort(
"Incompatible grids!")
5865 pw1%array = pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5868 CALL timestop(handle)
5870 END SUBROUTINE pw_multiply_with_c1d_r1d_gs
5885 FUNCTION pw_integral_ab_c1d_r1d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
5887 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
5888 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
5889 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
5890 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
5891 REAL(kind=dp) :: integral_value
5893 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_c1d_r1d_gs'
5895 INTEGER :: handle, loc_sumtype
5896 LOGICAL :: my_just_sum, my_local_only
5898 CALL timeset(routinen, handle)
5901 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
5903 my_local_only = .false.
5904 IF (
PRESENT(local_only)) my_local_only = local_only
5906 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
5907 cpabort(
"Grids incompatible")
5910 my_just_sum = .false.
5911 IF (
PRESENT(just_sum)) my_just_sum = just_sum
5918 integral_value = sum(real(pw1%array, kind=dp)*pw2%array)
5923 integral_value = accurate_sum(real(pw1%array, kind=dp)*pw2%array)
5927 IF (.NOT. my_just_sum)
THEN
5928 integral_value = integral_value*pw1%pw_grid%vol
5931 IF (pw1%pw_grid%grid_span == halfspace)
THEN
5932 integral_value = 2.0_dp*integral_value
5933 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
5934 REAL(conjg(pw1%array(1))*pw2%array(1), kind=dp)
5937 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
5938 CALL pw1%pw_grid%para%group%sum(integral_value)
5940 CALL timestop(handle)
5942 END FUNCTION pw_integral_ab_c1d_r1d_gs
5950 FUNCTION pw_integral_a2b_c1d_r1d (pw1, pw2)
RESULT(integral_value)
5952 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
5953 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
5954 REAL(kind=dp) :: integral_value
5956 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_a2b'
5960 CALL timeset(routinen, handle)
5962 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
5963 cpabort(
"Grids incompatible")
5966 integral_value = accurate_sum(real(conjg(pw1%array), kind=dp)*pw2%array*pw1%pw_grid%gsq)
5967 IF (pw1%pw_grid%grid_span == halfspace)
THEN
5968 integral_value = 2.0_dp*integral_value
5971 integral_value = integral_value*pw1%pw_grid%vol
5973 IF (pw1%pw_grid%para%mode == pw_mode_distributed) &
5974 CALL pw1%pw_grid%para%group%sum(integral_value)
5975 CALL timestop(handle)
5977 END FUNCTION pw_integral_a2b_c1d_r1d
5991 SUBROUTINE pw_copy_c1d_c1d_rs (pw1, pw2)
5993 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
5994 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw2
5996 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
5999 INTEGER :: i, j, ng, ng1, ng2, ns
6001 CALL timeset(routinen, handle)
6003 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
6004 cpabort(
"Both grids must be either spherical or non-spherical!")
6005 IF (pw1%pw_grid%spherical) &
6006 cpabort(
"Spherical grids only exist in reciprocal space!")
6008 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
6009 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
6010 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
6011 ng1 =
SIZE(pw1%array)
6012 ng2 =
SIZE(pw2%array)
6015 pw2%array(1:ng) = pw1%array(1:ng)
6019 pw2%array(ng + 1:ng2) = cmplx(0.0_dp, 0.0_dp, kind=dp)
6023 cpabort(
"Copies between spherical grids require compatible grids!")
6026 ng1 =
SIZE(pw1%array)
6027 ng2 =
SIZE(pw2%array)
6028 ns = 2*max(ng1, ng2)
6030 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
6031 IF (ng1 >= ng2)
THEN
6036 j = pw2%pw_grid%gidx(i)
6037 pw2%array(i) = pw1%array(j)
6046 j = pw2%pw_grid%gidx(i)
6047 pw2%array(j) = pw1%array(i)
6051 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
6052 IF (ng1 >= ng2)
THEN
6055 j = pw1%pw_grid%gidx(i)
6056 pw2%array(i) = pw1%array(j)
6063 j = pw1%pw_grid%gidx(i)
6064 pw2%array(j) = pw1%array(i)
6069 cpabort(
"Copy not implemented!")
6076 pw2%array = pw1%array
6080 CALL timestop(handle)
6082 END SUBROUTINE pw_copy_c1d_c1d_rs
6089 SUBROUTINE pw_copy_to_array_c1d_c1d_rs (pw, array)
6090 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw
6091 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(INOUT) :: array
6093 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
6097 CALL timeset(routinen, handle)
6100 array(:) = pw%array(:)
6103 CALL timestop(handle)
6104 END SUBROUTINE pw_copy_to_array_c1d_c1d_rs
6111 SUBROUTINE pw_copy_from_array_c1d_c1d_rs (pw, array)
6112 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw
6113 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(IN) :: array
6115 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
6119 CALL timeset(routinen, handle)
6125 CALL timestop(handle)
6126 END SUBROUTINE pw_copy_from_array_c1d_c1d_rs
6144 SUBROUTINE pw_axpy_c1d_c1d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
6146 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
6147 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw2
6148 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
6149 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
6151 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
6154 LOGICAL :: my_allow_noncompatible_grids
6155 REAL(KIND=dp) :: my_alpha, my_beta
6156 INTEGER :: i, j, ng, ng1, ng2
6158 CALL timeset(routinen, handle)
6161 IF (
PRESENT(alpha)) my_alpha = alpha
6164 IF (
PRESENT(beta)) my_beta = beta
6166 my_allow_noncompatible_grids = .false.
6167 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
6169 IF (my_beta /= 1.0_dp)
THEN
6170 IF (my_beta == 0.0_dp)
THEN
6174 pw2%array = pw2%array*my_beta
6179 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
6181 IF (my_alpha == 1.0_dp)
THEN
6183 pw2%array = pw2%array + pw1%array
6185 ELSE IF (my_alpha /= 0.0_dp)
THEN
6187 pw2%array = pw2%array + my_alpha* pw1%array
6191 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
6193 ng1 =
SIZE(pw1%array)
6194 ng2 =
SIZE(pw2%array)
6197 IF (pw1%pw_grid%spherical)
THEN
6198 IF (my_alpha == 1.0_dp)
THEN
6201 pw2%array(i) = pw2%array(i) + pw1%array(i)
6204 ELSE IF (my_alpha /= 0.0_dp)
THEN
6207 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(i)
6211 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
6212 IF (ng1 >= ng2)
THEN
6213 IF (my_alpha == 1.0_dp)
THEN
6216 j = pw2%pw_grid%gidx(i)
6217 pw2%array(i) = pw2%array(i) + pw1%array(j)
6220 ELSE IF (my_alpha /= 0.0_dp)
THEN
6223 j = pw2%pw_grid%gidx(i)
6224 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
6229 IF (my_alpha == 1.0_dp)
THEN
6232 j = pw2%pw_grid%gidx(i)
6233 pw2%array(j) = pw2%array(j) + pw1%array(i)
6236 ELSE IF (my_alpha /= 0.0_dp)
THEN
6239 j = pw2%pw_grid%gidx(i)
6240 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
6245 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
6246 IF (ng1 >= ng2)
THEN
6247 IF (my_alpha == 1.0_dp)
THEN
6250 j = pw1%pw_grid%gidx(i)
6251 pw2%array(i) = pw2%array(i) + pw1%array(j)
6254 ELSE IF (my_alpha /= 0.0_dp)
THEN
6257 j = pw1%pw_grid%gidx(i)
6258 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
6263 IF (my_alpha == 1.0_dp)
THEN
6266 j = pw1%pw_grid%gidx(i)
6267 pw2%array(j) = pw2%array(j) + pw1%array(i)
6270 ELSE IF (my_alpha /= 0.0_dp)
THEN
6273 j = pw1%pw_grid%gidx(i)
6274 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
6280 cpabort(
"Grids not compatible")
6285 cpabort(
"Grids not compatible")
6289 CALL timestop(handle)
6291 END SUBROUTINE pw_axpy_c1d_c1d_rs
6302 SUBROUTINE pw_multiply_c1d_c1d_rs (pw_out, pw1, pw2, alpha)
6304 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw_out
6305 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
6306 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw2
6307 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
6309 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
6312 REAL(KIND=dp) :: my_alpha
6314 CALL timeset(routinen, handle)
6317 IF (
PRESENT(alpha)) my_alpha = alpha
6319 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
6320 cpabort(
"pw_multiply not implemented for non-identical grids!")
6322#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
6323 IF (my_alpha == 1.0_dp)
THEN
6325 pw_out%array = pw_out%array + pw1%array* pw2%array
6329 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
6333 IF (my_alpha == 1.0_dp)
THEN
6334 pw_out%array = pw_out%array + pw1%array* pw2%array
6336 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
6340 CALL timestop(handle)
6342 END SUBROUTINE pw_multiply_c1d_c1d_rs
6349 SUBROUTINE pw_multiply_with_c1d_c1d_rs (pw1, pw2)
6350 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw1
6351 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw2
6353 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
6357 CALL timeset(routinen, handle)
6359 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
6360 cpabort(
"Incompatible grids!")
6363 pw1%array = pw1%array* pw2%array
6366 CALL timestop(handle)
6368 END SUBROUTINE pw_multiply_with_c1d_c1d_rs
6383 FUNCTION pw_integral_ab_c1d_c1d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
6385 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
6386 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw2
6387 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
6388 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
6389 REAL(kind=dp) :: integral_value
6391 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_c1d_c1d_rs'
6393 INTEGER :: handle, loc_sumtype
6394 LOGICAL :: my_just_sum, my_local_only
6396 CALL timeset(routinen, handle)
6399 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
6401 my_local_only = .false.
6402 IF (
PRESENT(local_only)) my_local_only = local_only
6404 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
6405 cpabort(
"Grids incompatible")
6408 my_just_sum = .false.
6409 IF (
PRESENT(just_sum)) my_just_sum = just_sum
6416 integral_value = sum(real(conjg(pw1%array) &
6417 *pw2%array, kind=dp))
6422 integral_value = accurate_sum(real(conjg(pw1%array)*pw2%array, kind=dp))
6426 IF (.NOT. my_just_sum)
THEN
6427 integral_value = integral_value*pw1%pw_grid%dvol
6430 IF (pw1%pw_grid%grid_span == halfspace)
THEN
6431 integral_value = 2.0_dp*integral_value
6432 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
6433 REAL(conjg(pw1%array(1))*pw2%array(1), kind=dp)
6436 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
6437 CALL pw1%pw_grid%para%group%sum(integral_value)
6439 CALL timestop(handle)
6441 END FUNCTION pw_integral_ab_c1d_c1d_rs
6455 SUBROUTINE pw_copy_c1d_c1d_gs (pw1, pw2)
6457 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
6458 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
6460 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
6463 INTEGER :: i, j, ng, ng1, ng2, ns
6465 CALL timeset(routinen, handle)
6467 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
6468 cpabort(
"Both grids must be either spherical or non-spherical!")
6470 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
6471 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
6472 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
6473 ng1 =
SIZE(pw1%array)
6474 ng2 =
SIZE(pw2%array)
6477 pw2%array(1:ng) = pw1%array(1:ng)
6481 pw2%array(ng + 1:ng2) = cmplx(0.0_dp, 0.0_dp, kind=dp)
6485 cpabort(
"Copies between spherical grids require compatible grids!")
6488 ng1 =
SIZE(pw1%array)
6489 ng2 =
SIZE(pw2%array)
6490 ns = 2*max(ng1, ng2)
6492 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
6493 IF (ng1 >= ng2)
THEN
6498 j = pw2%pw_grid%gidx(i)
6499 pw2%array(i) = pw1%array(j)
6508 j = pw2%pw_grid%gidx(i)
6509 pw2%array(j) = pw1%array(i)
6513 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
6514 IF (ng1 >= ng2)
THEN
6517 j = pw1%pw_grid%gidx(i)
6518 pw2%array(i) = pw1%array(j)
6525 j = pw1%pw_grid%gidx(i)
6526 pw2%array(j) = pw1%array(i)
6531 CALL pw_copy_match(pw1, pw2)
6538 pw2%array = pw1%array
6542 CALL timestop(handle)
6544 END SUBROUTINE pw_copy_c1d_c1d_gs
6551 SUBROUTINE pw_copy_to_array_c1d_c1d_gs (pw, array)
6552 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
6553 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(INOUT) :: array
6555 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
6559 CALL timeset(routinen, handle)
6562 array(:) = pw%array(:)
6565 CALL timestop(handle)
6566 END SUBROUTINE pw_copy_to_array_c1d_c1d_gs
6573 SUBROUTINE pw_copy_from_array_c1d_c1d_gs (pw, array)
6574 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
6575 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(IN) :: array
6577 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
6581 CALL timeset(routinen, handle)
6587 CALL timestop(handle)
6588 END SUBROUTINE pw_copy_from_array_c1d_c1d_gs
6606 SUBROUTINE pw_axpy_c1d_c1d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
6608 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
6609 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
6610 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
6611 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
6613 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
6616 LOGICAL :: my_allow_noncompatible_grids
6617 REAL(KIND=dp) :: my_alpha, my_beta
6618 INTEGER :: i, j, ng, ng1, ng2
6620 CALL timeset(routinen, handle)
6623 IF (
PRESENT(alpha)) my_alpha = alpha
6626 IF (
PRESENT(beta)) my_beta = beta
6628 my_allow_noncompatible_grids = .false.
6629 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
6631 IF (my_beta /= 1.0_dp)
THEN
6632 IF (my_beta == 0.0_dp)
THEN
6636 pw2%array = pw2%array*my_beta
6641 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
6643 IF (my_alpha == 1.0_dp)
THEN
6645 pw2%array = pw2%array + pw1%array
6647 ELSE IF (my_alpha /= 0.0_dp)
THEN
6649 pw2%array = pw2%array + my_alpha* pw1%array
6653 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
6655 ng1 =
SIZE(pw1%array)
6656 ng2 =
SIZE(pw2%array)
6659 IF (pw1%pw_grid%spherical)
THEN
6660 IF (my_alpha == 1.0_dp)
THEN
6663 pw2%array(i) = pw2%array(i) + pw1%array(i)
6666 ELSE IF (my_alpha /= 0.0_dp)
THEN
6669 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(i)
6673 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
6674 IF (ng1 >= ng2)
THEN
6675 IF (my_alpha == 1.0_dp)
THEN
6678 j = pw2%pw_grid%gidx(i)
6679 pw2%array(i) = pw2%array(i) + pw1%array(j)
6682 ELSE IF (my_alpha /= 0.0_dp)
THEN
6685 j = pw2%pw_grid%gidx(i)
6686 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
6691 IF (my_alpha == 1.0_dp)
THEN
6694 j = pw2%pw_grid%gidx(i)
6695 pw2%array(j) = pw2%array(j) + pw1%array(i)
6698 ELSE IF (my_alpha /= 0.0_dp)
THEN
6701 j = pw2%pw_grid%gidx(i)
6702 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
6707 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
6708 IF (ng1 >= ng2)
THEN
6709 IF (my_alpha == 1.0_dp)
THEN
6712 j = pw1%pw_grid%gidx(i)
6713 pw2%array(i) = pw2%array(i) + pw1%array(j)
6716 ELSE IF (my_alpha /= 0.0_dp)
THEN
6719 j = pw1%pw_grid%gidx(i)
6720 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
6725 IF (my_alpha == 1.0_dp)
THEN
6728 j = pw1%pw_grid%gidx(i)
6729 pw2%array(j) = pw2%array(j) + pw1%array(i)
6732 ELSE IF (my_alpha /= 0.0_dp)
THEN
6735 j = pw1%pw_grid%gidx(i)
6736 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
6742 cpabort(
"Grids not compatible")
6747 cpabort(
"Grids not compatible")
6751 CALL timestop(handle)
6753 END SUBROUTINE pw_axpy_c1d_c1d_gs
6764 SUBROUTINE pw_multiply_c1d_c1d_gs (pw_out, pw1, pw2, alpha)
6766 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw_out
6767 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
6768 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
6769 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
6771 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
6774 REAL(KIND=dp) :: my_alpha
6776 CALL timeset(routinen, handle)
6779 IF (
PRESENT(alpha)) my_alpha = alpha
6781 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
6782 cpabort(
"pw_multiply not implemented for non-identical grids!")
6784#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
6785 IF (my_alpha == 1.0_dp)
THEN
6787 pw_out%array = pw_out%array + pw1%array* pw2%array
6791 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
6795 IF (my_alpha == 1.0_dp)
THEN
6796 pw_out%array = pw_out%array + pw1%array* pw2%array
6798 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
6802 CALL timestop(handle)
6804 END SUBROUTINE pw_multiply_c1d_c1d_gs
6811 SUBROUTINE pw_multiply_with_c1d_c1d_gs (pw1, pw2)
6812 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw1
6813 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
6815 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
6819 CALL timeset(routinen, handle)
6821 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
6822 cpabort(
"Incompatible grids!")
6825 pw1%array = pw1%array* pw2%array
6828 CALL timestop(handle)
6830 END SUBROUTINE pw_multiply_with_c1d_c1d_gs
6845 FUNCTION pw_integral_ab_c1d_c1d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
6847 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
6848 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
6849 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
6850 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
6851 REAL(kind=dp) :: integral_value
6853 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_c1d_c1d_gs'
6855 INTEGER :: handle, loc_sumtype
6856 LOGICAL :: my_just_sum, my_local_only
6858 CALL timeset(routinen, handle)
6861 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
6863 my_local_only = .false.
6864 IF (
PRESENT(local_only)) my_local_only = local_only
6866 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
6867 cpabort(
"Grids incompatible")
6870 my_just_sum = .false.
6871 IF (
PRESENT(just_sum)) my_just_sum = just_sum
6878 integral_value = sum(real(conjg(pw1%array) &
6879 *pw2%array, kind=dp))
6884 integral_value = accurate_sum(real(conjg(pw1%array)*pw2%array, kind=dp))
6888 IF (.NOT. my_just_sum)
THEN
6889 integral_value = integral_value*pw1%pw_grid%vol
6892 IF (pw1%pw_grid%grid_span == halfspace)
THEN
6893 integral_value = 2.0_dp*integral_value
6894 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
6895 REAL(conjg(pw1%array(1))*pw2%array(1), kind=dp)
6898 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
6899 CALL pw1%pw_grid%para%group%sum(integral_value)
6901 CALL timestop(handle)
6903 END FUNCTION pw_integral_ab_c1d_c1d_gs
6911 FUNCTION pw_integral_a2b_c1d_c1d (pw1, pw2)
RESULT(integral_value)
6913 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
6914 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
6915 REAL(kind=dp) :: integral_value
6917 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_a2b'
6921 CALL timeset(routinen, handle)
6923 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
6924 cpabort(
"Grids incompatible")
6927 integral_value = accurate_sum(real(conjg(pw1%array)*pw2%array, kind=dp)*pw1%pw_grid%gsq)
6928 IF (pw1%pw_grid%grid_span == halfspace)
THEN
6929 integral_value = 2.0_dp*integral_value
6932 integral_value = integral_value*pw1%pw_grid%vol
6934 IF (pw1%pw_grid%para%mode == pw_mode_distributed) &
6935 CALL pw1%pw_grid%para%group%sum(integral_value)
6936 CALL timestop(handle)
6938 END FUNCTION pw_integral_a2b_c1d_c1d
6952 SUBROUTINE pw_copy_c3d_r3d_rs (pw1, pw2)
6954 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
6955 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw2
6957 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
6961 CALL timeset(routinen, handle)
6963 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
6964 cpabort(
"Both grids must be either spherical or non-spherical!")
6965 IF (pw1%pw_grid%spherical) &
6966 cpabort(
"Spherical grids only exist in reciprocal space!")
6968 IF (any(shape(pw2%array) /= shape(pw1%array))) &
6969 cpabort(
"3D grids must be compatible!")
6970 IF (pw1%pw_grid%spherical) &
6971 cpabort(
"3D grids must not be spherical!")
6973 pw2%array(:, :, :) = real(pw1%array(:, :, :), kind=dp)
6976 CALL timestop(handle)
6978 END SUBROUTINE pw_copy_c3d_r3d_rs
6985 SUBROUTINE pw_copy_to_array_c3d_r3d_rs (pw, array)
6986 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw
6987 REAL(KIND=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
6989 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
6993 CALL timeset(routinen, handle)
6996 array(:, :, :) = real(pw%array(:, :, :), kind=dp)
6999 CALL timestop(handle)
7000 END SUBROUTINE pw_copy_to_array_c3d_r3d_rs
7007 SUBROUTINE pw_copy_from_array_c3d_r3d_rs (pw, array)
7008 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw
7009 REAL(KIND=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
7011 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
7015 CALL timeset(routinen, handle)
7018 pw%array = cmplx(array, 0.0_dp, kind=dp)
7021 CALL timestop(handle)
7022 END SUBROUTINE pw_copy_from_array_c3d_r3d_rs
7040 SUBROUTINE pw_axpy_c3d_r3d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
7042 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
7043 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw2
7044 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
7045 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
7047 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
7050 LOGICAL :: my_allow_noncompatible_grids
7051 REAL(KIND=dp) :: my_alpha, my_beta
7053 CALL timeset(routinen, handle)
7056 IF (
PRESENT(alpha)) my_alpha = alpha
7059 IF (
PRESENT(beta)) my_beta = beta
7061 my_allow_noncompatible_grids = .false.
7062 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
7064 IF (my_beta /= 1.0_dp)
THEN
7065 IF (my_beta == 0.0_dp)
THEN
7069 pw2%array = pw2%array*my_beta
7074 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
7075 IF (my_alpha == 1.0_dp)
THEN
7077 pw2%array = pw2%array + real(pw1%array, kind=dp)
7079 ELSE IF (my_alpha /= 0.0_dp)
THEN
7081 pw2%array = pw2%array + my_alpha* real(pw1%array, kind=dp)
7085 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
7087 IF (any(shape(pw1%array) /= shape(pw2%array))) &
7088 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
7090 IF (my_alpha == 1.0_dp)
THEN
7092 pw2%array = pw2%array + real(pw1%array, kind=dp)
7096 pw2%array = pw2%array + my_alpha* real(pw1%array, kind=dp)
7102 cpabort(
"Grids not compatible")
7106 CALL timestop(handle)
7108 END SUBROUTINE pw_axpy_c3d_r3d_rs
7119 SUBROUTINE pw_multiply_c3d_r3d_rs (pw_out, pw1, pw2, alpha)
7121 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw_out
7122 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
7123 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw2
7124 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
7126 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
7129 REAL(KIND=dp) :: my_alpha
7131 CALL timeset(routinen, handle)
7134 IF (
PRESENT(alpha)) my_alpha = alpha
7136 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
7137 cpabort(
"pw_multiply not implemented for non-identical grids!")
7139#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
7140 IF (my_alpha == 1.0_dp)
THEN
7142 pw_out%array = pw_out%array + pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7146 pw_out%array = pw_out%array + my_alpha*pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7150 IF (my_alpha == 1.0_dp)
THEN
7151 pw_out%array = pw_out%array + pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7153 pw_out%array = pw_out%array + my_alpha*pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7157 CALL timestop(handle)
7159 END SUBROUTINE pw_multiply_c3d_r3d_rs
7166 SUBROUTINE pw_multiply_with_c3d_r3d_rs (pw1, pw2)
7167 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw1
7168 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw2
7170 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
7174 CALL timeset(routinen, handle)
7176 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
7177 cpabort(
"Incompatible grids!")
7180 pw1%array = pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7183 CALL timestop(handle)
7185 END SUBROUTINE pw_multiply_with_c3d_r3d_rs
7200 FUNCTION pw_integral_ab_c3d_r3d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
7202 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
7203 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw2
7204 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
7205 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
7206 REAL(kind=dp) :: integral_value
7208 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_c3d_r3d_rs'
7210 INTEGER :: handle, loc_sumtype
7211 LOGICAL :: my_just_sum, my_local_only
7213 CALL timeset(routinen, handle)
7216 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
7218 my_local_only = .false.
7219 IF (
PRESENT(local_only)) my_local_only = local_only
7221 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
7222 cpabort(
"Grids incompatible")
7225 my_just_sum = .false.
7226 IF (
PRESENT(just_sum)) my_just_sum = just_sum
7233 integral_value = sum(real(pw1%array, kind=dp)*pw2%array)
7238 integral_value = accurate_sum(real(pw1%array, kind=dp)*pw2%array)
7242 IF (.NOT. my_just_sum)
THEN
7243 integral_value = integral_value*pw1%pw_grid%dvol
7247 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
7248 CALL pw1%pw_grid%para%group%sum(integral_value)
7250 CALL timestop(handle)
7252 END FUNCTION pw_integral_ab_c3d_r3d_rs
7266 SUBROUTINE pw_copy_c3d_r3d_gs (pw1, pw2)
7268 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
7269 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw2
7271 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
7275 CALL timeset(routinen, handle)
7277 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
7278 cpabort(
"Both grids must be either spherical or non-spherical!")
7280 IF (any(shape(pw2%array) /= shape(pw1%array))) &
7281 cpabort(
"3D grids must be compatible!")
7282 IF (pw1%pw_grid%spherical) &
7283 cpabort(
"3D grids must not be spherical!")
7285 pw2%array(:, :, :) = real(pw1%array(:, :, :), kind=dp)
7288 CALL timestop(handle)
7290 END SUBROUTINE pw_copy_c3d_r3d_gs
7297 SUBROUTINE pw_copy_to_array_c3d_r3d_gs (pw, array)
7298 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw
7299 REAL(KIND=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
7301 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
7305 CALL timeset(routinen, handle)
7308 array(:, :, :) = real(pw%array(:, :, :), kind=dp)
7311 CALL timestop(handle)
7312 END SUBROUTINE pw_copy_to_array_c3d_r3d_gs
7319 SUBROUTINE pw_copy_from_array_c3d_r3d_gs (pw, array)
7320 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw
7321 REAL(KIND=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
7323 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
7327 CALL timeset(routinen, handle)
7330 pw%array = cmplx(array, 0.0_dp, kind=dp)
7333 CALL timestop(handle)
7334 END SUBROUTINE pw_copy_from_array_c3d_r3d_gs
7352 SUBROUTINE pw_axpy_c3d_r3d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
7354 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
7355 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw2
7356 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
7357 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
7359 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
7362 LOGICAL :: my_allow_noncompatible_grids
7363 REAL(KIND=dp) :: my_alpha, my_beta
7365 CALL timeset(routinen, handle)
7368 IF (
PRESENT(alpha)) my_alpha = alpha
7371 IF (
PRESENT(beta)) my_beta = beta
7373 my_allow_noncompatible_grids = .false.
7374 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
7376 IF (my_beta /= 1.0_dp)
THEN
7377 IF (my_beta == 0.0_dp)
THEN
7381 pw2%array = pw2%array*my_beta
7386 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
7387 IF (my_alpha == 1.0_dp)
THEN
7389 pw2%array = pw2%array + real(pw1%array, kind=dp)
7391 ELSE IF (my_alpha /= 0.0_dp)
THEN
7393 pw2%array = pw2%array + my_alpha* real(pw1%array, kind=dp)
7397 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
7399 IF (any(shape(pw1%array) /= shape(pw2%array))) &
7400 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
7402 IF (my_alpha == 1.0_dp)
THEN
7404 pw2%array = pw2%array + real(pw1%array, kind=dp)
7408 pw2%array = pw2%array + my_alpha* real(pw1%array, kind=dp)
7414 cpabort(
"Grids not compatible")
7418 CALL timestop(handle)
7420 END SUBROUTINE pw_axpy_c3d_r3d_gs
7431 SUBROUTINE pw_multiply_c3d_r3d_gs (pw_out, pw1, pw2, alpha)
7433 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw_out
7434 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
7435 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw2
7436 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
7438 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
7441 REAL(KIND=dp) :: my_alpha
7443 CALL timeset(routinen, handle)
7446 IF (
PRESENT(alpha)) my_alpha = alpha
7448 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
7449 cpabort(
"pw_multiply not implemented for non-identical grids!")
7451#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
7452 IF (my_alpha == 1.0_dp)
THEN
7454 pw_out%array = pw_out%array + pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7458 pw_out%array = pw_out%array + my_alpha*pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7462 IF (my_alpha == 1.0_dp)
THEN
7463 pw_out%array = pw_out%array + pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7465 pw_out%array = pw_out%array + my_alpha*pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7469 CALL timestop(handle)
7471 END SUBROUTINE pw_multiply_c3d_r3d_gs
7478 SUBROUTINE pw_multiply_with_c3d_r3d_gs (pw1, pw2)
7479 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw1
7480 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw2
7482 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
7486 CALL timeset(routinen, handle)
7488 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
7489 cpabort(
"Incompatible grids!")
7492 pw1%array = pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7495 CALL timestop(handle)
7497 END SUBROUTINE pw_multiply_with_c3d_r3d_gs
7512 FUNCTION pw_integral_ab_c3d_r3d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
7514 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
7515 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw2
7516 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
7517 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
7518 REAL(kind=dp) :: integral_value
7520 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_c3d_r3d_gs'
7522 INTEGER :: handle, loc_sumtype
7523 LOGICAL :: my_just_sum, my_local_only
7525 CALL timeset(routinen, handle)
7528 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
7530 my_local_only = .false.
7531 IF (
PRESENT(local_only)) my_local_only = local_only
7533 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
7534 cpabort(
"Grids incompatible")
7537 my_just_sum = .false.
7538 IF (
PRESENT(just_sum)) my_just_sum = just_sum
7545 integral_value = sum(real(pw1%array, kind=dp)*pw2%array)
7550 integral_value = accurate_sum(real(pw1%array, kind=dp)*pw2%array)
7554 IF (.NOT. my_just_sum)
THEN
7555 integral_value = integral_value*pw1%pw_grid%vol
7559 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
7560 CALL pw1%pw_grid%para%group%sum(integral_value)
7562 CALL timestop(handle)
7564 END FUNCTION pw_integral_ab_c3d_r3d_gs
7579 SUBROUTINE pw_copy_c3d_c3d_rs (pw1, pw2)
7581 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
7582 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw2
7584 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
7588 CALL timeset(routinen, handle)
7590 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
7591 cpabort(
"Both grids must be either spherical or non-spherical!")
7592 IF (pw1%pw_grid%spherical) &
7593 cpabort(
"Spherical grids only exist in reciprocal space!")
7595 IF (any(shape(pw2%array) /= shape(pw1%array))) &
7596 cpabort(
"3D grids must be compatible!")
7597 IF (pw1%pw_grid%spherical) &
7598 cpabort(
"3D grids must not be spherical!")
7600 pw2%array(:, :, :) = pw1%array(:, :, :)
7603 CALL timestop(handle)
7605 END SUBROUTINE pw_copy_c3d_c3d_rs
7612 SUBROUTINE pw_copy_to_array_c3d_c3d_rs (pw, array)
7613 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw
7614 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
7616 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
7620 CALL timeset(routinen, handle)
7623 array(:, :, :) = pw%array(:, :, :)
7626 CALL timestop(handle)
7627 END SUBROUTINE pw_copy_to_array_c3d_c3d_rs
7634 SUBROUTINE pw_copy_from_array_c3d_c3d_rs (pw, array)
7635 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw
7636 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
7638 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
7642 CALL timeset(routinen, handle)
7648 CALL timestop(handle)
7649 END SUBROUTINE pw_copy_from_array_c3d_c3d_rs
7667 SUBROUTINE pw_axpy_c3d_c3d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
7669 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
7670 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw2
7671 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
7672 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
7674 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
7677 LOGICAL :: my_allow_noncompatible_grids
7678 REAL(KIND=dp) :: my_alpha, my_beta
7680 CALL timeset(routinen, handle)
7683 IF (
PRESENT(alpha)) my_alpha = alpha
7686 IF (
PRESENT(beta)) my_beta = beta
7688 my_allow_noncompatible_grids = .false.
7689 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
7691 IF (my_beta /= 1.0_dp)
THEN
7692 IF (my_beta == 0.0_dp)
THEN
7696 pw2%array = pw2%array*my_beta
7701 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
7702 IF (my_alpha == 1.0_dp)
THEN
7704 pw2%array = pw2%array + pw1%array
7706 ELSE IF (my_alpha /= 0.0_dp)
THEN
7708 pw2%array = pw2%array + my_alpha* pw1%array
7712 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
7714 IF (any(shape(pw1%array) /= shape(pw2%array))) &
7715 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
7717 IF (my_alpha == 1.0_dp)
THEN
7719 pw2%array = pw2%array + pw1%array
7723 pw2%array = pw2%array + my_alpha* pw1%array
7729 cpabort(
"Grids not compatible")
7733 CALL timestop(handle)
7735 END SUBROUTINE pw_axpy_c3d_c3d_rs
7746 SUBROUTINE pw_multiply_c3d_c3d_rs (pw_out, pw1, pw2, alpha)
7748 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw_out
7749 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
7750 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw2
7751 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
7753 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
7756 REAL(KIND=dp) :: my_alpha
7758 CALL timeset(routinen, handle)
7761 IF (
PRESENT(alpha)) my_alpha = alpha
7763 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
7764 cpabort(
"pw_multiply not implemented for non-identical grids!")
7766#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
7767 IF (my_alpha == 1.0_dp)
THEN
7769 pw_out%array = pw_out%array + pw1%array* pw2%array
7773 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
7777 IF (my_alpha == 1.0_dp)
THEN
7778 pw_out%array = pw_out%array + pw1%array* pw2%array
7780 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
7784 CALL timestop(handle)
7786 END SUBROUTINE pw_multiply_c3d_c3d_rs
7793 SUBROUTINE pw_multiply_with_c3d_c3d_rs (pw1, pw2)
7794 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw1
7795 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw2
7797 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
7801 CALL timeset(routinen, handle)
7803 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
7804 cpabort(
"Incompatible grids!")
7807 pw1%array = pw1%array* pw2%array
7810 CALL timestop(handle)
7812 END SUBROUTINE pw_multiply_with_c3d_c3d_rs
7827 FUNCTION pw_integral_ab_c3d_c3d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
7829 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
7830 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw2
7831 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
7832 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
7833 REAL(kind=dp) :: integral_value
7835 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_c3d_c3d_rs'
7837 INTEGER :: handle, loc_sumtype
7838 LOGICAL :: my_just_sum, my_local_only
7840 CALL timeset(routinen, handle)
7843 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
7845 my_local_only = .false.
7846 IF (
PRESENT(local_only)) my_local_only = local_only
7848 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
7849 cpabort(
"Grids incompatible")
7852 my_just_sum = .false.
7853 IF (
PRESENT(just_sum)) my_just_sum = just_sum
7860 integral_value = sum(real(conjg(pw1%array) &
7861 *pw2%array, kind=dp))
7866 integral_value = accurate_sum(real(conjg(pw1%array)*pw2%array, kind=dp))
7870 IF (.NOT. my_just_sum)
THEN
7871 integral_value = integral_value*pw1%pw_grid%dvol
7875 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
7876 CALL pw1%pw_grid%para%group%sum(integral_value)
7878 CALL timestop(handle)
7880 END FUNCTION pw_integral_ab_c3d_c3d_rs
7894 SUBROUTINE pw_copy_c3d_c3d_gs (pw1, pw2)
7896 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
7897 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
7899 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
7903 CALL timeset(routinen, handle)
7905 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
7906 cpabort(
"Both grids must be either spherical or non-spherical!")
7908 IF (any(shape(pw2%array) /= shape(pw1%array))) &
7909 cpabort(
"3D grids must be compatible!")
7910 IF (pw1%pw_grid%spherical) &
7911 cpabort(
"3D grids must not be spherical!")
7913 pw2%array(:, :, :) = pw1%array(:, :, :)
7916 CALL timestop(handle)
7918 END SUBROUTINE pw_copy_c3d_c3d_gs
7925 SUBROUTINE pw_copy_to_array_c3d_c3d_gs (pw, array)
7926 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw
7927 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
7929 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
7933 CALL timeset(routinen, handle)
7936 array(:, :, :) = pw%array(:, :, :)
7939 CALL timestop(handle)
7940 END SUBROUTINE pw_copy_to_array_c3d_c3d_gs
7947 SUBROUTINE pw_copy_from_array_c3d_c3d_gs (pw, array)
7948 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw
7949 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
7951 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
7955 CALL timeset(routinen, handle)
7961 CALL timestop(handle)
7962 END SUBROUTINE pw_copy_from_array_c3d_c3d_gs
7980 SUBROUTINE pw_axpy_c3d_c3d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
7982 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
7983 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
7984 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
7985 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
7987 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
7990 LOGICAL :: my_allow_noncompatible_grids
7991 REAL(KIND=dp) :: my_alpha, my_beta
7993 CALL timeset(routinen, handle)
7996 IF (
PRESENT(alpha)) my_alpha = alpha
7999 IF (
PRESENT(beta)) my_beta = beta
8001 my_allow_noncompatible_grids = .false.
8002 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
8004 IF (my_beta /= 1.0_dp)
THEN
8005 IF (my_beta == 0.0_dp)
THEN
8009 pw2%array = pw2%array*my_beta
8014 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
8015 IF (my_alpha == 1.0_dp)
THEN
8017 pw2%array = pw2%array + pw1%array
8019 ELSE IF (my_alpha /= 0.0_dp)
THEN
8021 pw2%array = pw2%array + my_alpha* pw1%array
8025 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
8027 IF (any(shape(pw1%array) /= shape(pw2%array))) &
8028 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
8030 IF (my_alpha == 1.0_dp)
THEN
8032 pw2%array = pw2%array + pw1%array
8036 pw2%array = pw2%array + my_alpha* pw1%array
8042 cpabort(
"Grids not compatible")
8046 CALL timestop(handle)
8048 END SUBROUTINE pw_axpy_c3d_c3d_gs
8059 SUBROUTINE pw_multiply_c3d_c3d_gs (pw_out, pw1, pw2, alpha)
8061 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw_out
8062 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
8063 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw2
8064 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
8066 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
8069 REAL(KIND=dp) :: my_alpha
8071 CALL timeset(routinen, handle)
8074 IF (
PRESENT(alpha)) my_alpha = alpha
8076 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
8077 cpabort(
"pw_multiply not implemented for non-identical grids!")
8079#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
8080 IF (my_alpha == 1.0_dp)
THEN
8082 pw_out%array = pw_out%array + pw1%array* pw2%array
8086 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
8090 IF (my_alpha == 1.0_dp)
THEN
8091 pw_out%array = pw_out%array + pw1%array* pw2%array
8093 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
8097 CALL timestop(handle)
8099 END SUBROUTINE pw_multiply_c3d_c3d_gs
8106 SUBROUTINE pw_multiply_with_c3d_c3d_gs (pw1, pw2)
8107 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw1
8108 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw2
8110 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
8114 CALL timeset(routinen, handle)
8116 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
8117 cpabort(
"Incompatible grids!")
8120 pw1%array = pw1%array* pw2%array
8123 CALL timestop(handle)
8125 END SUBROUTINE pw_multiply_with_c3d_c3d_gs
8140 FUNCTION pw_integral_ab_c3d_c3d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
8142 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
8143 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw2
8144 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
8145 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
8146 REAL(kind=dp) :: integral_value
8148 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_c3d_c3d_gs'
8150 INTEGER :: handle, loc_sumtype
8151 LOGICAL :: my_just_sum, my_local_only
8153 CALL timeset(routinen, handle)
8156 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
8158 my_local_only = .false.
8159 IF (
PRESENT(local_only)) my_local_only = local_only
8161 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
8162 cpabort(
"Grids incompatible")
8165 my_just_sum = .false.
8166 IF (
PRESENT(just_sum)) my_just_sum = just_sum
8173 integral_value = sum(real(conjg(pw1%array) &
8174 *pw2%array, kind=dp))
8179 integral_value = accurate_sum(real(conjg(pw1%array)*pw2%array, kind=dp))
8183 IF (.NOT. my_just_sum)
THEN
8184 integral_value = integral_value*pw1%pw_grid%vol
8188 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
8189 CALL pw1%pw_grid%para%group%sum(integral_value)
8191 CALL timestop(handle)
8193 END FUNCTION pw_integral_ab_c3d_c3d_gs
8216 SUBROUTINE pw_gather_s_r1d_r3d_2(pw1, pw2, scale)
8218 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
8219 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw2
8220 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
8222 CALL pw_gather_s_r1d_r3d (pw2, pw1%array, scale)
8224 END SUBROUTINE pw_gather_s_r1d_r3d_2
8235 SUBROUTINE pw_gather_s_r1d_r3d (pw, c, scale)
8237 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw
8238 REAL(kind=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(IN) :: c
8239 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
8241 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_gather_s'
8243 INTEGER :: gpt, handle, l, m, n
8245 CALL timeset(routinen, handle)
8247 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
8248 ngpts =>
SIZE(pw%pw_grid%gsq), ghat => pw%pw_grid%g_hat)
8250 IF (
PRESENT(scale))
THEN
8253 l = mapl(ghat(1, gpt)) + 1
8254 m = mapm(ghat(2, gpt)) + 1
8255 n = mapn(ghat(3, gpt)) + 1
8256 pw%array(gpt) = scale* c(l, m, n)
8262 l = mapl(ghat(1, gpt)) + 1
8263 m = mapm(ghat(2, gpt)) + 1
8264 n = mapn(ghat(3, gpt)) + 1
8265 pw%array(gpt) = c(l, m, n)
8272 CALL timestop(handle)
8274 END SUBROUTINE pw_gather_s_r1d_r3d
8285 SUBROUTINE pw_scatter_s_r1d_r3d_2(pw1, pw2, scale)
8287 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
8288 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw2
8289 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
8291 CALL pw_scatter_s_r1d_r3d (pw1, pw2%array, scale)
8293 END SUBROUTINE pw_scatter_s_r1d_r3d_2
8304 SUBROUTINE pw_scatter_s_r1d_r3d (pw, c, scale)
8306 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw
8307 REAL(kind=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(INOUT) :: c
8308 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
8310 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_scatter_s'
8312 INTEGER :: gpt, handle, l, m, n
8314 CALL timeset(routinen, handle)
8316 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
8317 ghat => pw%pw_grid%g_hat, ngpts =>
SIZE(pw%pw_grid%gsq))
8320 IF (.NOT.
PRESENT(scale)) c = 0.0_dp
8322 IF (
PRESENT(scale))
THEN
8325 l = mapl(ghat(1, gpt)) + 1
8326 m = mapm(ghat(2, gpt)) + 1
8327 n = mapn(ghat(3, gpt)) + 1
8328 c(l, m, n) = scale* pw%array(gpt)
8334 l = mapl(ghat(1, gpt)) + 1
8335 m = mapm(ghat(2, gpt)) + 1
8336 n = mapn(ghat(3, gpt)) + 1
8337 c(l, m, n) = pw%array(gpt)
8344 IF (pw%pw_grid%grid_span == halfspace)
THEN
8346 associate(mapl => pw%pw_grid%mapl%neg, mapm => pw%pw_grid%mapm%neg, mapn => pw%pw_grid%mapn%neg, &
8347 ghat => pw%pw_grid%g_hat, ngpts =>
SIZE(pw%pw_grid%gsq))
8349 IF (
PRESENT(scale))
THEN
8352 l = mapl(ghat(1, gpt)) + 1
8353 m = mapm(ghat(2, gpt)) + 1
8354 n = mapn(ghat(3, gpt)) + 1
8355 c(l, m, n) = scale*( pw%array(gpt))
8361 l = mapl(ghat(1, gpt)) + 1
8362 m = mapm(ghat(2, gpt)) + 1
8363 n = mapn(ghat(3, gpt)) + 1
8364 c(l, m, n) = ( pw%array(gpt))
8373 CALL timestop(handle)
8375 END SUBROUTINE pw_scatter_s_r1d_r3d
8396 SUBROUTINE pw_gather_s_r1d_c3d_2(pw1, pw2, scale)
8398 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
8399 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw2
8400 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
8402 CALL pw_gather_s_r1d_c3d (pw2, pw1%array, scale)
8404 END SUBROUTINE pw_gather_s_r1d_c3d_2
8415 SUBROUTINE pw_gather_s_r1d_c3d (pw, c, scale)
8417 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw
8418 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(IN) :: c
8419 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
8421 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_gather_s'
8423 INTEGER :: gpt, handle, l, m, n
8425 CALL timeset(routinen, handle)
8427 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
8428 ngpts =>
SIZE(pw%pw_grid%gsq), ghat => pw%pw_grid%g_hat)
8430 IF (
PRESENT(scale))
THEN
8433 l = mapl(ghat(1, gpt)) + 1
8434 m = mapm(ghat(2, gpt)) + 1
8435 n = mapn(ghat(3, gpt)) + 1
8436 pw%array(gpt) = scale* real(c(l, m, n), kind=dp)
8442 l = mapl(ghat(1, gpt)) + 1
8443 m = mapm(ghat(2, gpt)) + 1
8444 n = mapn(ghat(3, gpt)) + 1
8445 pw%array(gpt) = real(c(l, m, n), kind=dp)
8452 CALL timestop(handle)
8454 END SUBROUTINE pw_gather_s_r1d_c3d
8465 SUBROUTINE pw_scatter_s_r1d_c3d_2(pw1, pw2, scale)
8467 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
8468 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
8469 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
8471 CALL pw_scatter_s_r1d_c3d (pw1, pw2%array, scale)
8473 END SUBROUTINE pw_scatter_s_r1d_c3d_2
8484 SUBROUTINE pw_scatter_s_r1d_c3d (pw, c, scale)
8486 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw
8487 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(INOUT) :: c
8488 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
8490 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_scatter_s'
8492 INTEGER :: gpt, handle, l, m, n
8494 CALL timeset(routinen, handle)
8496 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
8497 ghat => pw%pw_grid%g_hat, ngpts =>
SIZE(pw%pw_grid%gsq))
8500 IF (.NOT.
PRESENT(scale)) c = 0.0_dp
8502 IF (
PRESENT(scale))
THEN
8505 l = mapl(ghat(1, gpt)) + 1
8506 m = mapm(ghat(2, gpt)) + 1
8507 n = mapn(ghat(3, gpt)) + 1
8508 c(l, m, n) = scale* cmplx(pw%array(gpt), 0.0_dp, kind=dp)
8514 l = mapl(ghat(1, gpt)) + 1
8515 m = mapm(ghat(2, gpt)) + 1
8516 n = mapn(ghat(3, gpt)) + 1
8517 c(l, m, n) = cmplx(pw%array(gpt), 0.0_dp, kind=dp)
8524 IF (pw%pw_grid%grid_span == halfspace)
THEN
8526 associate(mapl => pw%pw_grid%mapl%neg, mapm => pw%pw_grid%mapm%neg, mapn => pw%pw_grid%mapn%neg, &
8527 ghat => pw%pw_grid%g_hat, ngpts =>
SIZE(pw%pw_grid%gsq))
8529 IF (
PRESENT(scale))
THEN
8532 l = mapl(ghat(1, gpt)) + 1
8533 m = mapm(ghat(2, gpt)) + 1
8534 n = mapn(ghat(3, gpt)) + 1
8535 c(l, m, n) = scale*( cmplx(pw%array(gpt), 0.0_dp, kind=dp))
8541 l = mapl(ghat(1, gpt)) + 1
8542 m = mapm(ghat(2, gpt)) + 1
8543 n = mapn(ghat(3, gpt)) + 1
8544 c(l, m, n) = ( cmplx(pw%array(gpt), 0.0_dp, kind=dp))
8553 CALL timestop(handle)
8555 END SUBROUTINE pw_scatter_s_r1d_c3d
8581 SUBROUTINE fft_wrap_pw1pw2_r3d_c1d_rs_gs (pw1, pw2, debug)
8583 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
8584 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
8585 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
8587 CHARACTER(len=*),
PARAMETER :: routineN =
'fft_wrap_pw1pw2'
8589 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
8590 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
8591 INTEGER :: handle, handle2, my_pos, nrays, &
8593 INTEGER,
DIMENSION(3) :: nloc
8594#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
8595 LOGICAL :: use_pw_gpu
8597 INTEGER,
DIMENSION(:),
POINTER :: n
8600 CALL timeset(routinen, handle2)
8601 out_unit = cp_logger_get_default_io_unit()
8602 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
8603 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
8608 IF (
PRESENT(debug))
THEN
8615 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
8616 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
8617 cpabort(
"PW grids not compatible")
8619 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
8620 cpabort(
"PW grids have not compatible MPI groups")
8624 n => pw1%pw_grid%npts
8626 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
8632 IF (test .AND. out_unit > 0)
THEN
8633 WRITE (out_unit,
'(A)')
" FFT Protocol "
8634 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
8635 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
8636 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
8639#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
8640 CALL pw_gpu_r3dc1d_3d(pw1, pw2)
8641#elif defined (__PW_FPGA)
8642 ALLOCATE (c_out(n(1), n(2), n(3)))
8645 IF (pw_fpga_init_bitstream(n) == 1)
THEN
8647#if (__PW_FPGA_SP && __PW_FPGA)
8648 CALL pw_fpga_r3dc1d_3d_sp(n, c_out)
8650 CALL pw_fpga_r3dc1d_3d_dp(n, c_out)
8652 CALL zdscal(n(1)*n(2)*n(3), 1.0_dp/pw1%pw_grid%ngpts, c_out, 1)
8653 CALL pw_gather_s_c1d_c3d(pw2, c_out)
8656 CALL fft3d(fwfft, n, c_out, debug=test)
8657 CALL pw_gather_s_c1d_c3d(pw2, c_out)
8661 ALLOCATE (c_out(n(1), n(2), n(3)))
8664 CALL fft3d(fwfft, n, c_out, debug=test)
8665 CALL pw_gather_s_c1d_c3d(pw2, c_out)
8669 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
8677 IF (test .AND. out_unit > 0)
THEN
8678 WRITE (out_unit,
'(A)')
" FFT Protocol "
8679 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
8680 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
8681 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
8684 my_pos = pw1%pw_grid%para%group%mepos
8685 nrays = pw1%pw_grid%para%nyzray(my_pos)
8686 grays => pw1%pw_grid%grays
8688#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
8690 use_pw_gpu = pw1%pw_grid%para%ray_distribution
8691 IF (use_pw_gpu)
THEN
8692 CALL pw_gpu_r3dc1d_3d_ps(pw1, pw2)
8696 nloc = pw1%pw_grid%npts_local
8697 ALLOCATE (c_in(nloc(1), nloc(2), nloc(3)))
8701 IF (pw1%pw_grid%para%ray_distribution)
THEN
8702 CALL fft3d(fwfft, n, c_in, grays, pw1%pw_grid%para%group, &
8703 pw1%pw_grid%para%yzp, pw1%pw_grid%para%nyzray, &
8704 pw1%pw_grid%para%bo, debug=test)
8706 CALL fft3d(fwfft, n, c_in, grays, pw1%pw_grid%para%group, &
8707 pw1%pw_grid%para%bo, debug=test)
8710 IF (test .AND. out_unit > 0) &
8711 WRITE (out_unit,
'(A)')
" PW_GATHER : 2d -> 1d "
8712 CALL pw_gather_p_c1d (pw2, grays)
8715#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
8720 IF (test .AND. out_unit > 0)
THEN
8721 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
8724 CALL timestop(handle)
8725 CALL timestop(handle2)
8727 END SUBROUTINE fft_wrap_pw1pw2_r3d_c1d_rs_gs
8746 SUBROUTINE fft_wrap_pw1pw2_r3d_c3d_rs_gs (pw1, pw2, debug)
8748 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
8749 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
8750 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
8752 CHARACTER(len=*),
PARAMETER :: routineN =
'fft_wrap_pw1pw2'
8754 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
8755 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
8756 INTEGER :: handle, handle2, my_pos, nrays, &
8758 INTEGER,
DIMENSION(:),
POINTER :: n
8761 CALL timeset(routinen, handle2)
8762 out_unit = cp_logger_get_default_io_unit()
8763 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
8764 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
8769 IF (
PRESENT(debug))
THEN
8776 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
8777 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
8778 cpabort(
"PW grids not compatible")
8780 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
8781 cpabort(
"PW grids have not compatible MPI groups")
8785 n => pw1%pw_grid%npts
8787 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
8793 IF (test .AND. out_unit > 0)
THEN
8794 WRITE (out_unit,
'(A)')
" FFT Protocol "
8795 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
8796 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
8797 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
8800 pw2%array = cmplx(pw1%array, 0.0_dp, kind=dp)
8802 CALL fft3d(fwfft, n, c_out, debug=test)
8804 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
8812 IF (test .AND. out_unit > 0)
THEN
8813 WRITE (out_unit,
'(A)')
" FFT Protocol "
8814 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
8815 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
8816 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
8819 my_pos = pw1%pw_grid%para%group%mepos
8820 nrays = pw1%pw_grid%para%nyzray(my_pos)
8821 grays => pw1%pw_grid%grays
8825 IF (test .AND. out_unit > 0)
THEN
8826 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
8829 CALL timestop(handle)
8830 CALL timestop(handle2)
8832 END SUBROUTINE fft_wrap_pw1pw2_r3d_c3d_rs_gs
8857 SUBROUTINE fft_wrap_pw1pw2_c1d_r3d_gs_rs (pw1, pw2, debug)
8859 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
8860 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw2
8861 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
8863 CHARACTER(len=*),
PARAMETER :: routineN =
'fft_wrap_pw1pw2'
8865 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
8866 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
8867 INTEGER :: handle, handle2, my_pos, nrays, &
8869 INTEGER,
DIMENSION(3) :: nloc
8870#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
8871 LOGICAL :: use_pw_gpu
8873 INTEGER,
DIMENSION(:),
POINTER :: n
8876 CALL timeset(routinen, handle2)
8877 out_unit = cp_logger_get_default_io_unit()
8878 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
8879 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
8884 IF (
PRESENT(debug))
THEN
8891 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
8892 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
8893 cpabort(
"PW grids not compatible")
8895 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
8896 cpabort(
"PW grids have not compatible MPI groups")
8900 n => pw1%pw_grid%npts
8902 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
8908 IF (test .AND. out_unit > 0)
THEN
8909 WRITE (out_unit,
'(A)')
" FFT Protocol "
8910 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
8911 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
8912 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
8915#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
8916 CALL pw_gpu_c1dr3d_3d(pw1, pw2)
8917#elif defined (__PW_FPGA)
8918 ALLOCATE (c_out(n(1), n(2), n(3)))
8921 IF (pw_fpga_init_bitstream(n) == 1)
THEN
8922 CALL pw_scatter_s_c1d_c3d(pw1, c_out)
8924#if (__PW_FPGA_SP && __PW_FPGA)
8925 CALL pw_fpga_c1dr3d_3d_sp(n, c_out)
8927 CALL pw_fpga_c1dr3d_3d_dp(n, c_out)
8929 CALL zdscal(n(1)*n(2)*n(3), 1.0_dp, c_out, 1)
8933 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" PW_SCATTER : 3d -> 1d "
8934 CALL pw_scatter_s_c1d_c3d(pw1, c_out)
8936 CALL fft3d(bwfft, n, c_out, debug=test)
8938 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" REAL part "
8943 ALLOCATE (c_out(n(1), n(2), n(3)))
8944 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" PW_SCATTER : 3d -> 1d "
8945 CALL pw_scatter_s_c1d_c3d(pw1, c_out)
8947 CALL fft3d(bwfft, n, c_out, debug=test)
8949 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" REAL part "
8954 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
8962 IF (test .AND. out_unit > 0)
THEN
8963 WRITE (out_unit,
'(A)')
" FFT Protocol "
8964 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
8965 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
8966 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
8969 my_pos = pw1%pw_grid%para%group%mepos
8970 nrays = pw1%pw_grid%para%nyzray(my_pos)
8971 grays => pw1%pw_grid%grays
8973#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
8975 use_pw_gpu = pw1%pw_grid%para%ray_distribution
8976 IF (use_pw_gpu)
THEN
8977 CALL pw_gpu_c1dr3d_3d_ps(pw1, pw2)
8981 IF (test .AND. out_unit > 0) &
8982 WRITE (out_unit,
'(A)')
" PW_SCATTER : 2d -> 1d "
8984 CALL pw_scatter_p_c1d (pw1, grays)
8985 nloc = pw2%pw_grid%npts_local
8986 ALLOCATE (c_in(nloc(1), nloc(2), nloc(3)))
8988 IF (pw1%pw_grid%para%ray_distribution)
THEN
8989 CALL fft3d(bwfft, n, c_in, grays, pw1%pw_grid%para%group, &
8990 pw1%pw_grid%para%yzp, pw1%pw_grid%para%nyzray, &
8991 pw1%pw_grid%para%bo, debug=test)
8993 CALL fft3d(bwfft, n, c_in, grays, pw1%pw_grid%para%group, &
8994 pw1%pw_grid%para%bo, debug=test)
8997 IF (test .AND. out_unit > 0) &
8998 WRITE (out_unit,
'(A)')
" Real part "
9001#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
9006 IF (test .AND. out_unit > 0)
THEN
9007 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9010 CALL timestop(handle)
9011 CALL timestop(handle2)
9013 END SUBROUTINE fft_wrap_pw1pw2_c1d_r3d_gs_rs
9026 SUBROUTINE pw_gather_s_c1d_r3d_2(pw1, pw2, scale)
9028 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
9029 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
9030 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
9032 CALL pw_gather_s_c1d_r3d (pw2, pw1%array, scale)
9034 END SUBROUTINE pw_gather_s_c1d_r3d_2
9045 SUBROUTINE pw_gather_s_c1d_r3d (pw, c, scale)
9047 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
9048 REAL(kind=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(IN) :: c
9049 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
9051 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_gather_s'
9053 INTEGER :: gpt, handle, l, m, n
9055 CALL timeset(routinen, handle)
9057 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
9058 ngpts =>
SIZE(pw%pw_grid%gsq), ghat => pw%pw_grid%g_hat)
9060 IF (
PRESENT(scale))
THEN
9063 l = mapl(ghat(1, gpt)) + 1
9064 m = mapm(ghat(2, gpt)) + 1
9065 n = mapn(ghat(3, gpt)) + 1
9066 pw%array(gpt) = scale* cmplx(c(l, m, n), 0.0_dp, kind=dp)
9072 l = mapl(ghat(1, gpt)) + 1
9073 m = mapm(ghat(2, gpt)) + 1
9074 n = mapn(ghat(3, gpt)) + 1
9075 pw%array(gpt) = cmplx(c(l, m, n), 0.0_dp, kind=dp)
9082 CALL timestop(handle)
9084 END SUBROUTINE pw_gather_s_c1d_r3d
9095 SUBROUTINE pw_scatter_s_c1d_r3d_2(pw1, pw2, scale)
9097 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
9098 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw2
9099 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
9101 CALL pw_scatter_s_c1d_r3d (pw1, pw2%array, scale)
9103 END SUBROUTINE pw_scatter_s_c1d_r3d_2
9114 SUBROUTINE pw_scatter_s_c1d_r3d (pw, c, scale)
9116 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
9117 REAL(kind=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(INOUT) :: c
9118 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
9120 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_scatter_s'
9122 INTEGER :: gpt, handle, l, m, n
9124 CALL timeset(routinen, handle)
9126 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
9127 ghat => pw%pw_grid%g_hat, ngpts =>
SIZE(pw%pw_grid%gsq))
9130 IF (.NOT.
PRESENT(scale)) c = 0.0_dp
9132 IF (
PRESENT(scale))
THEN
9135 l = mapl(ghat(1, gpt)) + 1
9136 m = mapm(ghat(2, gpt)) + 1
9137 n = mapn(ghat(3, gpt)) + 1
9138 c(l, m, n) = scale* real(pw%array(gpt), kind=dp)
9144 l = mapl(ghat(1, gpt)) + 1
9145 m = mapm(ghat(2, gpt)) + 1
9146 n = mapn(ghat(3, gpt)) + 1
9147 c(l, m, n) = real(pw%array(gpt), kind=dp)
9154 IF (pw%pw_grid%grid_span == halfspace)
THEN
9156 associate(mapl => pw%pw_grid%mapl%neg, mapm => pw%pw_grid%mapm%neg, mapn => pw%pw_grid%mapn%neg, &
9157 ghat => pw%pw_grid%g_hat, ngpts =>
SIZE(pw%pw_grid%gsq))
9159 IF (
PRESENT(scale))
THEN
9162 l = mapl(ghat(1, gpt)) + 1
9163 m = mapm(ghat(2, gpt)) + 1
9164 n = mapn(ghat(3, gpt)) + 1
9165 c(l, m, n) = scale*( real(pw%array(gpt), kind=dp))
9171 l = mapl(ghat(1, gpt)) + 1
9172 m = mapm(ghat(2, gpt)) + 1
9173 n = mapn(ghat(3, gpt)) + 1
9174 c(l, m, n) = ( real(pw%array(gpt), kind=dp))
9183 CALL timestop(handle)
9185 END SUBROUTINE pw_scatter_s_c1d_r3d
9207 SUBROUTINE fft_wrap_pw1pw2_c1d_c3d_gs_rs (pw1, pw2, debug)
9209 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
9210 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw2
9211 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
9213 CHARACTER(len=*),
PARAMETER :: routineN =
'fft_wrap_pw1pw2'
9215 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
9216 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
9217 INTEGER :: handle, handle2, my_pos, nrays, &
9219 INTEGER,
DIMENSION(:),
POINTER :: n
9222 CALL timeset(routinen, handle2)
9223 out_unit = cp_logger_get_default_io_unit()
9224 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
9225 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
9230 IF (
PRESENT(debug))
THEN
9237 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
9238 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
9239 cpabort(
"PW grids not compatible")
9241 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
9242 cpabort(
"PW grids have not compatible MPI groups")
9246 n => pw1%pw_grid%npts
9248 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
9254 IF (test .AND. out_unit > 0)
THEN
9255 WRITE (out_unit,
'(A)')
" FFT Protocol "
9256 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
9257 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
9258 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
9262 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" PW_SCATTER : 3d -> 1d "
9263 CALL pw_scatter_s_c1d_c3d(pw1, c_out)
9264 CALL fft3d(bwfft, n, c_out, debug=test)
9266 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9274 IF (test .AND. out_unit > 0)
THEN
9275 WRITE (out_unit,
'(A)')
" FFT Protocol "
9276 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
9277 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
9278 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
9281 my_pos = pw1%pw_grid%para%group%mepos
9282 nrays = pw1%pw_grid%para%nyzray(my_pos)
9283 grays => pw1%pw_grid%grays
9286 IF (test .AND. out_unit > 0) &
9287 WRITE (out_unit,
'(A)')
" PW_SCATTER : 2d -> 1d "
9289 CALL pw_scatter_p_c1d (pw1, grays)
9292 IF (pw1%pw_grid%para%ray_distribution)
THEN
9293 CALL fft3d(bwfft, n, c_in, grays, pw1%pw_grid%para%group, &
9294 pw1%pw_grid%para%yzp, pw1%pw_grid%para%nyzray, &
9295 pw1%pw_grid%para%bo, debug=test)
9297 CALL fft3d(bwfft, n, c_in, grays, pw1%pw_grid%para%group, &
9298 pw1%pw_grid%para%bo, debug=test)
9303 IF (test .AND. out_unit > 0)
THEN
9304 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9307 CALL timestop(handle)
9308 CALL timestop(handle2)
9310 END SUBROUTINE fft_wrap_pw1pw2_c1d_c3d_gs_rs
9323 SUBROUTINE pw_gather_s_c1d_c3d_2(pw1, pw2, scale)
9325 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
9326 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
9327 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
9329 CALL pw_gather_s_c1d_c3d (pw2, pw1%array, scale)
9331 END SUBROUTINE pw_gather_s_c1d_c3d_2
9342 SUBROUTINE pw_gather_s_c1d_c3d (pw, c, scale)
9344 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
9345 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(IN) :: c
9346 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
9348 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_gather_s'
9350 INTEGER :: gpt, handle, l, m, n
9352 CALL timeset(routinen, handle)
9354 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
9355 ngpts =>
SIZE(pw%pw_grid%gsq), ghat => pw%pw_grid%g_hat)
9357 IF (
PRESENT(scale))
THEN
9360 l = mapl(ghat(1, gpt)) + 1
9361 m = mapm(ghat(2, gpt)) + 1
9362 n = mapn(ghat(3, gpt)) + 1
9363 pw%array(gpt) = scale* c(l, m, n)
9369 l = mapl(ghat(1, gpt)) + 1
9370 m = mapm(ghat(2, gpt)) + 1
9371 n = mapn(ghat(3, gpt)) + 1
9372 pw%array(gpt) = c(l, m, n)
9379 CALL timestop(handle)
9381 END SUBROUTINE pw_gather_s_c1d_c3d
9392 SUBROUTINE pw_scatter_s_c1d_c3d_2(pw1, pw2, scale)
9394 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
9395 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
9396 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
9398 CALL pw_scatter_s_c1d_c3d (pw1, pw2%array, scale)
9400 END SUBROUTINE pw_scatter_s_c1d_c3d_2
9411 SUBROUTINE pw_scatter_s_c1d_c3d (pw, c, scale)
9413 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
9414 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(INOUT) :: c
9415 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
9417 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_scatter_s'
9419 INTEGER :: gpt, handle, l, m, n
9421 CALL timeset(routinen, handle)
9423 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
9424 ghat => pw%pw_grid%g_hat, ngpts =>
SIZE(pw%pw_grid%gsq))
9427 IF (.NOT.
PRESENT(scale)) c = 0.0_dp
9429 IF (
PRESENT(scale))
THEN
9432 l = mapl(ghat(1, gpt)) + 1
9433 m = mapm(ghat(2, gpt)) + 1
9434 n = mapn(ghat(3, gpt)) + 1
9435 c(l, m, n) = scale* pw%array(gpt)
9441 l = mapl(ghat(1, gpt)) + 1
9442 m = mapm(ghat(2, gpt)) + 1
9443 n = mapn(ghat(3, gpt)) + 1
9444 c(l, m, n) = pw%array(gpt)
9451 IF (pw%pw_grid%grid_span == halfspace)
THEN
9453 associate(mapl => pw%pw_grid%mapl%neg, mapm => pw%pw_grid%mapm%neg, mapn => pw%pw_grid%mapn%neg, &
9454 ghat => pw%pw_grid%g_hat, ngpts =>
SIZE(pw%pw_grid%gsq))
9456 IF (
PRESENT(scale))
THEN
9459 l = mapl(ghat(1, gpt)) + 1
9460 m = mapm(ghat(2, gpt)) + 1
9461 n = mapn(ghat(3, gpt)) + 1
9462 c(l, m, n) = scale*conjg( pw%array(gpt))
9468 l = mapl(ghat(1, gpt)) + 1
9469 m = mapm(ghat(2, gpt)) + 1
9470 n = mapn(ghat(3, gpt)) + 1
9471 c(l, m, n) = conjg( pw%array(gpt))
9480 CALL timestop(handle)
9482 END SUBROUTINE pw_scatter_s_c1d_c3d
9504 SUBROUTINE fft_wrap_pw1pw2_c3d_r3d_gs_rs (pw1, pw2, debug)
9506 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
9507 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw2
9508 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
9510 CHARACTER(len=*),
PARAMETER :: routineN =
'fft_wrap_pw1pw2'
9512 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
9513 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
9514 INTEGER :: handle, handle2, my_pos, nrays, &
9516 INTEGER,
DIMENSION(:),
POINTER :: n
9519 CALL timeset(routinen, handle2)
9520 out_unit = cp_logger_get_default_io_unit()
9521 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
9522 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
9527 IF (
PRESENT(debug))
THEN
9534 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
9535 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
9536 cpabort(
"PW grids not compatible")
9538 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
9539 cpabort(
"PW grids have not compatible MPI groups")
9543 n => pw1%pw_grid%npts
9545 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
9551 IF (test .AND. out_unit > 0)
THEN
9552 WRITE (out_unit,
'(A)')
" FFT Protocol "
9553 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
9554 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
9555 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
9559 ALLOCATE (c_out(n(1), n(2), n(3)))
9560 CALL fft3d(bwfft, n, c_in, c_out, debug=test)
9562 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" REAL part "
9563 pw2%array = real(c_out, kind=dp)
9566 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9574 IF (test .AND. out_unit > 0)
THEN
9575 WRITE (out_unit,
'(A)')
" FFT Protocol "
9576 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
9577 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
9578 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
9581 my_pos = pw1%pw_grid%para%group%mepos
9582 nrays = pw1%pw_grid%para%nyzray(my_pos)
9583 grays => pw1%pw_grid%grays
9587 IF (test .AND. out_unit > 0)
THEN
9588 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9591 CALL timestop(handle)
9592 CALL timestop(handle2)
9594 END SUBROUTINE fft_wrap_pw1pw2_c3d_r3d_gs_rs
9612 SUBROUTINE fft_wrap_pw1pw2_c3d_c1d_rs_gs (pw1, pw2, debug)
9614 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
9615 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
9616 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
9618 CHARACTER(len=*),
PARAMETER :: routineN =
'fft_wrap_pw1pw2'
9620 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
9621 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
9622 INTEGER :: handle, handle2, my_pos, nrays, &
9624 INTEGER,
DIMENSION(:),
POINTER :: n
9627 CALL timeset(routinen, handle2)
9628 out_unit = cp_logger_get_default_io_unit()
9629 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
9630 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
9635 IF (
PRESENT(debug))
THEN
9642 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
9643 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
9644 cpabort(
"PW grids not compatible")
9646 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
9647 cpabort(
"PW grids have not compatible MPI groups")
9651 n => pw1%pw_grid%npts
9653 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
9659 IF (test .AND. out_unit > 0)
THEN
9660 WRITE (out_unit,
'(A)')
" FFT Protocol "
9661 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
9662 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
9663 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
9667 ALLOCATE (c_out(n(1), n(2), n(3)))
9669 CALL fft3d(fwfft, n, c_in, c_out, debug=test)
9671 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" PW_GATHER : 3d -> 1d "
9672 CALL pw_gather_s_c1d_c3d(pw2, c_out)
9675 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9683 IF (test .AND. out_unit > 0)
THEN
9684 WRITE (out_unit,
'(A)')
" FFT Protocol "
9685 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
9686 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
9687 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
9690 my_pos = pw1%pw_grid%para%group%mepos
9691 nrays = pw1%pw_grid%para%nyzray(my_pos)
9692 grays => pw1%pw_grid%grays
9698 IF (pw1%pw_grid%para%ray_distribution)
THEN
9699 CALL fft3d(fwfft, n, c_in, grays, pw1%pw_grid%para%group, &
9700 pw1%pw_grid%para%yzp, pw1%pw_grid%para%nyzray, &
9701 pw1%pw_grid%para%bo, debug=test)
9703 CALL fft3d(fwfft, n, c_in, grays, pw1%pw_grid%para%group, &
9704 pw1%pw_grid%para%bo, debug=test)
9707 IF (test .AND. out_unit > 0) &
9708 WRITE (out_unit,
'(A)')
" PW_GATHER : 2d -> 1d "
9709 CALL pw_gather_p_c1d (pw2, grays)
9712 IF (test .AND. out_unit > 0)
THEN
9713 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9716 CALL timestop(handle)
9717 CALL timestop(handle2)
9719 END SUBROUTINE fft_wrap_pw1pw2_c3d_c1d_rs_gs
9738 SUBROUTINE fft_wrap_pw1pw2_c3d_c3d_rs_gs (pw1, pw2, debug)
9740 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
9741 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
9742 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
9744 CHARACTER(len=*),
PARAMETER :: routineN =
'fft_wrap_pw1pw2'
9746 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
9747 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
9748 INTEGER :: handle, handle2, my_pos, nrays, &
9750 INTEGER,
DIMENSION(:),
POINTER :: n
9753 CALL timeset(routinen, handle2)
9754 out_unit = cp_logger_get_default_io_unit()
9755 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
9756 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
9761 IF (
PRESENT(debug))
THEN
9768 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
9769 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
9770 cpabort(
"PW grids not compatible")
9772 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
9773 cpabort(
"PW grids have not compatible MPI groups")
9777 n => pw1%pw_grid%npts
9779 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
9785 IF (test .AND. out_unit > 0)
THEN
9786 WRITE (out_unit,
'(A)')
" FFT Protocol "
9787 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
9788 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
9789 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
9794 CALL fft3d(fwfft, n, c_in, c_out, debug=test)
9796 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9804 IF (test .AND. out_unit > 0)
THEN
9805 WRITE (out_unit,
'(A)')
" FFT Protocol "
9806 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
9807 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
9808 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
9811 my_pos = pw1%pw_grid%para%group%mepos
9812 nrays = pw1%pw_grid%para%nyzray(my_pos)
9813 grays => pw1%pw_grid%grays
9817 IF (test .AND. out_unit > 0)
THEN
9818 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9821 CALL timestop(handle)
9822 CALL timestop(handle2)
9824 END SUBROUTINE fft_wrap_pw1pw2_c3d_c3d_rs_gs
9839 SUBROUTINE fft_wrap_pw1pw2_c3d_c3d_gs_rs (pw1, pw2, debug)
9841 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
9842 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw2
9843 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
9845 CHARACTER(len=*),
PARAMETER :: routineN =
'fft_wrap_pw1pw2'
9847 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
9848 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
9849 INTEGER :: handle, handle2, my_pos, nrays, &
9851 INTEGER,
DIMENSION(:),
POINTER :: n
9854 CALL timeset(routinen, handle2)
9855 out_unit = cp_logger_get_default_io_unit()
9856 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
9857 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
9862 IF (
PRESENT(debug))
THEN
9869 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
9870 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
9871 cpabort(
"PW grids not compatible")
9873 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
9874 cpabort(
"PW grids have not compatible MPI groups")
9878 n => pw1%pw_grid%npts
9880 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
9886 IF (test .AND. out_unit > 0)
THEN
9887 WRITE (out_unit,
'(A)')
" FFT Protocol "
9888 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
9889 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
9890 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
9895 CALL fft3d(bwfft, n, c_in, c_out, debug=test)
9897 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9905 IF (test .AND. out_unit > 0)
THEN
9906 WRITE (out_unit,
'(A)')
" FFT Protocol "
9907 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
9908 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
9909 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
9912 my_pos = pw1%pw_grid%para%group%mepos
9913 nrays = pw1%pw_grid%para%nyzray(my_pos)
9914 grays => pw1%pw_grid%grays
9918 IF (test .AND. out_unit > 0)
THEN
9919 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9922 CALL timestop(handle)
9923 CALL timestop(handle2)
9925 END SUBROUTINE fft_wrap_pw1pw2_c3d_c3d_gs_rs
9944 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
9945 REAL(kind=dp),
INTENT(IN) :: omega
9947 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_gauss_damp'
9949 INTEGER :: handle, i
9950 REAL(kind=dp) :: omega_2, tmp
9952 CALL timeset(routinen, handle)
9953 cpassert(omega >= 0)
9955 omega_2 = omega*omega
9956 omega_2 = 0.25_dp/omega_2
9959 DO i = 1,
SIZE(pw%array)
9960 tmp = exp(-pw%pw_grid%gsq(i)*omega_2)
9961 pw%array(i) = pw%array(i)*tmp
9965 CALL timestop(handle)
9978 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
9979 REAL(kind=dp),
INTENT(IN) :: omega
9981 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_log_deriv_gauss'
9983 INTEGER :: handle, i
9984 REAL(kind=dp) :: omega_2, tmp
9986 CALL timeset(routinen, handle)
9987 cpassert(omega >= 0)
9989 omega_2 = omega*omega
9990 omega_2 = 0.25_dp/omega_2
9993 DO i = 1,
SIZE(pw%array)
9994 tmp = 1.0_dp + omega_2*pw%pw_grid%gsq(i)
9995 pw%array(i) = pw%array(i)*tmp
9999 CALL timestop(handle)
10017 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
10018 REAL(kind=dp),
INTENT(IN) :: omega
10020 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_compl_gauss_damp'
10022 INTEGER :: cnt, handle, i
10023 REAL(kind=dp) :: omega_2, tmp, tmp2
10025 CALL timeset(routinen, handle)
10027 omega_2 = omega*omega
10028 omega_2 = 0.25_dp/omega_2
10030 cnt =
SIZE(pw%array)
10034 tmp = -omega_2*pw%pw_grid%gsq(i)
10035 IF (abs(tmp) > 1.0e-5_dp)
THEN
10036 tmp2 = 1.0_dp - exp(tmp)
10038 tmp2 = tmp + 0.5_dp*tmp*(tmp + (1.0_dp/3.0_dp)*tmp**2)
10040 pw%array(i) = pw%array(i)*tmp2
10044 CALL timestop(handle)
10056 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
10057 REAL(kind=dp),
INTENT(IN) :: omega
10059 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_log_deriv_compl_gauss'
10061 INTEGER :: handle, i
10062 REAL(kind=dp) :: omega_2, tmp, tmp2
10064 CALL timeset(routinen, handle)
10066 omega_2 = omega*omega
10067 omega_2 = 0.25_dp/omega_2
10071 DO i = 1,
SIZE(pw%array)
10072 tmp = omega_2*pw%pw_grid%gsq(i)
10074 IF (abs(tmp) >= 0.003_dp)
THEN
10075 tmp2 = 1.0_dp - tmp*exp(-tmp)/(1.0_dp - exp(-tmp))
10077 tmp2 = 0.5_dp*tmp - tmp**2/12.0_dp
10079 pw%array(i) = pw%array(i)*tmp2
10083 CALL timestop(handle)
10104 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
10105 REAL(kind=dp),
INTENT(IN) :: omega, scale_coul, scale_long
10107 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_gauss_damp_mix'
10109 INTEGER :: handle, i
10110 REAL(kind=dp) :: omega_2, tmp
10112 CALL timeset(routinen, handle)
10114 omega_2 = omega*omega
10115 omega_2 = 0.25_dp/omega_2
10118 DO i = 1,
SIZE(pw%array)
10119 tmp = scale_coul + scale_long*exp(-pw%pw_grid%gsq(i)*omega_2)
10120 pw%array(i) = pw%array(i)*tmp
10124 CALL timestop(handle)
10139 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
10140 REAL(kind=dp),
INTENT(IN) :: omega, scale_coul, scale_long
10142 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_log_deriv_mix_cl'
10144 INTEGER :: handle, i
10145 REAL(kind=dp) :: omega_2, tmp, tmp2
10147 CALL timeset(routinen, handle)
10149 omega_2 = omega*omega
10150 omega_2 = 0.25_dp/omega_2
10154 DO i = 1,
SIZE(pw%array)
10155 tmp = omega_2*pw%pw_grid%gsq(i)
10156 tmp2 = 1.0_dp + scale_long*tmp*exp(-tmp)/(scale_coul + scale_long*exp(-tmp))
10157 pw%array(i) = pw%array(i)*tmp2
10161 CALL timestop(handle)
10180 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
10181 REAL(kind=dp),
INTENT(IN) :: rcutoff
10183 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_truncated'
10185 INTEGER :: handle, i
10186 REAL(kind=dp) :: tmp, tmp2
10188 CALL timeset(routinen, handle)
10189 cpassert(rcutoff >= 0)
10192 DO i = 1,
SIZE(pw%array)
10193 tmp = sqrt(pw%pw_grid%gsq(i))*rcutoff
10194 IF (tmp >= 0.005_dp)
THEN
10195 tmp2 = 1.0_dp - cos(tmp)
10197 tmp2 = tmp**2/2.0_dp*(1.0 - tmp**2/12.0_dp)
10199 pw%array(i) = pw%array(i)*tmp2
10203 CALL timestop(handle)
10216 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
10217 REAL(kind=dp),
INTENT(IN) :: rcutoff
10219 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_log_deriv_trunc'
10221 INTEGER :: handle, i
10222 REAL(kind=dp) :: rchalf, tmp, tmp2
10224 CALL timeset(routinen, handle)
10225 cpassert(rcutoff >= 0)
10227 rchalf = 0.5_dp*rcutoff
10230 DO i = 1,
SIZE(pw%array)
10231 tmp = rchalf*sqrt(pw%pw_grid%gsq(i))
10233 IF (abs(tmp) >= 0.0001_dp)
THEN
10234 tmp2 = 1.0_dp - tmp/tan(tmp)
10236 tmp2 = tmp**2*(1.0_dp + tmp**2/15.0_dp)/3.0_dp
10238 pw%array(i) = pw%array(i)*tmp2
10242 CALL timestop(handle)
10258 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
10259 INTEGER,
DIMENSION(3),
INTENT(IN) :: n
10261 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_derive'
10263 COMPLEX(KIND=dp) :: im
10264 INTEGER :: handle, m,
idx, idir
10266 CALL timeset(routinen, handle)
10269 cpabort(
"Nonnegative exponents are not supported!")
10275 IF (n(idir) == 1)
THEN
10277 DO idx = 1,
SIZE(pw%array)
10278 pw%array(
idx) = pw%array(
idx)*pw%pw_grid%g(idir,
idx)
10281 ELSE IF (n(idir) > 1)
THEN
10283 DO idx = 1,
SIZE(pw%array)
10284 pw%array(
idx) = pw%array(
idx)*pw%pw_grid%g(idir,
idx)**n(idir)
10292 IF (abs(real(im, kind=dp) - 1.0_dp) > 1.0e-10_dp)
THEN
10294 pw%array(:) = im*pw%array(:)
10298 CALL timestop(handle)
10313 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
10315 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_laplace'
10319 CALL timeset(routinen, handle)
10322 pw%array(:) = -pw%array(:)*pw%pw_grid%gsq(:)
10325 CALL timestop(handle)
10342 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw, pwdr2
10343 INTEGER,
INTENT(IN) :: i, j
10345 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_dr2'
10347 INTEGER :: cnt, handle, ig
10348 REAL(kind=dp) :: gg, o3
10350 CALL timeset(routinen, handle)
10354 cnt =
SIZE(pw%array)
10359 gg = pw%pw_grid%g(i, ig)**2 - o3*pw%pw_grid%gsq(ig)
10360 pwdr2%array(ig) = gg*pw%array(ig)
10366 pwdr2%array(ig) = pw%array(ig)*(pw%pw_grid%g(i, ig)*pw%pw_grid%g(j, ig))
10371 CALL timestop(handle)
10390 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw, pwdr2_gg
10391 INTEGER,
INTENT(IN) :: i, j
10393 INTEGER :: cnt, handle, ig
10394 REAL(kind=dp) :: gg, o3
10395 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_dr2_gg'
10397 CALL timeset(routinen, handle)
10401 cnt =
SIZE(pw%array)
10405 DO ig = pw%pw_grid%first_gne0, cnt
10406 gg = pw%pw_grid%g(i, ig)**2 - o3*pw%pw_grid%gsq(ig)
10407 pwdr2_gg%array(ig) = gg/pw%pw_grid%gsq(ig)*pw%array(ig)
10412 DO ig = pw%pw_grid%first_gne0, cnt
10413 pwdr2_gg%array(ig) = pw%array(ig)*((pw%pw_grid%g(i, ig)*pw%pw_grid%g(j, ig)) &
10414 /pw%pw_grid%gsq(ig))
10419 IF (pw%pw_grid%have_g0) pwdr2_gg%array(1) = 0.0_dp
10421 CALL timestop(handle)
10438 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
10439 REAL(kind=dp),
INTENT(IN) :: ecut, sigma
10441 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_smoothing'
10443 INTEGER :: cnt, handle, ig
10444 REAL(kind=dp) :: arg, f
10446 CALL timeset(routinen, handle)
10448 cnt =
SIZE(pw%array)
10452 arg = (ecut - pw%pw_grid%gsq(ig))/sigma
10453 f = exp(arg)/(1 + exp(arg))
10454 pw%array(ig) = f*pw%array(ig)
10458 CALL timestop(handle)
10468 ELEMENTAL FUNCTION pw_compatible(grida, gridb)
RESULT(compat)
10470 TYPE(pw_grid_type),
INTENT(IN) :: grida, gridb
10475 IF (grida%id_nr == gridb%id_nr)
THEN
10477 ELSE IF (grida%reference == gridb%id_nr)
THEN
10479 ELSE IF (gridb%reference == grida%id_nr)
THEN
10483 END FUNCTION pw_compatible
10496 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: sf
10497 REAL(kind=dp),
DIMENSION(:),
INTENT(IN) :: r
10499 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_structure_factor'
10501 INTEGER :: cnt, handle, ig
10502 REAL(kind=dp) :: arg
10504 CALL timeset(routinen, handle)
10506 cnt =
SIZE(sf%array)
10510 arg = dot_product(sf%pw_grid%g(:, ig), r)
10511 sf%array(ig) = cmplx(cos(arg), -sin(arg), kind=dp)
10515 CALL timestop(handle)
static GRID_HOST_DEVICE int idx(const orbital a)
Return coset index of given orbital angular momentum.
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Defines the basic variable types.
integer, parameter, public dp
Machine interface based on Fortran 2003 and POSIX.
subroutine, public m_memory(mem)
Returns the total amount of memory [bytes] in use, if known, zero otherwise.
Definition of mathematical constants and functions.
complex(kind=dp), parameter, public gaussi
complex(kind=dp), parameter, public z_zero
subroutine, public pw_copy_match(pw1, pw2)
copy a pw type variable
integer function, public pw_fpga_init_bitstream(n)
Invoke the pw_fpga_check_bitstream C function passing the path to the data dir.
subroutine, public pw_fpga_r3dc1d_3d_dp(n, c_out)
perform an in-place double precision fft3d on the FPGA
subroutine, public pw_fpga_r3dc1d_3d_sp(n, c_out)
perform an in-place single precision fft3d on the FPGA
subroutine, public pw_fpga_c1dr3d_3d_dp(n, c_out)
perform an in-place double precision inverse fft3d on the FPGA
subroutine, public pw_fpga_c1dr3d_3d_sp(n, c_out)
perform an in-place single precision inverse fft3d on the FPGA
subroutine, public pw_gpu_c1dr3d_3d_ps(pw1, pw2)
perform an parallel scatter followed by a fft on the gpu
subroutine, public pw_gpu_c1dr3d_3d(pw1, pw2)
perform an scatter followed by a fft on the gpu
subroutine, public pw_gpu_r3dc1d_3d(pw1, pw2)
perform an fft followed by a gather on the gpu
subroutine, public pw_gpu_r3dc1d_3d_ps(pw1, pw2)
perform an parallel fft followed by a gather on the gpu
integer, parameter, public halfspace
integer, parameter, public pw_mode_local
integer, parameter, public pw_mode_distributed
subroutine, public pw_smoothing(pw, ecut, sigma)
Multiplies a G-space function with a smoothing factor of the form f(|G|) = exp((ecut - G^2)/sigma)/(1...
subroutine, public pw_log_deriv_trunc(pw, rcutoff)
Multiply all data points with the logarithmic derivative of the complementary cosine This function is...
subroutine, public pw_gauss_damp(pw, omega)
Multiply all data points with a Gaussian damping factor Needed for longrange Coulomb potential V(\vec...
subroutine, public pw_log_deriv_compl_gauss(pw, omega)
Multiply all data points with the logarithmic derivative of the complementary Gaussian damping factor...
subroutine, public pw_dr2_gg(pw, pwdr2_gg, i, j)
Calculate the tensorial 2nd derivative of a plane wave vector and divides by |G|^2....
subroutine, public pw_laplace(pw)
Calculate the Laplacian of a plane wave vector.
subroutine, public pw_log_deriv_gauss(pw, omega)
Multiply all data points with the logarithmic derivative of a Gaussian.
subroutine, public pw_gauss_damp_mix(pw, omega, scale_coul, scale_long)
Multiply all data points with a Gaussian damping factor and mixes it with the original function Neede...
subroutine, public pw_truncated(pw, rcutoff)
Multiply all data points with a complementary cosine Needed for truncated Coulomb potential V(\vec r)...
integer, parameter, public do_accurate_sum
integer, parameter, public do_standard_sum
subroutine, public pw_derive(pw, n)
Calculate the derivative of a plane wave vector.
subroutine, public pw_structure_factor(sf, r)
Calculate the structure factor for point r.
subroutine, public pw_compl_gauss_damp(pw, omega)
Multiply all data points with a Gaussian damping factor Needed for longrange Coulomb potential V(\vec...
subroutine, public pw_log_deriv_mix_cl(pw, omega, scale_coul, scale_long)
Multiply all data points with the logarithmic derivative of the mixed longrange/Coulomb potential Nee...
subroutine, public pw_dr2(pw, pwdr2, i, j)
Calculate the tensorial 2nd derivative of a plane wave vector.