60#include "../base/base_uses.f90"
77 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'pw_methods'
78 LOGICAL,
PARAMETER,
PRIVATE :: debug_this_module = .false.
83 MODULE PROCEDURE pw_zero_r1d_rs
84 MODULE PROCEDURE pw_zero_r3d_rs
85 MODULE PROCEDURE pw_zero_c1d_rs
86 MODULE PROCEDURE pw_zero_c3d_rs
87 MODULE PROCEDURE pw_zero_r1d_gs
88 MODULE PROCEDURE pw_zero_r3d_gs
89 MODULE PROCEDURE pw_zero_c1d_gs
90 MODULE PROCEDURE pw_zero_c3d_gs
94 MODULE PROCEDURE pw_scale_r1d_rs
95 MODULE PROCEDURE pw_scale_r3d_rs
96 MODULE PROCEDURE pw_scale_c1d_rs
97 MODULE PROCEDURE pw_scale_c3d_rs
98 MODULE PROCEDURE pw_scale_r1d_gs
99 MODULE PROCEDURE pw_scale_r3d_gs
100 MODULE PROCEDURE pw_scale_c1d_gs
101 MODULE PROCEDURE pw_scale_c3d_gs
105 MODULE PROCEDURE pw_write_r1d_rs
106 MODULE PROCEDURE pw_write_r3d_rs
107 MODULE PROCEDURE pw_write_c1d_rs
108 MODULE PROCEDURE pw_write_c3d_rs
109 MODULE PROCEDURE pw_write_r1d_gs
110 MODULE PROCEDURE pw_write_r3d_gs
111 MODULE PROCEDURE pw_write_c1d_gs
112 MODULE PROCEDURE pw_write_c3d_gs
116 MODULE PROCEDURE pw_integrate_function_r1d_rs
117 MODULE PROCEDURE pw_integrate_function_r3d_rs
118 MODULE PROCEDURE pw_integrate_function_c1d_rs
119 MODULE PROCEDURE pw_integrate_function_c3d_rs
120 MODULE PROCEDURE pw_integrate_function_r1d_gs
121 MODULE PROCEDURE pw_integrate_function_r3d_gs
122 MODULE PROCEDURE pw_integrate_function_c1d_gs
123 MODULE PROCEDURE pw_integrate_function_c3d_gs
127 MODULE PROCEDURE pw_set_value_r1d_rs
128 MODULE PROCEDURE pw_zero_r1d_rs
129 MODULE PROCEDURE pw_set_value_r3d_rs
130 MODULE PROCEDURE pw_zero_r3d_rs
131 MODULE PROCEDURE pw_set_value_c1d_rs
132 MODULE PROCEDURE pw_zero_c1d_rs
133 MODULE PROCEDURE pw_set_value_c3d_rs
134 MODULE PROCEDURE pw_zero_c3d_rs
135 MODULE PROCEDURE pw_set_value_r1d_gs
136 MODULE PROCEDURE pw_zero_r1d_gs
137 MODULE PROCEDURE pw_set_value_r3d_gs
138 MODULE PROCEDURE pw_zero_r3d_gs
139 MODULE PROCEDURE pw_set_value_c1d_gs
140 MODULE PROCEDURE pw_zero_c1d_gs
141 MODULE PROCEDURE pw_set_value_c3d_gs
142 MODULE PROCEDURE pw_zero_c3d_gs
146 MODULE PROCEDURE pw_copy_r1d_r1d_rs
147 MODULE PROCEDURE pw_copy_r1d_c1d_rs
148 MODULE PROCEDURE pw_copy_r3d_r3d_rs
149 MODULE PROCEDURE pw_copy_r3d_c3d_rs
150 MODULE PROCEDURE pw_copy_c1d_r1d_rs
151 MODULE PROCEDURE pw_copy_c1d_c1d_rs
152 MODULE PROCEDURE pw_copy_c3d_r3d_rs
153 MODULE PROCEDURE pw_copy_c3d_c3d_rs
154 MODULE PROCEDURE pw_copy_r1d_r1d_gs
155 MODULE PROCEDURE pw_copy_r1d_c1d_gs
156 MODULE PROCEDURE pw_copy_r3d_r3d_gs
157 MODULE PROCEDURE pw_copy_r3d_c3d_gs
158 MODULE PROCEDURE pw_copy_c1d_r1d_gs
159 MODULE PROCEDURE pw_copy_c1d_c1d_gs
160 MODULE PROCEDURE pw_copy_c3d_r3d_gs
161 MODULE PROCEDURE pw_copy_c3d_c3d_gs
165 MODULE PROCEDURE pw_axpy_r1d_r1d_rs
166 MODULE PROCEDURE pw_axpy_r1d_c1d_rs
167 MODULE PROCEDURE pw_axpy_r3d_r3d_rs
168 MODULE PROCEDURE pw_axpy_r3d_c3d_rs
169 MODULE PROCEDURE pw_axpy_c1d_r1d_rs
170 MODULE PROCEDURE pw_axpy_c1d_c1d_rs
171 MODULE PROCEDURE pw_axpy_c3d_r3d_rs
172 MODULE PROCEDURE pw_axpy_c3d_c3d_rs
173 MODULE PROCEDURE pw_axpy_r1d_r1d_gs
174 MODULE PROCEDURE pw_axpy_r1d_c1d_gs
175 MODULE PROCEDURE pw_axpy_r3d_r3d_gs
176 MODULE PROCEDURE pw_axpy_r3d_c3d_gs
177 MODULE PROCEDURE pw_axpy_c1d_r1d_gs
178 MODULE PROCEDURE pw_axpy_c1d_c1d_gs
179 MODULE PROCEDURE pw_axpy_c3d_r3d_gs
180 MODULE PROCEDURE pw_axpy_c3d_c3d_gs
184 MODULE PROCEDURE pw_multiply_r1d_r1d_rs
185 MODULE PROCEDURE pw_multiply_r1d_c1d_rs
186 MODULE PROCEDURE pw_multiply_r3d_r3d_rs
187 MODULE PROCEDURE pw_multiply_r3d_c3d_rs
188 MODULE PROCEDURE pw_multiply_c1d_r1d_rs
189 MODULE PROCEDURE pw_multiply_c1d_c1d_rs
190 MODULE PROCEDURE pw_multiply_c3d_r3d_rs
191 MODULE PROCEDURE pw_multiply_c3d_c3d_rs
192 MODULE PROCEDURE pw_multiply_r1d_r1d_gs
193 MODULE PROCEDURE pw_multiply_r1d_c1d_gs
194 MODULE PROCEDURE pw_multiply_r3d_r3d_gs
195 MODULE PROCEDURE pw_multiply_r3d_c3d_gs
196 MODULE PROCEDURE pw_multiply_c1d_r1d_gs
197 MODULE PROCEDURE pw_multiply_c1d_c1d_gs
198 MODULE PROCEDURE pw_multiply_c3d_r3d_gs
199 MODULE PROCEDURE pw_multiply_c3d_c3d_gs
203 MODULE PROCEDURE pw_multiply_with_r1d_r1d_rs
204 MODULE PROCEDURE pw_multiply_with_r1d_c1d_rs
205 MODULE PROCEDURE pw_multiply_with_r3d_r3d_rs
206 MODULE PROCEDURE pw_multiply_with_r3d_c3d_rs
207 MODULE PROCEDURE pw_multiply_with_c1d_r1d_rs
208 MODULE PROCEDURE pw_multiply_with_c1d_c1d_rs
209 MODULE PROCEDURE pw_multiply_with_c3d_r3d_rs
210 MODULE PROCEDURE pw_multiply_with_c3d_c3d_rs
211 MODULE PROCEDURE pw_multiply_with_r1d_r1d_gs
212 MODULE PROCEDURE pw_multiply_with_r1d_c1d_gs
213 MODULE PROCEDURE pw_multiply_with_r3d_r3d_gs
214 MODULE PROCEDURE pw_multiply_with_r3d_c3d_gs
215 MODULE PROCEDURE pw_multiply_with_c1d_r1d_gs
216 MODULE PROCEDURE pw_multiply_with_c1d_c1d_gs
217 MODULE PROCEDURE pw_multiply_with_c3d_r3d_gs
218 MODULE PROCEDURE pw_multiply_with_c3d_c3d_gs
222 MODULE PROCEDURE pw_integral_ab_r1d_r1d_rs
223 MODULE PROCEDURE pw_integral_ab_r1d_c1d_rs
224 MODULE PROCEDURE pw_integral_ab_r3d_r3d_rs
225 MODULE PROCEDURE pw_integral_ab_r3d_c3d_rs
226 MODULE PROCEDURE pw_integral_ab_c1d_r1d_rs
227 MODULE PROCEDURE pw_integral_ab_c1d_c1d_rs
228 MODULE PROCEDURE pw_integral_ab_c3d_r3d_rs
229 MODULE PROCEDURE pw_integral_ab_c3d_c3d_rs
230 MODULE PROCEDURE pw_integral_ab_r1d_r1d_gs
231 MODULE PROCEDURE pw_integral_ab_r1d_c1d_gs
232 MODULE PROCEDURE pw_integral_ab_r3d_r3d_gs
233 MODULE PROCEDURE pw_integral_ab_r3d_c3d_gs
234 MODULE PROCEDURE pw_integral_ab_c1d_r1d_gs
235 MODULE PROCEDURE pw_integral_ab_c1d_c1d_gs
236 MODULE PROCEDURE pw_integral_ab_c3d_r3d_gs
237 MODULE PROCEDURE pw_integral_ab_c3d_c3d_gs
241 MODULE PROCEDURE pw_integral_a2b_r1d_r1d
242 MODULE PROCEDURE pw_integral_a2b_r1d_c1d
243 MODULE PROCEDURE pw_integral_a2b_c1d_r1d
244 MODULE PROCEDURE pw_integral_a2b_c1d_c1d
248 MODULE PROCEDURE pw_gather_p_r1d
249 MODULE PROCEDURE pw_gather_p_c1d
250 MODULE PROCEDURE pw_gather_s_r1d_r3d
251 MODULE PROCEDURE pw_gather_s_r1d_c3d
252 MODULE PROCEDURE pw_gather_s_c1d_r3d
253 MODULE PROCEDURE pw_gather_s_c1d_c3d
257 MODULE PROCEDURE pw_scatter_p_r1d
258 MODULE PROCEDURE pw_scatter_p_c1d
259 MODULE PROCEDURE pw_scatter_s_r1d_r3d
260 MODULE PROCEDURE pw_scatter_s_r1d_c3d
261 MODULE PROCEDURE pw_scatter_s_c1d_r3d
262 MODULE PROCEDURE pw_scatter_s_c1d_c3d
266 MODULE PROCEDURE pw_copy_to_array_r1d_r1d_rs
267 MODULE PROCEDURE pw_copy_to_array_r1d_c1d_rs
268 MODULE PROCEDURE pw_copy_to_array_r3d_r3d_rs
269 MODULE PROCEDURE pw_copy_to_array_r3d_c3d_rs
270 MODULE PROCEDURE pw_copy_to_array_c1d_r1d_rs
271 MODULE PROCEDURE pw_copy_to_array_c1d_c1d_rs
272 MODULE PROCEDURE pw_copy_to_array_c3d_r3d_rs
273 MODULE PROCEDURE pw_copy_to_array_c3d_c3d_rs
274 MODULE PROCEDURE pw_copy_to_array_r1d_r1d_gs
275 MODULE PROCEDURE pw_copy_to_array_r1d_c1d_gs
276 MODULE PROCEDURE pw_copy_to_array_r3d_r3d_gs
277 MODULE PROCEDURE pw_copy_to_array_r3d_c3d_gs
278 MODULE PROCEDURE pw_copy_to_array_c1d_r1d_gs
279 MODULE PROCEDURE pw_copy_to_array_c1d_c1d_gs
280 MODULE PROCEDURE pw_copy_to_array_c3d_r3d_gs
281 MODULE PROCEDURE pw_copy_to_array_c3d_c3d_gs
285 MODULE PROCEDURE pw_copy_from_array_r1d_r1d_rs
286 MODULE PROCEDURE pw_copy_from_array_r1d_c1d_rs
287 MODULE PROCEDURE pw_copy_from_array_r3d_r3d_rs
288 MODULE PROCEDURE pw_copy_from_array_r3d_c3d_rs
289 MODULE PROCEDURE pw_copy_from_array_c1d_r1d_rs
290 MODULE PROCEDURE pw_copy_from_array_c1d_c1d_rs
291 MODULE PROCEDURE pw_copy_from_array_c3d_r3d_rs
292 MODULE PROCEDURE pw_copy_from_array_c3d_c3d_rs
293 MODULE PROCEDURE pw_copy_from_array_r1d_r1d_gs
294 MODULE PROCEDURE pw_copy_from_array_r1d_c1d_gs
295 MODULE PROCEDURE pw_copy_from_array_r3d_r3d_gs
296 MODULE PROCEDURE pw_copy_from_array_r3d_c3d_gs
297 MODULE PROCEDURE pw_copy_from_array_c1d_r1d_gs
298 MODULE PROCEDURE pw_copy_from_array_c1d_c1d_gs
299 MODULE PROCEDURE pw_copy_from_array_c3d_r3d_gs
300 MODULE PROCEDURE pw_copy_from_array_c3d_c3d_gs
304 MODULE PROCEDURE pw_copy_r1d_r1d_rs
305 MODULE PROCEDURE pw_copy_r1d_r1d_gs
306 MODULE PROCEDURE pw_gather_s_r1d_r3d_2
307 MODULE PROCEDURE pw_scatter_s_r1d_r3d_2
308 MODULE PROCEDURE pw_copy_r1d_c1d_rs
309 MODULE PROCEDURE pw_copy_r1d_c1d_gs
310 MODULE PROCEDURE pw_gather_s_r1d_c3d_2
311 MODULE PROCEDURE pw_scatter_s_r1d_c3d_2
312 MODULE PROCEDURE pw_copy_r3d_r3d_rs
313 MODULE PROCEDURE pw_copy_r3d_r3d_gs
314 MODULE PROCEDURE fft_wrap_pw1pw2_r3d_c1d_rs_gs
315 MODULE PROCEDURE pw_copy_r3d_c3d_rs
316 MODULE PROCEDURE pw_copy_r3d_c3d_gs
317 MODULE PROCEDURE fft_wrap_pw1pw2_r3d_c3d_rs_gs
318 MODULE PROCEDURE pw_copy_c1d_r1d_rs
319 MODULE PROCEDURE pw_copy_c1d_r1d_gs
320 MODULE PROCEDURE pw_gather_s_c1d_r3d_2
321 MODULE PROCEDURE pw_scatter_s_c1d_r3d_2
322 MODULE PROCEDURE fft_wrap_pw1pw2_c1d_r3d_gs_rs
323 MODULE PROCEDURE pw_copy_c1d_c1d_rs
324 MODULE PROCEDURE pw_copy_c1d_c1d_gs
325 MODULE PROCEDURE pw_gather_s_c1d_c3d_2
326 MODULE PROCEDURE pw_scatter_s_c1d_c3d_2
327 MODULE PROCEDURE fft_wrap_pw1pw2_c1d_c3d_gs_rs
328 MODULE PROCEDURE pw_copy_c3d_r3d_rs
329 MODULE PROCEDURE pw_copy_c3d_r3d_gs
330 MODULE PROCEDURE fft_wrap_pw1pw2_c3d_r3d_gs_rs
331 MODULE PROCEDURE fft_wrap_pw1pw2_c3d_c1d_rs_gs
332 MODULE PROCEDURE pw_copy_c3d_c3d_rs
333 MODULE PROCEDURE pw_copy_c3d_c3d_gs
334 MODULE PROCEDURE fft_wrap_pw1pw2_c3d_c3d_rs_gs
335 MODULE PROCEDURE fft_wrap_pw1pw2_c3d_c3d_gs_rs
346 SUBROUTINE pw_zero_r1d_rs (pw)
350 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_zero'
354 CALL timeset(routinen, handle)
356#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
364 CALL timestop(handle)
366 END SUBROUTINE pw_zero_r1d_rs
375 SUBROUTINE pw_scale_r1d_rs (pw, a)
378 REAL(KIND=
dp),
INTENT(IN) :: a
380 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_scale'
384 CALL timeset(routinen, handle)
387 pw%array = a*pw%array
390 CALL timestop(handle)
392 END SUBROUTINE pw_scale_r1d_rs
404 SUBROUTINE pw_write_r1d_rs (pw, unit_nr)
407 INTEGER,
INTENT(in) :: unit_nr
411 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
413 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=r1d"
414 IF (
ASSOCIATED(pw%array))
THEN
415 WRITE (unit=unit_nr, fmt=
"(' array=<r(',i8,':',i8,')>')") &
416 lbound(pw%array, 1), ubound(pw%array, 1)
418 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
421 END SUBROUTINE pw_write_r1d_rs
430 FUNCTION pw_integrate_function_r1d_rs (fun, isign, oprt)
RESULT(total_fun)
433 INTEGER,
INTENT(IN),
OPTIONAL :: isign
434 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
435 REAL(kind=
dp) :: total_fun
441 IF (
PRESENT(oprt))
THEN
446 cpabort(
"Unknown operator")
454 total_fun = fun%pw_grid%dvol*
accurate_sum(abs(fun%array))
460 CALL fun%pw_grid%para%group%sum(total_fun)
463 IF (
PRESENT(isign))
THEN
464 total_fun = total_fun*sign(1._dp, real(isign,
dp))
467 END FUNCTION pw_integrate_function_r1d_rs
474 SUBROUTINE pw_set_value_r1d_rs (pw, value)
476 REAL(KIND=
dp),
INTENT(IN) ::
value
478 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_set_value'
482 CALL timeset(routinen, handle)
488 CALL timestop(handle)
490 END SUBROUTINE pw_set_value_r1d_rs
498 SUBROUTINE pw_zero_r1d_gs (pw)
502 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_zero'
506 CALL timeset(routinen, handle)
508#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
516 CALL timestop(handle)
518 END SUBROUTINE pw_zero_r1d_gs
527 SUBROUTINE pw_scale_r1d_gs (pw, a)
530 REAL(KIND=
dp),
INTENT(IN) :: a
532 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_scale'
536 CALL timeset(routinen, handle)
539 pw%array = a*pw%array
542 CALL timestop(handle)
544 END SUBROUTINE pw_scale_r1d_gs
556 SUBROUTINE pw_write_r1d_gs (pw, unit_nr)
559 INTEGER,
INTENT(in) :: unit_nr
563 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
565 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=r1d"
566 IF (
ASSOCIATED(pw%array))
THEN
567 WRITE (unit=unit_nr, fmt=
"(' array=<r(',i8,':',i8,')>')") &
568 lbound(pw%array, 1), ubound(pw%array, 1)
570 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
573 END SUBROUTINE pw_write_r1d_gs
582 FUNCTION pw_integrate_function_r1d_gs (fun, isign, oprt)
RESULT(total_fun)
585 INTEGER,
INTENT(IN),
OPTIONAL :: isign
586 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
587 REAL(kind=
dp) :: total_fun
593 IF (
PRESENT(oprt))
THEN
598 cpabort(
"Unknown operator")
605 cpabort(
"Operator ABS not implemented")
606 IF (fun%pw_grid%have_g0) total_fun = fun%pw_grid%vol* fun%array(1)
609 CALL fun%pw_grid%para%group%sum(total_fun)
612 IF (
PRESENT(isign))
THEN
613 total_fun = total_fun*sign(1._dp, real(isign,
dp))
616 END FUNCTION pw_integrate_function_r1d_gs
623 SUBROUTINE pw_set_value_r1d_gs (pw, value)
625 REAL(KIND=
dp),
INTENT(IN) ::
value
627 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_set_value'
631 CALL timeset(routinen, handle)
637 CALL timestop(handle)
639 END SUBROUTINE pw_set_value_r1d_gs
647 SUBROUTINE pw_gather_p_r1d (pw, c, scale)
650 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
INTENT(IN) :: c
651 REAL(KIND=
dp),
INTENT(IN),
OPTIONAL :: scale
653 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_gather_p'
655 INTEGER :: gpt, handle, l, m, mn, n
657 CALL timeset(routinen, handle)
660 cpabort(
"This grid type is not distributed")
663 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
664 ngpts =>
SIZE(pw%pw_grid%gsq), ghat => pw%pw_grid%g_hat, yzq => pw%pw_grid%para%yzq)
666 IF (
PRESENT(scale))
THEN
671 l = mapl(ghat(1, gpt)) + 1
672 m = mapm(ghat(2, gpt)) + 1
673 n = mapn(ghat(3, gpt)) + 1
675 pw%array(gpt) = scale* real(c(l, mn), kind=
dp)
683 l = mapl(ghat(1, gpt)) + 1
684 m = mapm(ghat(2, gpt)) + 1
685 n = mapn(ghat(3, gpt)) + 1
687 pw%array(gpt) = real(c(l, mn), kind=
dp)
694 CALL timestop(handle)
696 END SUBROUTINE pw_gather_p_r1d
704 SUBROUTINE pw_scatter_p_r1d (pw, c, scale)
705 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw
706 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
INTENT(INOUT) :: c
707 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
709 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_scatter_p'
711 INTEGER :: gpt, handle, l, m, mn, n
713 CALL timeset(routinen, handle)
715 IF (pw%pw_grid%para%mode /= pw_mode_distributed)
THEN
716 cpabort(
"This grid type is not distributed")
719 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
720 ghat => pw%pw_grid%g_hat, yzq => pw%pw_grid%para%yzq, ngpts =>
SIZE(pw%pw_grid%gsq))
722 IF (.NOT.
PRESENT(scale)) c = z_zero
724 IF (
PRESENT(scale))
THEN
729 l = mapl(ghat(1, gpt)) + 1
730 m = mapm(ghat(2, gpt)) + 1
731 n = mapn(ghat(3, gpt)) + 1
733 c(l, mn) = cmplx(scale*pw%array(gpt), 0.0_dp, kind=dp)
741 l = mapl(ghat(1, gpt)) + 1
742 m = mapm(ghat(2, gpt)) + 1
743 n = mapn(ghat(3, gpt)) + 1
745 c(l, mn) = cmplx(pw%array(gpt), 0.0_dp, kind=dp)
752 IF (pw%pw_grid%grid_span == halfspace)
THEN
754 associate(mapm => pw%pw_grid%mapm%neg, mapn => pw%pw_grid%mapn%neg, mapl => pw%pw_grid%mapl%neg, &
755 ghat => pw%pw_grid%g_hat, ngpts =>
SIZE(pw%pw_grid%gsq), yzq => pw%pw_grid%para%yzq)
757 IF (
PRESENT(scale))
THEN
762 l = mapl(ghat(1, gpt)) + 1
763 m = mapm(ghat(2, gpt)) + 1
764 n = mapn(ghat(3, gpt)) + 1
766 c(l, mn) = scale*( cmplx(pw%array(gpt), 0.0_dp, kind=dp))
774 l = mapl(ghat(1, gpt)) + 1
775 m = mapm(ghat(2, gpt)) + 1
776 n = mapn(ghat(3, gpt)) + 1
778 c(l, mn) = ( cmplx(pw%array(gpt), 0.0_dp, kind=dp))
785 CALL timestop(handle)
787 END SUBROUTINE pw_scatter_p_r1d
795 SUBROUTINE pw_zero_r3d_rs (pw)
797 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw
799 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_zero'
803 CALL timeset(routinen, handle)
805#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
813 CALL timestop(handle)
815 END SUBROUTINE pw_zero_r3d_rs
824 SUBROUTINE pw_scale_r3d_rs (pw, a)
826 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw
827 REAL(KIND=dp),
INTENT(IN) :: a
829 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_scale'
833 CALL timeset(routinen, handle)
836 pw%array = a*pw%array
839 CALL timestop(handle)
841 END SUBROUTINE pw_scale_r3d_rs
853 SUBROUTINE pw_write_r3d_rs (pw, unit_nr)
855 TYPE(pw_r3d_rs_type),
INTENT(in) :: pw
856 INTEGER,
INTENT(in) :: unit_nr
860 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
862 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=r3d"
863 IF (
ASSOCIATED(pw%array))
THEN
864 WRITE (unit=unit_nr, fmt=
"(' array=<r(',i8,':',i8,',',i8,':',i8,',',i8,':',i8,')>')") &
865 lbound(pw%array, 1), ubound(pw%array, 1), lbound(pw%array, 2), ubound(pw%array, 2), &
866 lbound(pw%array, 3), ubound(pw%array, 3)
868 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
871 END SUBROUTINE pw_write_r3d_rs
880 FUNCTION pw_integrate_function_r3d_rs (fun, isign, oprt)
RESULT(total_fun)
882 TYPE(pw_r3d_rs_type),
INTENT(IN) :: fun
883 INTEGER,
INTENT(IN),
OPTIONAL :: isign
884 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
885 REAL(kind=dp) :: total_fun
891 IF (
PRESENT(oprt))
THEN
896 cpabort(
"Unknown operator")
904 total_fun = fun%pw_grid%dvol*accurate_sum(abs(fun%array))
906 total_fun = fun%pw_grid%dvol*accurate_sum( fun%array)
909 IF (fun%pw_grid%para%mode /= pw_mode_local)
THEN
910 CALL fun%pw_grid%para%group%sum(total_fun)
913 IF (
PRESENT(isign))
THEN
914 total_fun = total_fun*sign(1._dp, real(isign, dp))
917 END FUNCTION pw_integrate_function_r3d_rs
924 SUBROUTINE pw_set_value_r3d_rs (pw, value)
925 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw
926 REAL(KIND=dp),
INTENT(IN) ::
value
928 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_set_value'
932 CALL timeset(routinen, handle)
938 CALL timestop(handle)
940 END SUBROUTINE pw_set_value_r3d_rs
948 SUBROUTINE pw_zero_r3d_gs (pw)
950 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw
952 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_zero'
956 CALL timeset(routinen, handle)
958#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
966 CALL timestop(handle)
968 END SUBROUTINE pw_zero_r3d_gs
977 SUBROUTINE pw_scale_r3d_gs (pw, a)
979 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw
980 REAL(KIND=dp),
INTENT(IN) :: a
982 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_scale'
986 CALL timeset(routinen, handle)
989 pw%array = a*pw%array
992 CALL timestop(handle)
994 END SUBROUTINE pw_scale_r3d_gs
1006 SUBROUTINE pw_write_r3d_gs (pw, unit_nr)
1008 TYPE(pw_r3d_gs_type),
INTENT(in) :: pw
1009 INTEGER,
INTENT(in) :: unit_nr
1013 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
1015 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=r3d"
1016 IF (
ASSOCIATED(pw%array))
THEN
1017 WRITE (unit=unit_nr, fmt=
"(' array=<r(',i8,':',i8,',',i8,':',i8,',',i8,':',i8,')>')") &
1018 lbound(pw%array, 1), ubound(pw%array, 1), lbound(pw%array, 2), ubound(pw%array, 2), &
1019 lbound(pw%array, 3), ubound(pw%array, 3)
1021 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
1024 END SUBROUTINE pw_write_r3d_gs
1033 FUNCTION pw_integrate_function_r3d_gs (fun, isign, oprt)
RESULT(total_fun)
1035 TYPE(pw_r3d_gs_type),
INTENT(IN) :: fun
1036 INTEGER,
INTENT(IN),
OPTIONAL :: isign
1037 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
1038 REAL(kind=dp) :: total_fun
1044 IF (
PRESENT(oprt))
THEN
1049 cpabort(
"Unknown operator")
1056 cpabort(
"Operator ABS not implemented")
1057 cpabort(
"Reciprocal space integration for 3D grids not implemented!")
1059 IF (fun%pw_grid%para%mode /= pw_mode_local)
THEN
1060 CALL fun%pw_grid%para%group%sum(total_fun)
1063 IF (
PRESENT(isign))
THEN
1064 total_fun = total_fun*sign(1._dp, real(isign, dp))
1067 END FUNCTION pw_integrate_function_r3d_gs
1074 SUBROUTINE pw_set_value_r3d_gs (pw, value)
1075 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw
1076 REAL(KIND=dp),
INTENT(IN) ::
value
1078 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_set_value'
1082 CALL timeset(routinen, handle)
1088 CALL timestop(handle)
1090 END SUBROUTINE pw_set_value_r3d_gs
1099 SUBROUTINE pw_zero_c1d_rs (pw)
1101 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw
1103 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_zero'
1107 CALL timeset(routinen, handle)
1109#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
1111 pw%array = cmplx(0.0_dp, 0.0_dp, kind=dp)
1114 pw%array = cmplx(0.0_dp, 0.0_dp, kind=dp)
1117 CALL timestop(handle)
1119 END SUBROUTINE pw_zero_c1d_rs
1128 SUBROUTINE pw_scale_c1d_rs (pw, a)
1130 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw
1131 REAL(KIND=dp),
INTENT(IN) :: a
1133 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_scale'
1137 CALL timeset(routinen, handle)
1140 pw%array = a*pw%array
1143 CALL timestop(handle)
1145 END SUBROUTINE pw_scale_c1d_rs
1157 SUBROUTINE pw_write_c1d_rs (pw, unit_nr)
1159 TYPE(pw_c1d_rs_type),
INTENT(in) :: pw
1160 INTEGER,
INTENT(in) :: unit_nr
1164 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
1166 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=c1d"
1167 IF (
ASSOCIATED(pw%array))
THEN
1168 WRITE (unit=unit_nr, fmt=
"(' array=<c(',i8,':',i8,')>')") &
1169 lbound(pw%array, 1), ubound(pw%array, 1)
1171 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
1174 END SUBROUTINE pw_write_c1d_rs
1183 FUNCTION pw_integrate_function_c1d_rs (fun, isign, oprt)
RESULT(total_fun)
1185 TYPE(pw_c1d_rs_type),
INTENT(IN) :: fun
1186 INTEGER,
INTENT(IN),
OPTIONAL :: isign
1187 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
1188 REAL(kind=dp) :: total_fun
1194 IF (
PRESENT(oprt))
THEN
1199 cpabort(
"Unknown operator")
1207 total_fun = fun%pw_grid%dvol*accurate_sum(abs(fun%array))
1209 total_fun = fun%pw_grid%dvol*accurate_sum( real(fun%array, kind=dp))
1212 IF (fun%pw_grid%para%mode /= pw_mode_local)
THEN
1213 CALL fun%pw_grid%para%group%sum(total_fun)
1216 IF (
PRESENT(isign))
THEN
1217 total_fun = total_fun*sign(1._dp, real(isign, dp))
1220 END FUNCTION pw_integrate_function_c1d_rs
1227 SUBROUTINE pw_set_value_c1d_rs (pw, value)
1228 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw
1229 REAL(KIND=dp),
INTENT(IN) ::
value
1231 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_set_value'
1235 CALL timeset(routinen, handle)
1238 pw%array = cmplx(
value, 0.0_dp, kind=dp)
1241 CALL timestop(handle)
1243 END SUBROUTINE pw_set_value_c1d_rs
1251 SUBROUTINE pw_zero_c1d_gs (pw)
1253 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
1255 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_zero'
1259 CALL timeset(routinen, handle)
1261#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
1263 pw%array = cmplx(0.0_dp, 0.0_dp, kind=dp)
1266 pw%array = cmplx(0.0_dp, 0.0_dp, kind=dp)
1269 CALL timestop(handle)
1271 END SUBROUTINE pw_zero_c1d_gs
1280 SUBROUTINE pw_scale_c1d_gs (pw, a)
1282 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
1283 REAL(KIND=dp),
INTENT(IN) :: a
1285 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_scale'
1289 CALL timeset(routinen, handle)
1292 pw%array = a*pw%array
1295 CALL timestop(handle)
1297 END SUBROUTINE pw_scale_c1d_gs
1309 SUBROUTINE pw_write_c1d_gs (pw, unit_nr)
1311 TYPE(pw_c1d_gs_type),
INTENT(in) :: pw
1312 INTEGER,
INTENT(in) :: unit_nr
1316 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
1318 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=c1d"
1319 IF (
ASSOCIATED(pw%array))
THEN
1320 WRITE (unit=unit_nr, fmt=
"(' array=<c(',i8,':',i8,')>')") &
1321 lbound(pw%array, 1), ubound(pw%array, 1)
1323 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
1326 END SUBROUTINE pw_write_c1d_gs
1335 FUNCTION pw_integrate_function_c1d_gs (fun, isign, oprt)
RESULT(total_fun)
1337 TYPE(pw_c1d_gs_type),
INTENT(IN) :: fun
1338 INTEGER,
INTENT(IN),
OPTIONAL :: isign
1339 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
1340 REAL(kind=dp) :: total_fun
1346 IF (
PRESENT(oprt))
THEN
1351 cpabort(
"Unknown operator")
1358 cpabort(
"Operator ABS not implemented")
1359 IF (fun%pw_grid%have_g0) total_fun = fun%pw_grid%vol* real(fun%array(1), kind=dp)
1361 IF (fun%pw_grid%para%mode /= pw_mode_local)
THEN
1362 CALL fun%pw_grid%para%group%sum(total_fun)
1365 IF (
PRESENT(isign))
THEN
1366 total_fun = total_fun*sign(1._dp, real(isign, dp))
1369 END FUNCTION pw_integrate_function_c1d_gs
1376 SUBROUTINE pw_set_value_c1d_gs (pw, value)
1377 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
1378 REAL(KIND=dp),
INTENT(IN) ::
value
1380 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_set_value'
1384 CALL timeset(routinen, handle)
1387 pw%array = cmplx(
value, 0.0_dp, kind=dp)
1390 CALL timestop(handle)
1392 END SUBROUTINE pw_set_value_c1d_gs
1400 SUBROUTINE pw_gather_p_c1d (pw, c, scale)
1402 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
1403 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
INTENT(IN) :: c
1404 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
1406 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_gather_p'
1408 INTEGER :: gpt, handle, l, m, mn, n
1410 CALL timeset(routinen, handle)
1412 IF (pw%pw_grid%para%mode /= pw_mode_distributed)
THEN
1413 cpabort(
"This grid type is not distributed")
1416 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
1417 ngpts =>
SIZE(pw%pw_grid%gsq), ghat => pw%pw_grid%g_hat, yzq => pw%pw_grid%para%yzq)
1419 IF (
PRESENT(scale))
THEN
1424 l = mapl(ghat(1, gpt)) + 1
1425 m = mapm(ghat(2, gpt)) + 1
1426 n = mapn(ghat(3, gpt)) + 1
1428 pw%array(gpt) = scale* c(l, mn)
1436 l = mapl(ghat(1, gpt)) + 1
1437 m = mapm(ghat(2, gpt)) + 1
1438 n = mapn(ghat(3, gpt)) + 1
1440 pw%array(gpt) = c(l, mn)
1447 CALL timestop(handle)
1449 END SUBROUTINE pw_gather_p_c1d
1457 SUBROUTINE pw_scatter_p_c1d (pw, c, scale)
1458 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
1459 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
INTENT(INOUT) :: c
1460 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
1462 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_scatter_p'
1464 INTEGER :: gpt, handle, l, m, mn, n
1466 CALL timeset(routinen, handle)
1468 IF (pw%pw_grid%para%mode /= pw_mode_distributed)
THEN
1469 cpabort(
"This grid type is not distributed")
1472 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
1473 ghat => pw%pw_grid%g_hat, yzq => pw%pw_grid%para%yzq, ngpts =>
SIZE(pw%pw_grid%gsq))
1475 IF (.NOT.
PRESENT(scale)) c = z_zero
1477 IF (
PRESENT(scale))
THEN
1482 l = mapl(ghat(1, gpt)) + 1
1483 m = mapm(ghat(2, gpt)) + 1
1484 n = mapn(ghat(3, gpt)) + 1
1486 c(l, mn) = scale*pw%array(gpt)
1494 l = mapl(ghat(1, gpt)) + 1
1495 m = mapm(ghat(2, gpt)) + 1
1496 n = mapn(ghat(3, gpt)) + 1
1498 c(l, mn) = pw%array(gpt)
1505 IF (pw%pw_grid%grid_span == halfspace)
THEN
1507 associate(mapm => pw%pw_grid%mapm%neg, mapn => pw%pw_grid%mapn%neg, mapl => pw%pw_grid%mapl%neg, &
1508 ghat => pw%pw_grid%g_hat, ngpts =>
SIZE(pw%pw_grid%gsq), yzq => pw%pw_grid%para%yzq)
1510 IF (
PRESENT(scale))
THEN
1515 l = mapl(ghat(1, gpt)) + 1
1516 m = mapm(ghat(2, gpt)) + 1
1517 n = mapn(ghat(3, gpt)) + 1
1519 c(l, mn) = scale*conjg( pw%array(gpt))
1527 l = mapl(ghat(1, gpt)) + 1
1528 m = mapm(ghat(2, gpt)) + 1
1529 n = mapn(ghat(3, gpt)) + 1
1531 c(l, mn) = conjg( pw%array(gpt))
1538 CALL timestop(handle)
1540 END SUBROUTINE pw_scatter_p_c1d
1548 SUBROUTINE pw_zero_c3d_rs (pw)
1550 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw
1552 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_zero'
1556 CALL timeset(routinen, handle)
1558#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
1560 pw%array = cmplx(0.0_dp, 0.0_dp, kind=dp)
1563 pw%array = cmplx(0.0_dp, 0.0_dp, kind=dp)
1566 CALL timestop(handle)
1568 END SUBROUTINE pw_zero_c3d_rs
1577 SUBROUTINE pw_scale_c3d_rs (pw, a)
1579 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw
1580 REAL(KIND=dp),
INTENT(IN) :: a
1582 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_scale'
1586 CALL timeset(routinen, handle)
1589 pw%array = a*pw%array
1592 CALL timestop(handle)
1594 END SUBROUTINE pw_scale_c3d_rs
1606 SUBROUTINE pw_write_c3d_rs (pw, unit_nr)
1608 TYPE(pw_c3d_rs_type),
INTENT(in) :: pw
1609 INTEGER,
INTENT(in) :: unit_nr
1613 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
1615 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=c3d"
1616 IF (
ASSOCIATED(pw%array))
THEN
1617 WRITE (unit=unit_nr, fmt=
"(' array=<c(',i8,':',i8,',',i8,':',i8,',',i8,':',i8,')>')") &
1618 lbound(pw%array, 1), ubound(pw%array, 1), lbound(pw%array, 2), ubound(pw%array, 2), &
1619 lbound(pw%array, 3), ubound(pw%array, 3)
1621 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
1624 END SUBROUTINE pw_write_c3d_rs
1633 FUNCTION pw_integrate_function_c3d_rs (fun, isign, oprt)
RESULT(total_fun)
1635 TYPE(pw_c3d_rs_type),
INTENT(IN) :: fun
1636 INTEGER,
INTENT(IN),
OPTIONAL :: isign
1637 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
1638 REAL(kind=dp) :: total_fun
1644 IF (
PRESENT(oprt))
THEN
1649 cpabort(
"Unknown operator")
1657 total_fun = fun%pw_grid%dvol*accurate_sum(abs(fun%array))
1659 total_fun = fun%pw_grid%dvol*accurate_sum( real(fun%array, kind=dp))
1662 IF (fun%pw_grid%para%mode /= pw_mode_local)
THEN
1663 CALL fun%pw_grid%para%group%sum(total_fun)
1666 IF (
PRESENT(isign))
THEN
1667 total_fun = total_fun*sign(1._dp, real(isign, dp))
1670 END FUNCTION pw_integrate_function_c3d_rs
1677 SUBROUTINE pw_set_value_c3d_rs (pw, value)
1678 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw
1679 REAL(KIND=dp),
INTENT(IN) ::
value
1681 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_set_value'
1685 CALL timeset(routinen, handle)
1688 pw%array = cmplx(
value, 0.0_dp, kind=dp)
1691 CALL timestop(handle)
1693 END SUBROUTINE pw_set_value_c3d_rs
1701 SUBROUTINE pw_zero_c3d_gs (pw)
1703 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw
1705 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_zero'
1709 CALL timeset(routinen, handle)
1711#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
1713 pw%array = cmplx(0.0_dp, 0.0_dp, kind=dp)
1716 pw%array = cmplx(0.0_dp, 0.0_dp, kind=dp)
1719 CALL timestop(handle)
1721 END SUBROUTINE pw_zero_c3d_gs
1730 SUBROUTINE pw_scale_c3d_gs (pw, a)
1732 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw
1733 REAL(KIND=dp),
INTENT(IN) :: a
1735 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_scale'
1739 CALL timeset(routinen, handle)
1742 pw%array = a*pw%array
1745 CALL timestop(handle)
1747 END SUBROUTINE pw_scale_c3d_gs
1759 SUBROUTINE pw_write_c3d_gs (pw, unit_nr)
1761 TYPE(pw_c3d_gs_type),
INTENT(in) :: pw
1762 INTEGER,
INTENT(in) :: unit_nr
1766 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
1768 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=c3d"
1769 IF (
ASSOCIATED(pw%array))
THEN
1770 WRITE (unit=unit_nr, fmt=
"(' array=<c(',i8,':',i8,',',i8,':',i8,',',i8,':',i8,')>')") &
1771 lbound(pw%array, 1), ubound(pw%array, 1), lbound(pw%array, 2), ubound(pw%array, 2), &
1772 lbound(pw%array, 3), ubound(pw%array, 3)
1774 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
1777 END SUBROUTINE pw_write_c3d_gs
1786 FUNCTION pw_integrate_function_c3d_gs (fun, isign, oprt)
RESULT(total_fun)
1788 TYPE(pw_c3d_gs_type),
INTENT(IN) :: fun
1789 INTEGER,
INTENT(IN),
OPTIONAL :: isign
1790 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
1791 REAL(kind=dp) :: total_fun
1797 IF (
PRESENT(oprt))
THEN
1802 cpabort(
"Unknown operator")
1809 cpabort(
"Operator ABS not implemented")
1810 cpabort(
"Reciprocal space integration for 3D grids not implemented!")
1812 IF (fun%pw_grid%para%mode /= pw_mode_local)
THEN
1813 CALL fun%pw_grid%para%group%sum(total_fun)
1816 IF (
PRESENT(isign))
THEN
1817 total_fun = total_fun*sign(1._dp, real(isign, dp))
1820 END FUNCTION pw_integrate_function_c3d_gs
1827 SUBROUTINE pw_set_value_c3d_gs (pw, value)
1828 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw
1829 REAL(KIND=dp),
INTENT(IN) ::
value
1831 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_set_value'
1835 CALL timeset(routinen, handle)
1838 pw%array = cmplx(
value, 0.0_dp, kind=dp)
1841 CALL timestop(handle)
1843 END SUBROUTINE pw_set_value_c3d_gs
1859 SUBROUTINE pw_copy_r1d_r1d_rs (pw1, pw2)
1861 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
1862 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw2
1864 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
1867 INTEGER :: i, j, ng, ng1, ng2, ns
1869 CALL timeset(routinen, handle)
1871 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
1872 cpabort(
"Both grids must be either spherical or non-spherical!")
1873 IF (pw1%pw_grid%spherical) &
1874 cpabort(
"Spherical grids only exist in reciprocal space!")
1876 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
1877 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
1878 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
1879 ng1 =
SIZE(pw1%array)
1880 ng2 =
SIZE(pw2%array)
1883 pw2%array(1:ng) = pw1%array(1:ng)
1887 pw2%array(ng + 1:ng2) = 0.0_dp
1891 cpabort(
"Copies between spherical grids require compatible grids!")
1894 ng1 =
SIZE(pw1%array)
1895 ng2 =
SIZE(pw2%array)
1896 ns = 2*max(ng1, ng2)
1898 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
1899 IF (ng1 >= ng2)
THEN
1904 j = pw2%pw_grid%gidx(i)
1905 pw2%array(i) = pw1%array(j)
1914 j = pw2%pw_grid%gidx(i)
1915 pw2%array(j) = pw1%array(i)
1919 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
1920 IF (ng1 >= ng2)
THEN
1923 j = pw1%pw_grid%gidx(i)
1924 pw2%array(i) = pw1%array(j)
1931 j = pw1%pw_grid%gidx(i)
1932 pw2%array(j) = pw1%array(i)
1937 cpabort(
"Copy not implemented!")
1944 pw2%array = pw1%array
1948 CALL timestop(handle)
1950 END SUBROUTINE pw_copy_r1d_r1d_rs
1957 SUBROUTINE pw_copy_to_array_r1d_r1d_rs (pw, array)
1958 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw
1959 REAL(KIND=dp),
DIMENSION(:),
INTENT(INOUT) :: array
1961 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
1965 CALL timeset(routinen, handle)
1968 array(:) = pw%array(:)
1971 CALL timestop(handle)
1972 END SUBROUTINE pw_copy_to_array_r1d_r1d_rs
1979 SUBROUTINE pw_copy_from_array_r1d_r1d_rs (pw, array)
1980 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw
1981 REAL(KIND=dp),
DIMENSION(:),
INTENT(IN) :: array
1983 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
1987 CALL timeset(routinen, handle)
1993 CALL timestop(handle)
1994 END SUBROUTINE pw_copy_from_array_r1d_r1d_rs
2012 SUBROUTINE pw_axpy_r1d_r1d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
2014 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
2015 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw2
2016 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
2017 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
2019 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
2022 LOGICAL :: my_allow_noncompatible_grids
2023 REAL(KIND=dp) :: my_alpha, my_beta
2024 INTEGER :: i, j, ng, ng1, ng2
2026 CALL timeset(routinen, handle)
2029 IF (
PRESENT(alpha)) my_alpha = alpha
2032 IF (
PRESENT(beta)) my_beta = beta
2034 my_allow_noncompatible_grids = .false.
2035 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
2037 IF (my_beta /= 1.0_dp)
THEN
2038 IF (my_beta == 0.0_dp)
THEN
2042 pw2%array = pw2%array*my_beta
2047 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
2049 IF (my_alpha == 1.0_dp)
THEN
2051 pw2%array = pw2%array + pw1%array
2053 ELSE IF (my_alpha /= 0.0_dp)
THEN
2055 pw2%array = pw2%array + my_alpha* pw1%array
2059 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
2061 ng1 =
SIZE(pw1%array)
2062 ng2 =
SIZE(pw2%array)
2065 IF (pw1%pw_grid%spherical)
THEN
2066 IF (my_alpha == 1.0_dp)
THEN
2069 pw2%array(i) = pw2%array(i) + pw1%array(i)
2072 ELSE IF (my_alpha /= 0.0_dp)
THEN
2075 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(i)
2079 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
2080 IF (ng1 >= ng2)
THEN
2081 IF (my_alpha == 1.0_dp)
THEN
2084 j = pw2%pw_grid%gidx(i)
2085 pw2%array(i) = pw2%array(i) + pw1%array(j)
2088 ELSE IF (my_alpha /= 0.0_dp)
THEN
2091 j = pw2%pw_grid%gidx(i)
2092 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
2097 IF (my_alpha == 1.0_dp)
THEN
2100 j = pw2%pw_grid%gidx(i)
2101 pw2%array(j) = pw2%array(j) + pw1%array(i)
2104 ELSE IF (my_alpha /= 0.0_dp)
THEN
2107 j = pw2%pw_grid%gidx(i)
2108 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
2113 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
2114 IF (ng1 >= ng2)
THEN
2115 IF (my_alpha == 1.0_dp)
THEN
2118 j = pw1%pw_grid%gidx(i)
2119 pw2%array(i) = pw2%array(i) + pw1%array(j)
2122 ELSE IF (my_alpha /= 0.0_dp)
THEN
2125 j = pw1%pw_grid%gidx(i)
2126 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
2131 IF (my_alpha == 1.0_dp)
THEN
2134 j = pw1%pw_grid%gidx(i)
2135 pw2%array(j) = pw2%array(j) + pw1%array(i)
2138 ELSE IF (my_alpha /= 0.0_dp)
THEN
2141 j = pw1%pw_grid%gidx(i)
2142 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
2148 cpabort(
"Grids not compatible")
2153 cpabort(
"Grids not compatible")
2157 CALL timestop(handle)
2159 END SUBROUTINE pw_axpy_r1d_r1d_rs
2170 SUBROUTINE pw_multiply_r1d_r1d_rs (pw_out, pw1, pw2, alpha)
2172 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw_out
2173 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
2174 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw2
2175 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
2177 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
2180 REAL(KIND=dp) :: my_alpha
2182 CALL timeset(routinen, handle)
2185 IF (
PRESENT(alpha)) my_alpha = alpha
2187 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
2188 cpabort(
"pw_multiply not implemented for non-identical grids!")
2190#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
2191 IF (my_alpha == 1.0_dp)
THEN
2193 pw_out%array = pw_out%array + pw1%array* pw2%array
2197 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
2201 IF (my_alpha == 1.0_dp)
THEN
2202 pw_out%array = pw_out%array + pw1%array* pw2%array
2204 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
2208 CALL timestop(handle)
2210 END SUBROUTINE pw_multiply_r1d_r1d_rs
2217 SUBROUTINE pw_multiply_with_r1d_r1d_rs (pw1, pw2)
2218 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw1
2219 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw2
2221 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
2225 CALL timeset(routinen, handle)
2227 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
2228 cpabort(
"Incompatible grids!")
2231 pw1%array = pw1%array* pw2%array
2234 CALL timestop(handle)
2236 END SUBROUTINE pw_multiply_with_r1d_r1d_rs
2251 FUNCTION pw_integral_ab_r1d_r1d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
2253 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
2254 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw2
2255 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
2256 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
2257 REAL(kind=dp) :: integral_value
2259 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_r1d_r1d_rs'
2261 INTEGER :: handle, loc_sumtype
2262 LOGICAL :: my_just_sum, my_local_only
2264 CALL timeset(routinen, handle)
2267 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
2269 my_local_only = .false.
2270 IF (
PRESENT(local_only)) my_local_only = local_only
2272 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
2273 cpabort(
"Grids incompatible")
2276 my_just_sum = .false.
2277 IF (
PRESENT(just_sum)) my_just_sum = just_sum
2284 integral_value = dot_product(pw1%array, pw2%array)
2289 integral_value = accurate_dot_product(pw1%array, pw2%array)
2293 IF (.NOT. my_just_sum)
THEN
2294 integral_value = integral_value*pw1%pw_grid%dvol
2297 IF (pw1%pw_grid%grid_span == halfspace)
THEN
2298 integral_value = 2.0_dp*integral_value
2299 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
2300 pw1%array(1)*pw2%array(1)
2303 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
2304 CALL pw1%pw_grid%para%group%sum(integral_value)
2306 CALL timestop(handle)
2308 END FUNCTION pw_integral_ab_r1d_r1d_rs
2322 SUBROUTINE pw_copy_r1d_r1d_gs (pw1, pw2)
2324 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
2325 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw2
2327 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
2330 INTEGER :: i, j, ng, ng1, ng2, ns
2332 CALL timeset(routinen, handle)
2334 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
2335 cpabort(
"Both grids must be either spherical or non-spherical!")
2337 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
2338 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
2339 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
2340 ng1 =
SIZE(pw1%array)
2341 ng2 =
SIZE(pw2%array)
2344 pw2%array(1:ng) = pw1%array(1:ng)
2348 pw2%array(ng + 1:ng2) = 0.0_dp
2352 cpabort(
"Copies between spherical grids require compatible grids!")
2355 ng1 =
SIZE(pw1%array)
2356 ng2 =
SIZE(pw2%array)
2357 ns = 2*max(ng1, ng2)
2359 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
2360 IF (ng1 >= ng2)
THEN
2365 j = pw2%pw_grid%gidx(i)
2366 pw2%array(i) = pw1%array(j)
2375 j = pw2%pw_grid%gidx(i)
2376 pw2%array(j) = pw1%array(i)
2380 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
2381 IF (ng1 >= ng2)
THEN
2384 j = pw1%pw_grid%gidx(i)
2385 pw2%array(i) = pw1%array(j)
2392 j = pw1%pw_grid%gidx(i)
2393 pw2%array(j) = pw1%array(i)
2398 cpabort(
"Copy not implemented!")
2405 pw2%array = pw1%array
2409 CALL timestop(handle)
2411 END SUBROUTINE pw_copy_r1d_r1d_gs
2418 SUBROUTINE pw_copy_to_array_r1d_r1d_gs (pw, array)
2419 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw
2420 REAL(KIND=dp),
DIMENSION(:),
INTENT(INOUT) :: array
2422 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
2426 CALL timeset(routinen, handle)
2429 array(:) = pw%array(:)
2432 CALL timestop(handle)
2433 END SUBROUTINE pw_copy_to_array_r1d_r1d_gs
2440 SUBROUTINE pw_copy_from_array_r1d_r1d_gs (pw, array)
2441 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw
2442 REAL(KIND=dp),
DIMENSION(:),
INTENT(IN) :: array
2444 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
2448 CALL timeset(routinen, handle)
2454 CALL timestop(handle)
2455 END SUBROUTINE pw_copy_from_array_r1d_r1d_gs
2473 SUBROUTINE pw_axpy_r1d_r1d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
2475 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
2476 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw2
2477 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
2478 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
2480 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
2483 LOGICAL :: my_allow_noncompatible_grids
2484 REAL(KIND=dp) :: my_alpha, my_beta
2485 INTEGER :: i, j, ng, ng1, ng2
2487 CALL timeset(routinen, handle)
2490 IF (
PRESENT(alpha)) my_alpha = alpha
2493 IF (
PRESENT(beta)) my_beta = beta
2495 my_allow_noncompatible_grids = .false.
2496 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
2498 IF (my_beta /= 1.0_dp)
THEN
2499 IF (my_beta == 0.0_dp)
THEN
2503 pw2%array = pw2%array*my_beta
2508 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
2510 IF (my_alpha == 1.0_dp)
THEN
2512 pw2%array = pw2%array + pw1%array
2514 ELSE IF (my_alpha /= 0.0_dp)
THEN
2516 pw2%array = pw2%array + my_alpha* pw1%array
2520 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
2522 ng1 =
SIZE(pw1%array)
2523 ng2 =
SIZE(pw2%array)
2526 IF (pw1%pw_grid%spherical)
THEN
2527 IF (my_alpha == 1.0_dp)
THEN
2530 pw2%array(i) = pw2%array(i) + pw1%array(i)
2533 ELSE IF (my_alpha /= 0.0_dp)
THEN
2536 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(i)
2540 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
2541 IF (ng1 >= ng2)
THEN
2542 IF (my_alpha == 1.0_dp)
THEN
2545 j = pw2%pw_grid%gidx(i)
2546 pw2%array(i) = pw2%array(i) + pw1%array(j)
2549 ELSE IF (my_alpha /= 0.0_dp)
THEN
2552 j = pw2%pw_grid%gidx(i)
2553 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
2558 IF (my_alpha == 1.0_dp)
THEN
2561 j = pw2%pw_grid%gidx(i)
2562 pw2%array(j) = pw2%array(j) + pw1%array(i)
2565 ELSE IF (my_alpha /= 0.0_dp)
THEN
2568 j = pw2%pw_grid%gidx(i)
2569 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
2574 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
2575 IF (ng1 >= ng2)
THEN
2576 IF (my_alpha == 1.0_dp)
THEN
2579 j = pw1%pw_grid%gidx(i)
2580 pw2%array(i) = pw2%array(i) + pw1%array(j)
2583 ELSE IF (my_alpha /= 0.0_dp)
THEN
2586 j = pw1%pw_grid%gidx(i)
2587 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
2592 IF (my_alpha == 1.0_dp)
THEN
2595 j = pw1%pw_grid%gidx(i)
2596 pw2%array(j) = pw2%array(j) + pw1%array(i)
2599 ELSE IF (my_alpha /= 0.0_dp)
THEN
2602 j = pw1%pw_grid%gidx(i)
2603 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
2609 cpabort(
"Grids not compatible")
2614 cpabort(
"Grids not compatible")
2618 CALL timestop(handle)
2620 END SUBROUTINE pw_axpy_r1d_r1d_gs
2631 SUBROUTINE pw_multiply_r1d_r1d_gs (pw_out, pw1, pw2, alpha)
2633 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw_out
2634 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
2635 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
2636 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
2638 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
2641 REAL(KIND=dp) :: my_alpha
2643 CALL timeset(routinen, handle)
2646 IF (
PRESENT(alpha)) my_alpha = alpha
2648 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
2649 cpabort(
"pw_multiply not implemented for non-identical grids!")
2651#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
2652 IF (my_alpha == 1.0_dp)
THEN
2654 pw_out%array = pw_out%array + pw1%array* pw2%array
2658 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
2662 IF (my_alpha == 1.0_dp)
THEN
2663 pw_out%array = pw_out%array + pw1%array* pw2%array
2665 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
2669 CALL timestop(handle)
2671 END SUBROUTINE pw_multiply_r1d_r1d_gs
2678 SUBROUTINE pw_multiply_with_r1d_r1d_gs (pw1, pw2)
2679 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw1
2680 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
2682 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
2686 CALL timeset(routinen, handle)
2688 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
2689 cpabort(
"Incompatible grids!")
2692 pw1%array = pw1%array* pw2%array
2695 CALL timestop(handle)
2697 END SUBROUTINE pw_multiply_with_r1d_r1d_gs
2712 FUNCTION pw_integral_ab_r1d_r1d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
2714 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
2715 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
2716 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
2717 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
2718 REAL(kind=dp) :: integral_value
2720 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_r1d_r1d_gs'
2722 INTEGER :: handle, loc_sumtype
2723 LOGICAL :: my_just_sum, my_local_only
2725 CALL timeset(routinen, handle)
2728 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
2730 my_local_only = .false.
2731 IF (
PRESENT(local_only)) my_local_only = local_only
2733 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
2734 cpabort(
"Grids incompatible")
2737 my_just_sum = .false.
2738 IF (
PRESENT(just_sum)) my_just_sum = just_sum
2745 integral_value = dot_product(pw1%array, pw2%array)
2750 integral_value = accurate_dot_product(pw1%array, pw2%array)
2754 IF (.NOT. my_just_sum)
THEN
2755 integral_value = integral_value*pw1%pw_grid%vol
2758 IF (pw1%pw_grid%grid_span == halfspace)
THEN
2759 integral_value = 2.0_dp*integral_value
2760 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
2761 pw1%array(1)*pw2%array(1)
2764 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
2765 CALL pw1%pw_grid%para%group%sum(integral_value)
2767 CALL timestop(handle)
2769 END FUNCTION pw_integral_ab_r1d_r1d_gs
2777 FUNCTION pw_integral_a2b_r1d_r1d (pw1, pw2)
RESULT(integral_value)
2779 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
2780 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
2781 REAL(kind=dp) :: integral_value
2783 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_a2b'
2787 CALL timeset(routinen, handle)
2789 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
2790 cpabort(
"Grids incompatible")
2793 integral_value = accurate_sum(pw1%array*pw2%array*pw1%pw_grid%gsq)
2794 IF (pw1%pw_grid%grid_span == halfspace)
THEN
2795 integral_value = 2.0_dp*integral_value
2798 integral_value = integral_value*pw1%pw_grid%vol
2800 IF (pw1%pw_grid%para%mode == pw_mode_distributed) &
2801 CALL pw1%pw_grid%para%group%sum(integral_value)
2802 CALL timestop(handle)
2804 END FUNCTION pw_integral_a2b_r1d_r1d
2818 SUBROUTINE pw_copy_r1d_c1d_rs (pw1, pw2)
2820 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
2821 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw2
2823 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
2826 INTEGER :: i, j, ng, ng1, ng2, ns
2828 CALL timeset(routinen, handle)
2830 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
2831 cpabort(
"Both grids must be either spherical or non-spherical!")
2832 IF (pw1%pw_grid%spherical) &
2833 cpabort(
"Spherical grids only exist in reciprocal space!")
2835 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
2836 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
2837 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
2838 ng1 =
SIZE(pw1%array)
2839 ng2 =
SIZE(pw2%array)
2842 pw2%array(1:ng) = cmplx(pw1%array(1:ng), 0.0_dp, kind=dp)
2846 pw2%array(ng + 1:ng2) = cmplx(0.0_dp, 0.0_dp, kind=dp)
2850 cpabort(
"Copies between spherical grids require compatible grids!")
2853 ng1 =
SIZE(pw1%array)
2854 ng2 =
SIZE(pw2%array)
2855 ns = 2*max(ng1, ng2)
2857 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
2858 IF (ng1 >= ng2)
THEN
2863 j = pw2%pw_grid%gidx(i)
2864 pw2%array(i) = cmplx(pw1%array(j), 0.0_dp, kind=dp)
2873 j = pw2%pw_grid%gidx(i)
2874 pw2%array(j) = cmplx(pw1%array(i), 0.0_dp, kind=dp)
2878 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
2879 IF (ng1 >= ng2)
THEN
2882 j = pw1%pw_grid%gidx(i)
2883 pw2%array(i) = cmplx(pw1%array(j), 0.0_dp, kind=dp)
2890 j = pw1%pw_grid%gidx(i)
2891 pw2%array(j) = cmplx(pw1%array(i), 0.0_dp, kind=dp)
2896 cpabort(
"Copy not implemented!")
2903 pw2%array = cmplx(pw1%array, 0.0_dp, kind=dp)
2907 CALL timestop(handle)
2909 END SUBROUTINE pw_copy_r1d_c1d_rs
2916 SUBROUTINE pw_copy_to_array_r1d_c1d_rs (pw, array)
2917 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw
2918 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(INOUT) :: array
2920 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
2924 CALL timeset(routinen, handle)
2927 array(:) = cmplx(pw%array(:), 0.0_dp, kind=dp)
2930 CALL timestop(handle)
2931 END SUBROUTINE pw_copy_to_array_r1d_c1d_rs
2938 SUBROUTINE pw_copy_from_array_r1d_c1d_rs (pw, array)
2939 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw
2940 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(IN) :: array
2942 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
2946 CALL timeset(routinen, handle)
2949 pw%array = real(array, kind=dp)
2952 CALL timestop(handle)
2953 END SUBROUTINE pw_copy_from_array_r1d_c1d_rs
2971 SUBROUTINE pw_axpy_r1d_c1d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
2973 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
2974 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw2
2975 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
2976 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
2978 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
2981 LOGICAL :: my_allow_noncompatible_grids
2982 REAL(KIND=dp) :: my_alpha, my_beta
2983 INTEGER :: i, j, ng, ng1, ng2
2985 CALL timeset(routinen, handle)
2988 IF (
PRESENT(alpha)) my_alpha = alpha
2991 IF (
PRESENT(beta)) my_beta = beta
2993 my_allow_noncompatible_grids = .false.
2994 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
2996 IF (my_beta /= 1.0_dp)
THEN
2997 IF (my_beta == 0.0_dp)
THEN
3001 pw2%array = pw2%array*my_beta
3006 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
3008 IF (my_alpha == 1.0_dp)
THEN
3010 pw2%array = pw2%array + cmplx(pw1%array, 0.0_dp, kind=dp)
3012 ELSE IF (my_alpha /= 0.0_dp)
THEN
3014 pw2%array = pw2%array + my_alpha* cmplx(pw1%array, 0.0_dp, kind=dp)
3018 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
3020 ng1 =
SIZE(pw1%array)
3021 ng2 =
SIZE(pw2%array)
3024 IF (pw1%pw_grid%spherical)
THEN
3025 IF (my_alpha == 1.0_dp)
THEN
3028 pw2%array(i) = pw2%array(i) + cmplx(pw1%array(i), 0.0_dp, kind=dp)
3031 ELSE IF (my_alpha /= 0.0_dp)
THEN
3034 pw2%array(i) = pw2%array(i) + my_alpha* cmplx(pw1%array(i), 0.0_dp, kind=dp)
3038 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
3039 IF (ng1 >= ng2)
THEN
3040 IF (my_alpha == 1.0_dp)
THEN
3043 j = pw2%pw_grid%gidx(i)
3044 pw2%array(i) = pw2%array(i) + cmplx(pw1%array(j), 0.0_dp, kind=dp)
3047 ELSE IF (my_alpha /= 0.0_dp)
THEN
3050 j = pw2%pw_grid%gidx(i)
3051 pw2%array(i) = pw2%array(i) + my_alpha* cmplx(pw1%array(j), 0.0_dp, kind=dp)
3056 IF (my_alpha == 1.0_dp)
THEN
3059 j = pw2%pw_grid%gidx(i)
3060 pw2%array(j) = pw2%array(j) + cmplx(pw1%array(i), 0.0_dp, kind=dp)
3063 ELSE IF (my_alpha /= 0.0_dp)
THEN
3066 j = pw2%pw_grid%gidx(i)
3067 pw2%array(j) = pw2%array(j) + my_alpha* cmplx(pw1%array(i), 0.0_dp, kind=dp)
3072 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
3073 IF (ng1 >= ng2)
THEN
3074 IF (my_alpha == 1.0_dp)
THEN
3077 j = pw1%pw_grid%gidx(i)
3078 pw2%array(i) = pw2%array(i) + cmplx(pw1%array(j), 0.0_dp, kind=dp)
3081 ELSE IF (my_alpha /= 0.0_dp)
THEN
3084 j = pw1%pw_grid%gidx(i)
3085 pw2%array(i) = pw2%array(i) + my_alpha* cmplx(pw1%array(j), 0.0_dp, kind=dp)
3090 IF (my_alpha == 1.0_dp)
THEN
3093 j = pw1%pw_grid%gidx(i)
3094 pw2%array(j) = pw2%array(j) + cmplx(pw1%array(i), 0.0_dp, kind=dp)
3097 ELSE IF (my_alpha /= 0.0_dp)
THEN
3100 j = pw1%pw_grid%gidx(i)
3101 pw2%array(j) = pw2%array(j) + my_alpha* cmplx(pw1%array(i), 0.0_dp, kind=dp)
3107 cpabort(
"Grids not compatible")
3112 cpabort(
"Grids not compatible")
3116 CALL timestop(handle)
3118 END SUBROUTINE pw_axpy_r1d_c1d_rs
3129 SUBROUTINE pw_multiply_r1d_c1d_rs (pw_out, pw1, pw2, alpha)
3131 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw_out
3132 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
3133 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw2
3134 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
3136 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
3139 REAL(KIND=dp) :: my_alpha
3141 CALL timeset(routinen, handle)
3144 IF (
PRESENT(alpha)) my_alpha = alpha
3146 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
3147 cpabort(
"pw_multiply not implemented for non-identical grids!")
3149#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
3150 IF (my_alpha == 1.0_dp)
THEN
3152 pw_out%array = pw_out%array + pw1%array* real(pw2%array, kind=dp)
3156 pw_out%array = pw_out%array + my_alpha*pw1%array* real(pw2%array, kind=dp)
3160 IF (my_alpha == 1.0_dp)
THEN
3161 pw_out%array = pw_out%array + pw1%array* real(pw2%array, kind=dp)
3163 pw_out%array = pw_out%array + my_alpha*pw1%array* real(pw2%array, kind=dp)
3167 CALL timestop(handle)
3169 END SUBROUTINE pw_multiply_r1d_c1d_rs
3176 SUBROUTINE pw_multiply_with_r1d_c1d_rs (pw1, pw2)
3177 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw1
3178 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw2
3180 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
3184 CALL timeset(routinen, handle)
3186 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
3187 cpabort(
"Incompatible grids!")
3190 pw1%array = pw1%array* real(pw2%array, kind=dp)
3193 CALL timestop(handle)
3195 END SUBROUTINE pw_multiply_with_r1d_c1d_rs
3210 FUNCTION pw_integral_ab_r1d_c1d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
3212 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
3213 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw2
3214 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
3215 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
3216 REAL(kind=dp) :: integral_value
3218 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_r1d_c1d_rs'
3220 INTEGER :: handle, loc_sumtype
3221 LOGICAL :: my_just_sum, my_local_only
3223 CALL timeset(routinen, handle)
3226 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
3228 my_local_only = .false.
3229 IF (
PRESENT(local_only)) my_local_only = local_only
3231 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
3232 cpabort(
"Grids incompatible")
3235 my_just_sum = .false.
3236 IF (
PRESENT(just_sum)) my_just_sum = just_sum
3243 integral_value = sum(pw1%array*real(pw2%array, kind=dp))
3248 integral_value = accurate_sum(pw1%array*real(pw2%array, kind=dp))
3252 IF (.NOT. my_just_sum)
THEN
3253 integral_value = integral_value*pw1%pw_grid%dvol
3256 IF (pw1%pw_grid%grid_span == halfspace)
THEN
3257 integral_value = 2.0_dp*integral_value
3258 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
3259 pw1%array(1)*real(pw2%array(1), kind=dp)
3262 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
3263 CALL pw1%pw_grid%para%group%sum(integral_value)
3265 CALL timestop(handle)
3267 END FUNCTION pw_integral_ab_r1d_c1d_rs
3281 SUBROUTINE pw_copy_r1d_c1d_gs (pw1, pw2)
3283 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
3284 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
3286 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
3289 INTEGER :: i, j, ng, ng1, ng2, ns
3291 CALL timeset(routinen, handle)
3293 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
3294 cpabort(
"Both grids must be either spherical or non-spherical!")
3296 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
3297 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
3298 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
3299 ng1 =
SIZE(pw1%array)
3300 ng2 =
SIZE(pw2%array)
3303 pw2%array(1:ng) = cmplx(pw1%array(1:ng), 0.0_dp, kind=dp)
3307 pw2%array(ng + 1:ng2) = cmplx(0.0_dp, 0.0_dp, kind=dp)
3311 cpabort(
"Copies between spherical grids require compatible grids!")
3314 ng1 =
SIZE(pw1%array)
3315 ng2 =
SIZE(pw2%array)
3316 ns = 2*max(ng1, ng2)
3318 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
3319 IF (ng1 >= ng2)
THEN
3324 j = pw2%pw_grid%gidx(i)
3325 pw2%array(i) = cmplx(pw1%array(j), 0.0_dp, kind=dp)
3334 j = pw2%pw_grid%gidx(i)
3335 pw2%array(j) = cmplx(pw1%array(i), 0.0_dp, kind=dp)
3339 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
3340 IF (ng1 >= ng2)
THEN
3343 j = pw1%pw_grid%gidx(i)
3344 pw2%array(i) = cmplx(pw1%array(j), 0.0_dp, kind=dp)
3351 j = pw1%pw_grid%gidx(i)
3352 pw2%array(j) = cmplx(pw1%array(i), 0.0_dp, kind=dp)
3357 cpabort(
"Copy not implemented!")
3364 pw2%array = cmplx(pw1%array, 0.0_dp, kind=dp)
3368 CALL timestop(handle)
3370 END SUBROUTINE pw_copy_r1d_c1d_gs
3377 SUBROUTINE pw_copy_to_array_r1d_c1d_gs (pw, array)
3378 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw
3379 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(INOUT) :: array
3381 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
3385 CALL timeset(routinen, handle)
3388 array(:) = cmplx(pw%array(:), 0.0_dp, kind=dp)
3391 CALL timestop(handle)
3392 END SUBROUTINE pw_copy_to_array_r1d_c1d_gs
3399 SUBROUTINE pw_copy_from_array_r1d_c1d_gs (pw, array)
3400 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw
3401 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(IN) :: array
3403 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
3407 CALL timeset(routinen, handle)
3410 pw%array = real(array, kind=dp)
3413 CALL timestop(handle)
3414 END SUBROUTINE pw_copy_from_array_r1d_c1d_gs
3432 SUBROUTINE pw_axpy_r1d_c1d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
3434 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
3435 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
3436 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
3437 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
3439 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
3442 LOGICAL :: my_allow_noncompatible_grids
3443 REAL(KIND=dp) :: my_alpha, my_beta
3444 INTEGER :: i, j, ng, ng1, ng2
3446 CALL timeset(routinen, handle)
3449 IF (
PRESENT(alpha)) my_alpha = alpha
3452 IF (
PRESENT(beta)) my_beta = beta
3454 my_allow_noncompatible_grids = .false.
3455 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
3457 IF (my_beta /= 1.0_dp)
THEN
3458 IF (my_beta == 0.0_dp)
THEN
3462 pw2%array = pw2%array*my_beta
3467 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
3469 IF (my_alpha == 1.0_dp)
THEN
3471 pw2%array = pw2%array + cmplx(pw1%array, 0.0_dp, kind=dp)
3473 ELSE IF (my_alpha /= 0.0_dp)
THEN
3475 pw2%array = pw2%array + my_alpha* cmplx(pw1%array, 0.0_dp, kind=dp)
3479 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
3481 ng1 =
SIZE(pw1%array)
3482 ng2 =
SIZE(pw2%array)
3485 IF (pw1%pw_grid%spherical)
THEN
3486 IF (my_alpha == 1.0_dp)
THEN
3489 pw2%array(i) = pw2%array(i) + cmplx(pw1%array(i), 0.0_dp, kind=dp)
3492 ELSE IF (my_alpha /= 0.0_dp)
THEN
3495 pw2%array(i) = pw2%array(i) + my_alpha* cmplx(pw1%array(i), 0.0_dp, kind=dp)
3499 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
3500 IF (ng1 >= ng2)
THEN
3501 IF (my_alpha == 1.0_dp)
THEN
3504 j = pw2%pw_grid%gidx(i)
3505 pw2%array(i) = pw2%array(i) + cmplx(pw1%array(j), 0.0_dp, kind=dp)
3508 ELSE IF (my_alpha /= 0.0_dp)
THEN
3511 j = pw2%pw_grid%gidx(i)
3512 pw2%array(i) = pw2%array(i) + my_alpha* cmplx(pw1%array(j), 0.0_dp, kind=dp)
3517 IF (my_alpha == 1.0_dp)
THEN
3520 j = pw2%pw_grid%gidx(i)
3521 pw2%array(j) = pw2%array(j) + cmplx(pw1%array(i), 0.0_dp, kind=dp)
3524 ELSE IF (my_alpha /= 0.0_dp)
THEN
3527 j = pw2%pw_grid%gidx(i)
3528 pw2%array(j) = pw2%array(j) + my_alpha* cmplx(pw1%array(i), 0.0_dp, kind=dp)
3533 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
3534 IF (ng1 >= ng2)
THEN
3535 IF (my_alpha == 1.0_dp)
THEN
3538 j = pw1%pw_grid%gidx(i)
3539 pw2%array(i) = pw2%array(i) + cmplx(pw1%array(j), 0.0_dp, kind=dp)
3542 ELSE IF (my_alpha /= 0.0_dp)
THEN
3545 j = pw1%pw_grid%gidx(i)
3546 pw2%array(i) = pw2%array(i) + my_alpha* cmplx(pw1%array(j), 0.0_dp, kind=dp)
3551 IF (my_alpha == 1.0_dp)
THEN
3554 j = pw1%pw_grid%gidx(i)
3555 pw2%array(j) = pw2%array(j) + cmplx(pw1%array(i), 0.0_dp, kind=dp)
3558 ELSE IF (my_alpha /= 0.0_dp)
THEN
3561 j = pw1%pw_grid%gidx(i)
3562 pw2%array(j) = pw2%array(j) + my_alpha* cmplx(pw1%array(i), 0.0_dp, kind=dp)
3568 cpabort(
"Grids not compatible")
3573 cpabort(
"Grids not compatible")
3577 CALL timestop(handle)
3579 END SUBROUTINE pw_axpy_r1d_c1d_gs
3590 SUBROUTINE pw_multiply_r1d_c1d_gs (pw_out, pw1, pw2, alpha)
3592 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw_out
3593 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
3594 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
3595 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
3597 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
3600 REAL(KIND=dp) :: my_alpha
3602 CALL timeset(routinen, handle)
3605 IF (
PRESENT(alpha)) my_alpha = alpha
3607 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
3608 cpabort(
"pw_multiply not implemented for non-identical grids!")
3610#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
3611 IF (my_alpha == 1.0_dp)
THEN
3613 pw_out%array = pw_out%array + pw1%array* real(pw2%array, kind=dp)
3617 pw_out%array = pw_out%array + my_alpha*pw1%array* real(pw2%array, kind=dp)
3621 IF (my_alpha == 1.0_dp)
THEN
3622 pw_out%array = pw_out%array + pw1%array* real(pw2%array, kind=dp)
3624 pw_out%array = pw_out%array + my_alpha*pw1%array* real(pw2%array, kind=dp)
3628 CALL timestop(handle)
3630 END SUBROUTINE pw_multiply_r1d_c1d_gs
3637 SUBROUTINE pw_multiply_with_r1d_c1d_gs (pw1, pw2)
3638 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw1
3639 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
3641 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
3645 CALL timeset(routinen, handle)
3647 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
3648 cpabort(
"Incompatible grids!")
3651 pw1%array = pw1%array* real(pw2%array, kind=dp)
3654 CALL timestop(handle)
3656 END SUBROUTINE pw_multiply_with_r1d_c1d_gs
3671 FUNCTION pw_integral_ab_r1d_c1d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
3673 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
3674 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
3675 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
3676 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
3677 REAL(kind=dp) :: integral_value
3679 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_r1d_c1d_gs'
3681 INTEGER :: handle, loc_sumtype
3682 LOGICAL :: my_just_sum, my_local_only
3684 CALL timeset(routinen, handle)
3687 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
3689 my_local_only = .false.
3690 IF (
PRESENT(local_only)) my_local_only = local_only
3692 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
3693 cpabort(
"Grids incompatible")
3696 my_just_sum = .false.
3697 IF (
PRESENT(just_sum)) my_just_sum = just_sum
3704 integral_value = sum(pw1%array*real(pw2%array, kind=dp))
3709 integral_value = accurate_sum(pw1%array*real(pw2%array, kind=dp))
3713 IF (.NOT. my_just_sum)
THEN
3714 integral_value = integral_value*pw1%pw_grid%vol
3717 IF (pw1%pw_grid%grid_span == halfspace)
THEN
3718 integral_value = 2.0_dp*integral_value
3719 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
3720 pw1%array(1)*real(pw2%array(1), kind=dp)
3723 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
3724 CALL pw1%pw_grid%para%group%sum(integral_value)
3726 CALL timestop(handle)
3728 END FUNCTION pw_integral_ab_r1d_c1d_gs
3736 FUNCTION pw_integral_a2b_r1d_c1d (pw1, pw2)
RESULT(integral_value)
3738 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
3739 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
3740 REAL(kind=dp) :: integral_value
3742 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_a2b'
3746 CALL timeset(routinen, handle)
3748 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
3749 cpabort(
"Grids incompatible")
3752 integral_value = accurate_sum(pw1%array*real(pw2%array, kind=dp)*pw1%pw_grid%gsq)
3753 IF (pw1%pw_grid%grid_span == halfspace)
THEN
3754 integral_value = 2.0_dp*integral_value
3757 integral_value = integral_value*pw1%pw_grid%vol
3759 IF (pw1%pw_grid%para%mode == pw_mode_distributed) &
3760 CALL pw1%pw_grid%para%group%sum(integral_value)
3761 CALL timestop(handle)
3763 END FUNCTION pw_integral_a2b_r1d_c1d
3777 SUBROUTINE pw_copy_r3d_r3d_rs (pw1, pw2)
3779 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
3780 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw2
3782 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
3786 CALL timeset(routinen, handle)
3788 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
3789 cpabort(
"Both grids must be either spherical or non-spherical!")
3790 IF (pw1%pw_grid%spherical) &
3791 cpabort(
"Spherical grids only exist in reciprocal space!")
3793 IF (any(shape(pw2%array) /= shape(pw1%array))) &
3794 cpabort(
"3D grids must be compatible!")
3795 IF (pw1%pw_grid%spherical) &
3796 cpabort(
"3D grids must not be spherical!")
3798 pw2%array(:, :, :) = pw1%array(:, :, :)
3801 CALL timestop(handle)
3803 END SUBROUTINE pw_copy_r3d_r3d_rs
3810 SUBROUTINE pw_copy_to_array_r3d_r3d_rs (pw, array)
3811 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw
3812 REAL(KIND=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
3814 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
3818 CALL timeset(routinen, handle)
3821 array(:, :, :) = pw%array(:, :, :)
3824 CALL timestop(handle)
3825 END SUBROUTINE pw_copy_to_array_r3d_r3d_rs
3832 SUBROUTINE pw_copy_from_array_r3d_r3d_rs (pw, array)
3833 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw
3834 REAL(KIND=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
3836 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
3840 CALL timeset(routinen, handle)
3846 CALL timestop(handle)
3847 END SUBROUTINE pw_copy_from_array_r3d_r3d_rs
3865 SUBROUTINE pw_axpy_r3d_r3d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
3867 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
3868 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw2
3869 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
3870 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
3872 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
3875 LOGICAL :: my_allow_noncompatible_grids
3876 REAL(KIND=dp) :: my_alpha, my_beta
3878 CALL timeset(routinen, handle)
3881 IF (
PRESENT(alpha)) my_alpha = alpha
3884 IF (
PRESENT(beta)) my_beta = beta
3886 my_allow_noncompatible_grids = .false.
3887 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
3889 IF (my_beta /= 1.0_dp)
THEN
3890 IF (my_beta == 0.0_dp)
THEN
3894 pw2%array = pw2%array*my_beta
3899 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
3900 IF (my_alpha == 1.0_dp)
THEN
3902 pw2%array = pw2%array + pw1%array
3904 ELSE IF (my_alpha /= 0.0_dp)
THEN
3906 pw2%array = pw2%array + my_alpha* pw1%array
3910 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
3912 IF (any(shape(pw1%array) /= shape(pw2%array))) &
3913 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
3915 IF (my_alpha == 1.0_dp)
THEN
3917 pw2%array = pw2%array + pw1%array
3921 pw2%array = pw2%array + my_alpha* pw1%array
3927 cpabort(
"Grids not compatible")
3931 CALL timestop(handle)
3933 END SUBROUTINE pw_axpy_r3d_r3d_rs
3944 SUBROUTINE pw_multiply_r3d_r3d_rs (pw_out, pw1, pw2, alpha)
3946 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw_out
3947 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
3948 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw2
3949 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
3951 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
3954 REAL(KIND=dp) :: my_alpha
3956 CALL timeset(routinen, handle)
3959 IF (
PRESENT(alpha)) my_alpha = alpha
3961 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
3962 cpabort(
"pw_multiply not implemented for non-identical grids!")
3964#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
3965 IF (my_alpha == 1.0_dp)
THEN
3967 pw_out%array = pw_out%array + pw1%array* pw2%array
3971 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
3975 IF (my_alpha == 1.0_dp)
THEN
3976 pw_out%array = pw_out%array + pw1%array* pw2%array
3978 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
3982 CALL timestop(handle)
3984 END SUBROUTINE pw_multiply_r3d_r3d_rs
3991 SUBROUTINE pw_multiply_with_r3d_r3d_rs (pw1, pw2)
3992 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw1
3993 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw2
3995 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
3999 CALL timeset(routinen, handle)
4001 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
4002 cpabort(
"Incompatible grids!")
4005 pw1%array = pw1%array* pw2%array
4008 CALL timestop(handle)
4010 END SUBROUTINE pw_multiply_with_r3d_r3d_rs
4025 FUNCTION pw_integral_ab_r3d_r3d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
4027 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
4028 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw2
4029 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
4030 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
4031 REAL(kind=dp) :: integral_value
4033 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_r3d_r3d_rs'
4035 INTEGER :: handle, loc_sumtype
4036 LOGICAL :: my_just_sum, my_local_only
4038 CALL timeset(routinen, handle)
4041 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
4043 my_local_only = .false.
4044 IF (
PRESENT(local_only)) my_local_only = local_only
4046 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
4047 cpabort(
"Grids incompatible")
4050 my_just_sum = .false.
4051 IF (
PRESENT(just_sum)) my_just_sum = just_sum
4058 integral_value = sum(pw1%array*pw2%array)
4063 integral_value = accurate_dot_product(pw1%array, pw2%array)
4067 IF (.NOT. my_just_sum)
THEN
4068 integral_value = integral_value*pw1%pw_grid%dvol
4072 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
4073 CALL pw1%pw_grid%para%group%sum(integral_value)
4075 CALL timestop(handle)
4077 END FUNCTION pw_integral_ab_r3d_r3d_rs
4091 SUBROUTINE pw_copy_r3d_r3d_gs (pw1, pw2)
4093 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4094 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw2
4096 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
4100 CALL timeset(routinen, handle)
4102 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
4103 cpabort(
"Both grids must be either spherical or non-spherical!")
4105 IF (any(shape(pw2%array) /= shape(pw1%array))) &
4106 cpabort(
"3D grids must be compatible!")
4107 IF (pw1%pw_grid%spherical) &
4108 cpabort(
"3D grids must not be spherical!")
4110 pw2%array(:, :, :) = pw1%array(:, :, :)
4113 CALL timestop(handle)
4115 END SUBROUTINE pw_copy_r3d_r3d_gs
4122 SUBROUTINE pw_copy_to_array_r3d_r3d_gs (pw, array)
4123 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw
4124 REAL(KIND=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
4126 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
4130 CALL timeset(routinen, handle)
4133 array(:, :, :) = pw%array(:, :, :)
4136 CALL timestop(handle)
4137 END SUBROUTINE pw_copy_to_array_r3d_r3d_gs
4144 SUBROUTINE pw_copy_from_array_r3d_r3d_gs (pw, array)
4145 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw
4146 REAL(KIND=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
4148 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
4152 CALL timeset(routinen, handle)
4158 CALL timestop(handle)
4159 END SUBROUTINE pw_copy_from_array_r3d_r3d_gs
4177 SUBROUTINE pw_axpy_r3d_r3d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
4179 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4180 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw2
4181 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
4182 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
4184 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
4187 LOGICAL :: my_allow_noncompatible_grids
4188 REAL(KIND=dp) :: my_alpha, my_beta
4190 CALL timeset(routinen, handle)
4193 IF (
PRESENT(alpha)) my_alpha = alpha
4196 IF (
PRESENT(beta)) my_beta = beta
4198 my_allow_noncompatible_grids = .false.
4199 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
4201 IF (my_beta /= 1.0_dp)
THEN
4202 IF (my_beta == 0.0_dp)
THEN
4206 pw2%array = pw2%array*my_beta
4211 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
4212 IF (my_alpha == 1.0_dp)
THEN
4214 pw2%array = pw2%array + pw1%array
4216 ELSE IF (my_alpha /= 0.0_dp)
THEN
4218 pw2%array = pw2%array + my_alpha* pw1%array
4222 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
4224 IF (any(shape(pw1%array) /= shape(pw2%array))) &
4225 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
4227 IF (my_alpha == 1.0_dp)
THEN
4229 pw2%array = pw2%array + pw1%array
4233 pw2%array = pw2%array + my_alpha* pw1%array
4239 cpabort(
"Grids not compatible")
4243 CALL timestop(handle)
4245 END SUBROUTINE pw_axpy_r3d_r3d_gs
4256 SUBROUTINE pw_multiply_r3d_r3d_gs (pw_out, pw1, pw2, alpha)
4258 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw_out
4259 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4260 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw2
4261 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
4263 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
4266 REAL(KIND=dp) :: my_alpha
4268 CALL timeset(routinen, handle)
4271 IF (
PRESENT(alpha)) my_alpha = alpha
4273 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
4274 cpabort(
"pw_multiply not implemented for non-identical grids!")
4276#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
4277 IF (my_alpha == 1.0_dp)
THEN
4279 pw_out%array = pw_out%array + pw1%array* pw2%array
4283 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
4287 IF (my_alpha == 1.0_dp)
THEN
4288 pw_out%array = pw_out%array + pw1%array* pw2%array
4290 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
4294 CALL timestop(handle)
4296 END SUBROUTINE pw_multiply_r3d_r3d_gs
4303 SUBROUTINE pw_multiply_with_r3d_r3d_gs (pw1, pw2)
4304 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw1
4305 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw2
4307 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
4311 CALL timeset(routinen, handle)
4313 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
4314 cpabort(
"Incompatible grids!")
4317 pw1%array = pw1%array* pw2%array
4320 CALL timestop(handle)
4322 END SUBROUTINE pw_multiply_with_r3d_r3d_gs
4337 FUNCTION pw_integral_ab_r3d_r3d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
4339 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4340 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw2
4341 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
4342 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
4343 REAL(kind=dp) :: integral_value
4345 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_r3d_r3d_gs'
4347 INTEGER :: handle, loc_sumtype
4348 LOGICAL :: my_just_sum, my_local_only
4350 CALL timeset(routinen, handle)
4353 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
4355 my_local_only = .false.
4356 IF (
PRESENT(local_only)) my_local_only = local_only
4358 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
4359 cpabort(
"Grids incompatible")
4362 my_just_sum = .false.
4363 IF (
PRESENT(just_sum)) my_just_sum = just_sum
4370 integral_value = sum(pw1%array*pw2%array)
4375 integral_value = accurate_dot_product(pw1%array, pw2%array)
4379 IF (.NOT. my_just_sum)
THEN
4380 integral_value = integral_value*pw1%pw_grid%vol
4384 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
4385 CALL pw1%pw_grid%para%group%sum(integral_value)
4387 CALL timestop(handle)
4389 END FUNCTION pw_integral_ab_r3d_r3d_gs
4404 SUBROUTINE pw_copy_r3d_c3d_rs (pw1, pw2)
4406 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
4407 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw2
4409 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
4413 CALL timeset(routinen, handle)
4415 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
4416 cpabort(
"Both grids must be either spherical or non-spherical!")
4417 IF (pw1%pw_grid%spherical) &
4418 cpabort(
"Spherical grids only exist in reciprocal space!")
4420 IF (any(shape(pw2%array) /= shape(pw1%array))) &
4421 cpabort(
"3D grids must be compatible!")
4422 IF (pw1%pw_grid%spherical) &
4423 cpabort(
"3D grids must not be spherical!")
4425 pw2%array(:, :, :) = cmplx(pw1%array(:, :, :), 0.0_dp, kind=dp)
4428 CALL timestop(handle)
4430 END SUBROUTINE pw_copy_r3d_c3d_rs
4437 SUBROUTINE pw_copy_to_array_r3d_c3d_rs (pw, array)
4438 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw
4439 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
4441 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
4445 CALL timeset(routinen, handle)
4448 array(:, :, :) = cmplx(pw%array(:, :, :), 0.0_dp, kind=dp)
4451 CALL timestop(handle)
4452 END SUBROUTINE pw_copy_to_array_r3d_c3d_rs
4459 SUBROUTINE pw_copy_from_array_r3d_c3d_rs (pw, array)
4460 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw
4461 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
4463 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
4467 CALL timeset(routinen, handle)
4470 pw%array = real(array, kind=dp)
4473 CALL timestop(handle)
4474 END SUBROUTINE pw_copy_from_array_r3d_c3d_rs
4492 SUBROUTINE pw_axpy_r3d_c3d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
4494 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
4495 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw2
4496 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
4497 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
4499 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
4502 LOGICAL :: my_allow_noncompatible_grids
4503 REAL(KIND=dp) :: my_alpha, my_beta
4505 CALL timeset(routinen, handle)
4508 IF (
PRESENT(alpha)) my_alpha = alpha
4511 IF (
PRESENT(beta)) my_beta = beta
4513 my_allow_noncompatible_grids = .false.
4514 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
4516 IF (my_beta /= 1.0_dp)
THEN
4517 IF (my_beta == 0.0_dp)
THEN
4521 pw2%array = pw2%array*my_beta
4526 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
4527 IF (my_alpha == 1.0_dp)
THEN
4529 pw2%array = pw2%array + cmplx(pw1%array, 0.0_dp, kind=dp)
4531 ELSE IF (my_alpha /= 0.0_dp)
THEN
4533 pw2%array = pw2%array + my_alpha* cmplx(pw1%array, 0.0_dp, kind=dp)
4537 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
4539 IF (any(shape(pw1%array) /= shape(pw2%array))) &
4540 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
4542 IF (my_alpha == 1.0_dp)
THEN
4544 pw2%array = pw2%array + cmplx(pw1%array, 0.0_dp, kind=dp)
4548 pw2%array = pw2%array + my_alpha* cmplx(pw1%array, 0.0_dp, kind=dp)
4554 cpabort(
"Grids not compatible")
4558 CALL timestop(handle)
4560 END SUBROUTINE pw_axpy_r3d_c3d_rs
4571 SUBROUTINE pw_multiply_r3d_c3d_rs (pw_out, pw1, pw2, alpha)
4573 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw_out
4574 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
4575 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw2
4576 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
4578 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
4581 REAL(KIND=dp) :: my_alpha
4583 CALL timeset(routinen, handle)
4586 IF (
PRESENT(alpha)) my_alpha = alpha
4588 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
4589 cpabort(
"pw_multiply not implemented for non-identical grids!")
4591#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
4592 IF (my_alpha == 1.0_dp)
THEN
4594 pw_out%array = pw_out%array + pw1%array* real(pw2%array, kind=dp)
4598 pw_out%array = pw_out%array + my_alpha*pw1%array* real(pw2%array, kind=dp)
4602 IF (my_alpha == 1.0_dp)
THEN
4603 pw_out%array = pw_out%array + pw1%array* real(pw2%array, kind=dp)
4605 pw_out%array = pw_out%array + my_alpha*pw1%array* real(pw2%array, kind=dp)
4609 CALL timestop(handle)
4611 END SUBROUTINE pw_multiply_r3d_c3d_rs
4618 SUBROUTINE pw_multiply_with_r3d_c3d_rs (pw1, pw2)
4619 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw1
4620 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw2
4622 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
4626 CALL timeset(routinen, handle)
4628 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
4629 cpabort(
"Incompatible grids!")
4632 pw1%array = pw1%array* real(pw2%array, kind=dp)
4635 CALL timestop(handle)
4637 END SUBROUTINE pw_multiply_with_r3d_c3d_rs
4652 FUNCTION pw_integral_ab_r3d_c3d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
4654 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
4655 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw2
4656 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
4657 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
4658 REAL(kind=dp) :: integral_value
4660 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_r3d_c3d_rs'
4662 INTEGER :: handle, loc_sumtype
4663 LOGICAL :: my_just_sum, my_local_only
4665 CALL timeset(routinen, handle)
4668 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
4670 my_local_only = .false.
4671 IF (
PRESENT(local_only)) my_local_only = local_only
4673 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
4674 cpabort(
"Grids incompatible")
4677 my_just_sum = .false.
4678 IF (
PRESENT(just_sum)) my_just_sum = just_sum
4685 integral_value = sum(pw1%array*real(pw2%array, kind=dp))
4690 integral_value = accurate_sum(pw1%array*real(pw2%array, kind=dp))
4694 IF (.NOT. my_just_sum)
THEN
4695 integral_value = integral_value*pw1%pw_grid%dvol
4699 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
4700 CALL pw1%pw_grid%para%group%sum(integral_value)
4702 CALL timestop(handle)
4704 END FUNCTION pw_integral_ab_r3d_c3d_rs
4718 SUBROUTINE pw_copy_r3d_c3d_gs (pw1, pw2)
4720 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4721 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
4723 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
4727 CALL timeset(routinen, handle)
4729 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
4730 cpabort(
"Both grids must be either spherical or non-spherical!")
4732 IF (any(shape(pw2%array) /= shape(pw1%array))) &
4733 cpabort(
"3D grids must be compatible!")
4734 IF (pw1%pw_grid%spherical) &
4735 cpabort(
"3D grids must not be spherical!")
4737 pw2%array(:, :, :) = cmplx(pw1%array(:, :, :), 0.0_dp, kind=dp)
4740 CALL timestop(handle)
4742 END SUBROUTINE pw_copy_r3d_c3d_gs
4749 SUBROUTINE pw_copy_to_array_r3d_c3d_gs (pw, array)
4750 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw
4751 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
4753 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
4757 CALL timeset(routinen, handle)
4760 array(:, :, :) = cmplx(pw%array(:, :, :), 0.0_dp, kind=dp)
4763 CALL timestop(handle)
4764 END SUBROUTINE pw_copy_to_array_r3d_c3d_gs
4771 SUBROUTINE pw_copy_from_array_r3d_c3d_gs (pw, array)
4772 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw
4773 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
4775 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
4779 CALL timeset(routinen, handle)
4782 pw%array = real(array, kind=dp)
4785 CALL timestop(handle)
4786 END SUBROUTINE pw_copy_from_array_r3d_c3d_gs
4804 SUBROUTINE pw_axpy_r3d_c3d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
4806 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4807 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
4808 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
4809 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
4811 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
4814 LOGICAL :: my_allow_noncompatible_grids
4815 REAL(KIND=dp) :: my_alpha, my_beta
4817 CALL timeset(routinen, handle)
4820 IF (
PRESENT(alpha)) my_alpha = alpha
4823 IF (
PRESENT(beta)) my_beta = beta
4825 my_allow_noncompatible_grids = .false.
4826 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
4828 IF (my_beta /= 1.0_dp)
THEN
4829 IF (my_beta == 0.0_dp)
THEN
4833 pw2%array = pw2%array*my_beta
4838 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
4839 IF (my_alpha == 1.0_dp)
THEN
4841 pw2%array = pw2%array + cmplx(pw1%array, 0.0_dp, kind=dp)
4843 ELSE IF (my_alpha /= 0.0_dp)
THEN
4845 pw2%array = pw2%array + my_alpha* cmplx(pw1%array, 0.0_dp, kind=dp)
4849 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
4851 IF (any(shape(pw1%array) /= shape(pw2%array))) &
4852 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
4854 IF (my_alpha == 1.0_dp)
THEN
4856 pw2%array = pw2%array + cmplx(pw1%array, 0.0_dp, kind=dp)
4860 pw2%array = pw2%array + my_alpha* cmplx(pw1%array, 0.0_dp, kind=dp)
4866 cpabort(
"Grids not compatible")
4870 CALL timestop(handle)
4872 END SUBROUTINE pw_axpy_r3d_c3d_gs
4883 SUBROUTINE pw_multiply_r3d_c3d_gs (pw_out, pw1, pw2, alpha)
4885 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw_out
4886 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4887 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw2
4888 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
4890 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
4893 REAL(KIND=dp) :: my_alpha
4895 CALL timeset(routinen, handle)
4898 IF (
PRESENT(alpha)) my_alpha = alpha
4900 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
4901 cpabort(
"pw_multiply not implemented for non-identical grids!")
4903#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
4904 IF (my_alpha == 1.0_dp)
THEN
4906 pw_out%array = pw_out%array + pw1%array* real(pw2%array, kind=dp)
4910 pw_out%array = pw_out%array + my_alpha*pw1%array* real(pw2%array, kind=dp)
4914 IF (my_alpha == 1.0_dp)
THEN
4915 pw_out%array = pw_out%array + pw1%array* real(pw2%array, kind=dp)
4917 pw_out%array = pw_out%array + my_alpha*pw1%array* real(pw2%array, kind=dp)
4921 CALL timestop(handle)
4923 END SUBROUTINE pw_multiply_r3d_c3d_gs
4930 SUBROUTINE pw_multiply_with_r3d_c3d_gs (pw1, pw2)
4931 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw1
4932 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw2
4934 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
4938 CALL timeset(routinen, handle)
4940 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
4941 cpabort(
"Incompatible grids!")
4944 pw1%array = pw1%array* real(pw2%array, kind=dp)
4947 CALL timestop(handle)
4949 END SUBROUTINE pw_multiply_with_r3d_c3d_gs
4964 FUNCTION pw_integral_ab_r3d_c3d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
4966 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4967 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw2
4968 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
4969 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
4970 REAL(kind=dp) :: integral_value
4972 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_r3d_c3d_gs'
4974 INTEGER :: handle, loc_sumtype
4975 LOGICAL :: my_just_sum, my_local_only
4977 CALL timeset(routinen, handle)
4980 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
4982 my_local_only = .false.
4983 IF (
PRESENT(local_only)) my_local_only = local_only
4985 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
4986 cpabort(
"Grids incompatible")
4989 my_just_sum = .false.
4990 IF (
PRESENT(just_sum)) my_just_sum = just_sum
4997 integral_value = sum(pw1%array*real(pw2%array, kind=dp))
5002 integral_value = accurate_sum(pw1%array*real(pw2%array, kind=dp))
5006 IF (.NOT. my_just_sum)
THEN
5007 integral_value = integral_value*pw1%pw_grid%vol
5011 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
5012 CALL pw1%pw_grid%para%group%sum(integral_value)
5014 CALL timestop(handle)
5016 END FUNCTION pw_integral_ab_r3d_c3d_gs
5031 SUBROUTINE pw_copy_c1d_r1d_rs (pw1, pw2)
5033 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
5034 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw2
5036 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
5039 INTEGER :: i, j, ng, ng1, ng2, ns
5041 CALL timeset(routinen, handle)
5043 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
5044 cpabort(
"Both grids must be either spherical or non-spherical!")
5045 IF (pw1%pw_grid%spherical) &
5046 cpabort(
"Spherical grids only exist in reciprocal space!")
5048 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
5049 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
5050 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
5051 ng1 =
SIZE(pw1%array)
5052 ng2 =
SIZE(pw2%array)
5055 pw2%array(1:ng) = real(pw1%array(1:ng), kind=dp)
5059 pw2%array(ng + 1:ng2) = 0.0_dp
5063 cpabort(
"Copies between spherical grids require compatible grids!")
5066 ng1 =
SIZE(pw1%array)
5067 ng2 =
SIZE(pw2%array)
5068 ns = 2*max(ng1, ng2)
5070 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
5071 IF (ng1 >= ng2)
THEN
5076 j = pw2%pw_grid%gidx(i)
5077 pw2%array(i) = real(pw1%array(j), kind=dp)
5086 j = pw2%pw_grid%gidx(i)
5087 pw2%array(j) = real(pw1%array(i), kind=dp)
5091 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
5092 IF (ng1 >= ng2)
THEN
5095 j = pw1%pw_grid%gidx(i)
5096 pw2%array(i) = real(pw1%array(j), kind=dp)
5103 j = pw1%pw_grid%gidx(i)
5104 pw2%array(j) = real(pw1%array(i), kind=dp)
5109 cpabort(
"Copy not implemented!")
5116 pw2%array = real(pw1%array, kind=dp)
5120 CALL timestop(handle)
5122 END SUBROUTINE pw_copy_c1d_r1d_rs
5129 SUBROUTINE pw_copy_to_array_c1d_r1d_rs (pw, array)
5130 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw
5131 REAL(KIND=dp),
DIMENSION(:),
INTENT(INOUT) :: array
5133 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
5137 CALL timeset(routinen, handle)
5140 array(:) = real(pw%array(:), kind=dp)
5143 CALL timestop(handle)
5144 END SUBROUTINE pw_copy_to_array_c1d_r1d_rs
5151 SUBROUTINE pw_copy_from_array_c1d_r1d_rs (pw, array)
5152 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw
5153 REAL(KIND=dp),
DIMENSION(:),
INTENT(IN) :: array
5155 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
5159 CALL timeset(routinen, handle)
5162 pw%array = cmplx(array, 0.0_dp, kind=dp)
5165 CALL timestop(handle)
5166 END SUBROUTINE pw_copy_from_array_c1d_r1d_rs
5184 SUBROUTINE pw_axpy_c1d_r1d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
5186 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
5187 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw2
5188 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
5189 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
5191 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
5194 LOGICAL :: my_allow_noncompatible_grids
5195 REAL(KIND=dp) :: my_alpha, my_beta
5196 INTEGER :: i, j, ng, ng1, ng2
5198 CALL timeset(routinen, handle)
5201 IF (
PRESENT(alpha)) my_alpha = alpha
5204 IF (
PRESENT(beta)) my_beta = beta
5206 my_allow_noncompatible_grids = .false.
5207 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
5209 IF (my_beta /= 1.0_dp)
THEN
5210 IF (my_beta == 0.0_dp)
THEN
5214 pw2%array = pw2%array*my_beta
5219 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
5221 IF (my_alpha == 1.0_dp)
THEN
5223 pw2%array = pw2%array + real(pw1%array, kind=dp)
5225 ELSE IF (my_alpha /= 0.0_dp)
THEN
5227 pw2%array = pw2%array + my_alpha* real(pw1%array, kind=dp)
5231 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
5233 ng1 =
SIZE(pw1%array)
5234 ng2 =
SIZE(pw2%array)
5237 IF (pw1%pw_grid%spherical)
THEN
5238 IF (my_alpha == 1.0_dp)
THEN
5241 pw2%array(i) = pw2%array(i) + real(pw1%array(i), kind=dp)
5244 ELSE IF (my_alpha /= 0.0_dp)
THEN
5247 pw2%array(i) = pw2%array(i) + my_alpha* real(pw1%array(i), kind=dp)
5251 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
5252 IF (ng1 >= ng2)
THEN
5253 IF (my_alpha == 1.0_dp)
THEN
5256 j = pw2%pw_grid%gidx(i)
5257 pw2%array(i) = pw2%array(i) + real(pw1%array(j), kind=dp)
5260 ELSE IF (my_alpha /= 0.0_dp)
THEN
5263 j = pw2%pw_grid%gidx(i)
5264 pw2%array(i) = pw2%array(i) + my_alpha* real(pw1%array(j), kind=dp)
5269 IF (my_alpha == 1.0_dp)
THEN
5272 j = pw2%pw_grid%gidx(i)
5273 pw2%array(j) = pw2%array(j) + real(pw1%array(i), kind=dp)
5276 ELSE IF (my_alpha /= 0.0_dp)
THEN
5279 j = pw2%pw_grid%gidx(i)
5280 pw2%array(j) = pw2%array(j) + my_alpha* real(pw1%array(i), kind=dp)
5285 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
5286 IF (ng1 >= ng2)
THEN
5287 IF (my_alpha == 1.0_dp)
THEN
5290 j = pw1%pw_grid%gidx(i)
5291 pw2%array(i) = pw2%array(i) + real(pw1%array(j), kind=dp)
5294 ELSE IF (my_alpha /= 0.0_dp)
THEN
5297 j = pw1%pw_grid%gidx(i)
5298 pw2%array(i) = pw2%array(i) + my_alpha* real(pw1%array(j), kind=dp)
5303 IF (my_alpha == 1.0_dp)
THEN
5306 j = pw1%pw_grid%gidx(i)
5307 pw2%array(j) = pw2%array(j) + real(pw1%array(i), kind=dp)
5310 ELSE IF (my_alpha /= 0.0_dp)
THEN
5313 j = pw1%pw_grid%gidx(i)
5314 pw2%array(j) = pw2%array(j) + my_alpha* real(pw1%array(i), kind=dp)
5320 cpabort(
"Grids not compatible")
5325 cpabort(
"Grids not compatible")
5329 CALL timestop(handle)
5331 END SUBROUTINE pw_axpy_c1d_r1d_rs
5342 SUBROUTINE pw_multiply_c1d_r1d_rs (pw_out, pw1, pw2, alpha)
5344 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw_out
5345 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
5346 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw2
5347 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
5349 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
5352 REAL(KIND=dp) :: my_alpha
5354 CALL timeset(routinen, handle)
5357 IF (
PRESENT(alpha)) my_alpha = alpha
5359 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
5360 cpabort(
"pw_multiply not implemented for non-identical grids!")
5362#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
5363 IF (my_alpha == 1.0_dp)
THEN
5365 pw_out%array = pw_out%array + pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5369 pw_out%array = pw_out%array + my_alpha*pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5373 IF (my_alpha == 1.0_dp)
THEN
5374 pw_out%array = pw_out%array + pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5376 pw_out%array = pw_out%array + my_alpha*pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5380 CALL timestop(handle)
5382 END SUBROUTINE pw_multiply_c1d_r1d_rs
5389 SUBROUTINE pw_multiply_with_c1d_r1d_rs (pw1, pw2)
5390 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw1
5391 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw2
5393 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
5397 CALL timeset(routinen, handle)
5399 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
5400 cpabort(
"Incompatible grids!")
5403 pw1%array = pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5406 CALL timestop(handle)
5408 END SUBROUTINE pw_multiply_with_c1d_r1d_rs
5423 FUNCTION pw_integral_ab_c1d_r1d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
5425 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
5426 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw2
5427 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
5428 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
5429 REAL(kind=dp) :: integral_value
5431 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_c1d_r1d_rs'
5433 INTEGER :: handle, loc_sumtype
5434 LOGICAL :: my_just_sum, my_local_only
5436 CALL timeset(routinen, handle)
5439 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
5441 my_local_only = .false.
5442 IF (
PRESENT(local_only)) my_local_only = local_only
5444 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
5445 cpabort(
"Grids incompatible")
5448 my_just_sum = .false.
5449 IF (
PRESENT(just_sum)) my_just_sum = just_sum
5456 integral_value = sum(real(pw1%array, kind=dp)*pw2%array)
5461 integral_value = accurate_sum(real(pw1%array, kind=dp)*pw2%array)
5465 IF (.NOT. my_just_sum)
THEN
5466 integral_value = integral_value*pw1%pw_grid%dvol
5469 IF (pw1%pw_grid%grid_span == halfspace)
THEN
5470 integral_value = 2.0_dp*integral_value
5471 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
5472 REAL(conjg(pw1%array(1))*pw2%array(1), kind=dp)
5475 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
5476 CALL pw1%pw_grid%para%group%sum(integral_value)
5478 CALL timestop(handle)
5480 END FUNCTION pw_integral_ab_c1d_r1d_rs
5494 SUBROUTINE pw_copy_c1d_r1d_gs (pw1, pw2)
5496 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
5497 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw2
5499 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
5502 INTEGER :: i, j, ng, ng1, ng2, ns
5504 CALL timeset(routinen, handle)
5506 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
5507 cpabort(
"Both grids must be either spherical or non-spherical!")
5509 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
5510 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
5511 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
5512 ng1 =
SIZE(pw1%array)
5513 ng2 =
SIZE(pw2%array)
5516 pw2%array(1:ng) = real(pw1%array(1:ng), kind=dp)
5520 pw2%array(ng + 1:ng2) = 0.0_dp
5524 cpabort(
"Copies between spherical grids require compatible grids!")
5527 ng1 =
SIZE(pw1%array)
5528 ng2 =
SIZE(pw2%array)
5529 ns = 2*max(ng1, ng2)
5531 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
5532 IF (ng1 >= ng2)
THEN
5537 j = pw2%pw_grid%gidx(i)
5538 pw2%array(i) = real(pw1%array(j), kind=dp)
5547 j = pw2%pw_grid%gidx(i)
5548 pw2%array(j) = real(pw1%array(i), kind=dp)
5552 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
5553 IF (ng1 >= ng2)
THEN
5556 j = pw1%pw_grid%gidx(i)
5557 pw2%array(i) = real(pw1%array(j), kind=dp)
5564 j = pw1%pw_grid%gidx(i)
5565 pw2%array(j) = real(pw1%array(i), kind=dp)
5570 cpabort(
"Copy not implemented!")
5577 pw2%array = real(pw1%array, kind=dp)
5581 CALL timestop(handle)
5583 END SUBROUTINE pw_copy_c1d_r1d_gs
5590 SUBROUTINE pw_copy_to_array_c1d_r1d_gs (pw, array)
5591 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
5592 REAL(KIND=dp),
DIMENSION(:),
INTENT(INOUT) :: array
5594 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
5598 CALL timeset(routinen, handle)
5601 array(:) = real(pw%array(:), kind=dp)
5604 CALL timestop(handle)
5605 END SUBROUTINE pw_copy_to_array_c1d_r1d_gs
5612 SUBROUTINE pw_copy_from_array_c1d_r1d_gs (pw, array)
5613 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
5614 REAL(KIND=dp),
DIMENSION(:),
INTENT(IN) :: array
5616 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
5620 CALL timeset(routinen, handle)
5623 pw%array = cmplx(array, 0.0_dp, kind=dp)
5626 CALL timestop(handle)
5627 END SUBROUTINE pw_copy_from_array_c1d_r1d_gs
5645 SUBROUTINE pw_axpy_c1d_r1d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
5647 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
5648 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw2
5649 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
5650 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
5652 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
5655 LOGICAL :: my_allow_noncompatible_grids
5656 REAL(KIND=dp) :: my_alpha, my_beta
5657 INTEGER :: i, j, ng, ng1, ng2
5659 CALL timeset(routinen, handle)
5662 IF (
PRESENT(alpha)) my_alpha = alpha
5665 IF (
PRESENT(beta)) my_beta = beta
5667 my_allow_noncompatible_grids = .false.
5668 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
5670 IF (my_beta /= 1.0_dp)
THEN
5671 IF (my_beta == 0.0_dp)
THEN
5675 pw2%array = pw2%array*my_beta
5680 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
5682 IF (my_alpha == 1.0_dp)
THEN
5684 pw2%array = pw2%array + real(pw1%array, kind=dp)
5686 ELSE IF (my_alpha /= 0.0_dp)
THEN
5688 pw2%array = pw2%array + my_alpha* real(pw1%array, kind=dp)
5692 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
5694 ng1 =
SIZE(pw1%array)
5695 ng2 =
SIZE(pw2%array)
5698 IF (pw1%pw_grid%spherical)
THEN
5699 IF (my_alpha == 1.0_dp)
THEN
5702 pw2%array(i) = pw2%array(i) + real(pw1%array(i), kind=dp)
5705 ELSE IF (my_alpha /= 0.0_dp)
THEN
5708 pw2%array(i) = pw2%array(i) + my_alpha* real(pw1%array(i), kind=dp)
5712 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
5713 IF (ng1 >= ng2)
THEN
5714 IF (my_alpha == 1.0_dp)
THEN
5717 j = pw2%pw_grid%gidx(i)
5718 pw2%array(i) = pw2%array(i) + real(pw1%array(j), kind=dp)
5721 ELSE IF (my_alpha /= 0.0_dp)
THEN
5724 j = pw2%pw_grid%gidx(i)
5725 pw2%array(i) = pw2%array(i) + my_alpha* real(pw1%array(j), kind=dp)
5730 IF (my_alpha == 1.0_dp)
THEN
5733 j = pw2%pw_grid%gidx(i)
5734 pw2%array(j) = pw2%array(j) + real(pw1%array(i), kind=dp)
5737 ELSE IF (my_alpha /= 0.0_dp)
THEN
5740 j = pw2%pw_grid%gidx(i)
5741 pw2%array(j) = pw2%array(j) + my_alpha* real(pw1%array(i), kind=dp)
5746 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
5747 IF (ng1 >= ng2)
THEN
5748 IF (my_alpha == 1.0_dp)
THEN
5751 j = pw1%pw_grid%gidx(i)
5752 pw2%array(i) = pw2%array(i) + real(pw1%array(j), kind=dp)
5755 ELSE IF (my_alpha /= 0.0_dp)
THEN
5758 j = pw1%pw_grid%gidx(i)
5759 pw2%array(i) = pw2%array(i) + my_alpha* real(pw1%array(j), kind=dp)
5764 IF (my_alpha == 1.0_dp)
THEN
5767 j = pw1%pw_grid%gidx(i)
5768 pw2%array(j) = pw2%array(j) + real(pw1%array(i), kind=dp)
5771 ELSE IF (my_alpha /= 0.0_dp)
THEN
5774 j = pw1%pw_grid%gidx(i)
5775 pw2%array(j) = pw2%array(j) + my_alpha* real(pw1%array(i), kind=dp)
5781 cpabort(
"Grids not compatible")
5786 cpabort(
"Grids not compatible")
5790 CALL timestop(handle)
5792 END SUBROUTINE pw_axpy_c1d_r1d_gs
5803 SUBROUTINE pw_multiply_c1d_r1d_gs (pw_out, pw1, pw2, alpha)
5805 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw_out
5806 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
5807 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
5808 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
5810 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
5813 REAL(KIND=dp) :: my_alpha
5815 CALL timeset(routinen, handle)
5818 IF (
PRESENT(alpha)) my_alpha = alpha
5820 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
5821 cpabort(
"pw_multiply not implemented for non-identical grids!")
5823#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
5824 IF (my_alpha == 1.0_dp)
THEN
5826 pw_out%array = pw_out%array + pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5830 pw_out%array = pw_out%array + my_alpha*pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5834 IF (my_alpha == 1.0_dp)
THEN
5835 pw_out%array = pw_out%array + pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5837 pw_out%array = pw_out%array + my_alpha*pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5841 CALL timestop(handle)
5843 END SUBROUTINE pw_multiply_c1d_r1d_gs
5850 SUBROUTINE pw_multiply_with_c1d_r1d_gs (pw1, pw2)
5851 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw1
5852 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
5854 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
5858 CALL timeset(routinen, handle)
5860 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
5861 cpabort(
"Incompatible grids!")
5864 pw1%array = pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5867 CALL timestop(handle)
5869 END SUBROUTINE pw_multiply_with_c1d_r1d_gs
5884 FUNCTION pw_integral_ab_c1d_r1d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
5886 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
5887 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
5888 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
5889 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
5890 REAL(kind=dp) :: integral_value
5892 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_c1d_r1d_gs'
5894 INTEGER :: handle, loc_sumtype
5895 LOGICAL :: my_just_sum, my_local_only
5897 CALL timeset(routinen, handle)
5900 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
5902 my_local_only = .false.
5903 IF (
PRESENT(local_only)) my_local_only = local_only
5905 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
5906 cpabort(
"Grids incompatible")
5909 my_just_sum = .false.
5910 IF (
PRESENT(just_sum)) my_just_sum = just_sum
5917 integral_value = sum(real(pw1%array, kind=dp)*pw2%array)
5922 integral_value = accurate_sum(real(pw1%array, kind=dp)*pw2%array)
5926 IF (.NOT. my_just_sum)
THEN
5927 integral_value = integral_value*pw1%pw_grid%vol
5930 IF (pw1%pw_grid%grid_span == halfspace)
THEN
5931 integral_value = 2.0_dp*integral_value
5932 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
5933 REAL(conjg(pw1%array(1))*pw2%array(1), kind=dp)
5936 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
5937 CALL pw1%pw_grid%para%group%sum(integral_value)
5939 CALL timestop(handle)
5941 END FUNCTION pw_integral_ab_c1d_r1d_gs
5949 FUNCTION pw_integral_a2b_c1d_r1d (pw1, pw2)
RESULT(integral_value)
5951 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
5952 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
5953 REAL(kind=dp) :: integral_value
5955 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_a2b'
5959 CALL timeset(routinen, handle)
5961 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
5962 cpabort(
"Grids incompatible")
5965 integral_value = accurate_sum(real(conjg(pw1%array), kind=dp)*pw2%array*pw1%pw_grid%gsq)
5966 IF (pw1%pw_grid%grid_span == halfspace)
THEN
5967 integral_value = 2.0_dp*integral_value
5970 integral_value = integral_value*pw1%pw_grid%vol
5972 IF (pw1%pw_grid%para%mode == pw_mode_distributed) &
5973 CALL pw1%pw_grid%para%group%sum(integral_value)
5974 CALL timestop(handle)
5976 END FUNCTION pw_integral_a2b_c1d_r1d
5990 SUBROUTINE pw_copy_c1d_c1d_rs (pw1, pw2)
5992 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
5993 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw2
5995 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
5998 INTEGER :: i, j, ng, ng1, ng2, ns
6000 CALL timeset(routinen, handle)
6002 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
6003 cpabort(
"Both grids must be either spherical or non-spherical!")
6004 IF (pw1%pw_grid%spherical) &
6005 cpabort(
"Spherical grids only exist in reciprocal space!")
6007 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
6008 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
6009 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
6010 ng1 =
SIZE(pw1%array)
6011 ng2 =
SIZE(pw2%array)
6014 pw2%array(1:ng) = pw1%array(1:ng)
6018 pw2%array(ng + 1:ng2) = cmplx(0.0_dp, 0.0_dp, kind=dp)
6022 cpabort(
"Copies between spherical grids require compatible grids!")
6025 ng1 =
SIZE(pw1%array)
6026 ng2 =
SIZE(pw2%array)
6027 ns = 2*max(ng1, ng2)
6029 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
6030 IF (ng1 >= ng2)
THEN
6035 j = pw2%pw_grid%gidx(i)
6036 pw2%array(i) = pw1%array(j)
6045 j = pw2%pw_grid%gidx(i)
6046 pw2%array(j) = pw1%array(i)
6050 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
6051 IF (ng1 >= ng2)
THEN
6054 j = pw1%pw_grid%gidx(i)
6055 pw2%array(i) = pw1%array(j)
6062 j = pw1%pw_grid%gidx(i)
6063 pw2%array(j) = pw1%array(i)
6068 cpabort(
"Copy not implemented!")
6075 pw2%array = pw1%array
6079 CALL timestop(handle)
6081 END SUBROUTINE pw_copy_c1d_c1d_rs
6088 SUBROUTINE pw_copy_to_array_c1d_c1d_rs (pw, array)
6089 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw
6090 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(INOUT) :: array
6092 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
6096 CALL timeset(routinen, handle)
6099 array(:) = pw%array(:)
6102 CALL timestop(handle)
6103 END SUBROUTINE pw_copy_to_array_c1d_c1d_rs
6110 SUBROUTINE pw_copy_from_array_c1d_c1d_rs (pw, array)
6111 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw
6112 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(IN) :: array
6114 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
6118 CALL timeset(routinen, handle)
6124 CALL timestop(handle)
6125 END SUBROUTINE pw_copy_from_array_c1d_c1d_rs
6143 SUBROUTINE pw_axpy_c1d_c1d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
6145 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
6146 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw2
6147 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
6148 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
6150 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
6153 LOGICAL :: my_allow_noncompatible_grids
6154 REAL(KIND=dp) :: my_alpha, my_beta
6155 INTEGER :: i, j, ng, ng1, ng2
6157 CALL timeset(routinen, handle)
6160 IF (
PRESENT(alpha)) my_alpha = alpha
6163 IF (
PRESENT(beta)) my_beta = beta
6165 my_allow_noncompatible_grids = .false.
6166 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
6168 IF (my_beta /= 1.0_dp)
THEN
6169 IF (my_beta == 0.0_dp)
THEN
6173 pw2%array = pw2%array*my_beta
6178 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
6180 IF (my_alpha == 1.0_dp)
THEN
6182 pw2%array = pw2%array + pw1%array
6184 ELSE IF (my_alpha /= 0.0_dp)
THEN
6186 pw2%array = pw2%array + my_alpha* pw1%array
6190 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
6192 ng1 =
SIZE(pw1%array)
6193 ng2 =
SIZE(pw2%array)
6196 IF (pw1%pw_grid%spherical)
THEN
6197 IF (my_alpha == 1.0_dp)
THEN
6200 pw2%array(i) = pw2%array(i) + pw1%array(i)
6203 ELSE IF (my_alpha /= 0.0_dp)
THEN
6206 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(i)
6210 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
6211 IF (ng1 >= ng2)
THEN
6212 IF (my_alpha == 1.0_dp)
THEN
6215 j = pw2%pw_grid%gidx(i)
6216 pw2%array(i) = pw2%array(i) + pw1%array(j)
6219 ELSE IF (my_alpha /= 0.0_dp)
THEN
6222 j = pw2%pw_grid%gidx(i)
6223 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
6228 IF (my_alpha == 1.0_dp)
THEN
6231 j = pw2%pw_grid%gidx(i)
6232 pw2%array(j) = pw2%array(j) + pw1%array(i)
6235 ELSE IF (my_alpha /= 0.0_dp)
THEN
6238 j = pw2%pw_grid%gidx(i)
6239 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
6244 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
6245 IF (ng1 >= ng2)
THEN
6246 IF (my_alpha == 1.0_dp)
THEN
6249 j = pw1%pw_grid%gidx(i)
6250 pw2%array(i) = pw2%array(i) + pw1%array(j)
6253 ELSE IF (my_alpha /= 0.0_dp)
THEN
6256 j = pw1%pw_grid%gidx(i)
6257 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
6262 IF (my_alpha == 1.0_dp)
THEN
6265 j = pw1%pw_grid%gidx(i)
6266 pw2%array(j) = pw2%array(j) + pw1%array(i)
6269 ELSE IF (my_alpha /= 0.0_dp)
THEN
6272 j = pw1%pw_grid%gidx(i)
6273 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
6279 cpabort(
"Grids not compatible")
6284 cpabort(
"Grids not compatible")
6288 CALL timestop(handle)
6290 END SUBROUTINE pw_axpy_c1d_c1d_rs
6301 SUBROUTINE pw_multiply_c1d_c1d_rs (pw_out, pw1, pw2, alpha)
6303 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw_out
6304 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
6305 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw2
6306 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
6308 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
6311 REAL(KIND=dp) :: my_alpha
6313 CALL timeset(routinen, handle)
6316 IF (
PRESENT(alpha)) my_alpha = alpha
6318 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
6319 cpabort(
"pw_multiply not implemented for non-identical grids!")
6321#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
6322 IF (my_alpha == 1.0_dp)
THEN
6324 pw_out%array = pw_out%array + pw1%array* pw2%array
6328 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
6332 IF (my_alpha == 1.0_dp)
THEN
6333 pw_out%array = pw_out%array + pw1%array* pw2%array
6335 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
6339 CALL timestop(handle)
6341 END SUBROUTINE pw_multiply_c1d_c1d_rs
6348 SUBROUTINE pw_multiply_with_c1d_c1d_rs (pw1, pw2)
6349 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw1
6350 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw2
6352 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
6356 CALL timeset(routinen, handle)
6358 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
6359 cpabort(
"Incompatible grids!")
6362 pw1%array = pw1%array* pw2%array
6365 CALL timestop(handle)
6367 END SUBROUTINE pw_multiply_with_c1d_c1d_rs
6382 FUNCTION pw_integral_ab_c1d_c1d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
6384 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
6385 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw2
6386 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
6387 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
6388 REAL(kind=dp) :: integral_value
6390 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_c1d_c1d_rs'
6392 INTEGER :: handle, loc_sumtype
6393 LOGICAL :: my_just_sum, my_local_only
6395 CALL timeset(routinen, handle)
6398 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
6400 my_local_only = .false.
6401 IF (
PRESENT(local_only)) my_local_only = local_only
6403 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
6404 cpabort(
"Grids incompatible")
6407 my_just_sum = .false.
6408 IF (
PRESENT(just_sum)) my_just_sum = just_sum
6415 integral_value = sum(real(conjg(pw1%array) &
6416 *pw2%array, kind=dp))
6421 integral_value = accurate_sum(real(conjg(pw1%array)*pw2%array, kind=dp))
6425 IF (.NOT. my_just_sum)
THEN
6426 integral_value = integral_value*pw1%pw_grid%dvol
6429 IF (pw1%pw_grid%grid_span == halfspace)
THEN
6430 integral_value = 2.0_dp*integral_value
6431 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
6432 REAL(conjg(pw1%array(1))*pw2%array(1), kind=dp)
6435 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
6436 CALL pw1%pw_grid%para%group%sum(integral_value)
6438 CALL timestop(handle)
6440 END FUNCTION pw_integral_ab_c1d_c1d_rs
6454 SUBROUTINE pw_copy_c1d_c1d_gs (pw1, pw2)
6456 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
6457 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
6459 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
6462 INTEGER :: i, j, ng, ng1, ng2, ns
6464 CALL timeset(routinen, handle)
6466 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
6467 cpabort(
"Both grids must be either spherical or non-spherical!")
6469 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
6470 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
6471 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
6472 ng1 =
SIZE(pw1%array)
6473 ng2 =
SIZE(pw2%array)
6476 pw2%array(1:ng) = pw1%array(1:ng)
6480 pw2%array(ng + 1:ng2) = cmplx(0.0_dp, 0.0_dp, kind=dp)
6484 cpabort(
"Copies between spherical grids require compatible grids!")
6487 ng1 =
SIZE(pw1%array)
6488 ng2 =
SIZE(pw2%array)
6489 ns = 2*max(ng1, ng2)
6491 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
6492 IF (ng1 >= ng2)
THEN
6497 j = pw2%pw_grid%gidx(i)
6498 pw2%array(i) = pw1%array(j)
6507 j = pw2%pw_grid%gidx(i)
6508 pw2%array(j) = pw1%array(i)
6512 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
6513 IF (ng1 >= ng2)
THEN
6516 j = pw1%pw_grid%gidx(i)
6517 pw2%array(i) = pw1%array(j)
6524 j = pw1%pw_grid%gidx(i)
6525 pw2%array(j) = pw1%array(i)
6530 CALL pw_copy_match(pw1, pw2)
6537 pw2%array = pw1%array
6541 CALL timestop(handle)
6543 END SUBROUTINE pw_copy_c1d_c1d_gs
6550 SUBROUTINE pw_copy_to_array_c1d_c1d_gs (pw, array)
6551 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
6552 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(INOUT) :: array
6554 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
6558 CALL timeset(routinen, handle)
6561 array(:) = pw%array(:)
6564 CALL timestop(handle)
6565 END SUBROUTINE pw_copy_to_array_c1d_c1d_gs
6572 SUBROUTINE pw_copy_from_array_c1d_c1d_gs (pw, array)
6573 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
6574 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(IN) :: array
6576 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
6580 CALL timeset(routinen, handle)
6586 CALL timestop(handle)
6587 END SUBROUTINE pw_copy_from_array_c1d_c1d_gs
6605 SUBROUTINE pw_axpy_c1d_c1d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
6607 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
6608 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
6609 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
6610 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
6612 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
6615 LOGICAL :: my_allow_noncompatible_grids
6616 REAL(KIND=dp) :: my_alpha, my_beta
6617 INTEGER :: i, j, ng, ng1, ng2
6619 CALL timeset(routinen, handle)
6622 IF (
PRESENT(alpha)) my_alpha = alpha
6625 IF (
PRESENT(beta)) my_beta = beta
6627 my_allow_noncompatible_grids = .false.
6628 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
6630 IF (my_beta /= 1.0_dp)
THEN
6631 IF (my_beta == 0.0_dp)
THEN
6635 pw2%array = pw2%array*my_beta
6640 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
6642 IF (my_alpha == 1.0_dp)
THEN
6644 pw2%array = pw2%array + pw1%array
6646 ELSE IF (my_alpha /= 0.0_dp)
THEN
6648 pw2%array = pw2%array + my_alpha* pw1%array
6652 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
6654 ng1 =
SIZE(pw1%array)
6655 ng2 =
SIZE(pw2%array)
6658 IF (pw1%pw_grid%spherical)
THEN
6659 IF (my_alpha == 1.0_dp)
THEN
6662 pw2%array(i) = pw2%array(i) + pw1%array(i)
6665 ELSE IF (my_alpha /= 0.0_dp)
THEN
6668 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(i)
6672 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
6673 IF (ng1 >= ng2)
THEN
6674 IF (my_alpha == 1.0_dp)
THEN
6677 j = pw2%pw_grid%gidx(i)
6678 pw2%array(i) = pw2%array(i) + pw1%array(j)
6681 ELSE IF (my_alpha /= 0.0_dp)
THEN
6684 j = pw2%pw_grid%gidx(i)
6685 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
6690 IF (my_alpha == 1.0_dp)
THEN
6693 j = pw2%pw_grid%gidx(i)
6694 pw2%array(j) = pw2%array(j) + pw1%array(i)
6697 ELSE IF (my_alpha /= 0.0_dp)
THEN
6700 j = pw2%pw_grid%gidx(i)
6701 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
6706 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
6707 IF (ng1 >= ng2)
THEN
6708 IF (my_alpha == 1.0_dp)
THEN
6711 j = pw1%pw_grid%gidx(i)
6712 pw2%array(i) = pw2%array(i) + pw1%array(j)
6715 ELSE IF (my_alpha /= 0.0_dp)
THEN
6718 j = pw1%pw_grid%gidx(i)
6719 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
6724 IF (my_alpha == 1.0_dp)
THEN
6727 j = pw1%pw_grid%gidx(i)
6728 pw2%array(j) = pw2%array(j) + pw1%array(i)
6731 ELSE IF (my_alpha /= 0.0_dp)
THEN
6734 j = pw1%pw_grid%gidx(i)
6735 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
6741 cpabort(
"Grids not compatible")
6746 cpabort(
"Grids not compatible")
6750 CALL timestop(handle)
6752 END SUBROUTINE pw_axpy_c1d_c1d_gs
6763 SUBROUTINE pw_multiply_c1d_c1d_gs (pw_out, pw1, pw2, alpha)
6765 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw_out
6766 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
6767 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
6768 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
6770 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
6773 REAL(KIND=dp) :: my_alpha
6775 CALL timeset(routinen, handle)
6778 IF (
PRESENT(alpha)) my_alpha = alpha
6780 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
6781 cpabort(
"pw_multiply not implemented for non-identical grids!")
6783#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
6784 IF (my_alpha == 1.0_dp)
THEN
6786 pw_out%array = pw_out%array + pw1%array* pw2%array
6790 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
6794 IF (my_alpha == 1.0_dp)
THEN
6795 pw_out%array = pw_out%array + pw1%array* pw2%array
6797 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
6801 CALL timestop(handle)
6803 END SUBROUTINE pw_multiply_c1d_c1d_gs
6810 SUBROUTINE pw_multiply_with_c1d_c1d_gs (pw1, pw2)
6811 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw1
6812 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
6814 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
6818 CALL timeset(routinen, handle)
6820 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
6821 cpabort(
"Incompatible grids!")
6824 pw1%array = pw1%array* pw2%array
6827 CALL timestop(handle)
6829 END SUBROUTINE pw_multiply_with_c1d_c1d_gs
6844 FUNCTION pw_integral_ab_c1d_c1d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
6846 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
6847 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
6848 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
6849 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
6850 REAL(kind=dp) :: integral_value
6852 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_c1d_c1d_gs'
6854 INTEGER :: handle, loc_sumtype
6855 LOGICAL :: my_just_sum, my_local_only
6857 CALL timeset(routinen, handle)
6860 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
6862 my_local_only = .false.
6863 IF (
PRESENT(local_only)) my_local_only = local_only
6865 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
6866 cpabort(
"Grids incompatible")
6869 my_just_sum = .false.
6870 IF (
PRESENT(just_sum)) my_just_sum = just_sum
6877 integral_value = sum(real(conjg(pw1%array) &
6878 *pw2%array, kind=dp))
6883 integral_value = accurate_sum(real(conjg(pw1%array)*pw2%array, kind=dp))
6887 IF (.NOT. my_just_sum)
THEN
6888 integral_value = integral_value*pw1%pw_grid%vol
6891 IF (pw1%pw_grid%grid_span == halfspace)
THEN
6892 integral_value = 2.0_dp*integral_value
6893 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
6894 REAL(conjg(pw1%array(1))*pw2%array(1), kind=dp)
6897 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
6898 CALL pw1%pw_grid%para%group%sum(integral_value)
6900 CALL timestop(handle)
6902 END FUNCTION pw_integral_ab_c1d_c1d_gs
6910 FUNCTION pw_integral_a2b_c1d_c1d (pw1, pw2)
RESULT(integral_value)
6912 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
6913 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
6914 REAL(kind=dp) :: integral_value
6916 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_a2b'
6920 CALL timeset(routinen, handle)
6922 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
6923 cpabort(
"Grids incompatible")
6926 integral_value = accurate_sum(real(conjg(pw1%array)*pw2%array, kind=dp)*pw1%pw_grid%gsq)
6927 IF (pw1%pw_grid%grid_span == halfspace)
THEN
6928 integral_value = 2.0_dp*integral_value
6931 integral_value = integral_value*pw1%pw_grid%vol
6933 IF (pw1%pw_grid%para%mode == pw_mode_distributed) &
6934 CALL pw1%pw_grid%para%group%sum(integral_value)
6935 CALL timestop(handle)
6937 END FUNCTION pw_integral_a2b_c1d_c1d
6951 SUBROUTINE pw_copy_c3d_r3d_rs (pw1, pw2)
6953 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
6954 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw2
6956 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
6960 CALL timeset(routinen, handle)
6962 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
6963 cpabort(
"Both grids must be either spherical or non-spherical!")
6964 IF (pw1%pw_grid%spherical) &
6965 cpabort(
"Spherical grids only exist in reciprocal space!")
6967 IF (any(shape(pw2%array) /= shape(pw1%array))) &
6968 cpabort(
"3D grids must be compatible!")
6969 IF (pw1%pw_grid%spherical) &
6970 cpabort(
"3D grids must not be spherical!")
6972 pw2%array(:, :, :) = real(pw1%array(:, :, :), kind=dp)
6975 CALL timestop(handle)
6977 END SUBROUTINE pw_copy_c3d_r3d_rs
6984 SUBROUTINE pw_copy_to_array_c3d_r3d_rs (pw, array)
6985 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw
6986 REAL(KIND=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
6988 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
6992 CALL timeset(routinen, handle)
6995 array(:, :, :) = real(pw%array(:, :, :), kind=dp)
6998 CALL timestop(handle)
6999 END SUBROUTINE pw_copy_to_array_c3d_r3d_rs
7006 SUBROUTINE pw_copy_from_array_c3d_r3d_rs (pw, array)
7007 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw
7008 REAL(KIND=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
7010 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
7014 CALL timeset(routinen, handle)
7017 pw%array = cmplx(array, 0.0_dp, kind=dp)
7020 CALL timestop(handle)
7021 END SUBROUTINE pw_copy_from_array_c3d_r3d_rs
7039 SUBROUTINE pw_axpy_c3d_r3d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
7041 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
7042 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw2
7043 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
7044 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
7046 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
7049 LOGICAL :: my_allow_noncompatible_grids
7050 REAL(KIND=dp) :: my_alpha, my_beta
7052 CALL timeset(routinen, handle)
7055 IF (
PRESENT(alpha)) my_alpha = alpha
7058 IF (
PRESENT(beta)) my_beta = beta
7060 my_allow_noncompatible_grids = .false.
7061 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
7063 IF (my_beta /= 1.0_dp)
THEN
7064 IF (my_beta == 0.0_dp)
THEN
7068 pw2%array = pw2%array*my_beta
7073 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
7074 IF (my_alpha == 1.0_dp)
THEN
7076 pw2%array = pw2%array + real(pw1%array, kind=dp)
7078 ELSE IF (my_alpha /= 0.0_dp)
THEN
7080 pw2%array = pw2%array + my_alpha* real(pw1%array, kind=dp)
7084 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
7086 IF (any(shape(pw1%array) /= shape(pw2%array))) &
7087 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
7089 IF (my_alpha == 1.0_dp)
THEN
7091 pw2%array = pw2%array + real(pw1%array, kind=dp)
7095 pw2%array = pw2%array + my_alpha* real(pw1%array, kind=dp)
7101 cpabort(
"Grids not compatible")
7105 CALL timestop(handle)
7107 END SUBROUTINE pw_axpy_c3d_r3d_rs
7118 SUBROUTINE pw_multiply_c3d_r3d_rs (pw_out, pw1, pw2, alpha)
7120 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw_out
7121 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
7122 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw2
7123 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
7125 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
7128 REAL(KIND=dp) :: my_alpha
7130 CALL timeset(routinen, handle)
7133 IF (
PRESENT(alpha)) my_alpha = alpha
7135 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
7136 cpabort(
"pw_multiply not implemented for non-identical grids!")
7138#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
7139 IF (my_alpha == 1.0_dp)
THEN
7141 pw_out%array = pw_out%array + pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7145 pw_out%array = pw_out%array + my_alpha*pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7149 IF (my_alpha == 1.0_dp)
THEN
7150 pw_out%array = pw_out%array + pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7152 pw_out%array = pw_out%array + my_alpha*pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7156 CALL timestop(handle)
7158 END SUBROUTINE pw_multiply_c3d_r3d_rs
7165 SUBROUTINE pw_multiply_with_c3d_r3d_rs (pw1, pw2)
7166 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw1
7167 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw2
7169 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
7173 CALL timeset(routinen, handle)
7175 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
7176 cpabort(
"Incompatible grids!")
7179 pw1%array = pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7182 CALL timestop(handle)
7184 END SUBROUTINE pw_multiply_with_c3d_r3d_rs
7199 FUNCTION pw_integral_ab_c3d_r3d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
7201 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
7202 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw2
7203 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
7204 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
7205 REAL(kind=dp) :: integral_value
7207 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_c3d_r3d_rs'
7209 INTEGER :: handle, loc_sumtype
7210 LOGICAL :: my_just_sum, my_local_only
7212 CALL timeset(routinen, handle)
7215 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
7217 my_local_only = .false.
7218 IF (
PRESENT(local_only)) my_local_only = local_only
7220 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
7221 cpabort(
"Grids incompatible")
7224 my_just_sum = .false.
7225 IF (
PRESENT(just_sum)) my_just_sum = just_sum
7232 integral_value = sum(real(pw1%array, kind=dp)*pw2%array)
7237 integral_value = accurate_sum(real(pw1%array, kind=dp)*pw2%array)
7241 IF (.NOT. my_just_sum)
THEN
7242 integral_value = integral_value*pw1%pw_grid%dvol
7246 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
7247 CALL pw1%pw_grid%para%group%sum(integral_value)
7249 CALL timestop(handle)
7251 END FUNCTION pw_integral_ab_c3d_r3d_rs
7265 SUBROUTINE pw_copy_c3d_r3d_gs (pw1, pw2)
7267 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
7268 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw2
7270 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
7274 CALL timeset(routinen, handle)
7276 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
7277 cpabort(
"Both grids must be either spherical or non-spherical!")
7279 IF (any(shape(pw2%array) /= shape(pw1%array))) &
7280 cpabort(
"3D grids must be compatible!")
7281 IF (pw1%pw_grid%spherical) &
7282 cpabort(
"3D grids must not be spherical!")
7284 pw2%array(:, :, :) = real(pw1%array(:, :, :), kind=dp)
7287 CALL timestop(handle)
7289 END SUBROUTINE pw_copy_c3d_r3d_gs
7296 SUBROUTINE pw_copy_to_array_c3d_r3d_gs (pw, array)
7297 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw
7298 REAL(KIND=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
7300 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
7304 CALL timeset(routinen, handle)
7307 array(:, :, :) = real(pw%array(:, :, :), kind=dp)
7310 CALL timestop(handle)
7311 END SUBROUTINE pw_copy_to_array_c3d_r3d_gs
7318 SUBROUTINE pw_copy_from_array_c3d_r3d_gs (pw, array)
7319 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw
7320 REAL(KIND=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
7322 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
7326 CALL timeset(routinen, handle)
7329 pw%array = cmplx(array, 0.0_dp, kind=dp)
7332 CALL timestop(handle)
7333 END SUBROUTINE pw_copy_from_array_c3d_r3d_gs
7351 SUBROUTINE pw_axpy_c3d_r3d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
7353 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
7354 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw2
7355 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
7356 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
7358 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
7361 LOGICAL :: my_allow_noncompatible_grids
7362 REAL(KIND=dp) :: my_alpha, my_beta
7364 CALL timeset(routinen, handle)
7367 IF (
PRESENT(alpha)) my_alpha = alpha
7370 IF (
PRESENT(beta)) my_beta = beta
7372 my_allow_noncompatible_grids = .false.
7373 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
7375 IF (my_beta /= 1.0_dp)
THEN
7376 IF (my_beta == 0.0_dp)
THEN
7380 pw2%array = pw2%array*my_beta
7385 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
7386 IF (my_alpha == 1.0_dp)
THEN
7388 pw2%array = pw2%array + real(pw1%array, kind=dp)
7390 ELSE IF (my_alpha /= 0.0_dp)
THEN
7392 pw2%array = pw2%array + my_alpha* real(pw1%array, kind=dp)
7396 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
7398 IF (any(shape(pw1%array) /= shape(pw2%array))) &
7399 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
7401 IF (my_alpha == 1.0_dp)
THEN
7403 pw2%array = pw2%array + real(pw1%array, kind=dp)
7407 pw2%array = pw2%array + my_alpha* real(pw1%array, kind=dp)
7413 cpabort(
"Grids not compatible")
7417 CALL timestop(handle)
7419 END SUBROUTINE pw_axpy_c3d_r3d_gs
7430 SUBROUTINE pw_multiply_c3d_r3d_gs (pw_out, pw1, pw2, alpha)
7432 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw_out
7433 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
7434 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw2
7435 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
7437 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
7440 REAL(KIND=dp) :: my_alpha
7442 CALL timeset(routinen, handle)
7445 IF (
PRESENT(alpha)) my_alpha = alpha
7447 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
7448 cpabort(
"pw_multiply not implemented for non-identical grids!")
7450#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
7451 IF (my_alpha == 1.0_dp)
THEN
7453 pw_out%array = pw_out%array + pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7457 pw_out%array = pw_out%array + my_alpha*pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7461 IF (my_alpha == 1.0_dp)
THEN
7462 pw_out%array = pw_out%array + pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7464 pw_out%array = pw_out%array + my_alpha*pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7468 CALL timestop(handle)
7470 END SUBROUTINE pw_multiply_c3d_r3d_gs
7477 SUBROUTINE pw_multiply_with_c3d_r3d_gs (pw1, pw2)
7478 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw1
7479 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw2
7481 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
7485 CALL timeset(routinen, handle)
7487 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
7488 cpabort(
"Incompatible grids!")
7491 pw1%array = pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7494 CALL timestop(handle)
7496 END SUBROUTINE pw_multiply_with_c3d_r3d_gs
7511 FUNCTION pw_integral_ab_c3d_r3d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
7513 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
7514 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw2
7515 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
7516 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
7517 REAL(kind=dp) :: integral_value
7519 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_c3d_r3d_gs'
7521 INTEGER :: handle, loc_sumtype
7522 LOGICAL :: my_just_sum, my_local_only
7524 CALL timeset(routinen, handle)
7527 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
7529 my_local_only = .false.
7530 IF (
PRESENT(local_only)) my_local_only = local_only
7532 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
7533 cpabort(
"Grids incompatible")
7536 my_just_sum = .false.
7537 IF (
PRESENT(just_sum)) my_just_sum = just_sum
7544 integral_value = sum(real(pw1%array, kind=dp)*pw2%array)
7549 integral_value = accurate_sum(real(pw1%array, kind=dp)*pw2%array)
7553 IF (.NOT. my_just_sum)
THEN
7554 integral_value = integral_value*pw1%pw_grid%vol
7558 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
7559 CALL pw1%pw_grid%para%group%sum(integral_value)
7561 CALL timestop(handle)
7563 END FUNCTION pw_integral_ab_c3d_r3d_gs
7578 SUBROUTINE pw_copy_c3d_c3d_rs (pw1, pw2)
7580 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
7581 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw2
7583 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
7587 CALL timeset(routinen, handle)
7589 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
7590 cpabort(
"Both grids must be either spherical or non-spherical!")
7591 IF (pw1%pw_grid%spherical) &
7592 cpabort(
"Spherical grids only exist in reciprocal space!")
7594 IF (any(shape(pw2%array) /= shape(pw1%array))) &
7595 cpabort(
"3D grids must be compatible!")
7596 IF (pw1%pw_grid%spherical) &
7597 cpabort(
"3D grids must not be spherical!")
7599 pw2%array(:, :, :) = pw1%array(:, :, :)
7602 CALL timestop(handle)
7604 END SUBROUTINE pw_copy_c3d_c3d_rs
7611 SUBROUTINE pw_copy_to_array_c3d_c3d_rs (pw, array)
7612 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw
7613 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
7615 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
7619 CALL timeset(routinen, handle)
7622 array(:, :, :) = pw%array(:, :, :)
7625 CALL timestop(handle)
7626 END SUBROUTINE pw_copy_to_array_c3d_c3d_rs
7633 SUBROUTINE pw_copy_from_array_c3d_c3d_rs (pw, array)
7634 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw
7635 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
7637 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
7641 CALL timeset(routinen, handle)
7647 CALL timestop(handle)
7648 END SUBROUTINE pw_copy_from_array_c3d_c3d_rs
7666 SUBROUTINE pw_axpy_c3d_c3d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
7668 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
7669 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw2
7670 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
7671 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
7673 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
7676 LOGICAL :: my_allow_noncompatible_grids
7677 REAL(KIND=dp) :: my_alpha, my_beta
7679 CALL timeset(routinen, handle)
7682 IF (
PRESENT(alpha)) my_alpha = alpha
7685 IF (
PRESENT(beta)) my_beta = beta
7687 my_allow_noncompatible_grids = .false.
7688 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
7690 IF (my_beta /= 1.0_dp)
THEN
7691 IF (my_beta == 0.0_dp)
THEN
7695 pw2%array = pw2%array*my_beta
7700 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
7701 IF (my_alpha == 1.0_dp)
THEN
7703 pw2%array = pw2%array + pw1%array
7705 ELSE IF (my_alpha /= 0.0_dp)
THEN
7707 pw2%array = pw2%array + my_alpha* pw1%array
7711 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
7713 IF (any(shape(pw1%array) /= shape(pw2%array))) &
7714 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
7716 IF (my_alpha == 1.0_dp)
THEN
7718 pw2%array = pw2%array + pw1%array
7722 pw2%array = pw2%array + my_alpha* pw1%array
7728 cpabort(
"Grids not compatible")
7732 CALL timestop(handle)
7734 END SUBROUTINE pw_axpy_c3d_c3d_rs
7745 SUBROUTINE pw_multiply_c3d_c3d_rs (pw_out, pw1, pw2, alpha)
7747 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw_out
7748 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
7749 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw2
7750 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
7752 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
7755 REAL(KIND=dp) :: my_alpha
7757 CALL timeset(routinen, handle)
7760 IF (
PRESENT(alpha)) my_alpha = alpha
7762 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
7763 cpabort(
"pw_multiply not implemented for non-identical grids!")
7765#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
7766 IF (my_alpha == 1.0_dp)
THEN
7768 pw_out%array = pw_out%array + pw1%array* pw2%array
7772 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
7776 IF (my_alpha == 1.0_dp)
THEN
7777 pw_out%array = pw_out%array + pw1%array* pw2%array
7779 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
7783 CALL timestop(handle)
7785 END SUBROUTINE pw_multiply_c3d_c3d_rs
7792 SUBROUTINE pw_multiply_with_c3d_c3d_rs (pw1, pw2)
7793 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw1
7794 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw2
7796 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
7800 CALL timeset(routinen, handle)
7802 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
7803 cpabort(
"Incompatible grids!")
7806 pw1%array = pw1%array* pw2%array
7809 CALL timestop(handle)
7811 END SUBROUTINE pw_multiply_with_c3d_c3d_rs
7826 FUNCTION pw_integral_ab_c3d_c3d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
7828 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
7829 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw2
7830 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
7831 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
7832 REAL(kind=dp) :: integral_value
7834 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_c3d_c3d_rs'
7836 INTEGER :: handle, loc_sumtype
7837 LOGICAL :: my_just_sum, my_local_only
7839 CALL timeset(routinen, handle)
7842 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
7844 my_local_only = .false.
7845 IF (
PRESENT(local_only)) my_local_only = local_only
7847 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
7848 cpabort(
"Grids incompatible")
7851 my_just_sum = .false.
7852 IF (
PRESENT(just_sum)) my_just_sum = just_sum
7859 integral_value = sum(real(conjg(pw1%array) &
7860 *pw2%array, kind=dp))
7865 integral_value = accurate_sum(real(conjg(pw1%array)*pw2%array, kind=dp))
7869 IF (.NOT. my_just_sum)
THEN
7870 integral_value = integral_value*pw1%pw_grid%dvol
7874 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
7875 CALL pw1%pw_grid%para%group%sum(integral_value)
7877 CALL timestop(handle)
7879 END FUNCTION pw_integral_ab_c3d_c3d_rs
7893 SUBROUTINE pw_copy_c3d_c3d_gs (pw1, pw2)
7895 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
7896 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
7898 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
7902 CALL timeset(routinen, handle)
7904 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
7905 cpabort(
"Both grids must be either spherical or non-spherical!")
7907 IF (any(shape(pw2%array) /= shape(pw1%array))) &
7908 cpabort(
"3D grids must be compatible!")
7909 IF (pw1%pw_grid%spherical) &
7910 cpabort(
"3D grids must not be spherical!")
7912 pw2%array(:, :, :) = pw1%array(:, :, :)
7915 CALL timestop(handle)
7917 END SUBROUTINE pw_copy_c3d_c3d_gs
7924 SUBROUTINE pw_copy_to_array_c3d_c3d_gs (pw, array)
7925 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw
7926 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
7928 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
7932 CALL timeset(routinen, handle)
7935 array(:, :, :) = pw%array(:, :, :)
7938 CALL timestop(handle)
7939 END SUBROUTINE pw_copy_to_array_c3d_c3d_gs
7946 SUBROUTINE pw_copy_from_array_c3d_c3d_gs (pw, array)
7947 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw
7948 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
7950 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
7954 CALL timeset(routinen, handle)
7960 CALL timestop(handle)
7961 END SUBROUTINE pw_copy_from_array_c3d_c3d_gs
7979 SUBROUTINE pw_axpy_c3d_c3d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
7981 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
7982 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
7983 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
7984 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
7986 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
7989 LOGICAL :: my_allow_noncompatible_grids
7990 REAL(KIND=dp) :: my_alpha, my_beta
7992 CALL timeset(routinen, handle)
7995 IF (
PRESENT(alpha)) my_alpha = alpha
7998 IF (
PRESENT(beta)) my_beta = beta
8000 my_allow_noncompatible_grids = .false.
8001 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
8003 IF (my_beta /= 1.0_dp)
THEN
8004 IF (my_beta == 0.0_dp)
THEN
8008 pw2%array = pw2%array*my_beta
8013 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
8014 IF (my_alpha == 1.0_dp)
THEN
8016 pw2%array = pw2%array + pw1%array
8018 ELSE IF (my_alpha /= 0.0_dp)
THEN
8020 pw2%array = pw2%array + my_alpha* pw1%array
8024 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
8026 IF (any(shape(pw1%array) /= shape(pw2%array))) &
8027 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
8029 IF (my_alpha == 1.0_dp)
THEN
8031 pw2%array = pw2%array + pw1%array
8035 pw2%array = pw2%array + my_alpha* pw1%array
8041 cpabort(
"Grids not compatible")
8045 CALL timestop(handle)
8047 END SUBROUTINE pw_axpy_c3d_c3d_gs
8058 SUBROUTINE pw_multiply_c3d_c3d_gs (pw_out, pw1, pw2, alpha)
8060 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw_out
8061 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
8062 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw2
8063 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
8065 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
8068 REAL(KIND=dp) :: my_alpha
8070 CALL timeset(routinen, handle)
8073 IF (
PRESENT(alpha)) my_alpha = alpha
8075 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
8076 cpabort(
"pw_multiply not implemented for non-identical grids!")
8078#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
8079 IF (my_alpha == 1.0_dp)
THEN
8081 pw_out%array = pw_out%array + pw1%array* pw2%array
8085 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
8089 IF (my_alpha == 1.0_dp)
THEN
8090 pw_out%array = pw_out%array + pw1%array* pw2%array
8092 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
8096 CALL timestop(handle)
8098 END SUBROUTINE pw_multiply_c3d_c3d_gs
8105 SUBROUTINE pw_multiply_with_c3d_c3d_gs (pw1, pw2)
8106 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw1
8107 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw2
8109 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
8113 CALL timeset(routinen, handle)
8115 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
8116 cpabort(
"Incompatible grids!")
8119 pw1%array = pw1%array* pw2%array
8122 CALL timestop(handle)
8124 END SUBROUTINE pw_multiply_with_c3d_c3d_gs
8139 FUNCTION pw_integral_ab_c3d_c3d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
8141 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
8142 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw2
8143 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
8144 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
8145 REAL(kind=dp) :: integral_value
8147 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_c3d_c3d_gs'
8149 INTEGER :: handle, loc_sumtype
8150 LOGICAL :: my_just_sum, my_local_only
8152 CALL timeset(routinen, handle)
8155 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
8157 my_local_only = .false.
8158 IF (
PRESENT(local_only)) my_local_only = local_only
8160 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
8161 cpabort(
"Grids incompatible")
8164 my_just_sum = .false.
8165 IF (
PRESENT(just_sum)) my_just_sum = just_sum
8172 integral_value = sum(real(conjg(pw1%array) &
8173 *pw2%array, kind=dp))
8178 integral_value = accurate_sum(real(conjg(pw1%array)*pw2%array, kind=dp))
8182 IF (.NOT. my_just_sum)
THEN
8183 integral_value = integral_value*pw1%pw_grid%vol
8187 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
8188 CALL pw1%pw_grid%para%group%sum(integral_value)
8190 CALL timestop(handle)
8192 END FUNCTION pw_integral_ab_c3d_c3d_gs
8215 SUBROUTINE pw_gather_s_r1d_r3d_2(pw1, pw2, scale)
8217 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
8218 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw2
8219 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
8221 CALL pw_gather_s_r1d_r3d (pw2, pw1%array, scale)
8223 END SUBROUTINE pw_gather_s_r1d_r3d_2
8234 SUBROUTINE pw_gather_s_r1d_r3d (pw, c, scale)
8236 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw
8237 REAL(kind=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(IN) :: c
8238 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
8240 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_gather_s'
8242 INTEGER :: gpt, handle, l, m, n
8244 CALL timeset(routinen, handle)
8246 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
8247 ngpts =>
SIZE(pw%pw_grid%gsq), ghat => pw%pw_grid%g_hat)
8249 IF (
PRESENT(scale))
THEN
8252 l = mapl(ghat(1, gpt)) + 1
8253 m = mapm(ghat(2, gpt)) + 1
8254 n = mapn(ghat(3, gpt)) + 1
8255 pw%array(gpt) = scale* c(l, m, n)
8261 l = mapl(ghat(1, gpt)) + 1
8262 m = mapm(ghat(2, gpt)) + 1
8263 n = mapn(ghat(3, gpt)) + 1
8264 pw%array(gpt) = c(l, m, n)
8271 CALL timestop(handle)
8273 END SUBROUTINE pw_gather_s_r1d_r3d
8284 SUBROUTINE pw_scatter_s_r1d_r3d_2(pw1, pw2, scale)
8286 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
8287 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw2
8288 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
8290 CALL pw_scatter_s_r1d_r3d (pw1, pw2%array, scale)
8292 END SUBROUTINE pw_scatter_s_r1d_r3d_2
8303 SUBROUTINE pw_scatter_s_r1d_r3d (pw, c, scale)
8305 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw
8306 REAL(kind=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(INOUT) :: c
8307 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
8309 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_scatter_s'
8311 INTEGER :: gpt, handle, l, m, n
8313 CALL timeset(routinen, handle)
8315 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
8316 ghat => pw%pw_grid%g_hat, ngpts =>
SIZE(pw%pw_grid%gsq))
8319 IF (.NOT.
PRESENT(scale)) c = 0.0_dp
8321 IF (
PRESENT(scale))
THEN
8324 l = mapl(ghat(1, gpt)) + 1
8325 m = mapm(ghat(2, gpt)) + 1
8326 n = mapn(ghat(3, gpt)) + 1
8327 c(l, m, n) = scale* pw%array(gpt)
8333 l = mapl(ghat(1, gpt)) + 1
8334 m = mapm(ghat(2, gpt)) + 1
8335 n = mapn(ghat(3, gpt)) + 1
8336 c(l, m, n) = pw%array(gpt)
8343 IF (pw%pw_grid%grid_span == halfspace)
THEN
8345 associate(mapl => pw%pw_grid%mapl%neg, mapm => pw%pw_grid%mapm%neg, mapn => pw%pw_grid%mapn%neg, &
8346 ghat => pw%pw_grid%g_hat, ngpts =>
SIZE(pw%pw_grid%gsq))
8348 IF (
PRESENT(scale))
THEN
8351 l = mapl(ghat(1, gpt)) + 1
8352 m = mapm(ghat(2, gpt)) + 1
8353 n = mapn(ghat(3, gpt)) + 1
8354 c(l, m, n) = scale*( pw%array(gpt))
8360 l = mapl(ghat(1, gpt)) + 1
8361 m = mapm(ghat(2, gpt)) + 1
8362 n = mapn(ghat(3, gpt)) + 1
8363 c(l, m, n) = ( pw%array(gpt))
8372 CALL timestop(handle)
8374 END SUBROUTINE pw_scatter_s_r1d_r3d
8395 SUBROUTINE pw_gather_s_r1d_c3d_2(pw1, pw2, scale)
8397 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
8398 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw2
8399 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
8401 CALL pw_gather_s_r1d_c3d (pw2, pw1%array, scale)
8403 END SUBROUTINE pw_gather_s_r1d_c3d_2
8414 SUBROUTINE pw_gather_s_r1d_c3d (pw, c, scale)
8416 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw
8417 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(IN) :: c
8418 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
8420 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_gather_s'
8422 INTEGER :: gpt, handle, l, m, n
8424 CALL timeset(routinen, handle)
8426 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
8427 ngpts =>
SIZE(pw%pw_grid%gsq), ghat => pw%pw_grid%g_hat)
8429 IF (
PRESENT(scale))
THEN
8432 l = mapl(ghat(1, gpt)) + 1
8433 m = mapm(ghat(2, gpt)) + 1
8434 n = mapn(ghat(3, gpt)) + 1
8435 pw%array(gpt) = scale* real(c(l, m, n), kind=dp)
8441 l = mapl(ghat(1, gpt)) + 1
8442 m = mapm(ghat(2, gpt)) + 1
8443 n = mapn(ghat(3, gpt)) + 1
8444 pw%array(gpt) = real(c(l, m, n), kind=dp)
8451 CALL timestop(handle)
8453 END SUBROUTINE pw_gather_s_r1d_c3d
8464 SUBROUTINE pw_scatter_s_r1d_c3d_2(pw1, pw2, scale)
8466 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
8467 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
8468 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
8470 CALL pw_scatter_s_r1d_c3d (pw1, pw2%array, scale)
8472 END SUBROUTINE pw_scatter_s_r1d_c3d_2
8483 SUBROUTINE pw_scatter_s_r1d_c3d (pw, c, scale)
8485 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw
8486 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(INOUT) :: c
8487 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
8489 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_scatter_s'
8491 INTEGER :: gpt, handle, l, m, n
8493 CALL timeset(routinen, handle)
8495 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
8496 ghat => pw%pw_grid%g_hat, ngpts =>
SIZE(pw%pw_grid%gsq))
8499 IF (.NOT.
PRESENT(scale)) c = 0.0_dp
8501 IF (
PRESENT(scale))
THEN
8504 l = mapl(ghat(1, gpt)) + 1
8505 m = mapm(ghat(2, gpt)) + 1
8506 n = mapn(ghat(3, gpt)) + 1
8507 c(l, m, n) = scale* cmplx(pw%array(gpt), 0.0_dp, kind=dp)
8513 l = mapl(ghat(1, gpt)) + 1
8514 m = mapm(ghat(2, gpt)) + 1
8515 n = mapn(ghat(3, gpt)) + 1
8516 c(l, m, n) = cmplx(pw%array(gpt), 0.0_dp, kind=dp)
8523 IF (pw%pw_grid%grid_span == halfspace)
THEN
8525 associate(mapl => pw%pw_grid%mapl%neg, mapm => pw%pw_grid%mapm%neg, mapn => pw%pw_grid%mapn%neg, &
8526 ghat => pw%pw_grid%g_hat, ngpts =>
SIZE(pw%pw_grid%gsq))
8528 IF (
PRESENT(scale))
THEN
8531 l = mapl(ghat(1, gpt)) + 1
8532 m = mapm(ghat(2, gpt)) + 1
8533 n = mapn(ghat(3, gpt)) + 1
8534 c(l, m, n) = scale*( cmplx(pw%array(gpt), 0.0_dp, kind=dp))
8540 l = mapl(ghat(1, gpt)) + 1
8541 m = mapm(ghat(2, gpt)) + 1
8542 n = mapn(ghat(3, gpt)) + 1
8543 c(l, m, n) = ( cmplx(pw%array(gpt), 0.0_dp, kind=dp))
8552 CALL timestop(handle)
8554 END SUBROUTINE pw_scatter_s_r1d_c3d
8580 SUBROUTINE fft_wrap_pw1pw2_r3d_c1d_rs_gs (pw1, pw2, debug)
8582 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
8583 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
8584 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
8586 CHARACTER(len=*),
PARAMETER :: routineN =
'fft_wrap_pw1pw2'
8588 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
8589 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
8590 INTEGER :: handle, handle2, my_pos, nrays, &
8592 INTEGER,
DIMENSION(3) :: nloc
8593#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
8594 LOGICAL :: use_pw_gpu
8596 INTEGER,
DIMENSION(:),
POINTER :: n
8599 CALL timeset(routinen, handle2)
8600 out_unit = cp_logger_get_default_io_unit()
8601 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
8602 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
8607 IF (
PRESENT(debug))
THEN
8614 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
8615 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
8616 cpabort(
"PW grids not compatible")
8618 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
8619 cpabort(
"PW grids have not compatible MPI groups")
8623 n => pw1%pw_grid%npts
8625 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
8631 IF (test .AND. out_unit > 0)
THEN
8632 WRITE (out_unit,
'(A)')
" FFT Protocol "
8633 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
8634 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
8635 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
8638#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
8639 CALL pw_gpu_r3dc1d_3d(pw1, pw2)
8640#elif defined (__PW_FPGA)
8641 ALLOCATE (c_out(n(1), n(2), n(3)))
8644 IF (pw_fpga_init_bitstream(n) == 1)
THEN
8646#if (__PW_FPGA_SP && __PW_FPGA)
8647 CALL pw_fpga_r3dc1d_3d_sp(n, c_out)
8649 CALL pw_fpga_r3dc1d_3d_dp(n, c_out)
8651 CALL zdscal(n(1)*n(2)*n(3), 1.0_dp/pw1%pw_grid%ngpts, c_out, 1)
8652 CALL pw_gather_s_c1d_c3d(pw2, c_out)
8655 CALL fft3d(fwfft, n, c_out, debug=test)
8656 CALL pw_gather_s_c1d_c3d(pw2, c_out)
8660 ALLOCATE (c_out(n(1), n(2), n(3)))
8663 CALL fft3d(fwfft, n, c_out, debug=test)
8664 CALL pw_gather_s_c1d_c3d(pw2, c_out)
8668 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
8676 IF (test .AND. out_unit > 0)
THEN
8677 WRITE (out_unit,
'(A)')
" FFT Protocol "
8678 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
8679 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
8680 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
8683 my_pos = pw1%pw_grid%para%group%mepos
8684 nrays = pw1%pw_grid%para%nyzray(my_pos)
8685 grays => pw1%pw_grid%grays
8687#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
8689 use_pw_gpu = pw1%pw_grid%para%ray_distribution
8690 IF (use_pw_gpu)
THEN
8691 CALL pw_gpu_r3dc1d_3d_ps(pw1, pw2)
8695 nloc = pw1%pw_grid%npts_local
8696 ALLOCATE (c_in(nloc(1), nloc(2), nloc(3)))
8700 IF (pw1%pw_grid%para%ray_distribution)
THEN
8701 CALL fft3d(fwfft, n, c_in, grays, pw1%pw_grid%para%group, &
8702 pw1%pw_grid%para%yzp, pw1%pw_grid%para%nyzray, &
8703 pw1%pw_grid%para%bo, debug=test)
8705 CALL fft3d(fwfft, n, c_in, grays, pw1%pw_grid%para%group, &
8706 pw1%pw_grid%para%bo, debug=test)
8709 IF (test .AND. out_unit > 0) &
8710 WRITE (out_unit,
'(A)')
" PW_GATHER : 2d -> 1d "
8711 CALL pw_gather_p_c1d (pw2, grays)
8714#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
8719 IF (test .AND. out_unit > 0)
THEN
8720 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
8723 CALL timestop(handle)
8724 CALL timestop(handle2)
8726 END SUBROUTINE fft_wrap_pw1pw2_r3d_c1d_rs_gs
8745 SUBROUTINE fft_wrap_pw1pw2_r3d_c3d_rs_gs (pw1, pw2, debug)
8747 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
8748 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
8749 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
8751 CHARACTER(len=*),
PARAMETER :: routineN =
'fft_wrap_pw1pw2'
8753 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
8754 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
8755 INTEGER :: handle, handle2, my_pos, nrays, &
8757 INTEGER,
DIMENSION(:),
POINTER :: n
8760 CALL timeset(routinen, handle2)
8761 out_unit = cp_logger_get_default_io_unit()
8762 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
8763 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
8768 IF (
PRESENT(debug))
THEN
8775 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
8776 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
8777 cpabort(
"PW grids not compatible")
8779 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
8780 cpabort(
"PW grids have not compatible MPI groups")
8784 n => pw1%pw_grid%npts
8786 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
8792 IF (test .AND. out_unit > 0)
THEN
8793 WRITE (out_unit,
'(A)')
" FFT Protocol "
8794 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
8795 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
8796 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
8799 pw2%array = cmplx(pw1%array, 0.0_dp, kind=dp)
8801 CALL fft3d(fwfft, n, c_out, debug=test)
8803 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
8811 IF (test .AND. out_unit > 0)
THEN
8812 WRITE (out_unit,
'(A)')
" FFT Protocol "
8813 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
8814 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
8815 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
8818 my_pos = pw1%pw_grid%para%group%mepos
8819 nrays = pw1%pw_grid%para%nyzray(my_pos)
8820 grays => pw1%pw_grid%grays
8824 IF (test .AND. out_unit > 0)
THEN
8825 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
8828 CALL timestop(handle)
8829 CALL timestop(handle2)
8831 END SUBROUTINE fft_wrap_pw1pw2_r3d_c3d_rs_gs
8856 SUBROUTINE fft_wrap_pw1pw2_c1d_r3d_gs_rs (pw1, pw2, debug)
8858 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
8859 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw2
8860 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
8862 CHARACTER(len=*),
PARAMETER :: routineN =
'fft_wrap_pw1pw2'
8864 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
8865 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
8866 INTEGER :: handle, handle2, my_pos, nrays, &
8868 INTEGER,
DIMENSION(3) :: nloc
8869#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
8870 LOGICAL :: use_pw_gpu
8872 INTEGER,
DIMENSION(:),
POINTER :: n
8875 CALL timeset(routinen, handle2)
8876 out_unit = cp_logger_get_default_io_unit()
8877 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
8878 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
8883 IF (
PRESENT(debug))
THEN
8890 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
8891 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
8892 cpabort(
"PW grids not compatible")
8894 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
8895 cpabort(
"PW grids have not compatible MPI groups")
8899 n => pw1%pw_grid%npts
8901 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
8907 IF (test .AND. out_unit > 0)
THEN
8908 WRITE (out_unit,
'(A)')
" FFT Protocol "
8909 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
8910 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
8911 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
8914#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
8915 CALL pw_gpu_c1dr3d_3d(pw1, pw2)
8916#elif defined (__PW_FPGA)
8917 ALLOCATE (c_out(n(1), n(2), n(3)))
8920 IF (pw_fpga_init_bitstream(n) == 1)
THEN
8921 CALL pw_scatter_s_c1d_c3d(pw1, c_out)
8923#if (__PW_FPGA_SP && __PW_FPGA)
8924 CALL pw_fpga_c1dr3d_3d_sp(n, c_out)
8926 CALL pw_fpga_c1dr3d_3d_dp(n, c_out)
8928 CALL zdscal(n(1)*n(2)*n(3), 1.0_dp, c_out, 1)
8932 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" PW_SCATTER : 3d -> 1d "
8933 CALL pw_scatter_s_c1d_c3d(pw1, c_out)
8935 CALL fft3d(bwfft, n, c_out, debug=test)
8937 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" REAL part "
8942 ALLOCATE (c_out(n(1), n(2), n(3)))
8943 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" PW_SCATTER : 3d -> 1d "
8944 CALL pw_scatter_s_c1d_c3d(pw1, c_out)
8946 CALL fft3d(bwfft, n, c_out, debug=test)
8948 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" REAL part "
8953 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
8961 IF (test .AND. out_unit > 0)
THEN
8962 WRITE (out_unit,
'(A)')
" FFT Protocol "
8963 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
8964 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
8965 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
8968 my_pos = pw1%pw_grid%para%group%mepos
8969 nrays = pw1%pw_grid%para%nyzray(my_pos)
8970 grays => pw1%pw_grid%grays
8972#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
8974 use_pw_gpu = pw1%pw_grid%para%ray_distribution
8975 IF (use_pw_gpu)
THEN
8976 CALL pw_gpu_c1dr3d_3d_ps(pw1, pw2)
8980 IF (test .AND. out_unit > 0) &
8981 WRITE (out_unit,
'(A)')
" PW_SCATTER : 2d -> 1d "
8983 CALL pw_scatter_p_c1d (pw1, grays)
8984 nloc = pw2%pw_grid%npts_local
8985 ALLOCATE (c_in(nloc(1), nloc(2), nloc(3)))
8987 IF (pw1%pw_grid%para%ray_distribution)
THEN
8988 CALL fft3d(bwfft, n, c_in, grays, pw1%pw_grid%para%group, &
8989 pw1%pw_grid%para%yzp, pw1%pw_grid%para%nyzray, &
8990 pw1%pw_grid%para%bo, debug=test)
8992 CALL fft3d(bwfft, n, c_in, grays, pw1%pw_grid%para%group, &
8993 pw1%pw_grid%para%bo, debug=test)
8996 IF (test .AND. out_unit > 0) &
8997 WRITE (out_unit,
'(A)')
" Real part "
9000#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
9005 IF (test .AND. out_unit > 0)
THEN
9006 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9009 CALL timestop(handle)
9010 CALL timestop(handle2)
9012 END SUBROUTINE fft_wrap_pw1pw2_c1d_r3d_gs_rs
9025 SUBROUTINE pw_gather_s_c1d_r3d_2(pw1, pw2, scale)
9027 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
9028 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
9029 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
9031 CALL pw_gather_s_c1d_r3d (pw2, pw1%array, scale)
9033 END SUBROUTINE pw_gather_s_c1d_r3d_2
9044 SUBROUTINE pw_gather_s_c1d_r3d (pw, c, scale)
9046 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
9047 REAL(kind=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(IN) :: c
9048 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
9050 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_gather_s'
9052 INTEGER :: gpt, handle, l, m, n
9054 CALL timeset(routinen, handle)
9056 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
9057 ngpts =>
SIZE(pw%pw_grid%gsq), ghat => pw%pw_grid%g_hat)
9059 IF (
PRESENT(scale))
THEN
9062 l = mapl(ghat(1, gpt)) + 1
9063 m = mapm(ghat(2, gpt)) + 1
9064 n = mapn(ghat(3, gpt)) + 1
9065 pw%array(gpt) = scale* cmplx(c(l, m, n), 0.0_dp, kind=dp)
9071 l = mapl(ghat(1, gpt)) + 1
9072 m = mapm(ghat(2, gpt)) + 1
9073 n = mapn(ghat(3, gpt)) + 1
9074 pw%array(gpt) = cmplx(c(l, m, n), 0.0_dp, kind=dp)
9081 CALL timestop(handle)
9083 END SUBROUTINE pw_gather_s_c1d_r3d
9094 SUBROUTINE pw_scatter_s_c1d_r3d_2(pw1, pw2, scale)
9096 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
9097 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw2
9098 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
9100 CALL pw_scatter_s_c1d_r3d (pw1, pw2%array, scale)
9102 END SUBROUTINE pw_scatter_s_c1d_r3d_2
9113 SUBROUTINE pw_scatter_s_c1d_r3d (pw, c, scale)
9115 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
9116 REAL(kind=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(INOUT) :: c
9117 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
9119 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_scatter_s'
9121 INTEGER :: gpt, handle, l, m, n
9123 CALL timeset(routinen, handle)
9125 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
9126 ghat => pw%pw_grid%g_hat, ngpts =>
SIZE(pw%pw_grid%gsq))
9129 IF (.NOT.
PRESENT(scale)) c = 0.0_dp
9131 IF (
PRESENT(scale))
THEN
9134 l = mapl(ghat(1, gpt)) + 1
9135 m = mapm(ghat(2, gpt)) + 1
9136 n = mapn(ghat(3, gpt)) + 1
9137 c(l, m, n) = scale* real(pw%array(gpt), kind=dp)
9143 l = mapl(ghat(1, gpt)) + 1
9144 m = mapm(ghat(2, gpt)) + 1
9145 n = mapn(ghat(3, gpt)) + 1
9146 c(l, m, n) = real(pw%array(gpt), kind=dp)
9153 IF (pw%pw_grid%grid_span == halfspace)
THEN
9155 associate(mapl => pw%pw_grid%mapl%neg, mapm => pw%pw_grid%mapm%neg, mapn => pw%pw_grid%mapn%neg, &
9156 ghat => pw%pw_grid%g_hat, ngpts =>
SIZE(pw%pw_grid%gsq))
9158 IF (
PRESENT(scale))
THEN
9161 l = mapl(ghat(1, gpt)) + 1
9162 m = mapm(ghat(2, gpt)) + 1
9163 n = mapn(ghat(3, gpt)) + 1
9164 c(l, m, n) = scale*( real(pw%array(gpt), kind=dp))
9170 l = mapl(ghat(1, gpt)) + 1
9171 m = mapm(ghat(2, gpt)) + 1
9172 n = mapn(ghat(3, gpt)) + 1
9173 c(l, m, n) = ( real(pw%array(gpt), kind=dp))
9182 CALL timestop(handle)
9184 END SUBROUTINE pw_scatter_s_c1d_r3d
9206 SUBROUTINE fft_wrap_pw1pw2_c1d_c3d_gs_rs (pw1, pw2, debug)
9208 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
9209 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw2
9210 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
9212 CHARACTER(len=*),
PARAMETER :: routineN =
'fft_wrap_pw1pw2'
9214 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
9215 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
9216 INTEGER :: handle, handle2, my_pos, nrays, &
9218 INTEGER,
DIMENSION(:),
POINTER :: n
9221 CALL timeset(routinen, handle2)
9222 out_unit = cp_logger_get_default_io_unit()
9223 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
9224 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
9229 IF (
PRESENT(debug))
THEN
9236 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
9237 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
9238 cpabort(
"PW grids not compatible")
9240 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
9241 cpabort(
"PW grids have not compatible MPI groups")
9245 n => pw1%pw_grid%npts
9247 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
9253 IF (test .AND. out_unit > 0)
THEN
9254 WRITE (out_unit,
'(A)')
" FFT Protocol "
9255 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
9256 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
9257 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
9261 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" PW_SCATTER : 3d -> 1d "
9262 CALL pw_scatter_s_c1d_c3d(pw1, c_out)
9263 CALL fft3d(bwfft, n, c_out, debug=test)
9265 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9273 IF (test .AND. out_unit > 0)
THEN
9274 WRITE (out_unit,
'(A)')
" FFT Protocol "
9275 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
9276 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
9277 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
9280 my_pos = pw1%pw_grid%para%group%mepos
9281 nrays = pw1%pw_grid%para%nyzray(my_pos)
9282 grays => pw1%pw_grid%grays
9285 IF (test .AND. out_unit > 0) &
9286 WRITE (out_unit,
'(A)')
" PW_SCATTER : 2d -> 1d "
9288 CALL pw_scatter_p_c1d (pw1, grays)
9291 IF (pw1%pw_grid%para%ray_distribution)
THEN
9292 CALL fft3d(bwfft, n, c_in, grays, pw1%pw_grid%para%group, &
9293 pw1%pw_grid%para%yzp, pw1%pw_grid%para%nyzray, &
9294 pw1%pw_grid%para%bo, debug=test)
9296 CALL fft3d(bwfft, n, c_in, grays, pw1%pw_grid%para%group, &
9297 pw1%pw_grid%para%bo, debug=test)
9302 IF (test .AND. out_unit > 0)
THEN
9303 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9306 CALL timestop(handle)
9307 CALL timestop(handle2)
9309 END SUBROUTINE fft_wrap_pw1pw2_c1d_c3d_gs_rs
9322 SUBROUTINE pw_gather_s_c1d_c3d_2(pw1, pw2, scale)
9324 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
9325 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
9326 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
9328 CALL pw_gather_s_c1d_c3d (pw2, pw1%array, scale)
9330 END SUBROUTINE pw_gather_s_c1d_c3d_2
9341 SUBROUTINE pw_gather_s_c1d_c3d (pw, c, scale)
9343 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
9344 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(IN) :: c
9345 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
9347 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_gather_s'
9349 INTEGER :: gpt, handle, l, m, n
9351 CALL timeset(routinen, handle)
9353 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
9354 ngpts =>
SIZE(pw%pw_grid%gsq), ghat => pw%pw_grid%g_hat)
9356 IF (
PRESENT(scale))
THEN
9359 l = mapl(ghat(1, gpt)) + 1
9360 m = mapm(ghat(2, gpt)) + 1
9361 n = mapn(ghat(3, gpt)) + 1
9362 pw%array(gpt) = scale* c(l, m, n)
9368 l = mapl(ghat(1, gpt)) + 1
9369 m = mapm(ghat(2, gpt)) + 1
9370 n = mapn(ghat(3, gpt)) + 1
9371 pw%array(gpt) = c(l, m, n)
9378 CALL timestop(handle)
9380 END SUBROUTINE pw_gather_s_c1d_c3d
9391 SUBROUTINE pw_scatter_s_c1d_c3d_2(pw1, pw2, scale)
9393 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
9394 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
9395 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
9397 CALL pw_scatter_s_c1d_c3d (pw1, pw2%array, scale)
9399 END SUBROUTINE pw_scatter_s_c1d_c3d_2
9410 SUBROUTINE pw_scatter_s_c1d_c3d (pw, c, scale)
9412 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
9413 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(INOUT) :: c
9414 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
9416 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_scatter_s'
9418 INTEGER :: gpt, handle, l, m, n
9420 CALL timeset(routinen, handle)
9422 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
9423 ghat => pw%pw_grid%g_hat, ngpts =>
SIZE(pw%pw_grid%gsq))
9426 IF (.NOT.
PRESENT(scale)) c = 0.0_dp
9428 IF (
PRESENT(scale))
THEN
9431 l = mapl(ghat(1, gpt)) + 1
9432 m = mapm(ghat(2, gpt)) + 1
9433 n = mapn(ghat(3, gpt)) + 1
9434 c(l, m, n) = scale* pw%array(gpt)
9440 l = mapl(ghat(1, gpt)) + 1
9441 m = mapm(ghat(2, gpt)) + 1
9442 n = mapn(ghat(3, gpt)) + 1
9443 c(l, m, n) = pw%array(gpt)
9450 IF (pw%pw_grid%grid_span == halfspace)
THEN
9452 associate(mapl => pw%pw_grid%mapl%neg, mapm => pw%pw_grid%mapm%neg, mapn => pw%pw_grid%mapn%neg, &
9453 ghat => pw%pw_grid%g_hat, ngpts =>
SIZE(pw%pw_grid%gsq))
9455 IF (
PRESENT(scale))
THEN
9458 l = mapl(ghat(1, gpt)) + 1
9459 m = mapm(ghat(2, gpt)) + 1
9460 n = mapn(ghat(3, gpt)) + 1
9461 c(l, m, n) = scale*conjg( pw%array(gpt))
9467 l = mapl(ghat(1, gpt)) + 1
9468 m = mapm(ghat(2, gpt)) + 1
9469 n = mapn(ghat(3, gpt)) + 1
9470 c(l, m, n) = conjg( pw%array(gpt))
9479 CALL timestop(handle)
9481 END SUBROUTINE pw_scatter_s_c1d_c3d
9503 SUBROUTINE fft_wrap_pw1pw2_c3d_r3d_gs_rs (pw1, pw2, debug)
9505 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
9506 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw2
9507 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
9509 CHARACTER(len=*),
PARAMETER :: routineN =
'fft_wrap_pw1pw2'
9511 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
9512 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
9513 INTEGER :: handle, handle2, my_pos, nrays, &
9515 INTEGER,
DIMENSION(:),
POINTER :: n
9518 CALL timeset(routinen, handle2)
9519 out_unit = cp_logger_get_default_io_unit()
9520 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
9521 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
9526 IF (
PRESENT(debug))
THEN
9533 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
9534 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
9535 cpabort(
"PW grids not compatible")
9537 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
9538 cpabort(
"PW grids have not compatible MPI groups")
9542 n => pw1%pw_grid%npts
9544 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
9550 IF (test .AND. out_unit > 0)
THEN
9551 WRITE (out_unit,
'(A)')
" FFT Protocol "
9552 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
9553 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
9554 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
9558 ALLOCATE (c_out(n(1), n(2), n(3)))
9559 CALL fft3d(bwfft, n, c_in, c_out, debug=test)
9561 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" REAL part "
9562 pw2%array = real(c_out, kind=dp)
9565 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9573 IF (test .AND. out_unit > 0)
THEN
9574 WRITE (out_unit,
'(A)')
" FFT Protocol "
9575 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
9576 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
9577 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
9580 my_pos = pw1%pw_grid%para%group%mepos
9581 nrays = pw1%pw_grid%para%nyzray(my_pos)
9582 grays => pw1%pw_grid%grays
9586 IF (test .AND. out_unit > 0)
THEN
9587 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9590 CALL timestop(handle)
9591 CALL timestop(handle2)
9593 END SUBROUTINE fft_wrap_pw1pw2_c3d_r3d_gs_rs
9611 SUBROUTINE fft_wrap_pw1pw2_c3d_c1d_rs_gs (pw1, pw2, debug)
9613 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
9614 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
9615 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
9617 CHARACTER(len=*),
PARAMETER :: routineN =
'fft_wrap_pw1pw2'
9619 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
9620 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
9621 INTEGER :: handle, handle2, my_pos, nrays, &
9623 INTEGER,
DIMENSION(:),
POINTER :: n
9626 CALL timeset(routinen, handle2)
9627 out_unit = cp_logger_get_default_io_unit()
9628 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
9629 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
9634 IF (
PRESENT(debug))
THEN
9641 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
9642 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
9643 cpabort(
"PW grids not compatible")
9645 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
9646 cpabort(
"PW grids have not compatible MPI groups")
9650 n => pw1%pw_grid%npts
9652 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
9658 IF (test .AND. out_unit > 0)
THEN
9659 WRITE (out_unit,
'(A)')
" FFT Protocol "
9660 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
9661 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
9662 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
9666 ALLOCATE (c_out(n(1), n(2), n(3)))
9668 CALL fft3d(fwfft, n, c_in, c_out, debug=test)
9670 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" PW_GATHER : 3d -> 1d "
9671 CALL pw_gather_s_c1d_c3d(pw2, c_out)
9674 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9682 IF (test .AND. out_unit > 0)
THEN
9683 WRITE (out_unit,
'(A)')
" FFT Protocol "
9684 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
9685 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
9686 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
9689 my_pos = pw1%pw_grid%para%group%mepos
9690 nrays = pw1%pw_grid%para%nyzray(my_pos)
9691 grays => pw1%pw_grid%grays
9697 IF (pw1%pw_grid%para%ray_distribution)
THEN
9698 CALL fft3d(fwfft, n, c_in, grays, pw1%pw_grid%para%group, &
9699 pw1%pw_grid%para%yzp, pw1%pw_grid%para%nyzray, &
9700 pw1%pw_grid%para%bo, debug=test)
9702 CALL fft3d(fwfft, n, c_in, grays, pw1%pw_grid%para%group, &
9703 pw1%pw_grid%para%bo, debug=test)
9706 IF (test .AND. out_unit > 0) &
9707 WRITE (out_unit,
'(A)')
" PW_GATHER : 2d -> 1d "
9708 CALL pw_gather_p_c1d (pw2, grays)
9711 IF (test .AND. out_unit > 0)
THEN
9712 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9715 CALL timestop(handle)
9716 CALL timestop(handle2)
9718 END SUBROUTINE fft_wrap_pw1pw2_c3d_c1d_rs_gs
9737 SUBROUTINE fft_wrap_pw1pw2_c3d_c3d_rs_gs (pw1, pw2, debug)
9739 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
9740 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
9741 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
9743 CHARACTER(len=*),
PARAMETER :: routineN =
'fft_wrap_pw1pw2'
9745 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
9746 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
9747 INTEGER :: handle, handle2, my_pos, nrays, &
9749 INTEGER,
DIMENSION(:),
POINTER :: n
9752 CALL timeset(routinen, handle2)
9753 out_unit = cp_logger_get_default_io_unit()
9754 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
9755 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
9760 IF (
PRESENT(debug))
THEN
9767 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
9768 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
9769 cpabort(
"PW grids not compatible")
9771 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
9772 cpabort(
"PW grids have not compatible MPI groups")
9776 n => pw1%pw_grid%npts
9778 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
9784 IF (test .AND. out_unit > 0)
THEN
9785 WRITE (out_unit,
'(A)')
" FFT Protocol "
9786 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
9787 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
9788 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
9793 CALL fft3d(fwfft, n, c_in, c_out, debug=test)
9795 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9803 IF (test .AND. out_unit > 0)
THEN
9804 WRITE (out_unit,
'(A)')
" FFT Protocol "
9805 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
9806 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
9807 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
9810 my_pos = pw1%pw_grid%para%group%mepos
9811 nrays = pw1%pw_grid%para%nyzray(my_pos)
9812 grays => pw1%pw_grid%grays
9816 IF (test .AND. out_unit > 0)
THEN
9817 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9820 CALL timestop(handle)
9821 CALL timestop(handle2)
9823 END SUBROUTINE fft_wrap_pw1pw2_c3d_c3d_rs_gs
9838 SUBROUTINE fft_wrap_pw1pw2_c3d_c3d_gs_rs (pw1, pw2, debug)
9840 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
9841 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw2
9842 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
9844 CHARACTER(len=*),
PARAMETER :: routineN =
'fft_wrap_pw1pw2'
9846 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
9847 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
9848 INTEGER :: handle, handle2, my_pos, nrays, &
9850 INTEGER,
DIMENSION(:),
POINTER :: n
9853 CALL timeset(routinen, handle2)
9854 out_unit = cp_logger_get_default_io_unit()
9855 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
9856 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
9861 IF (
PRESENT(debug))
THEN
9868 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
9869 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
9870 cpabort(
"PW grids not compatible")
9872 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
9873 cpabort(
"PW grids have not compatible MPI groups")
9877 n => pw1%pw_grid%npts
9879 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
9885 IF (test .AND. out_unit > 0)
THEN
9886 WRITE (out_unit,
'(A)')
" FFT Protocol "
9887 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
9888 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
9889 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
9894 CALL fft3d(bwfft, n, c_in, c_out, debug=test)
9896 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9904 IF (test .AND. out_unit > 0)
THEN
9905 WRITE (out_unit,
'(A)')
" FFT Protocol "
9906 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
9907 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
9908 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
9911 my_pos = pw1%pw_grid%para%group%mepos
9912 nrays = pw1%pw_grid%para%nyzray(my_pos)
9913 grays => pw1%pw_grid%grays
9917 IF (test .AND. out_unit > 0)
THEN
9918 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9921 CALL timestop(handle)
9922 CALL timestop(handle2)
9924 END SUBROUTINE fft_wrap_pw1pw2_c3d_c3d_gs_rs
9943 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
9944 REAL(kind=dp),
INTENT(IN) :: omega
9946 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_gauss_damp'
9948 INTEGER :: handle, i
9949 REAL(kind=dp) :: omega_2, tmp
9951 CALL timeset(routinen, handle)
9952 cpassert(omega >= 0)
9954 omega_2 = omega*omega
9955 omega_2 = 0.25_dp/omega_2
9958 DO i = 1,
SIZE(pw%array)
9959 tmp = exp(-pw%pw_grid%gsq(i)*omega_2)
9960 pw%array(i) = pw%array(i)*tmp
9964 CALL timestop(handle)
9977 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
9978 REAL(kind=dp),
INTENT(IN) :: omega
9980 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_log_deriv_gauss'
9982 INTEGER :: handle, i
9983 REAL(kind=dp) :: omega_2, tmp
9985 CALL timeset(routinen, handle)
9986 cpassert(omega >= 0)
9988 omega_2 = omega*omega
9989 omega_2 = 0.25_dp/omega_2
9992 DO i = 1,
SIZE(pw%array)
9993 tmp = 1.0_dp + omega_2*pw%pw_grid%gsq(i)
9994 pw%array(i) = pw%array(i)*tmp
9998 CALL timestop(handle)
10016 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
10017 REAL(kind=dp),
INTENT(IN) :: omega
10019 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_compl_gauss_damp'
10021 INTEGER :: cnt, handle, i
10022 REAL(kind=dp) :: omega_2, tmp, tmp2
10024 CALL timeset(routinen, handle)
10026 omega_2 = omega*omega
10027 omega_2 = 0.25_dp/omega_2
10029 cnt =
SIZE(pw%array)
10033 tmp = -omega_2*pw%pw_grid%gsq(i)
10034 IF (abs(tmp) > 1.0e-5_dp)
THEN
10035 tmp2 = 1.0_dp - exp(tmp)
10037 tmp2 = tmp + 0.5_dp*tmp*(tmp + (1.0_dp/3.0_dp)*tmp**2)
10039 pw%array(i) = pw%array(i)*tmp2
10043 CALL timestop(handle)
10055 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
10056 REAL(kind=dp),
INTENT(IN) :: omega
10058 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_log_deriv_compl_gauss'
10060 INTEGER :: handle, i
10061 REAL(kind=dp) :: omega_2, tmp, tmp2
10063 CALL timeset(routinen, handle)
10065 omega_2 = omega*omega
10066 omega_2 = 0.25_dp/omega_2
10070 DO i = 1,
SIZE(pw%array)
10071 tmp = omega_2*pw%pw_grid%gsq(i)
10073 IF (abs(tmp) >= 0.003_dp)
THEN
10074 tmp2 = 1.0_dp - tmp*exp(-tmp)/(1.0_dp - exp(-tmp))
10076 tmp2 = 0.5_dp*tmp - tmp**2/12.0_dp
10078 pw%array(i) = pw%array(i)*tmp2
10082 CALL timestop(handle)
10103 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
10104 REAL(kind=dp),
INTENT(IN) :: omega, scale_coul, scale_long
10106 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_gauss_damp_mix'
10108 INTEGER :: handle, i
10109 REAL(kind=dp) :: omega_2, tmp
10111 CALL timeset(routinen, handle)
10113 omega_2 = omega*omega
10114 omega_2 = 0.25_dp/omega_2
10117 DO i = 1,
SIZE(pw%array)
10118 tmp = scale_coul + scale_long*exp(-pw%pw_grid%gsq(i)*omega_2)
10119 pw%array(i) = pw%array(i)*tmp
10123 CALL timestop(handle)
10138 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
10139 REAL(kind=dp),
INTENT(IN) :: omega, scale_coul, scale_long
10141 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_log_deriv_mix_cl'
10143 INTEGER :: handle, i
10144 REAL(kind=dp) :: omega_2, tmp, tmp2
10146 CALL timeset(routinen, handle)
10148 omega_2 = omega*omega
10149 omega_2 = 0.25_dp/omega_2
10153 DO i = 1,
SIZE(pw%array)
10154 tmp = omega_2*pw%pw_grid%gsq(i)
10155 tmp2 = 1.0_dp + scale_long*tmp*exp(-tmp)/(scale_coul + scale_long*exp(-tmp))
10156 pw%array(i) = pw%array(i)*tmp2
10160 CALL timestop(handle)
10179 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
10180 REAL(kind=dp),
INTENT(IN) :: rcutoff
10182 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_truncated'
10184 INTEGER :: handle, i
10185 REAL(kind=dp) :: tmp, tmp2
10187 CALL timeset(routinen, handle)
10188 cpassert(rcutoff >= 0)
10191 DO i = 1,
SIZE(pw%array)
10192 tmp = sqrt(pw%pw_grid%gsq(i))*rcutoff
10193 IF (tmp >= 0.005_dp)
THEN
10194 tmp2 = 1.0_dp - cos(tmp)
10196 tmp2 = tmp**2/2.0_dp*(1.0 - tmp**2/12.0_dp)
10198 pw%array(i) = pw%array(i)*tmp2
10202 CALL timestop(handle)
10215 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
10216 REAL(kind=dp),
INTENT(IN) :: rcutoff
10218 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_log_deriv_trunc'
10220 INTEGER :: handle, i
10221 REAL(kind=dp) :: rchalf, tmp, tmp2
10223 CALL timeset(routinen, handle)
10224 cpassert(rcutoff >= 0)
10226 rchalf = 0.5_dp*rcutoff
10229 DO i = 1,
SIZE(pw%array)
10230 tmp = rchalf*sqrt(pw%pw_grid%gsq(i))
10232 IF (abs(tmp) >= 0.0001_dp)
THEN
10233 tmp2 = 1.0_dp - tmp/tan(tmp)
10235 tmp2 = tmp**2*(1.0_dp + tmp**2/15.0_dp)/3.0_dp
10237 pw%array(i) = pw%array(i)*tmp2
10241 CALL timestop(handle)
10257 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
10258 INTEGER,
DIMENSION(3),
INTENT(IN) :: n
10260 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_derive'
10262 COMPLEX(KIND=dp) :: im
10263 INTEGER :: handle, m,
idx, idir
10265 CALL timeset(routinen, handle)
10268 cpabort(
"Nonnegative exponents are not supported!")
10271 im = cmplx(0.0_dp, 1.0_dp, kind=dp)**m
10274 IF (n(idir) == 1)
THEN
10276 DO idx = 1,
SIZE(pw%array)
10277 pw%array(
idx) = pw%array(
idx)*pw%pw_grid%g(idir,
idx)
10280 ELSE IF (n(idir) > 1)
THEN
10282 DO idx = 1,
SIZE(pw%array)
10283 pw%array(
idx) = pw%array(
idx)*pw%pw_grid%g(idir,
idx)**n(idir)
10291 IF (abs(real(im, kind=dp) - 1.0_dp) > 1.0e-10_dp)
THEN
10293 pw%array(:) = im*pw%array(:)
10297 CALL timestop(handle)
10312 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
10314 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_laplace'
10318 CALL timeset(routinen, handle)
10321 pw%array(:) = -pw%array(:)*pw%pw_grid%gsq(:)
10324 CALL timestop(handle)
10341 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw, pwdr2
10342 INTEGER,
INTENT(IN) :: i, j
10344 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_dr2'
10346 INTEGER :: cnt, handle, ig
10347 REAL(kind=dp) :: gg, o3
10349 CALL timeset(routinen, handle)
10353 cnt =
SIZE(pw%array)
10358 gg = pw%pw_grid%g(i, ig)**2 - o3*pw%pw_grid%gsq(ig)
10359 pwdr2%array(ig) = gg*pw%array(ig)
10365 pwdr2%array(ig) = pw%array(ig)*(pw%pw_grid%g(i, ig)*pw%pw_grid%g(j, ig))
10370 CALL timestop(handle)
10389 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw, pwdr2_gg
10390 INTEGER,
INTENT(IN) :: i, j
10392 INTEGER :: cnt, handle, ig
10393 REAL(kind=dp) :: gg, o3
10394 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_dr2_gg'
10396 CALL timeset(routinen, handle)
10400 cnt =
SIZE(pw%array)
10404 DO ig = pw%pw_grid%first_gne0, cnt
10405 gg = pw%pw_grid%g(i, ig)**2 - o3*pw%pw_grid%gsq(ig)
10406 pwdr2_gg%array(ig) = gg/pw%pw_grid%gsq(ig)*pw%array(ig)
10411 DO ig = pw%pw_grid%first_gne0, cnt
10412 pwdr2_gg%array(ig) = pw%array(ig)*((pw%pw_grid%g(i, ig)*pw%pw_grid%g(j, ig)) &
10413 /pw%pw_grid%gsq(ig))
10418 IF (pw%pw_grid%have_g0) pwdr2_gg%array(1) = 0.0_dp
10420 CALL timestop(handle)
10437 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
10438 REAL(kind=dp),
INTENT(IN) :: ecut, sigma
10440 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_smoothing'
10442 INTEGER :: cnt, handle, ig
10443 REAL(kind=dp) :: arg, f
10445 CALL timeset(routinen, handle)
10447 cnt =
SIZE(pw%array)
10451 arg = (ecut - pw%pw_grid%gsq(ig))/sigma
10452 f = exp(arg)/(1 + exp(arg))
10453 pw%array(ig) = f*pw%array(ig)
10457 CALL timestop(handle)
10467 ELEMENTAL FUNCTION pw_compatible(grida, gridb)
RESULT(compat)
10469 TYPE(pw_grid_type),
INTENT(IN) :: grida, gridb
10474 IF (grida%id_nr == gridb%id_nr)
THEN
10476 ELSE IF (grida%reference == gridb%id_nr)
THEN
10478 ELSE IF (gridb%reference == grida%id_nr)
THEN
10482 END FUNCTION pw_compatible
10495 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: sf
10496 REAL(kind=dp),
DIMENSION(:),
INTENT(IN) :: r
10498 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_structure_factor'
10500 INTEGER :: cnt, handle, ig
10501 REAL(kind=dp) :: arg
10503 CALL timeset(routinen, handle)
10505 cnt =
SIZE(sf%array)
10509 arg = dot_product(sf%pw_grid%g(:, ig), r)
10510 sf%array(ig) = cmplx(cos(arg), -sin(arg), kind=dp)
10514 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 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.