61#include "../base/base_uses.f90"
78 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'pw_methods'
79 LOGICAL,
PARAMETER,
PRIVATE :: debug_this_module = .false.
84 MODULE PROCEDURE pw_zero_r1d_rs
85 MODULE PROCEDURE pw_zero_r3d_rs
86 MODULE PROCEDURE pw_zero_c1d_rs
87 MODULE PROCEDURE pw_zero_c3d_rs
88 MODULE PROCEDURE pw_zero_r1d_gs
89 MODULE PROCEDURE pw_zero_r3d_gs
90 MODULE PROCEDURE pw_zero_c1d_gs
91 MODULE PROCEDURE pw_zero_c3d_gs
95 MODULE PROCEDURE pw_scale_r1d_rs
96 MODULE PROCEDURE pw_scale_r3d_rs
97 MODULE PROCEDURE pw_scale_c1d_rs
98 MODULE PROCEDURE pw_scale_c3d_rs
99 MODULE PROCEDURE pw_scale_r1d_gs
100 MODULE PROCEDURE pw_scale_r3d_gs
101 MODULE PROCEDURE pw_scale_c1d_gs
102 MODULE PROCEDURE pw_scale_c3d_gs
106 MODULE PROCEDURE pw_write_r1d_rs
107 MODULE PROCEDURE pw_write_r3d_rs
108 MODULE PROCEDURE pw_write_c1d_rs
109 MODULE PROCEDURE pw_write_c3d_rs
110 MODULE PROCEDURE pw_write_r1d_gs
111 MODULE PROCEDURE pw_write_r3d_gs
112 MODULE PROCEDURE pw_write_c1d_gs
113 MODULE PROCEDURE pw_write_c3d_gs
117 MODULE PROCEDURE pw_integrate_function_r1d_rs
118 MODULE PROCEDURE pw_integrate_function_r3d_rs
119 MODULE PROCEDURE pw_integrate_function_c1d_rs
120 MODULE PROCEDURE pw_integrate_function_c3d_rs
121 MODULE PROCEDURE pw_integrate_function_r1d_gs
122 MODULE PROCEDURE pw_integrate_function_r3d_gs
123 MODULE PROCEDURE pw_integrate_function_c1d_gs
124 MODULE PROCEDURE pw_integrate_function_c3d_gs
128 MODULE PROCEDURE pw_set_value_r1d_rs
129 MODULE PROCEDURE pw_zero_r1d_rs
130 MODULE PROCEDURE pw_set_value_r3d_rs
131 MODULE PROCEDURE pw_zero_r3d_rs
132 MODULE PROCEDURE pw_set_value_c1d_rs
133 MODULE PROCEDURE pw_zero_c1d_rs
134 MODULE PROCEDURE pw_set_value_c3d_rs
135 MODULE PROCEDURE pw_zero_c3d_rs
136 MODULE PROCEDURE pw_set_value_r1d_gs
137 MODULE PROCEDURE pw_zero_r1d_gs
138 MODULE PROCEDURE pw_set_value_r3d_gs
139 MODULE PROCEDURE pw_zero_r3d_gs
140 MODULE PROCEDURE pw_set_value_c1d_gs
141 MODULE PROCEDURE pw_zero_c1d_gs
142 MODULE PROCEDURE pw_set_value_c3d_gs
143 MODULE PROCEDURE pw_zero_c3d_gs
147 MODULE PROCEDURE pw_copy_r1d_r1d_rs
148 MODULE PROCEDURE pw_copy_r1d_c1d_rs
149 MODULE PROCEDURE pw_copy_r3d_r3d_rs
150 MODULE PROCEDURE pw_copy_r3d_c3d_rs
151 MODULE PROCEDURE pw_copy_c1d_r1d_rs
152 MODULE PROCEDURE pw_copy_c1d_c1d_rs
153 MODULE PROCEDURE pw_copy_c3d_r3d_rs
154 MODULE PROCEDURE pw_copy_c3d_c3d_rs
155 MODULE PROCEDURE pw_copy_r1d_r1d_gs
156 MODULE PROCEDURE pw_copy_r1d_c1d_gs
157 MODULE PROCEDURE pw_copy_r3d_r3d_gs
158 MODULE PROCEDURE pw_copy_r3d_c3d_gs
159 MODULE PROCEDURE pw_copy_c1d_r1d_gs
160 MODULE PROCEDURE pw_copy_c1d_c1d_gs
161 MODULE PROCEDURE pw_copy_c3d_r3d_gs
162 MODULE PROCEDURE pw_copy_c3d_c3d_gs
166 MODULE PROCEDURE pw_axpy_r1d_r1d_rs
167 MODULE PROCEDURE pw_axpy_r1d_c1d_rs
168 MODULE PROCEDURE pw_axpy_r3d_r3d_rs
169 MODULE PROCEDURE pw_axpy_r3d_c3d_rs
170 MODULE PROCEDURE pw_axpy_c1d_r1d_rs
171 MODULE PROCEDURE pw_axpy_c1d_c1d_rs
172 MODULE PROCEDURE pw_axpy_c3d_r3d_rs
173 MODULE PROCEDURE pw_axpy_c3d_c3d_rs
174 MODULE PROCEDURE pw_axpy_r1d_r1d_gs
175 MODULE PROCEDURE pw_axpy_r1d_c1d_gs
176 MODULE PROCEDURE pw_axpy_r3d_r3d_gs
177 MODULE PROCEDURE pw_axpy_r3d_c3d_gs
178 MODULE PROCEDURE pw_axpy_c1d_r1d_gs
179 MODULE PROCEDURE pw_axpy_c1d_c1d_gs
180 MODULE PROCEDURE pw_axpy_c3d_r3d_gs
181 MODULE PROCEDURE pw_axpy_c3d_c3d_gs
185 MODULE PROCEDURE pw_multiply_r1d_r1d_rs
186 MODULE PROCEDURE pw_multiply_r1d_c1d_rs
187 MODULE PROCEDURE pw_multiply_r3d_r3d_rs
188 MODULE PROCEDURE pw_multiply_r3d_c3d_rs
189 MODULE PROCEDURE pw_multiply_c1d_r1d_rs
190 MODULE PROCEDURE pw_multiply_c1d_c1d_rs
191 MODULE PROCEDURE pw_multiply_c3d_r3d_rs
192 MODULE PROCEDURE pw_multiply_c3d_c3d_rs
193 MODULE PROCEDURE pw_multiply_r1d_r1d_gs
194 MODULE PROCEDURE pw_multiply_r1d_c1d_gs
195 MODULE PROCEDURE pw_multiply_r3d_r3d_gs
196 MODULE PROCEDURE pw_multiply_r3d_c3d_gs
197 MODULE PROCEDURE pw_multiply_c1d_r1d_gs
198 MODULE PROCEDURE pw_multiply_c1d_c1d_gs
199 MODULE PROCEDURE pw_multiply_c3d_r3d_gs
200 MODULE PROCEDURE pw_multiply_c3d_c3d_gs
204 MODULE PROCEDURE pw_multiply_with_r1d_r1d_rs
205 MODULE PROCEDURE pw_multiply_with_r1d_c1d_rs
206 MODULE PROCEDURE pw_multiply_with_r3d_r3d_rs
207 MODULE PROCEDURE pw_multiply_with_r3d_c3d_rs
208 MODULE PROCEDURE pw_multiply_with_c1d_r1d_rs
209 MODULE PROCEDURE pw_multiply_with_c1d_c1d_rs
210 MODULE PROCEDURE pw_multiply_with_c3d_r3d_rs
211 MODULE PROCEDURE pw_multiply_with_c3d_c3d_rs
212 MODULE PROCEDURE pw_multiply_with_r1d_r1d_gs
213 MODULE PROCEDURE pw_multiply_with_r1d_c1d_gs
214 MODULE PROCEDURE pw_multiply_with_r3d_r3d_gs
215 MODULE PROCEDURE pw_multiply_with_r3d_c3d_gs
216 MODULE PROCEDURE pw_multiply_with_c1d_r1d_gs
217 MODULE PROCEDURE pw_multiply_with_c1d_c1d_gs
218 MODULE PROCEDURE pw_multiply_with_c3d_r3d_gs
219 MODULE PROCEDURE pw_multiply_with_c3d_c3d_gs
223 MODULE PROCEDURE pw_integral_ab_r1d_r1d_rs
224 MODULE PROCEDURE pw_integral_ab_r1d_c1d_rs
225 MODULE PROCEDURE pw_integral_ab_r3d_r3d_rs
226 MODULE PROCEDURE pw_integral_ab_r3d_c3d_rs
227 MODULE PROCEDURE pw_integral_ab_c1d_r1d_rs
228 MODULE PROCEDURE pw_integral_ab_c1d_c1d_rs
229 MODULE PROCEDURE pw_integral_ab_c3d_r3d_rs
230 MODULE PROCEDURE pw_integral_ab_c3d_c3d_rs
231 MODULE PROCEDURE pw_integral_ab_r1d_r1d_gs
232 MODULE PROCEDURE pw_integral_ab_r1d_c1d_gs
233 MODULE PROCEDURE pw_integral_ab_r3d_r3d_gs
234 MODULE PROCEDURE pw_integral_ab_r3d_c3d_gs
235 MODULE PROCEDURE pw_integral_ab_c1d_r1d_gs
236 MODULE PROCEDURE pw_integral_ab_c1d_c1d_gs
237 MODULE PROCEDURE pw_integral_ab_c3d_r3d_gs
238 MODULE PROCEDURE pw_integral_ab_c3d_c3d_gs
242 MODULE PROCEDURE pw_integral_a2b_r1d_r1d
243 MODULE PROCEDURE pw_integral_a2b_r1d_c1d
244 MODULE PROCEDURE pw_integral_a2b_c1d_r1d
245 MODULE PROCEDURE pw_integral_a2b_c1d_c1d
249 MODULE PROCEDURE pw_gather_p_r1d
250 MODULE PROCEDURE pw_gather_p_c1d
251 MODULE PROCEDURE pw_gather_s_r1d_r3d
252 MODULE PROCEDURE pw_gather_s_r1d_c3d
253 MODULE PROCEDURE pw_gather_s_c1d_r3d
254 MODULE PROCEDURE pw_gather_s_c1d_c3d
258 MODULE PROCEDURE pw_scatter_p_r1d
259 MODULE PROCEDURE pw_scatter_p_c1d
260 MODULE PROCEDURE pw_scatter_s_r1d_r3d
261 MODULE PROCEDURE pw_scatter_s_r1d_c3d
262 MODULE PROCEDURE pw_scatter_s_c1d_r3d
263 MODULE PROCEDURE pw_scatter_s_c1d_c3d
267 MODULE PROCEDURE pw_copy_to_array_r1d_r1d_rs
268 MODULE PROCEDURE pw_copy_to_array_r1d_c1d_rs
269 MODULE PROCEDURE pw_copy_to_array_r3d_r3d_rs
270 MODULE PROCEDURE pw_copy_to_array_r3d_c3d_rs
271 MODULE PROCEDURE pw_copy_to_array_c1d_r1d_rs
272 MODULE PROCEDURE pw_copy_to_array_c1d_c1d_rs
273 MODULE PROCEDURE pw_copy_to_array_c3d_r3d_rs
274 MODULE PROCEDURE pw_copy_to_array_c3d_c3d_rs
275 MODULE PROCEDURE pw_copy_to_array_r1d_r1d_gs
276 MODULE PROCEDURE pw_copy_to_array_r1d_c1d_gs
277 MODULE PROCEDURE pw_copy_to_array_r3d_r3d_gs
278 MODULE PROCEDURE pw_copy_to_array_r3d_c3d_gs
279 MODULE PROCEDURE pw_copy_to_array_c1d_r1d_gs
280 MODULE PROCEDURE pw_copy_to_array_c1d_c1d_gs
281 MODULE PROCEDURE pw_copy_to_array_c3d_r3d_gs
282 MODULE PROCEDURE pw_copy_to_array_c3d_c3d_gs
286 MODULE PROCEDURE pw_copy_from_array_r1d_r1d_rs
287 MODULE PROCEDURE pw_copy_from_array_r1d_c1d_rs
288 MODULE PROCEDURE pw_copy_from_array_r3d_r3d_rs
289 MODULE PROCEDURE pw_copy_from_array_r3d_c3d_rs
290 MODULE PROCEDURE pw_copy_from_array_c1d_r1d_rs
291 MODULE PROCEDURE pw_copy_from_array_c1d_c1d_rs
292 MODULE PROCEDURE pw_copy_from_array_c3d_r3d_rs
293 MODULE PROCEDURE pw_copy_from_array_c3d_c3d_rs
294 MODULE PROCEDURE pw_copy_from_array_r1d_r1d_gs
295 MODULE PROCEDURE pw_copy_from_array_r1d_c1d_gs
296 MODULE PROCEDURE pw_copy_from_array_r3d_r3d_gs
297 MODULE PROCEDURE pw_copy_from_array_r3d_c3d_gs
298 MODULE PROCEDURE pw_copy_from_array_c1d_r1d_gs
299 MODULE PROCEDURE pw_copy_from_array_c1d_c1d_gs
300 MODULE PROCEDURE pw_copy_from_array_c3d_r3d_gs
301 MODULE PROCEDURE pw_copy_from_array_c3d_c3d_gs
305 MODULE PROCEDURE pw_copy_r1d_r1d_rs
306 MODULE PROCEDURE pw_copy_r1d_r1d_gs
307 MODULE PROCEDURE pw_gather_s_r1d_r3d_2
308 MODULE PROCEDURE pw_scatter_s_r1d_r3d_2
309 MODULE PROCEDURE pw_copy_r1d_c1d_rs
310 MODULE PROCEDURE pw_copy_r1d_c1d_gs
311 MODULE PROCEDURE pw_gather_s_r1d_c3d_2
312 MODULE PROCEDURE pw_scatter_s_r1d_c3d_2
313 MODULE PROCEDURE pw_copy_r3d_r3d_rs
314 MODULE PROCEDURE pw_copy_r3d_r3d_gs
315 MODULE PROCEDURE fft_wrap_pw1pw2_r3d_c1d_rs_gs
316 MODULE PROCEDURE pw_copy_r3d_c3d_rs
317 MODULE PROCEDURE pw_copy_r3d_c3d_gs
318 MODULE PROCEDURE fft_wrap_pw1pw2_r3d_c3d_rs_gs
319 MODULE PROCEDURE pw_copy_c1d_r1d_rs
320 MODULE PROCEDURE pw_copy_c1d_r1d_gs
321 MODULE PROCEDURE pw_gather_s_c1d_r3d_2
322 MODULE PROCEDURE pw_scatter_s_c1d_r3d_2
323 MODULE PROCEDURE fft_wrap_pw1pw2_c1d_r3d_gs_rs
324 MODULE PROCEDURE pw_copy_c1d_c1d_rs
325 MODULE PROCEDURE pw_copy_c1d_c1d_gs
326 MODULE PROCEDURE pw_gather_s_c1d_c3d_2
327 MODULE PROCEDURE pw_scatter_s_c1d_c3d_2
328 MODULE PROCEDURE fft_wrap_pw1pw2_c1d_c3d_gs_rs
329 MODULE PROCEDURE pw_copy_c3d_r3d_rs
330 MODULE PROCEDURE pw_copy_c3d_r3d_gs
331 MODULE PROCEDURE fft_wrap_pw1pw2_c3d_r3d_gs_rs
332 MODULE PROCEDURE fft_wrap_pw1pw2_c3d_c1d_rs_gs
333 MODULE PROCEDURE pw_copy_c3d_c3d_rs
334 MODULE PROCEDURE pw_copy_c3d_c3d_gs
335 MODULE PROCEDURE fft_wrap_pw1pw2_c3d_c3d_rs_gs
336 MODULE PROCEDURE fft_wrap_pw1pw2_c3d_c3d_gs_rs
347 SUBROUTINE pw_zero_r1d_rs (pw)
351 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_zero'
355 CALL timeset(routinen, handle)
357#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
365 CALL timestop(handle)
367 END SUBROUTINE pw_zero_r1d_rs
376 SUBROUTINE pw_scale_r1d_rs (pw, a)
379 REAL(KIND=
dp),
INTENT(IN) :: a
381 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_scale'
385 CALL timeset(routinen, handle)
388 pw%array = a*pw%array
391 CALL timestop(handle)
393 END SUBROUTINE pw_scale_r1d_rs
405 SUBROUTINE pw_write_r1d_rs (pw, unit_nr)
408 INTEGER,
INTENT(in) :: unit_nr
412 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
414 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=r1d"
415 IF (
ASSOCIATED(pw%array))
THEN
416 WRITE (unit=unit_nr, fmt=
"(' array=<r(',i8,':',i8,')>')") &
417 lbound(pw%array, 1), ubound(pw%array, 1)
419 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
422 END SUBROUTINE pw_write_r1d_rs
431 FUNCTION pw_integrate_function_r1d_rs (fun, isign, oprt)
RESULT(total_fun)
434 INTEGER,
INTENT(IN),
OPTIONAL :: isign
435 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
436 REAL(kind=
dp) :: total_fun
442 IF (
PRESENT(oprt))
THEN
447 cpabort(
"Unknown operator")
455 total_fun = fun%pw_grid%dvol*
accurate_sum(abs(fun%array))
461 CALL fun%pw_grid%para%group%sum(total_fun)
464 IF (
PRESENT(isign))
THEN
465 total_fun = total_fun*sign(1._dp, real(isign,
dp))
468 END FUNCTION pw_integrate_function_r1d_rs
475 SUBROUTINE pw_set_value_r1d_rs (pw, value)
477 REAL(KIND=
dp),
INTENT(IN) ::
value
479 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_set_value'
483 CALL timeset(routinen, handle)
489 CALL timestop(handle)
491 END SUBROUTINE pw_set_value_r1d_rs
499 SUBROUTINE pw_zero_r1d_gs (pw)
503 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_zero'
507 CALL timeset(routinen, handle)
509#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
517 CALL timestop(handle)
519 END SUBROUTINE pw_zero_r1d_gs
528 SUBROUTINE pw_scale_r1d_gs (pw, a)
531 REAL(KIND=
dp),
INTENT(IN) :: a
533 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_scale'
537 CALL timeset(routinen, handle)
540 pw%array = a*pw%array
543 CALL timestop(handle)
545 END SUBROUTINE pw_scale_r1d_gs
557 SUBROUTINE pw_write_r1d_gs (pw, unit_nr)
560 INTEGER,
INTENT(in) :: unit_nr
564 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
566 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=r1d"
567 IF (
ASSOCIATED(pw%array))
THEN
568 WRITE (unit=unit_nr, fmt=
"(' array=<r(',i8,':',i8,')>')") &
569 lbound(pw%array, 1), ubound(pw%array, 1)
571 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
574 END SUBROUTINE pw_write_r1d_gs
583 FUNCTION pw_integrate_function_r1d_gs (fun, isign, oprt)
RESULT(total_fun)
586 INTEGER,
INTENT(IN),
OPTIONAL :: isign
587 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
588 REAL(kind=
dp) :: total_fun
594 IF (
PRESENT(oprt))
THEN
599 cpabort(
"Unknown operator")
606 cpabort(
"Operator ABS not implemented")
607 IF (fun%pw_grid%have_g0) total_fun = fun%pw_grid%vol* fun%array(1)
610 CALL fun%pw_grid%para%group%sum(total_fun)
613 IF (
PRESENT(isign))
THEN
614 total_fun = total_fun*sign(1._dp, real(isign,
dp))
617 END FUNCTION pw_integrate_function_r1d_gs
624 SUBROUTINE pw_set_value_r1d_gs (pw, value)
626 REAL(KIND=
dp),
INTENT(IN) ::
value
628 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_set_value'
632 CALL timeset(routinen, handle)
638 CALL timestop(handle)
640 END SUBROUTINE pw_set_value_r1d_gs
648 SUBROUTINE pw_gather_p_r1d (pw, c, scale)
651 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
INTENT(IN),
TARGET :: c
652 REAL(KIND=
dp),
INTENT(IN),
OPTIONAL :: scale
654 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_gather_p'
656 INTEGER :: gpt, handle
657 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:),
POINTER :: c_1d
659 CALL timeset(routinen, handle)
662 cpabort(
"This grid type is not distributed")
667 associate(ngpts =>
SIZE(pw%pw_grid%gsq), g_hatmap => pw%pw_grid%g_hatmap)
669 IF (
PRESENT(scale))
THEN
674 pw%array(gpt) = scale* real(c_1d(g_hatmap(gpt, 1) + 1), kind=
dp)
682 pw%array(gpt) = real(c_1d(g_hatmap(gpt, 1) + 1), kind=
dp)
689 CALL timestop(handle)
691 END SUBROUTINE pw_gather_p_r1d
699 SUBROUTINE pw_scatter_p_r1d (pw, c, scale)
700 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw
701 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
INTENT(INOUT),
TARGET :: c
702 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
704 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_scatter_p'
706 INTEGER :: gpt, handle
707 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:),
POINTER :: c_1d
709 CALL timeset(routinen, handle)
711 IF (pw%pw_grid%para%mode /= pw_mode_distributed)
THEN
712 cpabort(
"This grid type is not distributed")
717 associate(g_hatmap => pw%pw_grid%g_hatmap, ngpts =>
SIZE(pw%pw_grid%gsq))
719 IF (.NOT.
PRESENT(scale)) c = z_zero
721 IF (
PRESENT(scale))
THEN
726 c_1d(g_hatmap(gpt, 1) + 1) = cmplx(scale*pw%array(gpt), 0.0_dp, kind=dp)
734 c_1d(g_hatmap(gpt, 1) + 1) = cmplx(pw%array(gpt), 0.0_dp, kind=dp)
739 IF (pw%pw_grid%grid_span == halfspace)
THEN
741 IF (
PRESENT(scale))
THEN
746 c_1d(g_hatmap(gpt, 2) + 1) = scale*( cmplx(pw%array(gpt), 0.0_dp, kind=dp))
754 c_1d(g_hatmap(gpt, 2) + 1) = ( cmplx(pw%array(gpt), 0.0_dp, kind=dp))
762 CALL timestop(handle)
764 END SUBROUTINE pw_scatter_p_r1d
772 SUBROUTINE pw_zero_r3d_rs (pw)
774 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw
776 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_zero'
780 CALL timeset(routinen, handle)
782#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
790 CALL timestop(handle)
792 END SUBROUTINE pw_zero_r3d_rs
801 SUBROUTINE pw_scale_r3d_rs (pw, a)
803 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw
804 REAL(KIND=dp),
INTENT(IN) :: a
806 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_scale'
810 CALL timeset(routinen, handle)
813 pw%array = a*pw%array
816 CALL timestop(handle)
818 END SUBROUTINE pw_scale_r3d_rs
830 SUBROUTINE pw_write_r3d_rs (pw, unit_nr)
832 TYPE(pw_r3d_rs_type),
INTENT(in) :: pw
833 INTEGER,
INTENT(in) :: unit_nr
837 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
839 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=r3d"
840 IF (
ASSOCIATED(pw%array))
THEN
841 WRITE (unit=unit_nr, fmt=
"(' array=<r(',i8,':',i8,',',i8,':',i8,',',i8,':',i8,')>')") &
842 lbound(pw%array, 1), ubound(pw%array, 1), lbound(pw%array, 2), ubound(pw%array, 2), &
843 lbound(pw%array, 3), ubound(pw%array, 3)
845 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
848 END SUBROUTINE pw_write_r3d_rs
857 FUNCTION pw_integrate_function_r3d_rs (fun, isign, oprt)
RESULT(total_fun)
859 TYPE(pw_r3d_rs_type),
INTENT(IN) :: fun
860 INTEGER,
INTENT(IN),
OPTIONAL :: isign
861 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
862 REAL(kind=dp) :: total_fun
868 IF (
PRESENT(oprt))
THEN
873 cpabort(
"Unknown operator")
881 total_fun = fun%pw_grid%dvol*accurate_sum(abs(fun%array))
883 total_fun = fun%pw_grid%dvol*accurate_sum( fun%array)
886 IF (fun%pw_grid%para%mode /= pw_mode_local)
THEN
887 CALL fun%pw_grid%para%group%sum(total_fun)
890 IF (
PRESENT(isign))
THEN
891 total_fun = total_fun*sign(1._dp, real(isign, dp))
894 END FUNCTION pw_integrate_function_r3d_rs
901 SUBROUTINE pw_set_value_r3d_rs (pw, value)
902 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw
903 REAL(KIND=dp),
INTENT(IN) ::
value
905 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_set_value'
909 CALL timeset(routinen, handle)
915 CALL timestop(handle)
917 END SUBROUTINE pw_set_value_r3d_rs
925 SUBROUTINE pw_zero_r3d_gs (pw)
927 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw
929 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_zero'
933 CALL timeset(routinen, handle)
935#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
943 CALL timestop(handle)
945 END SUBROUTINE pw_zero_r3d_gs
954 SUBROUTINE pw_scale_r3d_gs (pw, a)
956 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw
957 REAL(KIND=dp),
INTENT(IN) :: a
959 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_scale'
963 CALL timeset(routinen, handle)
966 pw%array = a*pw%array
969 CALL timestop(handle)
971 END SUBROUTINE pw_scale_r3d_gs
983 SUBROUTINE pw_write_r3d_gs (pw, unit_nr)
985 TYPE(pw_r3d_gs_type),
INTENT(in) :: pw
986 INTEGER,
INTENT(in) :: unit_nr
990 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
992 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=r3d"
993 IF (
ASSOCIATED(pw%array))
THEN
994 WRITE (unit=unit_nr, fmt=
"(' array=<r(',i8,':',i8,',',i8,':',i8,',',i8,':',i8,')>')") &
995 lbound(pw%array, 1), ubound(pw%array, 1), lbound(pw%array, 2), ubound(pw%array, 2), &
996 lbound(pw%array, 3), ubound(pw%array, 3)
998 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
1001 END SUBROUTINE pw_write_r3d_gs
1010 FUNCTION pw_integrate_function_r3d_gs (fun, isign, oprt)
RESULT(total_fun)
1012 TYPE(pw_r3d_gs_type),
INTENT(IN) :: fun
1013 INTEGER,
INTENT(IN),
OPTIONAL :: isign
1014 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
1015 REAL(kind=dp) :: total_fun
1021 IF (
PRESENT(oprt))
THEN
1026 cpabort(
"Unknown operator")
1033 cpabort(
"Operator ABS not implemented")
1034 cpabort(
"Reciprocal space integration for 3D grids not implemented!")
1036 IF (fun%pw_grid%para%mode /= pw_mode_local)
THEN
1037 CALL fun%pw_grid%para%group%sum(total_fun)
1040 IF (
PRESENT(isign))
THEN
1041 total_fun = total_fun*sign(1._dp, real(isign, dp))
1044 END FUNCTION pw_integrate_function_r3d_gs
1051 SUBROUTINE pw_set_value_r3d_gs (pw, value)
1052 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw
1053 REAL(KIND=dp),
INTENT(IN) ::
value
1055 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_set_value'
1059 CALL timeset(routinen, handle)
1065 CALL timestop(handle)
1067 END SUBROUTINE pw_set_value_r3d_gs
1076 SUBROUTINE pw_zero_c1d_rs (pw)
1078 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw
1080 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_zero'
1084 CALL timeset(routinen, handle)
1086#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
1088 pw%array = cmplx(0.0_dp, 0.0_dp, kind=dp)
1091 pw%array = cmplx(0.0_dp, 0.0_dp, kind=dp)
1094 CALL timestop(handle)
1096 END SUBROUTINE pw_zero_c1d_rs
1105 SUBROUTINE pw_scale_c1d_rs (pw, a)
1107 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw
1108 REAL(KIND=dp),
INTENT(IN) :: a
1110 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_scale'
1114 CALL timeset(routinen, handle)
1117 pw%array = a*pw%array
1120 CALL timestop(handle)
1122 END SUBROUTINE pw_scale_c1d_rs
1134 SUBROUTINE pw_write_c1d_rs (pw, unit_nr)
1136 TYPE(pw_c1d_rs_type),
INTENT(in) :: pw
1137 INTEGER,
INTENT(in) :: unit_nr
1141 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
1143 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=c1d"
1144 IF (
ASSOCIATED(pw%array))
THEN
1145 WRITE (unit=unit_nr, fmt=
"(' array=<c(',i8,':',i8,')>')") &
1146 lbound(pw%array, 1), ubound(pw%array, 1)
1148 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
1151 END SUBROUTINE pw_write_c1d_rs
1160 FUNCTION pw_integrate_function_c1d_rs (fun, isign, oprt)
RESULT(total_fun)
1162 TYPE(pw_c1d_rs_type),
INTENT(IN) :: fun
1163 INTEGER,
INTENT(IN),
OPTIONAL :: isign
1164 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
1165 REAL(kind=dp) :: total_fun
1171 IF (
PRESENT(oprt))
THEN
1176 cpabort(
"Unknown operator")
1184 total_fun = fun%pw_grid%dvol*accurate_sum(abs(fun%array))
1186 total_fun = fun%pw_grid%dvol*accurate_sum( real(fun%array, kind=dp))
1189 IF (fun%pw_grid%para%mode /= pw_mode_local)
THEN
1190 CALL fun%pw_grid%para%group%sum(total_fun)
1193 IF (
PRESENT(isign))
THEN
1194 total_fun = total_fun*sign(1._dp, real(isign, dp))
1197 END FUNCTION pw_integrate_function_c1d_rs
1204 SUBROUTINE pw_set_value_c1d_rs (pw, value)
1205 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw
1206 REAL(KIND=dp),
INTENT(IN) ::
value
1208 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_set_value'
1212 CALL timeset(routinen, handle)
1215 pw%array = cmplx(
value, 0.0_dp, kind=dp)
1218 CALL timestop(handle)
1220 END SUBROUTINE pw_set_value_c1d_rs
1228 SUBROUTINE pw_zero_c1d_gs (pw)
1230 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
1232 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_zero'
1236 CALL timeset(routinen, handle)
1238#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
1240 pw%array = cmplx(0.0_dp, 0.0_dp, kind=dp)
1243 pw%array = cmplx(0.0_dp, 0.0_dp, kind=dp)
1246 CALL timestop(handle)
1248 END SUBROUTINE pw_zero_c1d_gs
1257 SUBROUTINE pw_scale_c1d_gs (pw, a)
1259 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
1260 REAL(KIND=dp),
INTENT(IN) :: a
1262 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_scale'
1266 CALL timeset(routinen, handle)
1269 pw%array = a*pw%array
1272 CALL timestop(handle)
1274 END SUBROUTINE pw_scale_c1d_gs
1286 SUBROUTINE pw_write_c1d_gs (pw, unit_nr)
1288 TYPE(pw_c1d_gs_type),
INTENT(in) :: pw
1289 INTEGER,
INTENT(in) :: unit_nr
1293 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
1295 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=c1d"
1296 IF (
ASSOCIATED(pw%array))
THEN
1297 WRITE (unit=unit_nr, fmt=
"(' array=<c(',i8,':',i8,')>')") &
1298 lbound(pw%array, 1), ubound(pw%array, 1)
1300 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
1303 END SUBROUTINE pw_write_c1d_gs
1312 FUNCTION pw_integrate_function_c1d_gs (fun, isign, oprt)
RESULT(total_fun)
1314 TYPE(pw_c1d_gs_type),
INTENT(IN) :: fun
1315 INTEGER,
INTENT(IN),
OPTIONAL :: isign
1316 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
1317 REAL(kind=dp) :: total_fun
1323 IF (
PRESENT(oprt))
THEN
1328 cpabort(
"Unknown operator")
1335 cpabort(
"Operator ABS not implemented")
1336 IF (fun%pw_grid%have_g0) total_fun = fun%pw_grid%vol* real(fun%array(1), kind=dp)
1338 IF (fun%pw_grid%para%mode /= pw_mode_local)
THEN
1339 CALL fun%pw_grid%para%group%sum(total_fun)
1342 IF (
PRESENT(isign))
THEN
1343 total_fun = total_fun*sign(1._dp, real(isign, dp))
1346 END FUNCTION pw_integrate_function_c1d_gs
1353 SUBROUTINE pw_set_value_c1d_gs (pw, value)
1354 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
1355 REAL(KIND=dp),
INTENT(IN) ::
value
1357 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_set_value'
1361 CALL timeset(routinen, handle)
1364 pw%array = cmplx(
value, 0.0_dp, kind=dp)
1367 CALL timestop(handle)
1369 END SUBROUTINE pw_set_value_c1d_gs
1377 SUBROUTINE pw_gather_p_c1d (pw, c, scale)
1379 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
1380 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
INTENT(IN),
TARGET :: c
1381 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
1383 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_gather_p'
1385 INTEGER :: gpt, handle
1386 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:),
POINTER :: c_1d
1388 CALL timeset(routinen, handle)
1390 IF (pw%pw_grid%para%mode /= pw_mode_distributed)
THEN
1391 cpabort(
"This grid type is not distributed")
1395 c_1d(1:
SIZE(c)) => c
1396 associate(ngpts =>
SIZE(pw%pw_grid%gsq), g_hatmap => pw%pw_grid%g_hatmap)
1398 IF (
PRESENT(scale))
THEN
1403 pw%array(gpt) = scale* c_1d(g_hatmap(gpt, 1) + 1)
1411 pw%array(gpt) = c_1d(g_hatmap(gpt, 1) + 1)
1418 CALL timestop(handle)
1420 END SUBROUTINE pw_gather_p_c1d
1428 SUBROUTINE pw_scatter_p_c1d (pw, c, scale)
1429 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
1430 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
INTENT(INOUT),
TARGET :: c
1431 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
1433 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_scatter_p'
1435 INTEGER :: gpt, handle
1436 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:),
POINTER :: c_1d
1438 CALL timeset(routinen, handle)
1440 IF (pw%pw_grid%para%mode /= pw_mode_distributed)
THEN
1441 cpabort(
"This grid type is not distributed")
1445 c_1d(1:
SIZE(c)) => c
1446 associate(g_hatmap => pw%pw_grid%g_hatmap, ngpts =>
SIZE(pw%pw_grid%gsq))
1448 IF (.NOT.
PRESENT(scale)) c = z_zero
1450 IF (
PRESENT(scale))
THEN
1455 c_1d(g_hatmap(gpt, 1) + 1) = scale*pw%array(gpt)
1463 c_1d(g_hatmap(gpt, 1) + 1) = pw%array(gpt)
1468 IF (pw%pw_grid%grid_span == halfspace)
THEN
1470 IF (
PRESENT(scale))
THEN
1475 c_1d(g_hatmap(gpt, 2) + 1) = scale*conjg( pw%array(gpt))
1483 c_1d(g_hatmap(gpt, 2) + 1) = conjg( pw%array(gpt))
1491 CALL timestop(handle)
1493 END SUBROUTINE pw_scatter_p_c1d
1501 SUBROUTINE pw_zero_c3d_rs (pw)
1503 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw
1505 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_zero'
1509 CALL timeset(routinen, handle)
1511#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
1513 pw%array = cmplx(0.0_dp, 0.0_dp, kind=dp)
1516 pw%array = cmplx(0.0_dp, 0.0_dp, kind=dp)
1519 CALL timestop(handle)
1521 END SUBROUTINE pw_zero_c3d_rs
1530 SUBROUTINE pw_scale_c3d_rs (pw, a)
1532 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw
1533 REAL(KIND=dp),
INTENT(IN) :: a
1535 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_scale'
1539 CALL timeset(routinen, handle)
1542 pw%array = a*pw%array
1545 CALL timestop(handle)
1547 END SUBROUTINE pw_scale_c3d_rs
1559 SUBROUTINE pw_write_c3d_rs (pw, unit_nr)
1561 TYPE(pw_c3d_rs_type),
INTENT(in) :: pw
1562 INTEGER,
INTENT(in) :: unit_nr
1566 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
1568 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=c3d"
1569 IF (
ASSOCIATED(pw%array))
THEN
1570 WRITE (unit=unit_nr, fmt=
"(' array=<c(',i8,':',i8,',',i8,':',i8,',',i8,':',i8,')>')") &
1571 lbound(pw%array, 1), ubound(pw%array, 1), lbound(pw%array, 2), ubound(pw%array, 2), &
1572 lbound(pw%array, 3), ubound(pw%array, 3)
1574 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
1577 END SUBROUTINE pw_write_c3d_rs
1586 FUNCTION pw_integrate_function_c3d_rs (fun, isign, oprt)
RESULT(total_fun)
1588 TYPE(pw_c3d_rs_type),
INTENT(IN) :: fun
1589 INTEGER,
INTENT(IN),
OPTIONAL :: isign
1590 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
1591 REAL(kind=dp) :: total_fun
1597 IF (
PRESENT(oprt))
THEN
1602 cpabort(
"Unknown operator")
1610 total_fun = fun%pw_grid%dvol*accurate_sum(abs(fun%array))
1612 total_fun = fun%pw_grid%dvol*accurate_sum( real(fun%array, kind=dp))
1615 IF (fun%pw_grid%para%mode /= pw_mode_local)
THEN
1616 CALL fun%pw_grid%para%group%sum(total_fun)
1619 IF (
PRESENT(isign))
THEN
1620 total_fun = total_fun*sign(1._dp, real(isign, dp))
1623 END FUNCTION pw_integrate_function_c3d_rs
1630 SUBROUTINE pw_set_value_c3d_rs (pw, value)
1631 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw
1632 REAL(KIND=dp),
INTENT(IN) ::
value
1634 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_set_value'
1638 CALL timeset(routinen, handle)
1641 pw%array = cmplx(
value, 0.0_dp, kind=dp)
1644 CALL timestop(handle)
1646 END SUBROUTINE pw_set_value_c3d_rs
1654 SUBROUTINE pw_zero_c3d_gs (pw)
1656 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw
1658 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_zero'
1662 CALL timeset(routinen, handle)
1664#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
1666 pw%array = cmplx(0.0_dp, 0.0_dp, kind=dp)
1669 pw%array = cmplx(0.0_dp, 0.0_dp, kind=dp)
1672 CALL timestop(handle)
1674 END SUBROUTINE pw_zero_c3d_gs
1683 SUBROUTINE pw_scale_c3d_gs (pw, a)
1685 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw
1686 REAL(KIND=dp),
INTENT(IN) :: a
1688 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_scale'
1692 CALL timeset(routinen, handle)
1695 pw%array = a*pw%array
1698 CALL timestop(handle)
1700 END SUBROUTINE pw_scale_c3d_gs
1712 SUBROUTINE pw_write_c3d_gs (pw, unit_nr)
1714 TYPE(pw_c3d_gs_type),
INTENT(in) :: pw
1715 INTEGER,
INTENT(in) :: unit_nr
1719 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
1721 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=c3d"
1722 IF (
ASSOCIATED(pw%array))
THEN
1723 WRITE (unit=unit_nr, fmt=
"(' array=<c(',i8,':',i8,',',i8,':',i8,',',i8,':',i8,')>')") &
1724 lbound(pw%array, 1), ubound(pw%array, 1), lbound(pw%array, 2), ubound(pw%array, 2), &
1725 lbound(pw%array, 3), ubound(pw%array, 3)
1727 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
1730 END SUBROUTINE pw_write_c3d_gs
1739 FUNCTION pw_integrate_function_c3d_gs (fun, isign, oprt)
RESULT(total_fun)
1741 TYPE(pw_c3d_gs_type),
INTENT(IN) :: fun
1742 INTEGER,
INTENT(IN),
OPTIONAL :: isign
1743 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
1744 REAL(kind=dp) :: total_fun
1750 IF (
PRESENT(oprt))
THEN
1755 cpabort(
"Unknown operator")
1762 cpabort(
"Operator ABS not implemented")
1763 cpabort(
"Reciprocal space integration for 3D grids not implemented!")
1765 IF (fun%pw_grid%para%mode /= pw_mode_local)
THEN
1766 CALL fun%pw_grid%para%group%sum(total_fun)
1769 IF (
PRESENT(isign))
THEN
1770 total_fun = total_fun*sign(1._dp, real(isign, dp))
1773 END FUNCTION pw_integrate_function_c3d_gs
1780 SUBROUTINE pw_set_value_c3d_gs (pw, value)
1781 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw
1782 REAL(KIND=dp),
INTENT(IN) ::
value
1784 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_set_value'
1788 CALL timeset(routinen, handle)
1791 pw%array = cmplx(
value, 0.0_dp, kind=dp)
1794 CALL timestop(handle)
1796 END SUBROUTINE pw_set_value_c3d_gs
1812 SUBROUTINE pw_copy_r1d_r1d_rs (pw1, pw2)
1814 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
1815 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw2
1817 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
1820 INTEGER :: i, j, ng, ng1, ng2, ns
1822 CALL timeset(routinen, handle)
1824 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
1825 cpabort(
"Both grids must be either spherical or non-spherical!")
1826 IF (pw1%pw_grid%spherical) &
1827 cpabort(
"Spherical grids only exist in reciprocal space!")
1829 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
1830 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
1831 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
1832 ng1 =
SIZE(pw1%array)
1833 ng2 =
SIZE(pw2%array)
1836 pw2%array(1:ng) = pw1%array(1:ng)
1840 pw2%array(ng + 1:ng2) = 0.0_dp
1844 cpabort(
"Copies between spherical grids require compatible grids!")
1847 ng1 =
SIZE(pw1%array)
1848 ng2 =
SIZE(pw2%array)
1849 ns = 2*max(ng1, ng2)
1851 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
1852 IF (ng1 >= ng2)
THEN
1857 j = pw2%pw_grid%gidx(i)
1858 pw2%array(i) = pw1%array(j)
1867 j = pw2%pw_grid%gidx(i)
1868 pw2%array(j) = pw1%array(i)
1872 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
1873 IF (ng1 >= ng2)
THEN
1876 j = pw1%pw_grid%gidx(i)
1877 pw2%array(i) = pw1%array(j)
1884 j = pw1%pw_grid%gidx(i)
1885 pw2%array(j) = pw1%array(i)
1890 cpabort(
"Copy not implemented!")
1897 pw2%array = pw1%array
1901 CALL timestop(handle)
1903 END SUBROUTINE pw_copy_r1d_r1d_rs
1910 SUBROUTINE pw_copy_to_array_r1d_r1d_rs (pw, array)
1911 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw
1912 REAL(KIND=dp),
DIMENSION(:),
INTENT(INOUT) :: array
1914 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
1918 CALL timeset(routinen, handle)
1921 array(:) = pw%array(:)
1924 CALL timestop(handle)
1925 END SUBROUTINE pw_copy_to_array_r1d_r1d_rs
1932 SUBROUTINE pw_copy_from_array_r1d_r1d_rs (pw, array)
1933 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw
1934 REAL(KIND=dp),
DIMENSION(:),
INTENT(IN) :: array
1936 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
1940 CALL timeset(routinen, handle)
1946 CALL timestop(handle)
1947 END SUBROUTINE pw_copy_from_array_r1d_r1d_rs
1965 SUBROUTINE pw_axpy_r1d_r1d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
1967 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
1968 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw2
1969 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
1970 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
1972 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
1975 LOGICAL :: my_allow_noncompatible_grids
1976 REAL(KIND=dp) :: my_alpha, my_beta
1977 INTEGER :: i, j, ng, ng1, ng2
1979 CALL timeset(routinen, handle)
1982 IF (
PRESENT(alpha)) my_alpha = alpha
1985 IF (
PRESENT(beta)) my_beta = beta
1987 my_allow_noncompatible_grids = .false.
1988 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
1990 IF (my_beta /= 1.0_dp)
THEN
1991 IF (my_beta == 0.0_dp)
THEN
1995 pw2%array = pw2%array*my_beta
2000 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
2002 IF (my_alpha == 1.0_dp)
THEN
2004 pw2%array = pw2%array + pw1%array
2006 ELSE IF (my_alpha /= 0.0_dp)
THEN
2008 pw2%array = pw2%array + my_alpha* pw1%array
2012 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
2014 ng1 =
SIZE(pw1%array)
2015 ng2 =
SIZE(pw2%array)
2018 IF (pw1%pw_grid%spherical)
THEN
2019 IF (my_alpha == 1.0_dp)
THEN
2022 pw2%array(i) = pw2%array(i) + pw1%array(i)
2025 ELSE IF (my_alpha /= 0.0_dp)
THEN
2028 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(i)
2032 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
2033 IF (ng1 >= ng2)
THEN
2034 IF (my_alpha == 1.0_dp)
THEN
2037 j = pw2%pw_grid%gidx(i)
2038 pw2%array(i) = pw2%array(i) + pw1%array(j)
2041 ELSE IF (my_alpha /= 0.0_dp)
THEN
2044 j = pw2%pw_grid%gidx(i)
2045 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
2050 IF (my_alpha == 1.0_dp)
THEN
2053 j = pw2%pw_grid%gidx(i)
2054 pw2%array(j) = pw2%array(j) + pw1%array(i)
2057 ELSE IF (my_alpha /= 0.0_dp)
THEN
2060 j = pw2%pw_grid%gidx(i)
2061 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
2066 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
2067 IF (ng1 >= ng2)
THEN
2068 IF (my_alpha == 1.0_dp)
THEN
2071 j = pw1%pw_grid%gidx(i)
2072 pw2%array(i) = pw2%array(i) + pw1%array(j)
2075 ELSE IF (my_alpha /= 0.0_dp)
THEN
2078 j = pw1%pw_grid%gidx(i)
2079 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
2084 IF (my_alpha == 1.0_dp)
THEN
2087 j = pw1%pw_grid%gidx(i)
2088 pw2%array(j) = pw2%array(j) + pw1%array(i)
2091 ELSE IF (my_alpha /= 0.0_dp)
THEN
2094 j = pw1%pw_grid%gidx(i)
2095 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
2101 cpabort(
"Grids not compatible")
2106 cpabort(
"Grids not compatible")
2110 CALL timestop(handle)
2112 END SUBROUTINE pw_axpy_r1d_r1d_rs
2123 SUBROUTINE pw_multiply_r1d_r1d_rs (pw_out, pw1, pw2, alpha)
2125 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw_out
2126 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
2127 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw2
2128 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
2130 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
2133 REAL(KIND=dp) :: my_alpha
2135 CALL timeset(routinen, handle)
2138 IF (
PRESENT(alpha)) my_alpha = alpha
2140 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
2141 cpabort(
"pw_multiply not implemented for non-identical grids!")
2143#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
2144 IF (my_alpha == 1.0_dp)
THEN
2146 pw_out%array = pw_out%array + pw1%array* pw2%array
2150 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
2154 IF (my_alpha == 1.0_dp)
THEN
2155 pw_out%array = pw_out%array + pw1%array* pw2%array
2157 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
2161 CALL timestop(handle)
2163 END SUBROUTINE pw_multiply_r1d_r1d_rs
2170 SUBROUTINE pw_multiply_with_r1d_r1d_rs (pw1, pw2)
2171 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw1
2172 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw2
2174 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
2178 CALL timeset(routinen, handle)
2180 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
2181 cpabort(
"Incompatible grids!")
2184 pw1%array = pw1%array* pw2%array
2187 CALL timestop(handle)
2189 END SUBROUTINE pw_multiply_with_r1d_r1d_rs
2204 FUNCTION pw_integral_ab_r1d_r1d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
2206 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
2207 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw2
2208 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
2209 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
2210 REAL(kind=dp) :: integral_value
2212 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_r1d_r1d_rs'
2214 INTEGER :: handle, loc_sumtype
2215 LOGICAL :: my_just_sum, my_local_only
2217 CALL timeset(routinen, handle)
2220 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
2222 my_local_only = .false.
2223 IF (
PRESENT(local_only)) my_local_only = local_only
2225 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
2226 cpabort(
"Grids incompatible")
2229 my_just_sum = .false.
2230 IF (
PRESENT(just_sum)) my_just_sum = just_sum
2237 integral_value = dot_product(pw1%array, pw2%array)
2242 integral_value = accurate_dot_product(pw1%array, pw2%array)
2248 IF (pw1%pw_grid%grid_span == halfspace)
THEN
2249 integral_value = 2.0_dp*integral_value
2250 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
2251 pw1%array(1)*pw2%array(1)
2254 IF (.NOT. my_just_sum)
THEN
2255 integral_value = integral_value*pw1%pw_grid%dvol
2258 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
2259 CALL pw1%pw_grid%para%group%sum(integral_value)
2261 CALL timestop(handle)
2263 END FUNCTION pw_integral_ab_r1d_r1d_rs
2277 SUBROUTINE pw_copy_r1d_r1d_gs (pw1, pw2)
2279 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
2280 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw2
2282 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
2285 INTEGER :: i, j, ng, ng1, ng2, ns
2287 CALL timeset(routinen, handle)
2289 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
2290 cpabort(
"Both grids must be either spherical or non-spherical!")
2292 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
2293 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
2294 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
2295 ng1 =
SIZE(pw1%array)
2296 ng2 =
SIZE(pw2%array)
2299 pw2%array(1:ng) = pw1%array(1:ng)
2303 pw2%array(ng + 1:ng2) = 0.0_dp
2307 cpabort(
"Copies between spherical grids require compatible grids!")
2310 ng1 =
SIZE(pw1%array)
2311 ng2 =
SIZE(pw2%array)
2312 ns = 2*max(ng1, ng2)
2314 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
2315 IF (ng1 >= ng2)
THEN
2320 j = pw2%pw_grid%gidx(i)
2321 pw2%array(i) = pw1%array(j)
2330 j = pw2%pw_grid%gidx(i)
2331 pw2%array(j) = pw1%array(i)
2335 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
2336 IF (ng1 >= ng2)
THEN
2339 j = pw1%pw_grid%gidx(i)
2340 pw2%array(i) = pw1%array(j)
2347 j = pw1%pw_grid%gidx(i)
2348 pw2%array(j) = pw1%array(i)
2353 cpabort(
"Copy not implemented!")
2360 pw2%array = pw1%array
2364 CALL timestop(handle)
2366 END SUBROUTINE pw_copy_r1d_r1d_gs
2373 SUBROUTINE pw_copy_to_array_r1d_r1d_gs (pw, array)
2374 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw
2375 REAL(KIND=dp),
DIMENSION(:),
INTENT(INOUT) :: array
2377 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
2381 CALL timeset(routinen, handle)
2384 array(:) = pw%array(:)
2387 CALL timestop(handle)
2388 END SUBROUTINE pw_copy_to_array_r1d_r1d_gs
2395 SUBROUTINE pw_copy_from_array_r1d_r1d_gs (pw, array)
2396 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw
2397 REAL(KIND=dp),
DIMENSION(:),
INTENT(IN) :: array
2399 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
2403 CALL timeset(routinen, handle)
2409 CALL timestop(handle)
2410 END SUBROUTINE pw_copy_from_array_r1d_r1d_gs
2428 SUBROUTINE pw_axpy_r1d_r1d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
2430 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
2431 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw2
2432 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
2433 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
2435 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
2438 LOGICAL :: my_allow_noncompatible_grids
2439 REAL(KIND=dp) :: my_alpha, my_beta
2440 INTEGER :: i, j, ng, ng1, ng2
2442 CALL timeset(routinen, handle)
2445 IF (
PRESENT(alpha)) my_alpha = alpha
2448 IF (
PRESENT(beta)) my_beta = beta
2450 my_allow_noncompatible_grids = .false.
2451 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
2453 IF (my_beta /= 1.0_dp)
THEN
2454 IF (my_beta == 0.0_dp)
THEN
2458 pw2%array = pw2%array*my_beta
2463 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
2465 IF (my_alpha == 1.0_dp)
THEN
2467 pw2%array = pw2%array + pw1%array
2469 ELSE IF (my_alpha /= 0.0_dp)
THEN
2471 pw2%array = pw2%array + my_alpha* pw1%array
2475 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
2477 ng1 =
SIZE(pw1%array)
2478 ng2 =
SIZE(pw2%array)
2481 IF (pw1%pw_grid%spherical)
THEN
2482 IF (my_alpha == 1.0_dp)
THEN
2485 pw2%array(i) = pw2%array(i) + pw1%array(i)
2488 ELSE IF (my_alpha /= 0.0_dp)
THEN
2491 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(i)
2495 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
2496 IF (ng1 >= ng2)
THEN
2497 IF (my_alpha == 1.0_dp)
THEN
2500 j = pw2%pw_grid%gidx(i)
2501 pw2%array(i) = pw2%array(i) + pw1%array(j)
2504 ELSE IF (my_alpha /= 0.0_dp)
THEN
2507 j = pw2%pw_grid%gidx(i)
2508 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
2513 IF (my_alpha == 1.0_dp)
THEN
2516 j = pw2%pw_grid%gidx(i)
2517 pw2%array(j) = pw2%array(j) + pw1%array(i)
2520 ELSE IF (my_alpha /= 0.0_dp)
THEN
2523 j = pw2%pw_grid%gidx(i)
2524 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
2529 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
2530 IF (ng1 >= ng2)
THEN
2531 IF (my_alpha == 1.0_dp)
THEN
2534 j = pw1%pw_grid%gidx(i)
2535 pw2%array(i) = pw2%array(i) + pw1%array(j)
2538 ELSE IF (my_alpha /= 0.0_dp)
THEN
2541 j = pw1%pw_grid%gidx(i)
2542 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
2547 IF (my_alpha == 1.0_dp)
THEN
2550 j = pw1%pw_grid%gidx(i)
2551 pw2%array(j) = pw2%array(j) + pw1%array(i)
2554 ELSE IF (my_alpha /= 0.0_dp)
THEN
2557 j = pw1%pw_grid%gidx(i)
2558 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
2564 cpabort(
"Grids not compatible")
2569 cpabort(
"Grids not compatible")
2573 CALL timestop(handle)
2575 END SUBROUTINE pw_axpy_r1d_r1d_gs
2586 SUBROUTINE pw_multiply_r1d_r1d_gs (pw_out, pw1, pw2, alpha)
2588 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw_out
2589 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
2590 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
2591 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
2593 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
2596 REAL(KIND=dp) :: my_alpha
2598 CALL timeset(routinen, handle)
2601 IF (
PRESENT(alpha)) my_alpha = alpha
2603 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
2604 cpabort(
"pw_multiply not implemented for non-identical grids!")
2606#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
2607 IF (my_alpha == 1.0_dp)
THEN
2609 pw_out%array = pw_out%array + pw1%array* pw2%array
2613 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
2617 IF (my_alpha == 1.0_dp)
THEN
2618 pw_out%array = pw_out%array + pw1%array* pw2%array
2620 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
2624 CALL timestop(handle)
2626 END SUBROUTINE pw_multiply_r1d_r1d_gs
2633 SUBROUTINE pw_multiply_with_r1d_r1d_gs (pw1, pw2)
2634 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw1
2635 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
2637 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
2641 CALL timeset(routinen, handle)
2643 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
2644 cpabort(
"Incompatible grids!")
2647 pw1%array = pw1%array* pw2%array
2650 CALL timestop(handle)
2652 END SUBROUTINE pw_multiply_with_r1d_r1d_gs
2667 FUNCTION pw_integral_ab_r1d_r1d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
2669 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
2670 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
2671 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
2672 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
2673 REAL(kind=dp) :: integral_value
2675 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_r1d_r1d_gs'
2677 INTEGER :: handle, loc_sumtype
2678 LOGICAL :: my_just_sum, my_local_only
2680 CALL timeset(routinen, handle)
2683 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
2685 my_local_only = .false.
2686 IF (
PRESENT(local_only)) my_local_only = local_only
2688 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
2689 cpabort(
"Grids incompatible")
2692 my_just_sum = .false.
2693 IF (
PRESENT(just_sum)) my_just_sum = just_sum
2700 integral_value = dot_product(pw1%array, pw2%array)
2705 integral_value = accurate_dot_product(pw1%array, pw2%array)
2711 IF (pw1%pw_grid%grid_span == halfspace)
THEN
2712 integral_value = 2.0_dp*integral_value
2713 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
2714 pw1%array(1)*pw2%array(1)
2717 IF (.NOT. my_just_sum)
THEN
2718 integral_value = integral_value*pw1%pw_grid%vol
2721 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
2722 CALL pw1%pw_grid%para%group%sum(integral_value)
2724 CALL timestop(handle)
2726 END FUNCTION pw_integral_ab_r1d_r1d_gs
2734 FUNCTION pw_integral_a2b_r1d_r1d (pw1, pw2)
RESULT(integral_value)
2736 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
2737 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
2738 REAL(kind=dp) :: integral_value
2740 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_a2b'
2744 CALL timeset(routinen, handle)
2746 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
2747 cpabort(
"Grids incompatible")
2750 integral_value = accurate_sum(pw1%array*pw2%array*pw1%pw_grid%gsq)
2751 IF (pw1%pw_grid%grid_span == halfspace)
THEN
2752 integral_value = 2.0_dp*integral_value
2755 integral_value = integral_value*pw1%pw_grid%vol
2757 IF (pw1%pw_grid%para%mode == pw_mode_distributed) &
2758 CALL pw1%pw_grid%para%group%sum(integral_value)
2759 CALL timestop(handle)
2761 END FUNCTION pw_integral_a2b_r1d_r1d
2775 SUBROUTINE pw_copy_r1d_c1d_rs (pw1, pw2)
2777 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
2778 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw2
2780 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
2783 INTEGER :: i, j, ng, ng1, ng2, ns
2785 CALL timeset(routinen, handle)
2787 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
2788 cpabort(
"Both grids must be either spherical or non-spherical!")
2789 IF (pw1%pw_grid%spherical) &
2790 cpabort(
"Spherical grids only exist in reciprocal space!")
2792 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
2793 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
2794 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
2795 ng1 =
SIZE(pw1%array)
2796 ng2 =
SIZE(pw2%array)
2799 pw2%array(1:ng) = cmplx(pw1%array(1:ng), 0.0_dp, kind=dp)
2803 pw2%array(ng + 1:ng2) = cmplx(0.0_dp, 0.0_dp, kind=dp)
2807 cpabort(
"Copies between spherical grids require compatible grids!")
2810 ng1 =
SIZE(pw1%array)
2811 ng2 =
SIZE(pw2%array)
2812 ns = 2*max(ng1, ng2)
2814 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
2815 IF (ng1 >= ng2)
THEN
2820 j = pw2%pw_grid%gidx(i)
2821 pw2%array(i) = cmplx(pw1%array(j), 0.0_dp, kind=dp)
2830 j = pw2%pw_grid%gidx(i)
2831 pw2%array(j) = cmplx(pw1%array(i), 0.0_dp, kind=dp)
2835 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
2836 IF (ng1 >= ng2)
THEN
2839 j = pw1%pw_grid%gidx(i)
2840 pw2%array(i) = cmplx(pw1%array(j), 0.0_dp, kind=dp)
2847 j = pw1%pw_grid%gidx(i)
2848 pw2%array(j) = cmplx(pw1%array(i), 0.0_dp, kind=dp)
2853 cpabort(
"Copy not implemented!")
2860 pw2%array = cmplx(pw1%array, 0.0_dp, kind=dp)
2864 CALL timestop(handle)
2866 END SUBROUTINE pw_copy_r1d_c1d_rs
2873 SUBROUTINE pw_copy_to_array_r1d_c1d_rs (pw, array)
2874 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw
2875 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(INOUT) :: array
2877 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
2881 CALL timeset(routinen, handle)
2884 array(:) = cmplx(pw%array(:), 0.0_dp, kind=dp)
2887 CALL timestop(handle)
2888 END SUBROUTINE pw_copy_to_array_r1d_c1d_rs
2895 SUBROUTINE pw_copy_from_array_r1d_c1d_rs (pw, array)
2896 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw
2897 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(IN) :: array
2899 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
2903 CALL timeset(routinen, handle)
2906 pw%array = real(array, kind=dp)
2909 CALL timestop(handle)
2910 END SUBROUTINE pw_copy_from_array_r1d_c1d_rs
2928 SUBROUTINE pw_axpy_r1d_c1d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
2930 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
2931 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw2
2932 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
2933 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
2935 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
2938 LOGICAL :: my_allow_noncompatible_grids
2939 REAL(KIND=dp) :: my_alpha, my_beta
2940 INTEGER :: i, j, ng, ng1, ng2
2942 CALL timeset(routinen, handle)
2945 IF (
PRESENT(alpha)) my_alpha = alpha
2948 IF (
PRESENT(beta)) my_beta = beta
2950 my_allow_noncompatible_grids = .false.
2951 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
2953 IF (my_beta /= 1.0_dp)
THEN
2954 IF (my_beta == 0.0_dp)
THEN
2958 pw2%array = pw2%array*my_beta
2963 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
2965 IF (my_alpha == 1.0_dp)
THEN
2967 pw2%array = pw2%array + cmplx(pw1%array, 0.0_dp, kind=dp)
2969 ELSE IF (my_alpha /= 0.0_dp)
THEN
2971 pw2%array = pw2%array + my_alpha* cmplx(pw1%array, 0.0_dp, kind=dp)
2975 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
2977 ng1 =
SIZE(pw1%array)
2978 ng2 =
SIZE(pw2%array)
2981 IF (pw1%pw_grid%spherical)
THEN
2982 IF (my_alpha == 1.0_dp)
THEN
2985 pw2%array(i) = pw2%array(i) + cmplx(pw1%array(i), 0.0_dp, kind=dp)
2988 ELSE IF (my_alpha /= 0.0_dp)
THEN
2991 pw2%array(i) = pw2%array(i) + my_alpha* cmplx(pw1%array(i), 0.0_dp, kind=dp)
2995 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
2996 IF (ng1 >= ng2)
THEN
2997 IF (my_alpha == 1.0_dp)
THEN
3000 j = pw2%pw_grid%gidx(i)
3001 pw2%array(i) = pw2%array(i) + cmplx(pw1%array(j), 0.0_dp, kind=dp)
3004 ELSE IF (my_alpha /= 0.0_dp)
THEN
3007 j = pw2%pw_grid%gidx(i)
3008 pw2%array(i) = pw2%array(i) + my_alpha* cmplx(pw1%array(j), 0.0_dp, kind=dp)
3013 IF (my_alpha == 1.0_dp)
THEN
3016 j = pw2%pw_grid%gidx(i)
3017 pw2%array(j) = pw2%array(j) + cmplx(pw1%array(i), 0.0_dp, kind=dp)
3020 ELSE IF (my_alpha /= 0.0_dp)
THEN
3023 j = pw2%pw_grid%gidx(i)
3024 pw2%array(j) = pw2%array(j) + my_alpha* cmplx(pw1%array(i), 0.0_dp, kind=dp)
3029 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
3030 IF (ng1 >= ng2)
THEN
3031 IF (my_alpha == 1.0_dp)
THEN
3034 j = pw1%pw_grid%gidx(i)
3035 pw2%array(i) = pw2%array(i) + cmplx(pw1%array(j), 0.0_dp, kind=dp)
3038 ELSE IF (my_alpha /= 0.0_dp)
THEN
3041 j = pw1%pw_grid%gidx(i)
3042 pw2%array(i) = pw2%array(i) + my_alpha* cmplx(pw1%array(j), 0.0_dp, kind=dp)
3047 IF (my_alpha == 1.0_dp)
THEN
3050 j = pw1%pw_grid%gidx(i)
3051 pw2%array(j) = pw2%array(j) + cmplx(pw1%array(i), 0.0_dp, kind=dp)
3054 ELSE IF (my_alpha /= 0.0_dp)
THEN
3057 j = pw1%pw_grid%gidx(i)
3058 pw2%array(j) = pw2%array(j) + my_alpha* cmplx(pw1%array(i), 0.0_dp, kind=dp)
3064 cpabort(
"Grids not compatible")
3069 cpabort(
"Grids not compatible")
3073 CALL timestop(handle)
3075 END SUBROUTINE pw_axpy_r1d_c1d_rs
3086 SUBROUTINE pw_multiply_r1d_c1d_rs (pw_out, pw1, pw2, alpha)
3088 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw_out
3089 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
3090 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw2
3091 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
3093 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
3096 REAL(KIND=dp) :: my_alpha
3098 CALL timeset(routinen, handle)
3101 IF (
PRESENT(alpha)) my_alpha = alpha
3103 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
3104 cpabort(
"pw_multiply not implemented for non-identical grids!")
3106#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
3107 IF (my_alpha == 1.0_dp)
THEN
3109 pw_out%array = pw_out%array + pw1%array* real(pw2%array, kind=dp)
3113 pw_out%array = pw_out%array + my_alpha*pw1%array* real(pw2%array, kind=dp)
3117 IF (my_alpha == 1.0_dp)
THEN
3118 pw_out%array = pw_out%array + pw1%array* real(pw2%array, kind=dp)
3120 pw_out%array = pw_out%array + my_alpha*pw1%array* real(pw2%array, kind=dp)
3124 CALL timestop(handle)
3126 END SUBROUTINE pw_multiply_r1d_c1d_rs
3133 SUBROUTINE pw_multiply_with_r1d_c1d_rs (pw1, pw2)
3134 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw1
3135 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw2
3137 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
3141 CALL timeset(routinen, handle)
3143 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
3144 cpabort(
"Incompatible grids!")
3147 pw1%array = pw1%array* real(pw2%array, kind=dp)
3150 CALL timestop(handle)
3152 END SUBROUTINE pw_multiply_with_r1d_c1d_rs
3167 FUNCTION pw_integral_ab_r1d_c1d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
3169 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
3170 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw2
3171 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
3172 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
3173 REAL(kind=dp) :: integral_value
3175 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_r1d_c1d_rs'
3177 INTEGER :: handle, loc_sumtype
3178 LOGICAL :: my_just_sum, my_local_only
3180 CALL timeset(routinen, handle)
3183 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
3185 my_local_only = .false.
3186 IF (
PRESENT(local_only)) my_local_only = local_only
3188 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
3189 cpabort(
"Grids incompatible")
3192 my_just_sum = .false.
3193 IF (
PRESENT(just_sum)) my_just_sum = just_sum
3200 integral_value = sum(pw1%array*real(pw2%array, kind=dp))
3205 integral_value = accurate_sum(pw1%array*real(pw2%array, kind=dp))
3211 IF (pw1%pw_grid%grid_span == halfspace)
THEN
3212 integral_value = 2.0_dp*integral_value
3213 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
3214 pw1%array(1)*real(pw2%array(1), kind=dp)
3217 IF (.NOT. my_just_sum)
THEN
3218 integral_value = integral_value*pw1%pw_grid%dvol
3221 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
3222 CALL pw1%pw_grid%para%group%sum(integral_value)
3224 CALL timestop(handle)
3226 END FUNCTION pw_integral_ab_r1d_c1d_rs
3240 SUBROUTINE pw_copy_r1d_c1d_gs (pw1, pw2)
3242 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
3243 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
3245 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
3248 INTEGER :: i, j, ng, ng1, ng2, ns
3250 CALL timeset(routinen, handle)
3252 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
3253 cpabort(
"Both grids must be either spherical or non-spherical!")
3255 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
3256 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
3257 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
3258 ng1 =
SIZE(pw1%array)
3259 ng2 =
SIZE(pw2%array)
3262 pw2%array(1:ng) = cmplx(pw1%array(1:ng), 0.0_dp, kind=dp)
3266 pw2%array(ng + 1:ng2) = cmplx(0.0_dp, 0.0_dp, kind=dp)
3270 cpabort(
"Copies between spherical grids require compatible grids!")
3273 ng1 =
SIZE(pw1%array)
3274 ng2 =
SIZE(pw2%array)
3275 ns = 2*max(ng1, ng2)
3277 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
3278 IF (ng1 >= ng2)
THEN
3283 j = pw2%pw_grid%gidx(i)
3284 pw2%array(i) = cmplx(pw1%array(j), 0.0_dp, kind=dp)
3293 j = pw2%pw_grid%gidx(i)
3294 pw2%array(j) = cmplx(pw1%array(i), 0.0_dp, kind=dp)
3298 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
3299 IF (ng1 >= ng2)
THEN
3302 j = pw1%pw_grid%gidx(i)
3303 pw2%array(i) = cmplx(pw1%array(j), 0.0_dp, kind=dp)
3310 j = pw1%pw_grid%gidx(i)
3311 pw2%array(j) = cmplx(pw1%array(i), 0.0_dp, kind=dp)
3316 cpabort(
"Copy not implemented!")
3323 pw2%array = cmplx(pw1%array, 0.0_dp, kind=dp)
3327 CALL timestop(handle)
3329 END SUBROUTINE pw_copy_r1d_c1d_gs
3336 SUBROUTINE pw_copy_to_array_r1d_c1d_gs (pw, array)
3337 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw
3338 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(INOUT) :: array
3340 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
3344 CALL timeset(routinen, handle)
3347 array(:) = cmplx(pw%array(:), 0.0_dp, kind=dp)
3350 CALL timestop(handle)
3351 END SUBROUTINE pw_copy_to_array_r1d_c1d_gs
3358 SUBROUTINE pw_copy_from_array_r1d_c1d_gs (pw, array)
3359 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw
3360 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(IN) :: array
3362 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
3366 CALL timeset(routinen, handle)
3369 pw%array = real(array, kind=dp)
3372 CALL timestop(handle)
3373 END SUBROUTINE pw_copy_from_array_r1d_c1d_gs
3391 SUBROUTINE pw_axpy_r1d_c1d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
3393 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
3394 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
3395 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
3396 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
3398 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
3401 LOGICAL :: my_allow_noncompatible_grids
3402 REAL(KIND=dp) :: my_alpha, my_beta
3403 INTEGER :: i, j, ng, ng1, ng2
3405 CALL timeset(routinen, handle)
3408 IF (
PRESENT(alpha)) my_alpha = alpha
3411 IF (
PRESENT(beta)) my_beta = beta
3413 my_allow_noncompatible_grids = .false.
3414 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
3416 IF (my_beta /= 1.0_dp)
THEN
3417 IF (my_beta == 0.0_dp)
THEN
3421 pw2%array = pw2%array*my_beta
3426 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
3428 IF (my_alpha == 1.0_dp)
THEN
3430 pw2%array = pw2%array + cmplx(pw1%array, 0.0_dp, kind=dp)
3432 ELSE IF (my_alpha /= 0.0_dp)
THEN
3434 pw2%array = pw2%array + my_alpha* cmplx(pw1%array, 0.0_dp, kind=dp)
3438 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
3440 ng1 =
SIZE(pw1%array)
3441 ng2 =
SIZE(pw2%array)
3444 IF (pw1%pw_grid%spherical)
THEN
3445 IF (my_alpha == 1.0_dp)
THEN
3448 pw2%array(i) = pw2%array(i) + cmplx(pw1%array(i), 0.0_dp, kind=dp)
3451 ELSE IF (my_alpha /= 0.0_dp)
THEN
3454 pw2%array(i) = pw2%array(i) + my_alpha* cmplx(pw1%array(i), 0.0_dp, kind=dp)
3458 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
3459 IF (ng1 >= ng2)
THEN
3460 IF (my_alpha == 1.0_dp)
THEN
3463 j = pw2%pw_grid%gidx(i)
3464 pw2%array(i) = pw2%array(i) + cmplx(pw1%array(j), 0.0_dp, kind=dp)
3467 ELSE IF (my_alpha /= 0.0_dp)
THEN
3470 j = pw2%pw_grid%gidx(i)
3471 pw2%array(i) = pw2%array(i) + my_alpha* cmplx(pw1%array(j), 0.0_dp, kind=dp)
3476 IF (my_alpha == 1.0_dp)
THEN
3479 j = pw2%pw_grid%gidx(i)
3480 pw2%array(j) = pw2%array(j) + cmplx(pw1%array(i), 0.0_dp, kind=dp)
3483 ELSE IF (my_alpha /= 0.0_dp)
THEN
3486 j = pw2%pw_grid%gidx(i)
3487 pw2%array(j) = pw2%array(j) + my_alpha* cmplx(pw1%array(i), 0.0_dp, kind=dp)
3492 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
3493 IF (ng1 >= ng2)
THEN
3494 IF (my_alpha == 1.0_dp)
THEN
3497 j = pw1%pw_grid%gidx(i)
3498 pw2%array(i) = pw2%array(i) + cmplx(pw1%array(j), 0.0_dp, kind=dp)
3501 ELSE IF (my_alpha /= 0.0_dp)
THEN
3504 j = pw1%pw_grid%gidx(i)
3505 pw2%array(i) = pw2%array(i) + my_alpha* cmplx(pw1%array(j), 0.0_dp, kind=dp)
3510 IF (my_alpha == 1.0_dp)
THEN
3513 j = pw1%pw_grid%gidx(i)
3514 pw2%array(j) = pw2%array(j) + cmplx(pw1%array(i), 0.0_dp, kind=dp)
3517 ELSE IF (my_alpha /= 0.0_dp)
THEN
3520 j = pw1%pw_grid%gidx(i)
3521 pw2%array(j) = pw2%array(j) + my_alpha* cmplx(pw1%array(i), 0.0_dp, kind=dp)
3527 cpabort(
"Grids not compatible")
3532 cpabort(
"Grids not compatible")
3536 CALL timestop(handle)
3538 END SUBROUTINE pw_axpy_r1d_c1d_gs
3549 SUBROUTINE pw_multiply_r1d_c1d_gs (pw_out, pw1, pw2, alpha)
3551 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw_out
3552 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
3553 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
3554 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
3556 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
3559 REAL(KIND=dp) :: my_alpha
3561 CALL timeset(routinen, handle)
3564 IF (
PRESENT(alpha)) my_alpha = alpha
3566 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
3567 cpabort(
"pw_multiply not implemented for non-identical grids!")
3569#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
3570 IF (my_alpha == 1.0_dp)
THEN
3572 pw_out%array = pw_out%array + pw1%array* real(pw2%array, kind=dp)
3576 pw_out%array = pw_out%array + my_alpha*pw1%array* real(pw2%array, kind=dp)
3580 IF (my_alpha == 1.0_dp)
THEN
3581 pw_out%array = pw_out%array + pw1%array* real(pw2%array, kind=dp)
3583 pw_out%array = pw_out%array + my_alpha*pw1%array* real(pw2%array, kind=dp)
3587 CALL timestop(handle)
3589 END SUBROUTINE pw_multiply_r1d_c1d_gs
3596 SUBROUTINE pw_multiply_with_r1d_c1d_gs (pw1, pw2)
3597 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw1
3598 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
3600 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
3604 CALL timeset(routinen, handle)
3606 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
3607 cpabort(
"Incompatible grids!")
3610 pw1%array = pw1%array* real(pw2%array, kind=dp)
3613 CALL timestop(handle)
3615 END SUBROUTINE pw_multiply_with_r1d_c1d_gs
3630 FUNCTION pw_integral_ab_r1d_c1d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
3632 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
3633 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
3634 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
3635 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
3636 REAL(kind=dp) :: integral_value
3638 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_r1d_c1d_gs'
3640 INTEGER :: handle, loc_sumtype
3641 LOGICAL :: my_just_sum, my_local_only
3643 CALL timeset(routinen, handle)
3646 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
3648 my_local_only = .false.
3649 IF (
PRESENT(local_only)) my_local_only = local_only
3651 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
3652 cpabort(
"Grids incompatible")
3655 my_just_sum = .false.
3656 IF (
PRESENT(just_sum)) my_just_sum = just_sum
3663 integral_value = sum(pw1%array*real(pw2%array, kind=dp))
3668 integral_value = accurate_sum(pw1%array*real(pw2%array, kind=dp))
3674 IF (pw1%pw_grid%grid_span == halfspace)
THEN
3675 integral_value = 2.0_dp*integral_value
3676 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
3677 pw1%array(1)*real(pw2%array(1), kind=dp)
3680 IF (.NOT. my_just_sum)
THEN
3681 integral_value = integral_value*pw1%pw_grid%vol
3684 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
3685 CALL pw1%pw_grid%para%group%sum(integral_value)
3687 CALL timestop(handle)
3689 END FUNCTION pw_integral_ab_r1d_c1d_gs
3697 FUNCTION pw_integral_a2b_r1d_c1d (pw1, pw2)
RESULT(integral_value)
3699 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
3700 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
3701 REAL(kind=dp) :: integral_value
3703 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_a2b'
3707 CALL timeset(routinen, handle)
3709 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
3710 cpabort(
"Grids incompatible")
3713 integral_value = accurate_sum(pw1%array*real(pw2%array, kind=dp)*pw1%pw_grid%gsq)
3714 IF (pw1%pw_grid%grid_span == halfspace)
THEN
3715 integral_value = 2.0_dp*integral_value
3718 integral_value = integral_value*pw1%pw_grid%vol
3720 IF (pw1%pw_grid%para%mode == pw_mode_distributed) &
3721 CALL pw1%pw_grid%para%group%sum(integral_value)
3722 CALL timestop(handle)
3724 END FUNCTION pw_integral_a2b_r1d_c1d
3738 SUBROUTINE pw_copy_r3d_r3d_rs (pw1, pw2)
3740 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
3741 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw2
3743 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
3747 CALL timeset(routinen, handle)
3749 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
3750 cpabort(
"Both grids must be either spherical or non-spherical!")
3751 IF (pw1%pw_grid%spherical) &
3752 cpabort(
"Spherical grids only exist in reciprocal space!")
3754 IF (any(shape(pw2%array) /= shape(pw1%array))) &
3755 cpabort(
"3D grids must be compatible!")
3756 IF (pw1%pw_grid%spherical) &
3757 cpabort(
"3D grids must not be spherical!")
3759 pw2%array(:, :, :) = pw1%array(:, :, :)
3762 CALL timestop(handle)
3764 END SUBROUTINE pw_copy_r3d_r3d_rs
3771 SUBROUTINE pw_copy_to_array_r3d_r3d_rs (pw, array)
3772 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw
3773 REAL(KIND=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
3775 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
3779 CALL timeset(routinen, handle)
3782 array(:, :, :) = pw%array(:, :, :)
3785 CALL timestop(handle)
3786 END SUBROUTINE pw_copy_to_array_r3d_r3d_rs
3793 SUBROUTINE pw_copy_from_array_r3d_r3d_rs (pw, array)
3794 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw
3795 REAL(KIND=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
3797 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
3801 CALL timeset(routinen, handle)
3807 CALL timestop(handle)
3808 END SUBROUTINE pw_copy_from_array_r3d_r3d_rs
3826 SUBROUTINE pw_axpy_r3d_r3d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
3828 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
3829 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw2
3830 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
3831 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
3833 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
3836 LOGICAL :: my_allow_noncompatible_grids
3837 REAL(KIND=dp) :: my_alpha, my_beta
3839 CALL timeset(routinen, handle)
3842 IF (
PRESENT(alpha)) my_alpha = alpha
3845 IF (
PRESENT(beta)) my_beta = beta
3847 my_allow_noncompatible_grids = .false.
3848 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
3850 IF (my_beta /= 1.0_dp)
THEN
3851 IF (my_beta == 0.0_dp)
THEN
3855 pw2%array = pw2%array*my_beta
3860 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
3861 IF (my_alpha == 1.0_dp)
THEN
3863 pw2%array = pw2%array + pw1%array
3865 ELSE IF (my_alpha /= 0.0_dp)
THEN
3867 pw2%array = pw2%array + my_alpha* pw1%array
3871 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
3873 IF (any(shape(pw1%array) /= shape(pw2%array))) &
3874 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
3876 IF (my_alpha == 1.0_dp)
THEN
3878 pw2%array = pw2%array + pw1%array
3882 pw2%array = pw2%array + my_alpha* pw1%array
3888 cpabort(
"Grids not compatible")
3892 CALL timestop(handle)
3894 END SUBROUTINE pw_axpy_r3d_r3d_rs
3905 SUBROUTINE pw_multiply_r3d_r3d_rs (pw_out, pw1, pw2, alpha)
3907 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw_out
3908 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
3909 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw2
3910 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
3912 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
3915 REAL(KIND=dp) :: my_alpha
3917 CALL timeset(routinen, handle)
3920 IF (
PRESENT(alpha)) my_alpha = alpha
3922 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
3923 cpabort(
"pw_multiply not implemented for non-identical grids!")
3925#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
3926 IF (my_alpha == 1.0_dp)
THEN
3928 pw_out%array = pw_out%array + pw1%array* pw2%array
3932 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
3936 IF (my_alpha == 1.0_dp)
THEN
3937 pw_out%array = pw_out%array + pw1%array* pw2%array
3939 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
3943 CALL timestop(handle)
3945 END SUBROUTINE pw_multiply_r3d_r3d_rs
3952 SUBROUTINE pw_multiply_with_r3d_r3d_rs (pw1, pw2)
3953 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw1
3954 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw2
3956 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
3960 CALL timeset(routinen, handle)
3962 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
3963 cpabort(
"Incompatible grids!")
3966 pw1%array = pw1%array* pw2%array
3969 CALL timestop(handle)
3971 END SUBROUTINE pw_multiply_with_r3d_r3d_rs
3986 FUNCTION pw_integral_ab_r3d_r3d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
3988 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
3989 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw2
3990 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
3991 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
3992 REAL(kind=dp) :: integral_value
3994 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_r3d_r3d_rs'
3996 INTEGER :: handle, loc_sumtype
3997 LOGICAL :: my_just_sum, my_local_only
3999 CALL timeset(routinen, handle)
4002 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
4004 my_local_only = .false.
4005 IF (
PRESENT(local_only)) my_local_only = local_only
4007 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
4008 cpabort(
"Grids incompatible")
4011 my_just_sum = .false.
4012 IF (
PRESENT(just_sum)) my_just_sum = just_sum
4019 integral_value = sum(pw1%array*pw2%array)
4024 integral_value = accurate_dot_product(pw1%array, pw2%array)
4029 IF (.NOT. my_just_sum)
THEN
4030 integral_value = integral_value*pw1%pw_grid%dvol
4033 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
4034 CALL pw1%pw_grid%para%group%sum(integral_value)
4036 CALL timestop(handle)
4038 END FUNCTION pw_integral_ab_r3d_r3d_rs
4052 SUBROUTINE pw_copy_r3d_r3d_gs (pw1, pw2)
4054 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4055 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw2
4057 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
4061 CALL timeset(routinen, handle)
4063 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
4064 cpabort(
"Both grids must be either spherical or non-spherical!")
4066 IF (any(shape(pw2%array) /= shape(pw1%array))) &
4067 cpabort(
"3D grids must be compatible!")
4068 IF (pw1%pw_grid%spherical) &
4069 cpabort(
"3D grids must not be spherical!")
4071 pw2%array(:, :, :) = pw1%array(:, :, :)
4074 CALL timestop(handle)
4076 END SUBROUTINE pw_copy_r3d_r3d_gs
4083 SUBROUTINE pw_copy_to_array_r3d_r3d_gs (pw, array)
4084 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw
4085 REAL(KIND=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
4087 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
4091 CALL timeset(routinen, handle)
4094 array(:, :, :) = pw%array(:, :, :)
4097 CALL timestop(handle)
4098 END SUBROUTINE pw_copy_to_array_r3d_r3d_gs
4105 SUBROUTINE pw_copy_from_array_r3d_r3d_gs (pw, array)
4106 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw
4107 REAL(KIND=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
4109 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
4113 CALL timeset(routinen, handle)
4119 CALL timestop(handle)
4120 END SUBROUTINE pw_copy_from_array_r3d_r3d_gs
4138 SUBROUTINE pw_axpy_r3d_r3d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
4140 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4141 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw2
4142 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
4143 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
4145 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
4148 LOGICAL :: my_allow_noncompatible_grids
4149 REAL(KIND=dp) :: my_alpha, my_beta
4151 CALL timeset(routinen, handle)
4154 IF (
PRESENT(alpha)) my_alpha = alpha
4157 IF (
PRESENT(beta)) my_beta = beta
4159 my_allow_noncompatible_grids = .false.
4160 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
4162 IF (my_beta /= 1.0_dp)
THEN
4163 IF (my_beta == 0.0_dp)
THEN
4167 pw2%array = pw2%array*my_beta
4172 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
4173 IF (my_alpha == 1.0_dp)
THEN
4175 pw2%array = pw2%array + pw1%array
4177 ELSE IF (my_alpha /= 0.0_dp)
THEN
4179 pw2%array = pw2%array + my_alpha* pw1%array
4183 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
4185 IF (any(shape(pw1%array) /= shape(pw2%array))) &
4186 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
4188 IF (my_alpha == 1.0_dp)
THEN
4190 pw2%array = pw2%array + pw1%array
4194 pw2%array = pw2%array + my_alpha* pw1%array
4200 cpabort(
"Grids not compatible")
4204 CALL timestop(handle)
4206 END SUBROUTINE pw_axpy_r3d_r3d_gs
4217 SUBROUTINE pw_multiply_r3d_r3d_gs (pw_out, pw1, pw2, alpha)
4219 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw_out
4220 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4221 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw2
4222 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
4224 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
4227 REAL(KIND=dp) :: my_alpha
4229 CALL timeset(routinen, handle)
4232 IF (
PRESENT(alpha)) my_alpha = alpha
4234 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
4235 cpabort(
"pw_multiply not implemented for non-identical grids!")
4237#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
4238 IF (my_alpha == 1.0_dp)
THEN
4240 pw_out%array = pw_out%array + pw1%array* pw2%array
4244 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
4248 IF (my_alpha == 1.0_dp)
THEN
4249 pw_out%array = pw_out%array + pw1%array* pw2%array
4251 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
4255 CALL timestop(handle)
4257 END SUBROUTINE pw_multiply_r3d_r3d_gs
4264 SUBROUTINE pw_multiply_with_r3d_r3d_gs (pw1, pw2)
4265 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw1
4266 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw2
4268 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
4272 CALL timeset(routinen, handle)
4274 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
4275 cpabort(
"Incompatible grids!")
4278 pw1%array = pw1%array* pw2%array
4281 CALL timestop(handle)
4283 END SUBROUTINE pw_multiply_with_r3d_r3d_gs
4298 FUNCTION pw_integral_ab_r3d_r3d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
4300 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4301 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw2
4302 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
4303 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
4304 REAL(kind=dp) :: integral_value
4306 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_r3d_r3d_gs'
4308 INTEGER :: handle, loc_sumtype
4309 LOGICAL :: my_just_sum, my_local_only
4311 CALL timeset(routinen, handle)
4314 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
4316 my_local_only = .false.
4317 IF (
PRESENT(local_only)) my_local_only = local_only
4319 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
4320 cpabort(
"Grids incompatible")
4323 my_just_sum = .false.
4324 IF (
PRESENT(just_sum)) my_just_sum = just_sum
4331 integral_value = sum(pw1%array*pw2%array)
4336 integral_value = accurate_dot_product(pw1%array, pw2%array)
4341 IF (.NOT. my_just_sum)
THEN
4342 integral_value = integral_value*pw1%pw_grid%vol
4345 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
4346 CALL pw1%pw_grid%para%group%sum(integral_value)
4348 CALL timestop(handle)
4350 END FUNCTION pw_integral_ab_r3d_r3d_gs
4365 SUBROUTINE pw_copy_r3d_c3d_rs (pw1, pw2)
4367 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
4368 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw2
4370 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
4374 CALL timeset(routinen, handle)
4376 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
4377 cpabort(
"Both grids must be either spherical or non-spherical!")
4378 IF (pw1%pw_grid%spherical) &
4379 cpabort(
"Spherical grids only exist in reciprocal space!")
4381 IF (any(shape(pw2%array) /= shape(pw1%array))) &
4382 cpabort(
"3D grids must be compatible!")
4383 IF (pw1%pw_grid%spherical) &
4384 cpabort(
"3D grids must not be spherical!")
4386 pw2%array(:, :, :) = cmplx(pw1%array(:, :, :), 0.0_dp, kind=dp)
4389 CALL timestop(handle)
4391 END SUBROUTINE pw_copy_r3d_c3d_rs
4398 SUBROUTINE pw_copy_to_array_r3d_c3d_rs (pw, array)
4399 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw
4400 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
4402 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
4406 CALL timeset(routinen, handle)
4409 array(:, :, :) = cmplx(pw%array(:, :, :), 0.0_dp, kind=dp)
4412 CALL timestop(handle)
4413 END SUBROUTINE pw_copy_to_array_r3d_c3d_rs
4420 SUBROUTINE pw_copy_from_array_r3d_c3d_rs (pw, array)
4421 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw
4422 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
4424 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
4428 CALL timeset(routinen, handle)
4431 pw%array = real(array, kind=dp)
4434 CALL timestop(handle)
4435 END SUBROUTINE pw_copy_from_array_r3d_c3d_rs
4453 SUBROUTINE pw_axpy_r3d_c3d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
4455 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
4456 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw2
4457 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
4458 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
4460 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
4463 LOGICAL :: my_allow_noncompatible_grids
4464 REAL(KIND=dp) :: my_alpha, my_beta
4466 CALL timeset(routinen, handle)
4469 IF (
PRESENT(alpha)) my_alpha = alpha
4472 IF (
PRESENT(beta)) my_beta = beta
4474 my_allow_noncompatible_grids = .false.
4475 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
4477 IF (my_beta /= 1.0_dp)
THEN
4478 IF (my_beta == 0.0_dp)
THEN
4482 pw2%array = pw2%array*my_beta
4487 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
4488 IF (my_alpha == 1.0_dp)
THEN
4490 pw2%array = pw2%array + cmplx(pw1%array, 0.0_dp, kind=dp)
4492 ELSE IF (my_alpha /= 0.0_dp)
THEN
4494 pw2%array = pw2%array + my_alpha* cmplx(pw1%array, 0.0_dp, kind=dp)
4498 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
4500 IF (any(shape(pw1%array) /= shape(pw2%array))) &
4501 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
4503 IF (my_alpha == 1.0_dp)
THEN
4505 pw2%array = pw2%array + cmplx(pw1%array, 0.0_dp, kind=dp)
4509 pw2%array = pw2%array + my_alpha* cmplx(pw1%array, 0.0_dp, kind=dp)
4515 cpabort(
"Grids not compatible")
4519 CALL timestop(handle)
4521 END SUBROUTINE pw_axpy_r3d_c3d_rs
4532 SUBROUTINE pw_multiply_r3d_c3d_rs (pw_out, pw1, pw2, alpha)
4534 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw_out
4535 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
4536 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw2
4537 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
4539 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
4542 REAL(KIND=dp) :: my_alpha
4544 CALL timeset(routinen, handle)
4547 IF (
PRESENT(alpha)) my_alpha = alpha
4549 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
4550 cpabort(
"pw_multiply not implemented for non-identical grids!")
4552#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
4553 IF (my_alpha == 1.0_dp)
THEN
4555 pw_out%array = pw_out%array + pw1%array* real(pw2%array, kind=dp)
4559 pw_out%array = pw_out%array + my_alpha*pw1%array* real(pw2%array, kind=dp)
4563 IF (my_alpha == 1.0_dp)
THEN
4564 pw_out%array = pw_out%array + pw1%array* real(pw2%array, kind=dp)
4566 pw_out%array = pw_out%array + my_alpha*pw1%array* real(pw2%array, kind=dp)
4570 CALL timestop(handle)
4572 END SUBROUTINE pw_multiply_r3d_c3d_rs
4579 SUBROUTINE pw_multiply_with_r3d_c3d_rs (pw1, pw2)
4580 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw1
4581 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw2
4583 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
4587 CALL timeset(routinen, handle)
4589 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
4590 cpabort(
"Incompatible grids!")
4593 pw1%array = pw1%array* real(pw2%array, kind=dp)
4596 CALL timestop(handle)
4598 END SUBROUTINE pw_multiply_with_r3d_c3d_rs
4613 FUNCTION pw_integral_ab_r3d_c3d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
4615 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
4616 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw2
4617 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
4618 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
4619 REAL(kind=dp) :: integral_value
4621 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_r3d_c3d_rs'
4623 INTEGER :: handle, loc_sumtype
4624 LOGICAL :: my_just_sum, my_local_only
4626 CALL timeset(routinen, handle)
4629 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
4631 my_local_only = .false.
4632 IF (
PRESENT(local_only)) my_local_only = local_only
4634 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
4635 cpabort(
"Grids incompatible")
4638 my_just_sum = .false.
4639 IF (
PRESENT(just_sum)) my_just_sum = just_sum
4646 integral_value = sum(pw1%array*real(pw2%array, kind=dp))
4651 integral_value = accurate_sum(pw1%array*real(pw2%array, kind=dp))
4656 IF (.NOT. my_just_sum)
THEN
4657 integral_value = integral_value*pw1%pw_grid%dvol
4660 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
4661 CALL pw1%pw_grid%para%group%sum(integral_value)
4663 CALL timestop(handle)
4665 END FUNCTION pw_integral_ab_r3d_c3d_rs
4679 SUBROUTINE pw_copy_r3d_c3d_gs (pw1, pw2)
4681 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4682 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
4684 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
4688 CALL timeset(routinen, handle)
4690 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
4691 cpabort(
"Both grids must be either spherical or non-spherical!")
4693 IF (any(shape(pw2%array) /= shape(pw1%array))) &
4694 cpabort(
"3D grids must be compatible!")
4695 IF (pw1%pw_grid%spherical) &
4696 cpabort(
"3D grids must not be spherical!")
4698 pw2%array(:, :, :) = cmplx(pw1%array(:, :, :), 0.0_dp, kind=dp)
4701 CALL timestop(handle)
4703 END SUBROUTINE pw_copy_r3d_c3d_gs
4710 SUBROUTINE pw_copy_to_array_r3d_c3d_gs (pw, array)
4711 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw
4712 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
4714 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
4718 CALL timeset(routinen, handle)
4721 array(:, :, :) = cmplx(pw%array(:, :, :), 0.0_dp, kind=dp)
4724 CALL timestop(handle)
4725 END SUBROUTINE pw_copy_to_array_r3d_c3d_gs
4732 SUBROUTINE pw_copy_from_array_r3d_c3d_gs (pw, array)
4733 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw
4734 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
4736 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
4740 CALL timeset(routinen, handle)
4743 pw%array = real(array, kind=dp)
4746 CALL timestop(handle)
4747 END SUBROUTINE pw_copy_from_array_r3d_c3d_gs
4765 SUBROUTINE pw_axpy_r3d_c3d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
4767 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4768 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
4769 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
4770 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
4772 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
4775 LOGICAL :: my_allow_noncompatible_grids
4776 REAL(KIND=dp) :: my_alpha, my_beta
4778 CALL timeset(routinen, handle)
4781 IF (
PRESENT(alpha)) my_alpha = alpha
4784 IF (
PRESENT(beta)) my_beta = beta
4786 my_allow_noncompatible_grids = .false.
4787 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
4789 IF (my_beta /= 1.0_dp)
THEN
4790 IF (my_beta == 0.0_dp)
THEN
4794 pw2%array = pw2%array*my_beta
4799 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
4800 IF (my_alpha == 1.0_dp)
THEN
4802 pw2%array = pw2%array + cmplx(pw1%array, 0.0_dp, kind=dp)
4804 ELSE IF (my_alpha /= 0.0_dp)
THEN
4806 pw2%array = pw2%array + my_alpha* cmplx(pw1%array, 0.0_dp, kind=dp)
4810 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
4812 IF (any(shape(pw1%array) /= shape(pw2%array))) &
4813 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
4815 IF (my_alpha == 1.0_dp)
THEN
4817 pw2%array = pw2%array + cmplx(pw1%array, 0.0_dp, kind=dp)
4821 pw2%array = pw2%array + my_alpha* cmplx(pw1%array, 0.0_dp, kind=dp)
4827 cpabort(
"Grids not compatible")
4831 CALL timestop(handle)
4833 END SUBROUTINE pw_axpy_r3d_c3d_gs
4844 SUBROUTINE pw_multiply_r3d_c3d_gs (pw_out, pw1, pw2, alpha)
4846 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw_out
4847 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4848 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw2
4849 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
4851 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
4854 REAL(KIND=dp) :: my_alpha
4856 CALL timeset(routinen, handle)
4859 IF (
PRESENT(alpha)) my_alpha = alpha
4861 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
4862 cpabort(
"pw_multiply not implemented for non-identical grids!")
4864#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
4865 IF (my_alpha == 1.0_dp)
THEN
4867 pw_out%array = pw_out%array + pw1%array* real(pw2%array, kind=dp)
4871 pw_out%array = pw_out%array + my_alpha*pw1%array* real(pw2%array, kind=dp)
4875 IF (my_alpha == 1.0_dp)
THEN
4876 pw_out%array = pw_out%array + pw1%array* real(pw2%array, kind=dp)
4878 pw_out%array = pw_out%array + my_alpha*pw1%array* real(pw2%array, kind=dp)
4882 CALL timestop(handle)
4884 END SUBROUTINE pw_multiply_r3d_c3d_gs
4891 SUBROUTINE pw_multiply_with_r3d_c3d_gs (pw1, pw2)
4892 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw1
4893 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw2
4895 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
4899 CALL timeset(routinen, handle)
4901 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
4902 cpabort(
"Incompatible grids!")
4905 pw1%array = pw1%array* real(pw2%array, kind=dp)
4908 CALL timestop(handle)
4910 END SUBROUTINE pw_multiply_with_r3d_c3d_gs
4925 FUNCTION pw_integral_ab_r3d_c3d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
4927 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4928 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw2
4929 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
4930 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
4931 REAL(kind=dp) :: integral_value
4933 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_r3d_c3d_gs'
4935 INTEGER :: handle, loc_sumtype
4936 LOGICAL :: my_just_sum, my_local_only
4938 CALL timeset(routinen, handle)
4941 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
4943 my_local_only = .false.
4944 IF (
PRESENT(local_only)) my_local_only = local_only
4946 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
4947 cpabort(
"Grids incompatible")
4950 my_just_sum = .false.
4951 IF (
PRESENT(just_sum)) my_just_sum = just_sum
4958 integral_value = sum(pw1%array*real(pw2%array, kind=dp))
4963 integral_value = accurate_sum(pw1%array*real(pw2%array, kind=dp))
4968 IF (.NOT. my_just_sum)
THEN
4969 integral_value = integral_value*pw1%pw_grid%vol
4972 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
4973 CALL pw1%pw_grid%para%group%sum(integral_value)
4975 CALL timestop(handle)
4977 END FUNCTION pw_integral_ab_r3d_c3d_gs
4992 SUBROUTINE pw_copy_c1d_r1d_rs (pw1, pw2)
4994 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
4995 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw2
4997 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
5000 INTEGER :: i, j, ng, ng1, ng2, ns
5002 CALL timeset(routinen, handle)
5004 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
5005 cpabort(
"Both grids must be either spherical or non-spherical!")
5006 IF (pw1%pw_grid%spherical) &
5007 cpabort(
"Spherical grids only exist in reciprocal space!")
5009 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
5010 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
5011 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
5012 ng1 =
SIZE(pw1%array)
5013 ng2 =
SIZE(pw2%array)
5016 pw2%array(1:ng) = real(pw1%array(1:ng), kind=dp)
5020 pw2%array(ng + 1:ng2) = 0.0_dp
5024 cpabort(
"Copies between spherical grids require compatible grids!")
5027 ng1 =
SIZE(pw1%array)
5028 ng2 =
SIZE(pw2%array)
5029 ns = 2*max(ng1, ng2)
5031 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
5032 IF (ng1 >= ng2)
THEN
5037 j = pw2%pw_grid%gidx(i)
5038 pw2%array(i) = real(pw1%array(j), kind=dp)
5047 j = pw2%pw_grid%gidx(i)
5048 pw2%array(j) = real(pw1%array(i), kind=dp)
5052 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
5053 IF (ng1 >= ng2)
THEN
5056 j = pw1%pw_grid%gidx(i)
5057 pw2%array(i) = real(pw1%array(j), kind=dp)
5064 j = pw1%pw_grid%gidx(i)
5065 pw2%array(j) = real(pw1%array(i), kind=dp)
5070 cpabort(
"Copy not implemented!")
5077 pw2%array = real(pw1%array, kind=dp)
5081 CALL timestop(handle)
5083 END SUBROUTINE pw_copy_c1d_r1d_rs
5090 SUBROUTINE pw_copy_to_array_c1d_r1d_rs (pw, array)
5091 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw
5092 REAL(KIND=dp),
DIMENSION(:),
INTENT(INOUT) :: array
5094 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
5098 CALL timeset(routinen, handle)
5101 array(:) = real(pw%array(:), kind=dp)
5104 CALL timestop(handle)
5105 END SUBROUTINE pw_copy_to_array_c1d_r1d_rs
5112 SUBROUTINE pw_copy_from_array_c1d_r1d_rs (pw, array)
5113 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw
5114 REAL(KIND=dp),
DIMENSION(:),
INTENT(IN) :: array
5116 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
5120 CALL timeset(routinen, handle)
5123 pw%array = cmplx(array, 0.0_dp, kind=dp)
5126 CALL timestop(handle)
5127 END SUBROUTINE pw_copy_from_array_c1d_r1d_rs
5145 SUBROUTINE pw_axpy_c1d_r1d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
5147 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
5148 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw2
5149 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
5150 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
5152 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
5155 LOGICAL :: my_allow_noncompatible_grids
5156 REAL(KIND=dp) :: my_alpha, my_beta
5157 INTEGER :: i, j, ng, ng1, ng2
5159 CALL timeset(routinen, handle)
5162 IF (
PRESENT(alpha)) my_alpha = alpha
5165 IF (
PRESENT(beta)) my_beta = beta
5167 my_allow_noncompatible_grids = .false.
5168 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
5170 IF (my_beta /= 1.0_dp)
THEN
5171 IF (my_beta == 0.0_dp)
THEN
5175 pw2%array = pw2%array*my_beta
5180 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
5182 IF (my_alpha == 1.0_dp)
THEN
5184 pw2%array = pw2%array + real(pw1%array, kind=dp)
5186 ELSE IF (my_alpha /= 0.0_dp)
THEN
5188 pw2%array = pw2%array + my_alpha* real(pw1%array, kind=dp)
5192 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
5194 ng1 =
SIZE(pw1%array)
5195 ng2 =
SIZE(pw2%array)
5198 IF (pw1%pw_grid%spherical)
THEN
5199 IF (my_alpha == 1.0_dp)
THEN
5202 pw2%array(i) = pw2%array(i) + real(pw1%array(i), kind=dp)
5205 ELSE IF (my_alpha /= 0.0_dp)
THEN
5208 pw2%array(i) = pw2%array(i) + my_alpha* real(pw1%array(i), kind=dp)
5212 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
5213 IF (ng1 >= ng2)
THEN
5214 IF (my_alpha == 1.0_dp)
THEN
5217 j = pw2%pw_grid%gidx(i)
5218 pw2%array(i) = pw2%array(i) + real(pw1%array(j), kind=dp)
5221 ELSE IF (my_alpha /= 0.0_dp)
THEN
5224 j = pw2%pw_grid%gidx(i)
5225 pw2%array(i) = pw2%array(i) + my_alpha* real(pw1%array(j), kind=dp)
5230 IF (my_alpha == 1.0_dp)
THEN
5233 j = pw2%pw_grid%gidx(i)
5234 pw2%array(j) = pw2%array(j) + real(pw1%array(i), kind=dp)
5237 ELSE IF (my_alpha /= 0.0_dp)
THEN
5240 j = pw2%pw_grid%gidx(i)
5241 pw2%array(j) = pw2%array(j) + my_alpha* real(pw1%array(i), kind=dp)
5246 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
5247 IF (ng1 >= ng2)
THEN
5248 IF (my_alpha == 1.0_dp)
THEN
5251 j = pw1%pw_grid%gidx(i)
5252 pw2%array(i) = pw2%array(i) + real(pw1%array(j), kind=dp)
5255 ELSE IF (my_alpha /= 0.0_dp)
THEN
5258 j = pw1%pw_grid%gidx(i)
5259 pw2%array(i) = pw2%array(i) + my_alpha* real(pw1%array(j), kind=dp)
5264 IF (my_alpha == 1.0_dp)
THEN
5267 j = pw1%pw_grid%gidx(i)
5268 pw2%array(j) = pw2%array(j) + real(pw1%array(i), kind=dp)
5271 ELSE IF (my_alpha /= 0.0_dp)
THEN
5274 j = pw1%pw_grid%gidx(i)
5275 pw2%array(j) = pw2%array(j) + my_alpha* real(pw1%array(i), kind=dp)
5281 cpabort(
"Grids not compatible")
5286 cpabort(
"Grids not compatible")
5290 CALL timestop(handle)
5292 END SUBROUTINE pw_axpy_c1d_r1d_rs
5303 SUBROUTINE pw_multiply_c1d_r1d_rs (pw_out, pw1, pw2, alpha)
5305 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw_out
5306 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
5307 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw2
5308 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
5310 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
5313 REAL(KIND=dp) :: my_alpha
5315 CALL timeset(routinen, handle)
5318 IF (
PRESENT(alpha)) my_alpha = alpha
5320 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
5321 cpabort(
"pw_multiply not implemented for non-identical grids!")
5323#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
5324 IF (my_alpha == 1.0_dp)
THEN
5326 pw_out%array = pw_out%array + pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5330 pw_out%array = pw_out%array + my_alpha*pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5334 IF (my_alpha == 1.0_dp)
THEN
5335 pw_out%array = pw_out%array + pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5337 pw_out%array = pw_out%array + my_alpha*pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5341 CALL timestop(handle)
5343 END SUBROUTINE pw_multiply_c1d_r1d_rs
5350 SUBROUTINE pw_multiply_with_c1d_r1d_rs (pw1, pw2)
5351 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw1
5352 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw2
5354 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
5358 CALL timeset(routinen, handle)
5360 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
5361 cpabort(
"Incompatible grids!")
5364 pw1%array = pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5367 CALL timestop(handle)
5369 END SUBROUTINE pw_multiply_with_c1d_r1d_rs
5384 FUNCTION pw_integral_ab_c1d_r1d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
5386 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
5387 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw2
5388 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
5389 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
5390 REAL(kind=dp) :: integral_value
5392 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_c1d_r1d_rs'
5394 INTEGER :: handle, loc_sumtype
5395 LOGICAL :: my_just_sum, my_local_only
5397 CALL timeset(routinen, handle)
5400 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
5402 my_local_only = .false.
5403 IF (
PRESENT(local_only)) my_local_only = local_only
5405 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
5406 cpabort(
"Grids incompatible")
5409 my_just_sum = .false.
5410 IF (
PRESENT(just_sum)) my_just_sum = just_sum
5417 integral_value = sum(real(pw1%array, kind=dp)*pw2%array)
5422 integral_value = accurate_sum(real(pw1%array, kind=dp)*pw2%array)
5428 IF (pw1%pw_grid%grid_span == halfspace)
THEN
5429 integral_value = 2.0_dp*integral_value
5430 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
5431 REAL(conjg(pw1%array(1))*pw2%array(1), kind=dp)
5434 IF (.NOT. my_just_sum)
THEN
5435 integral_value = integral_value*pw1%pw_grid%dvol
5438 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
5439 CALL pw1%pw_grid%para%group%sum(integral_value)
5441 CALL timestop(handle)
5443 END FUNCTION pw_integral_ab_c1d_r1d_rs
5457 SUBROUTINE pw_copy_c1d_r1d_gs (pw1, pw2)
5459 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
5460 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw2
5462 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
5465 INTEGER :: i, j, ng, ng1, ng2, ns
5467 CALL timeset(routinen, handle)
5469 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
5470 cpabort(
"Both grids must be either spherical or non-spherical!")
5472 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
5473 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
5474 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
5475 ng1 =
SIZE(pw1%array)
5476 ng2 =
SIZE(pw2%array)
5479 pw2%array(1:ng) = real(pw1%array(1:ng), kind=dp)
5483 pw2%array(ng + 1:ng2) = 0.0_dp
5487 cpabort(
"Copies between spherical grids require compatible grids!")
5490 ng1 =
SIZE(pw1%array)
5491 ng2 =
SIZE(pw2%array)
5492 ns = 2*max(ng1, ng2)
5494 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
5495 IF (ng1 >= ng2)
THEN
5500 j = pw2%pw_grid%gidx(i)
5501 pw2%array(i) = real(pw1%array(j), kind=dp)
5510 j = pw2%pw_grid%gidx(i)
5511 pw2%array(j) = real(pw1%array(i), kind=dp)
5515 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
5516 IF (ng1 >= ng2)
THEN
5519 j = pw1%pw_grid%gidx(i)
5520 pw2%array(i) = real(pw1%array(j), kind=dp)
5527 j = pw1%pw_grid%gidx(i)
5528 pw2%array(j) = real(pw1%array(i), kind=dp)
5533 cpabort(
"Copy not implemented!")
5540 pw2%array = real(pw1%array, kind=dp)
5544 CALL timestop(handle)
5546 END SUBROUTINE pw_copy_c1d_r1d_gs
5553 SUBROUTINE pw_copy_to_array_c1d_r1d_gs (pw, array)
5554 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
5555 REAL(KIND=dp),
DIMENSION(:),
INTENT(INOUT) :: array
5557 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
5561 CALL timeset(routinen, handle)
5564 array(:) = real(pw%array(:), kind=dp)
5567 CALL timestop(handle)
5568 END SUBROUTINE pw_copy_to_array_c1d_r1d_gs
5575 SUBROUTINE pw_copy_from_array_c1d_r1d_gs (pw, array)
5576 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
5577 REAL(KIND=dp),
DIMENSION(:),
INTENT(IN) :: array
5579 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
5583 CALL timeset(routinen, handle)
5586 pw%array = cmplx(array, 0.0_dp, kind=dp)
5589 CALL timestop(handle)
5590 END SUBROUTINE pw_copy_from_array_c1d_r1d_gs
5608 SUBROUTINE pw_axpy_c1d_r1d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
5610 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
5611 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw2
5612 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
5613 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
5615 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
5618 LOGICAL :: my_allow_noncompatible_grids
5619 REAL(KIND=dp) :: my_alpha, my_beta
5620 INTEGER :: i, j, ng, ng1, ng2
5622 CALL timeset(routinen, handle)
5625 IF (
PRESENT(alpha)) my_alpha = alpha
5628 IF (
PRESENT(beta)) my_beta = beta
5630 my_allow_noncompatible_grids = .false.
5631 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
5633 IF (my_beta /= 1.0_dp)
THEN
5634 IF (my_beta == 0.0_dp)
THEN
5638 pw2%array = pw2%array*my_beta
5643 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
5645 IF (my_alpha == 1.0_dp)
THEN
5647 pw2%array = pw2%array + real(pw1%array, kind=dp)
5649 ELSE IF (my_alpha /= 0.0_dp)
THEN
5651 pw2%array = pw2%array + my_alpha* real(pw1%array, kind=dp)
5655 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
5657 ng1 =
SIZE(pw1%array)
5658 ng2 =
SIZE(pw2%array)
5661 IF (pw1%pw_grid%spherical)
THEN
5662 IF (my_alpha == 1.0_dp)
THEN
5665 pw2%array(i) = pw2%array(i) + real(pw1%array(i), kind=dp)
5668 ELSE IF (my_alpha /= 0.0_dp)
THEN
5671 pw2%array(i) = pw2%array(i) + my_alpha* real(pw1%array(i), kind=dp)
5675 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
5676 IF (ng1 >= ng2)
THEN
5677 IF (my_alpha == 1.0_dp)
THEN
5680 j = pw2%pw_grid%gidx(i)
5681 pw2%array(i) = pw2%array(i) + real(pw1%array(j), kind=dp)
5684 ELSE IF (my_alpha /= 0.0_dp)
THEN
5687 j = pw2%pw_grid%gidx(i)
5688 pw2%array(i) = pw2%array(i) + my_alpha* real(pw1%array(j), kind=dp)
5693 IF (my_alpha == 1.0_dp)
THEN
5696 j = pw2%pw_grid%gidx(i)
5697 pw2%array(j) = pw2%array(j) + real(pw1%array(i), kind=dp)
5700 ELSE IF (my_alpha /= 0.0_dp)
THEN
5703 j = pw2%pw_grid%gidx(i)
5704 pw2%array(j) = pw2%array(j) + my_alpha* real(pw1%array(i), kind=dp)
5709 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
5710 IF (ng1 >= ng2)
THEN
5711 IF (my_alpha == 1.0_dp)
THEN
5714 j = pw1%pw_grid%gidx(i)
5715 pw2%array(i) = pw2%array(i) + real(pw1%array(j), kind=dp)
5718 ELSE IF (my_alpha /= 0.0_dp)
THEN
5721 j = pw1%pw_grid%gidx(i)
5722 pw2%array(i) = pw2%array(i) + my_alpha* real(pw1%array(j), kind=dp)
5727 IF (my_alpha == 1.0_dp)
THEN
5730 j = pw1%pw_grid%gidx(i)
5731 pw2%array(j) = pw2%array(j) + real(pw1%array(i), kind=dp)
5734 ELSE IF (my_alpha /= 0.0_dp)
THEN
5737 j = pw1%pw_grid%gidx(i)
5738 pw2%array(j) = pw2%array(j) + my_alpha* real(pw1%array(i), kind=dp)
5744 cpabort(
"Grids not compatible")
5749 cpabort(
"Grids not compatible")
5753 CALL timestop(handle)
5755 END SUBROUTINE pw_axpy_c1d_r1d_gs
5766 SUBROUTINE pw_multiply_c1d_r1d_gs (pw_out, pw1, pw2, alpha)
5768 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw_out
5769 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
5770 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
5771 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
5773 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
5776 REAL(KIND=dp) :: my_alpha
5778 CALL timeset(routinen, handle)
5781 IF (
PRESENT(alpha)) my_alpha = alpha
5783 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
5784 cpabort(
"pw_multiply not implemented for non-identical grids!")
5786#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
5787 IF (my_alpha == 1.0_dp)
THEN
5789 pw_out%array = pw_out%array + pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5793 pw_out%array = pw_out%array + my_alpha*pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5797 IF (my_alpha == 1.0_dp)
THEN
5798 pw_out%array = pw_out%array + pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5800 pw_out%array = pw_out%array + my_alpha*pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5804 CALL timestop(handle)
5806 END SUBROUTINE pw_multiply_c1d_r1d_gs
5813 SUBROUTINE pw_multiply_with_c1d_r1d_gs (pw1, pw2)
5814 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw1
5815 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
5817 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
5821 CALL timeset(routinen, handle)
5823 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
5824 cpabort(
"Incompatible grids!")
5827 pw1%array = pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5830 CALL timestop(handle)
5832 END SUBROUTINE pw_multiply_with_c1d_r1d_gs
5847 FUNCTION pw_integral_ab_c1d_r1d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
5849 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
5850 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
5851 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
5852 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
5853 REAL(kind=dp) :: integral_value
5855 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_c1d_r1d_gs'
5857 INTEGER :: handle, loc_sumtype
5858 LOGICAL :: my_just_sum, my_local_only
5860 CALL timeset(routinen, handle)
5863 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
5865 my_local_only = .false.
5866 IF (
PRESENT(local_only)) my_local_only = local_only
5868 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
5869 cpabort(
"Grids incompatible")
5872 my_just_sum = .false.
5873 IF (
PRESENT(just_sum)) my_just_sum = just_sum
5880 integral_value = sum(real(pw1%array, kind=dp)*pw2%array)
5885 integral_value = accurate_sum(real(pw1%array, kind=dp)*pw2%array)
5891 IF (pw1%pw_grid%grid_span == halfspace)
THEN
5892 integral_value = 2.0_dp*integral_value
5893 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
5894 REAL(conjg(pw1%array(1))*pw2%array(1), kind=dp)
5897 IF (.NOT. my_just_sum)
THEN
5898 integral_value = integral_value*pw1%pw_grid%vol
5901 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
5902 CALL pw1%pw_grid%para%group%sum(integral_value)
5904 CALL timestop(handle)
5906 END FUNCTION pw_integral_ab_c1d_r1d_gs
5914 FUNCTION pw_integral_a2b_c1d_r1d (pw1, pw2)
RESULT(integral_value)
5916 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
5917 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
5918 REAL(kind=dp) :: integral_value
5920 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_a2b'
5924 CALL timeset(routinen, handle)
5926 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
5927 cpabort(
"Grids incompatible")
5930 integral_value = accurate_sum(real(conjg(pw1%array), kind=dp)*pw2%array*pw1%pw_grid%gsq)
5931 IF (pw1%pw_grid%grid_span == halfspace)
THEN
5932 integral_value = 2.0_dp*integral_value
5935 integral_value = integral_value*pw1%pw_grid%vol
5937 IF (pw1%pw_grid%para%mode == pw_mode_distributed) &
5938 CALL pw1%pw_grid%para%group%sum(integral_value)
5939 CALL timestop(handle)
5941 END FUNCTION pw_integral_a2b_c1d_r1d
5955 SUBROUTINE pw_copy_c1d_c1d_rs (pw1, pw2)
5957 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
5958 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw2
5960 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
5963 INTEGER :: i, j, ng, ng1, ng2, ns
5965 CALL timeset(routinen, handle)
5967 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
5968 cpabort(
"Both grids must be either spherical or non-spherical!")
5969 IF (pw1%pw_grid%spherical) &
5970 cpabort(
"Spherical grids only exist in reciprocal space!")
5972 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
5973 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
5974 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
5975 ng1 =
SIZE(pw1%array)
5976 ng2 =
SIZE(pw2%array)
5979 pw2%array(1:ng) = pw1%array(1:ng)
5983 pw2%array(ng + 1:ng2) = cmplx(0.0_dp, 0.0_dp, kind=dp)
5987 cpabort(
"Copies between spherical grids require compatible grids!")
5990 ng1 =
SIZE(pw1%array)
5991 ng2 =
SIZE(pw2%array)
5992 ns = 2*max(ng1, ng2)
5994 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
5995 IF (ng1 >= ng2)
THEN
6000 j = pw2%pw_grid%gidx(i)
6001 pw2%array(i) = pw1%array(j)
6010 j = pw2%pw_grid%gidx(i)
6011 pw2%array(j) = pw1%array(i)
6015 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
6016 IF (ng1 >= ng2)
THEN
6019 j = pw1%pw_grid%gidx(i)
6020 pw2%array(i) = pw1%array(j)
6027 j = pw1%pw_grid%gidx(i)
6028 pw2%array(j) = pw1%array(i)
6033 cpabort(
"Copy not implemented!")
6040 pw2%array = pw1%array
6044 CALL timestop(handle)
6046 END SUBROUTINE pw_copy_c1d_c1d_rs
6053 SUBROUTINE pw_copy_to_array_c1d_c1d_rs (pw, array)
6054 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw
6055 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(INOUT) :: array
6057 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
6061 CALL timeset(routinen, handle)
6064 array(:) = pw%array(:)
6067 CALL timestop(handle)
6068 END SUBROUTINE pw_copy_to_array_c1d_c1d_rs
6075 SUBROUTINE pw_copy_from_array_c1d_c1d_rs (pw, array)
6076 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw
6077 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(IN) :: array
6079 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
6083 CALL timeset(routinen, handle)
6089 CALL timestop(handle)
6090 END SUBROUTINE pw_copy_from_array_c1d_c1d_rs
6108 SUBROUTINE pw_axpy_c1d_c1d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
6110 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
6111 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw2
6112 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
6113 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
6115 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
6118 LOGICAL :: my_allow_noncompatible_grids
6119 REAL(KIND=dp) :: my_alpha, my_beta
6120 INTEGER :: i, j, ng, ng1, ng2
6122 CALL timeset(routinen, handle)
6125 IF (
PRESENT(alpha)) my_alpha = alpha
6128 IF (
PRESENT(beta)) my_beta = beta
6130 my_allow_noncompatible_grids = .false.
6131 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
6133 IF (my_beta /= 1.0_dp)
THEN
6134 IF (my_beta == 0.0_dp)
THEN
6138 pw2%array = pw2%array*my_beta
6143 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
6145 IF (my_alpha == 1.0_dp)
THEN
6147 pw2%array = pw2%array + pw1%array
6149 ELSE IF (my_alpha /= 0.0_dp)
THEN
6151 pw2%array = pw2%array + my_alpha* pw1%array
6155 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
6157 ng1 =
SIZE(pw1%array)
6158 ng2 =
SIZE(pw2%array)
6161 IF (pw1%pw_grid%spherical)
THEN
6162 IF (my_alpha == 1.0_dp)
THEN
6165 pw2%array(i) = pw2%array(i) + pw1%array(i)
6168 ELSE IF (my_alpha /= 0.0_dp)
THEN
6171 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(i)
6175 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
6176 IF (ng1 >= ng2)
THEN
6177 IF (my_alpha == 1.0_dp)
THEN
6180 j = pw2%pw_grid%gidx(i)
6181 pw2%array(i) = pw2%array(i) + pw1%array(j)
6184 ELSE IF (my_alpha /= 0.0_dp)
THEN
6187 j = pw2%pw_grid%gidx(i)
6188 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
6193 IF (my_alpha == 1.0_dp)
THEN
6196 j = pw2%pw_grid%gidx(i)
6197 pw2%array(j) = pw2%array(j) + pw1%array(i)
6200 ELSE IF (my_alpha /= 0.0_dp)
THEN
6203 j = pw2%pw_grid%gidx(i)
6204 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
6209 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
6210 IF (ng1 >= ng2)
THEN
6211 IF (my_alpha == 1.0_dp)
THEN
6214 j = pw1%pw_grid%gidx(i)
6215 pw2%array(i) = pw2%array(i) + pw1%array(j)
6218 ELSE IF (my_alpha /= 0.0_dp)
THEN
6221 j = pw1%pw_grid%gidx(i)
6222 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
6227 IF (my_alpha == 1.0_dp)
THEN
6230 j = pw1%pw_grid%gidx(i)
6231 pw2%array(j) = pw2%array(j) + pw1%array(i)
6234 ELSE IF (my_alpha /= 0.0_dp)
THEN
6237 j = pw1%pw_grid%gidx(i)
6238 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
6244 cpabort(
"Grids not compatible")
6249 cpabort(
"Grids not compatible")
6253 CALL timestop(handle)
6255 END SUBROUTINE pw_axpy_c1d_c1d_rs
6266 SUBROUTINE pw_multiply_c1d_c1d_rs (pw_out, pw1, pw2, alpha)
6268 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw_out
6269 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
6270 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw2
6271 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
6273 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
6276 REAL(KIND=dp) :: my_alpha
6278 CALL timeset(routinen, handle)
6281 IF (
PRESENT(alpha)) my_alpha = alpha
6283 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
6284 cpabort(
"pw_multiply not implemented for non-identical grids!")
6286#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
6287 IF (my_alpha == 1.0_dp)
THEN
6289 pw_out%array = pw_out%array + pw1%array* pw2%array
6293 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
6297 IF (my_alpha == 1.0_dp)
THEN
6298 pw_out%array = pw_out%array + pw1%array* pw2%array
6300 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
6304 CALL timestop(handle)
6306 END SUBROUTINE pw_multiply_c1d_c1d_rs
6313 SUBROUTINE pw_multiply_with_c1d_c1d_rs (pw1, pw2)
6314 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw1
6315 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw2
6317 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
6321 CALL timeset(routinen, handle)
6323 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
6324 cpabort(
"Incompatible grids!")
6327 pw1%array = pw1%array* pw2%array
6330 CALL timestop(handle)
6332 END SUBROUTINE pw_multiply_with_c1d_c1d_rs
6347 FUNCTION pw_integral_ab_c1d_c1d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
6349 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
6350 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw2
6351 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
6352 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
6353 REAL(kind=dp) :: integral_value
6355 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_c1d_c1d_rs'
6357 INTEGER :: handle, loc_sumtype
6358 LOGICAL :: my_just_sum, my_local_only
6360 CALL timeset(routinen, handle)
6363 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
6365 my_local_only = .false.
6366 IF (
PRESENT(local_only)) my_local_only = local_only
6368 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
6369 cpabort(
"Grids incompatible")
6372 my_just_sum = .false.
6373 IF (
PRESENT(just_sum)) my_just_sum = just_sum
6380 integral_value = sum(real(conjg(pw1%array) &
6381 *pw2%array, kind=dp))
6386 integral_value = accurate_sum(real(conjg(pw1%array)*pw2%array, kind=dp))
6392 IF (pw1%pw_grid%grid_span == halfspace)
THEN
6393 integral_value = 2.0_dp*integral_value
6394 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
6395 REAL(conjg(pw1%array(1))*pw2%array(1), kind=dp)
6398 IF (.NOT. my_just_sum)
THEN
6399 integral_value = integral_value*pw1%pw_grid%dvol
6402 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
6403 CALL pw1%pw_grid%para%group%sum(integral_value)
6405 CALL timestop(handle)
6407 END FUNCTION pw_integral_ab_c1d_c1d_rs
6421 SUBROUTINE pw_copy_c1d_c1d_gs (pw1, pw2)
6423 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
6424 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
6426 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
6429 INTEGER :: i, j, ng, ng1, ng2, ns
6431 CALL timeset(routinen, handle)
6433 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
6434 cpabort(
"Both grids must be either spherical or non-spherical!")
6436 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
6437 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
6438 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
6439 ng1 =
SIZE(pw1%array)
6440 ng2 =
SIZE(pw2%array)
6443 pw2%array(1:ng) = pw1%array(1:ng)
6447 pw2%array(ng + 1:ng2) = cmplx(0.0_dp, 0.0_dp, kind=dp)
6451 cpabort(
"Copies between spherical grids require compatible grids!")
6454 ng1 =
SIZE(pw1%array)
6455 ng2 =
SIZE(pw2%array)
6456 ns = 2*max(ng1, ng2)
6458 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
6459 IF (ng1 >= ng2)
THEN
6464 j = pw2%pw_grid%gidx(i)
6465 pw2%array(i) = pw1%array(j)
6474 j = pw2%pw_grid%gidx(i)
6475 pw2%array(j) = pw1%array(i)
6479 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
6480 IF (ng1 >= ng2)
THEN
6483 j = pw1%pw_grid%gidx(i)
6484 pw2%array(i) = pw1%array(j)
6491 j = pw1%pw_grid%gidx(i)
6492 pw2%array(j) = pw1%array(i)
6497 CALL pw_copy_match(pw1, pw2)
6504 pw2%array = pw1%array
6508 CALL timestop(handle)
6510 END SUBROUTINE pw_copy_c1d_c1d_gs
6517 SUBROUTINE pw_copy_to_array_c1d_c1d_gs (pw, array)
6518 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
6519 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(INOUT) :: array
6521 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
6525 CALL timeset(routinen, handle)
6528 array(:) = pw%array(:)
6531 CALL timestop(handle)
6532 END SUBROUTINE pw_copy_to_array_c1d_c1d_gs
6539 SUBROUTINE pw_copy_from_array_c1d_c1d_gs (pw, array)
6540 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
6541 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(IN) :: array
6543 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
6547 CALL timeset(routinen, handle)
6553 CALL timestop(handle)
6554 END SUBROUTINE pw_copy_from_array_c1d_c1d_gs
6572 SUBROUTINE pw_axpy_c1d_c1d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
6574 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
6575 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
6576 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
6577 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
6579 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
6582 LOGICAL :: my_allow_noncompatible_grids
6583 REAL(KIND=dp) :: my_alpha, my_beta
6584 INTEGER :: i, j, ng, ng1, ng2
6586 CALL timeset(routinen, handle)
6589 IF (
PRESENT(alpha)) my_alpha = alpha
6592 IF (
PRESENT(beta)) my_beta = beta
6594 my_allow_noncompatible_grids = .false.
6595 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
6597 IF (my_beta /= 1.0_dp)
THEN
6598 IF (my_beta == 0.0_dp)
THEN
6602 pw2%array = pw2%array*my_beta
6607 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
6609 IF (my_alpha == 1.0_dp)
THEN
6611 pw2%array = pw2%array + pw1%array
6613 ELSE IF (my_alpha /= 0.0_dp)
THEN
6615 pw2%array = pw2%array + my_alpha* pw1%array
6619 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
6621 ng1 =
SIZE(pw1%array)
6622 ng2 =
SIZE(pw2%array)
6625 IF (pw1%pw_grid%spherical)
THEN
6626 IF (my_alpha == 1.0_dp)
THEN
6629 pw2%array(i) = pw2%array(i) + pw1%array(i)
6632 ELSE IF (my_alpha /= 0.0_dp)
THEN
6635 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(i)
6639 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
6640 IF (ng1 >= ng2)
THEN
6641 IF (my_alpha == 1.0_dp)
THEN
6644 j = pw2%pw_grid%gidx(i)
6645 pw2%array(i) = pw2%array(i) + pw1%array(j)
6648 ELSE IF (my_alpha /= 0.0_dp)
THEN
6651 j = pw2%pw_grid%gidx(i)
6652 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
6657 IF (my_alpha == 1.0_dp)
THEN
6660 j = pw2%pw_grid%gidx(i)
6661 pw2%array(j) = pw2%array(j) + pw1%array(i)
6664 ELSE IF (my_alpha /= 0.0_dp)
THEN
6667 j = pw2%pw_grid%gidx(i)
6668 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
6673 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
6674 IF (ng1 >= ng2)
THEN
6675 IF (my_alpha == 1.0_dp)
THEN
6678 j = pw1%pw_grid%gidx(i)
6679 pw2%array(i) = pw2%array(i) + pw1%array(j)
6682 ELSE IF (my_alpha /= 0.0_dp)
THEN
6685 j = pw1%pw_grid%gidx(i)
6686 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
6691 IF (my_alpha == 1.0_dp)
THEN
6694 j = pw1%pw_grid%gidx(i)
6695 pw2%array(j) = pw2%array(j) + pw1%array(i)
6698 ELSE IF (my_alpha /= 0.0_dp)
THEN
6701 j = pw1%pw_grid%gidx(i)
6702 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
6708 cpabort(
"Grids not compatible")
6713 cpabort(
"Grids not compatible")
6717 CALL timestop(handle)
6719 END SUBROUTINE pw_axpy_c1d_c1d_gs
6730 SUBROUTINE pw_multiply_c1d_c1d_gs (pw_out, pw1, pw2, alpha)
6732 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw_out
6733 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
6734 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
6735 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
6737 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
6740 REAL(KIND=dp) :: my_alpha
6742 CALL timeset(routinen, handle)
6745 IF (
PRESENT(alpha)) my_alpha = alpha
6747 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
6748 cpabort(
"pw_multiply not implemented for non-identical grids!")
6750#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
6751 IF (my_alpha == 1.0_dp)
THEN
6753 pw_out%array = pw_out%array + pw1%array* pw2%array
6757 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
6761 IF (my_alpha == 1.0_dp)
THEN
6762 pw_out%array = pw_out%array + pw1%array* pw2%array
6764 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
6768 CALL timestop(handle)
6770 END SUBROUTINE pw_multiply_c1d_c1d_gs
6777 SUBROUTINE pw_multiply_with_c1d_c1d_gs (pw1, pw2)
6778 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw1
6779 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
6781 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
6785 CALL timeset(routinen, handle)
6787 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
6788 cpabort(
"Incompatible grids!")
6791 pw1%array = pw1%array* pw2%array
6794 CALL timestop(handle)
6796 END SUBROUTINE pw_multiply_with_c1d_c1d_gs
6811 FUNCTION pw_integral_ab_c1d_c1d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
6813 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
6814 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
6815 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
6816 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
6817 REAL(kind=dp) :: integral_value
6819 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_c1d_c1d_gs'
6821 INTEGER :: handle, loc_sumtype
6822 LOGICAL :: my_just_sum, my_local_only
6824 CALL timeset(routinen, handle)
6827 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
6829 my_local_only = .false.
6830 IF (
PRESENT(local_only)) my_local_only = local_only
6832 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
6833 cpabort(
"Grids incompatible")
6836 my_just_sum = .false.
6837 IF (
PRESENT(just_sum)) my_just_sum = just_sum
6844 integral_value = sum(real(conjg(pw1%array) &
6845 *pw2%array, kind=dp))
6850 integral_value = accurate_sum(real(conjg(pw1%array)*pw2%array, kind=dp))
6856 IF (pw1%pw_grid%grid_span == halfspace)
THEN
6857 integral_value = 2.0_dp*integral_value
6858 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
6859 REAL(conjg(pw1%array(1))*pw2%array(1), kind=dp)
6862 IF (.NOT. my_just_sum)
THEN
6863 integral_value = integral_value*pw1%pw_grid%vol
6866 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
6867 CALL pw1%pw_grid%para%group%sum(integral_value)
6869 CALL timestop(handle)
6871 END FUNCTION pw_integral_ab_c1d_c1d_gs
6879 FUNCTION pw_integral_a2b_c1d_c1d (pw1, pw2)
RESULT(integral_value)
6881 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
6882 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
6883 REAL(kind=dp) :: integral_value
6885 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_a2b'
6889 CALL timeset(routinen, handle)
6891 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
6892 cpabort(
"Grids incompatible")
6895 integral_value = accurate_sum(real(conjg(pw1%array)*pw2%array, kind=dp)*pw1%pw_grid%gsq)
6896 IF (pw1%pw_grid%grid_span == halfspace)
THEN
6897 integral_value = 2.0_dp*integral_value
6900 integral_value = integral_value*pw1%pw_grid%vol
6902 IF (pw1%pw_grid%para%mode == pw_mode_distributed) &
6903 CALL pw1%pw_grid%para%group%sum(integral_value)
6904 CALL timestop(handle)
6906 END FUNCTION pw_integral_a2b_c1d_c1d
6920 SUBROUTINE pw_copy_c3d_r3d_rs (pw1, pw2)
6922 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
6923 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw2
6925 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
6929 CALL timeset(routinen, handle)
6931 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
6932 cpabort(
"Both grids must be either spherical or non-spherical!")
6933 IF (pw1%pw_grid%spherical) &
6934 cpabort(
"Spherical grids only exist in reciprocal space!")
6936 IF (any(shape(pw2%array) /= shape(pw1%array))) &
6937 cpabort(
"3D grids must be compatible!")
6938 IF (pw1%pw_grid%spherical) &
6939 cpabort(
"3D grids must not be spherical!")
6941 pw2%array(:, :, :) = real(pw1%array(:, :, :), kind=dp)
6944 CALL timestop(handle)
6946 END SUBROUTINE pw_copy_c3d_r3d_rs
6953 SUBROUTINE pw_copy_to_array_c3d_r3d_rs (pw, array)
6954 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw
6955 REAL(KIND=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
6957 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
6961 CALL timeset(routinen, handle)
6964 array(:, :, :) = real(pw%array(:, :, :), kind=dp)
6967 CALL timestop(handle)
6968 END SUBROUTINE pw_copy_to_array_c3d_r3d_rs
6975 SUBROUTINE pw_copy_from_array_c3d_r3d_rs (pw, array)
6976 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw
6977 REAL(KIND=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
6979 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
6983 CALL timeset(routinen, handle)
6986 pw%array = cmplx(array, 0.0_dp, kind=dp)
6989 CALL timestop(handle)
6990 END SUBROUTINE pw_copy_from_array_c3d_r3d_rs
7008 SUBROUTINE pw_axpy_c3d_r3d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
7010 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
7011 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw2
7012 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
7013 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
7015 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
7018 LOGICAL :: my_allow_noncompatible_grids
7019 REAL(KIND=dp) :: my_alpha, my_beta
7021 CALL timeset(routinen, handle)
7024 IF (
PRESENT(alpha)) my_alpha = alpha
7027 IF (
PRESENT(beta)) my_beta = beta
7029 my_allow_noncompatible_grids = .false.
7030 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
7032 IF (my_beta /= 1.0_dp)
THEN
7033 IF (my_beta == 0.0_dp)
THEN
7037 pw2%array = pw2%array*my_beta
7042 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
7043 IF (my_alpha == 1.0_dp)
THEN
7045 pw2%array = pw2%array + real(pw1%array, kind=dp)
7047 ELSE IF (my_alpha /= 0.0_dp)
THEN
7049 pw2%array = pw2%array + my_alpha* real(pw1%array, kind=dp)
7053 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
7055 IF (any(shape(pw1%array) /= shape(pw2%array))) &
7056 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
7058 IF (my_alpha == 1.0_dp)
THEN
7060 pw2%array = pw2%array + real(pw1%array, kind=dp)
7064 pw2%array = pw2%array + my_alpha* real(pw1%array, kind=dp)
7070 cpabort(
"Grids not compatible")
7074 CALL timestop(handle)
7076 END SUBROUTINE pw_axpy_c3d_r3d_rs
7087 SUBROUTINE pw_multiply_c3d_r3d_rs (pw_out, pw1, pw2, alpha)
7089 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw_out
7090 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
7091 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw2
7092 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
7094 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
7097 REAL(KIND=dp) :: my_alpha
7099 CALL timeset(routinen, handle)
7102 IF (
PRESENT(alpha)) my_alpha = alpha
7104 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
7105 cpabort(
"pw_multiply not implemented for non-identical grids!")
7107#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
7108 IF (my_alpha == 1.0_dp)
THEN
7110 pw_out%array = pw_out%array + pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7114 pw_out%array = pw_out%array + my_alpha*pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7118 IF (my_alpha == 1.0_dp)
THEN
7119 pw_out%array = pw_out%array + pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7121 pw_out%array = pw_out%array + my_alpha*pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7125 CALL timestop(handle)
7127 END SUBROUTINE pw_multiply_c3d_r3d_rs
7134 SUBROUTINE pw_multiply_with_c3d_r3d_rs (pw1, pw2)
7135 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw1
7136 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw2
7138 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
7142 CALL timeset(routinen, handle)
7144 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
7145 cpabort(
"Incompatible grids!")
7148 pw1%array = pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7151 CALL timestop(handle)
7153 END SUBROUTINE pw_multiply_with_c3d_r3d_rs
7168 FUNCTION pw_integral_ab_c3d_r3d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
7170 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
7171 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw2
7172 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
7173 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
7174 REAL(kind=dp) :: integral_value
7176 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_c3d_r3d_rs'
7178 INTEGER :: handle, loc_sumtype
7179 LOGICAL :: my_just_sum, my_local_only
7181 CALL timeset(routinen, handle)
7184 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
7186 my_local_only = .false.
7187 IF (
PRESENT(local_only)) my_local_only = local_only
7189 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
7190 cpabort(
"Grids incompatible")
7193 my_just_sum = .false.
7194 IF (
PRESENT(just_sum)) my_just_sum = just_sum
7201 integral_value = sum(real(pw1%array, kind=dp)*pw2%array)
7206 integral_value = accurate_sum(real(pw1%array, kind=dp)*pw2%array)
7211 IF (.NOT. my_just_sum)
THEN
7212 integral_value = integral_value*pw1%pw_grid%dvol
7215 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
7216 CALL pw1%pw_grid%para%group%sum(integral_value)
7218 CALL timestop(handle)
7220 END FUNCTION pw_integral_ab_c3d_r3d_rs
7234 SUBROUTINE pw_copy_c3d_r3d_gs (pw1, pw2)
7236 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
7237 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw2
7239 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
7243 CALL timeset(routinen, handle)
7245 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
7246 cpabort(
"Both grids must be either spherical or non-spherical!")
7248 IF (any(shape(pw2%array) /= shape(pw1%array))) &
7249 cpabort(
"3D grids must be compatible!")
7250 IF (pw1%pw_grid%spherical) &
7251 cpabort(
"3D grids must not be spherical!")
7253 pw2%array(:, :, :) = real(pw1%array(:, :, :), kind=dp)
7256 CALL timestop(handle)
7258 END SUBROUTINE pw_copy_c3d_r3d_gs
7265 SUBROUTINE pw_copy_to_array_c3d_r3d_gs (pw, array)
7266 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw
7267 REAL(KIND=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
7269 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
7273 CALL timeset(routinen, handle)
7276 array(:, :, :) = real(pw%array(:, :, :), kind=dp)
7279 CALL timestop(handle)
7280 END SUBROUTINE pw_copy_to_array_c3d_r3d_gs
7287 SUBROUTINE pw_copy_from_array_c3d_r3d_gs (pw, array)
7288 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw
7289 REAL(KIND=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
7291 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
7295 CALL timeset(routinen, handle)
7298 pw%array = cmplx(array, 0.0_dp, kind=dp)
7301 CALL timestop(handle)
7302 END SUBROUTINE pw_copy_from_array_c3d_r3d_gs
7320 SUBROUTINE pw_axpy_c3d_r3d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
7322 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
7323 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw2
7324 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
7325 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
7327 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
7330 LOGICAL :: my_allow_noncompatible_grids
7331 REAL(KIND=dp) :: my_alpha, my_beta
7333 CALL timeset(routinen, handle)
7336 IF (
PRESENT(alpha)) my_alpha = alpha
7339 IF (
PRESENT(beta)) my_beta = beta
7341 my_allow_noncompatible_grids = .false.
7342 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
7344 IF (my_beta /= 1.0_dp)
THEN
7345 IF (my_beta == 0.0_dp)
THEN
7349 pw2%array = pw2%array*my_beta
7354 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
7355 IF (my_alpha == 1.0_dp)
THEN
7357 pw2%array = pw2%array + real(pw1%array, kind=dp)
7359 ELSE IF (my_alpha /= 0.0_dp)
THEN
7361 pw2%array = pw2%array + my_alpha* real(pw1%array, kind=dp)
7365 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
7367 IF (any(shape(pw1%array) /= shape(pw2%array))) &
7368 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
7370 IF (my_alpha == 1.0_dp)
THEN
7372 pw2%array = pw2%array + real(pw1%array, kind=dp)
7376 pw2%array = pw2%array + my_alpha* real(pw1%array, kind=dp)
7382 cpabort(
"Grids not compatible")
7386 CALL timestop(handle)
7388 END SUBROUTINE pw_axpy_c3d_r3d_gs
7399 SUBROUTINE pw_multiply_c3d_r3d_gs (pw_out, pw1, pw2, alpha)
7401 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw_out
7402 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
7403 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw2
7404 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
7406 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
7409 REAL(KIND=dp) :: my_alpha
7411 CALL timeset(routinen, handle)
7414 IF (
PRESENT(alpha)) my_alpha = alpha
7416 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
7417 cpabort(
"pw_multiply not implemented for non-identical grids!")
7419#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
7420 IF (my_alpha == 1.0_dp)
THEN
7422 pw_out%array = pw_out%array + pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7426 pw_out%array = pw_out%array + my_alpha*pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7430 IF (my_alpha == 1.0_dp)
THEN
7431 pw_out%array = pw_out%array + pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7433 pw_out%array = pw_out%array + my_alpha*pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7437 CALL timestop(handle)
7439 END SUBROUTINE pw_multiply_c3d_r3d_gs
7446 SUBROUTINE pw_multiply_with_c3d_r3d_gs (pw1, pw2)
7447 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw1
7448 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw2
7450 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
7454 CALL timeset(routinen, handle)
7456 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
7457 cpabort(
"Incompatible grids!")
7460 pw1%array = pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7463 CALL timestop(handle)
7465 END SUBROUTINE pw_multiply_with_c3d_r3d_gs
7480 FUNCTION pw_integral_ab_c3d_r3d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
7482 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
7483 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw2
7484 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
7485 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
7486 REAL(kind=dp) :: integral_value
7488 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_c3d_r3d_gs'
7490 INTEGER :: handle, loc_sumtype
7491 LOGICAL :: my_just_sum, my_local_only
7493 CALL timeset(routinen, handle)
7496 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
7498 my_local_only = .false.
7499 IF (
PRESENT(local_only)) my_local_only = local_only
7501 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
7502 cpabort(
"Grids incompatible")
7505 my_just_sum = .false.
7506 IF (
PRESENT(just_sum)) my_just_sum = just_sum
7513 integral_value = sum(real(pw1%array, kind=dp)*pw2%array)
7518 integral_value = accurate_sum(real(pw1%array, kind=dp)*pw2%array)
7523 IF (.NOT. my_just_sum)
THEN
7524 integral_value = integral_value*pw1%pw_grid%vol
7527 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
7528 CALL pw1%pw_grid%para%group%sum(integral_value)
7530 CALL timestop(handle)
7532 END FUNCTION pw_integral_ab_c3d_r3d_gs
7547 SUBROUTINE pw_copy_c3d_c3d_rs (pw1, pw2)
7549 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
7550 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw2
7552 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
7556 CALL timeset(routinen, handle)
7558 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
7559 cpabort(
"Both grids must be either spherical or non-spherical!")
7560 IF (pw1%pw_grid%spherical) &
7561 cpabort(
"Spherical grids only exist in reciprocal space!")
7563 IF (any(shape(pw2%array) /= shape(pw1%array))) &
7564 cpabort(
"3D grids must be compatible!")
7565 IF (pw1%pw_grid%spherical) &
7566 cpabort(
"3D grids must not be spherical!")
7568 pw2%array(:, :, :) = pw1%array(:, :, :)
7571 CALL timestop(handle)
7573 END SUBROUTINE pw_copy_c3d_c3d_rs
7580 SUBROUTINE pw_copy_to_array_c3d_c3d_rs (pw, array)
7581 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw
7582 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
7584 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
7588 CALL timeset(routinen, handle)
7591 array(:, :, :) = pw%array(:, :, :)
7594 CALL timestop(handle)
7595 END SUBROUTINE pw_copy_to_array_c3d_c3d_rs
7602 SUBROUTINE pw_copy_from_array_c3d_c3d_rs (pw, array)
7603 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw
7604 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
7606 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
7610 CALL timeset(routinen, handle)
7616 CALL timestop(handle)
7617 END SUBROUTINE pw_copy_from_array_c3d_c3d_rs
7635 SUBROUTINE pw_axpy_c3d_c3d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
7637 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
7638 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw2
7639 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
7640 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
7642 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
7645 LOGICAL :: my_allow_noncompatible_grids
7646 REAL(KIND=dp) :: my_alpha, my_beta
7648 CALL timeset(routinen, handle)
7651 IF (
PRESENT(alpha)) my_alpha = alpha
7654 IF (
PRESENT(beta)) my_beta = beta
7656 my_allow_noncompatible_grids = .false.
7657 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
7659 IF (my_beta /= 1.0_dp)
THEN
7660 IF (my_beta == 0.0_dp)
THEN
7664 pw2%array = pw2%array*my_beta
7669 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
7670 IF (my_alpha == 1.0_dp)
THEN
7672 pw2%array = pw2%array + pw1%array
7674 ELSE IF (my_alpha /= 0.0_dp)
THEN
7676 pw2%array = pw2%array + my_alpha* pw1%array
7680 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
7682 IF (any(shape(pw1%array) /= shape(pw2%array))) &
7683 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
7685 IF (my_alpha == 1.0_dp)
THEN
7687 pw2%array = pw2%array + pw1%array
7691 pw2%array = pw2%array + my_alpha* pw1%array
7697 cpabort(
"Grids not compatible")
7701 CALL timestop(handle)
7703 END SUBROUTINE pw_axpy_c3d_c3d_rs
7714 SUBROUTINE pw_multiply_c3d_c3d_rs (pw_out, pw1, pw2, alpha)
7716 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw_out
7717 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
7718 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw2
7719 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
7721 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
7724 REAL(KIND=dp) :: my_alpha
7726 CALL timeset(routinen, handle)
7729 IF (
PRESENT(alpha)) my_alpha = alpha
7731 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
7732 cpabort(
"pw_multiply not implemented for non-identical grids!")
7734#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
7735 IF (my_alpha == 1.0_dp)
THEN
7737 pw_out%array = pw_out%array + pw1%array* pw2%array
7741 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
7745 IF (my_alpha == 1.0_dp)
THEN
7746 pw_out%array = pw_out%array + pw1%array* pw2%array
7748 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
7752 CALL timestop(handle)
7754 END SUBROUTINE pw_multiply_c3d_c3d_rs
7761 SUBROUTINE pw_multiply_with_c3d_c3d_rs (pw1, pw2)
7762 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw1
7763 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw2
7765 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
7769 CALL timeset(routinen, handle)
7771 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
7772 cpabort(
"Incompatible grids!")
7775 pw1%array = pw1%array* pw2%array
7778 CALL timestop(handle)
7780 END SUBROUTINE pw_multiply_with_c3d_c3d_rs
7795 FUNCTION pw_integral_ab_c3d_c3d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
7797 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
7798 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw2
7799 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
7800 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
7801 REAL(kind=dp) :: integral_value
7803 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_c3d_c3d_rs'
7805 INTEGER :: handle, loc_sumtype
7806 LOGICAL :: my_just_sum, my_local_only
7808 CALL timeset(routinen, handle)
7811 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
7813 my_local_only = .false.
7814 IF (
PRESENT(local_only)) my_local_only = local_only
7816 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
7817 cpabort(
"Grids incompatible")
7820 my_just_sum = .false.
7821 IF (
PRESENT(just_sum)) my_just_sum = just_sum
7828 integral_value = sum(real(conjg(pw1%array) &
7829 *pw2%array, kind=dp))
7834 integral_value = accurate_sum(real(conjg(pw1%array)*pw2%array, kind=dp))
7839 IF (.NOT. my_just_sum)
THEN
7840 integral_value = integral_value*pw1%pw_grid%dvol
7843 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
7844 CALL pw1%pw_grid%para%group%sum(integral_value)
7846 CALL timestop(handle)
7848 END FUNCTION pw_integral_ab_c3d_c3d_rs
7862 SUBROUTINE pw_copy_c3d_c3d_gs (pw1, pw2)
7864 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
7865 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
7867 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy'
7871 CALL timeset(routinen, handle)
7873 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
7874 cpabort(
"Both grids must be either spherical or non-spherical!")
7876 IF (any(shape(pw2%array) /= shape(pw1%array))) &
7877 cpabort(
"3D grids must be compatible!")
7878 IF (pw1%pw_grid%spherical) &
7879 cpabort(
"3D grids must not be spherical!")
7881 pw2%array(:, :, :) = pw1%array(:, :, :)
7884 CALL timestop(handle)
7886 END SUBROUTINE pw_copy_c3d_c3d_gs
7893 SUBROUTINE pw_copy_to_array_c3d_c3d_gs (pw, array)
7894 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw
7895 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
7897 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_to_array'
7901 CALL timeset(routinen, handle)
7904 array(:, :, :) = pw%array(:, :, :)
7907 CALL timestop(handle)
7908 END SUBROUTINE pw_copy_to_array_c3d_c3d_gs
7915 SUBROUTINE pw_copy_from_array_c3d_c3d_gs (pw, array)
7916 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw
7917 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
7919 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_copy_from_array'
7923 CALL timeset(routinen, handle)
7929 CALL timestop(handle)
7930 END SUBROUTINE pw_copy_from_array_c3d_c3d_gs
7948 SUBROUTINE pw_axpy_c3d_c3d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
7950 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
7951 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
7952 REAL(KIND=dp),
INTENT(in),
OPTIONAL :: alpha, beta
7953 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
7955 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_axpy'
7958 LOGICAL :: my_allow_noncompatible_grids
7959 REAL(KIND=dp) :: my_alpha, my_beta
7961 CALL timeset(routinen, handle)
7964 IF (
PRESENT(alpha)) my_alpha = alpha
7967 IF (
PRESENT(beta)) my_beta = beta
7969 my_allow_noncompatible_grids = .false.
7970 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
7972 IF (my_beta /= 1.0_dp)
THEN
7973 IF (my_beta == 0.0_dp)
THEN
7977 pw2%array = pw2%array*my_beta
7982 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
7983 IF (my_alpha == 1.0_dp)
THEN
7985 pw2%array = pw2%array + pw1%array
7987 ELSE IF (my_alpha /= 0.0_dp)
THEN
7989 pw2%array = pw2%array + my_alpha* pw1%array
7993 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
7995 IF (any(shape(pw1%array) /= shape(pw2%array))) &
7996 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
7998 IF (my_alpha == 1.0_dp)
THEN
8000 pw2%array = pw2%array + pw1%array
8004 pw2%array = pw2%array + my_alpha* pw1%array
8010 cpabort(
"Grids not compatible")
8014 CALL timestop(handle)
8016 END SUBROUTINE pw_axpy_c3d_c3d_gs
8027 SUBROUTINE pw_multiply_c3d_c3d_gs (pw_out, pw1, pw2, alpha)
8029 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw_out
8030 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
8031 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw2
8032 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: alpha
8034 CHARACTER(len=*),
PARAMETER :: routineN =
'pw_multiply'
8037 REAL(KIND=dp) :: my_alpha
8039 CALL timeset(routinen, handle)
8042 IF (
PRESENT(alpha)) my_alpha = alpha
8044 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
8045 cpabort(
"pw_multiply not implemented for non-identical grids!")
8047#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
8048 IF (my_alpha == 1.0_dp)
THEN
8050 pw_out%array = pw_out%array + pw1%array* pw2%array
8054 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
8058 IF (my_alpha == 1.0_dp)
THEN
8059 pw_out%array = pw_out%array + pw1%array* pw2%array
8061 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
8065 CALL timestop(handle)
8067 END SUBROUTINE pw_multiply_c3d_c3d_gs
8074 SUBROUTINE pw_multiply_with_c3d_c3d_gs (pw1, pw2)
8075 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw1
8076 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw2
8078 CHARACTER(LEN=*),
PARAMETER :: routineN =
'pw_multiply_with'
8082 CALL timeset(routinen, handle)
8084 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
8085 cpabort(
"Incompatible grids!")
8088 pw1%array = pw1%array* pw2%array
8091 CALL timestop(handle)
8093 END SUBROUTINE pw_multiply_with_c3d_c3d_gs
8108 FUNCTION pw_integral_ab_c3d_c3d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
8110 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
8111 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw2
8112 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
8113 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
8114 REAL(kind=dp) :: integral_value
8116 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab_c3d_c3d_gs'
8118 INTEGER :: handle, loc_sumtype
8119 LOGICAL :: my_just_sum, my_local_only
8121 CALL timeset(routinen, handle)
8124 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
8126 my_local_only = .false.
8127 IF (
PRESENT(local_only)) my_local_only = local_only
8129 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
8130 cpabort(
"Grids incompatible")
8133 my_just_sum = .false.
8134 IF (
PRESENT(just_sum)) my_just_sum = just_sum
8141 integral_value = sum(real(conjg(pw1%array) &
8142 *pw2%array, kind=dp))
8147 integral_value = accurate_sum(real(conjg(pw1%array)*pw2%array, kind=dp))
8152 IF (.NOT. my_just_sum)
THEN
8153 integral_value = integral_value*pw1%pw_grid%vol
8156 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
8157 CALL pw1%pw_grid%para%group%sum(integral_value)
8159 CALL timestop(handle)
8161 END FUNCTION pw_integral_ab_c3d_c3d_gs
8184 SUBROUTINE pw_gather_s_r1d_r3d_2(pw1, pw2, scale)
8186 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
8187 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw2
8188 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
8190 CALL pw_gather_s_r1d_r3d (pw2, pw1%array, scale)
8192 END SUBROUTINE pw_gather_s_r1d_r3d_2
8203 SUBROUTINE pw_gather_s_r1d_r3d (pw, c, scale)
8205 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw
8206 REAL(kind=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(IN),
TARGET :: c
8207 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
8209 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_gather_s'
8211 INTEGER :: gpt, handle
8212 REAL(kind=dp),
CONTIGUOUS,
DIMENSION(:),
POINTER :: c_1d
8214 CALL timeset(routinen, handle)
8217 c_1d(1:
SIZE(c)) => c
8218 associate(ngpts =>
SIZE(pw%pw_grid%gsq), g_hatmap => pw%pw_grid%g_hatmap)
8220 IF (
PRESENT(scale))
THEN
8223 pw%array(gpt) = scale* c_1d(g_hatmap(gpt, 1) + 1)
8229 pw%array(gpt) = c_1d(g_hatmap(gpt, 1) + 1)
8236 CALL timestop(handle)
8238 END SUBROUTINE pw_gather_s_r1d_r3d
8249 SUBROUTINE pw_scatter_s_r1d_r3d_2(pw1, pw2, scale)
8251 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
8252 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw2
8253 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
8255 CALL pw_scatter_s_r1d_r3d (pw1, pw2%array, scale)
8257 END SUBROUTINE pw_scatter_s_r1d_r3d_2
8268 SUBROUTINE pw_scatter_s_r1d_r3d (pw, c, scale)
8270 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw
8271 REAL(kind=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(INOUT),
TARGET :: c
8272 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
8274 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_scatter_s'
8276 INTEGER :: gpt, handle
8277 REAL(kind=dp),
CONTIGUOUS,
DIMENSION(:),
POINTER :: c_1d
8279 CALL timeset(routinen, handle)
8282 c_1d(1:
SIZE(c)) => c
8283 associate(g_hatmap => pw%pw_grid%g_hatmap, ngpts =>
SIZE(pw%pw_grid%gsq))
8286 IF (.NOT.
PRESENT(scale)) c = 0.0_dp
8288 IF (
PRESENT(scale))
THEN
8291 c_1d(g_hatmap(gpt, 1) + 1) = scale* pw%array(gpt)
8297 c_1d(g_hatmap(gpt, 1) + 1) = pw%array(gpt)
8302 IF (pw%pw_grid%grid_span == halfspace)
THEN
8304 IF (
PRESENT(scale))
THEN
8307 c_1d(g_hatmap(gpt, 2) + 1) = scale*( pw%array(gpt))
8313 c_1d(g_hatmap(gpt, 2) + 1) = ( pw%array(gpt))
8322 CALL timestop(handle)
8324 END SUBROUTINE pw_scatter_s_r1d_r3d
8345 SUBROUTINE pw_gather_s_r1d_c3d_2(pw1, pw2, scale)
8347 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
8348 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw2
8349 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
8351 CALL pw_gather_s_r1d_c3d (pw2, pw1%array, scale)
8353 END SUBROUTINE pw_gather_s_r1d_c3d_2
8364 SUBROUTINE pw_gather_s_r1d_c3d (pw, c, scale)
8366 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw
8367 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(IN),
TARGET :: c
8368 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
8370 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_gather_s'
8372 INTEGER :: gpt, handle
8373 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:),
POINTER :: c_1d
8375 CALL timeset(routinen, handle)
8378 c_1d(1:
SIZE(c)) => c
8379 associate(ngpts =>
SIZE(pw%pw_grid%gsq), g_hatmap => pw%pw_grid%g_hatmap)
8381 IF (
PRESENT(scale))
THEN
8384 pw%array(gpt) = scale* real(c_1d(g_hatmap(gpt, 1) + 1), kind=dp)
8390 pw%array(gpt) = real(c_1d(g_hatmap(gpt, 1) + 1), kind=dp)
8397 CALL timestop(handle)
8399 END SUBROUTINE pw_gather_s_r1d_c3d
8410 SUBROUTINE pw_scatter_s_r1d_c3d_2(pw1, pw2, scale)
8412 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
8413 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
8414 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
8416 CALL pw_scatter_s_r1d_c3d (pw1, pw2%array, scale)
8418 END SUBROUTINE pw_scatter_s_r1d_c3d_2
8429 SUBROUTINE pw_scatter_s_r1d_c3d (pw, c, scale)
8431 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw
8432 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(INOUT),
TARGET :: c
8433 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
8435 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_scatter_s'
8437 INTEGER :: gpt, handle
8438 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:),
POINTER :: c_1d
8440 CALL timeset(routinen, handle)
8443 c_1d(1:
SIZE(c)) => c
8444 associate(g_hatmap => pw%pw_grid%g_hatmap, ngpts =>
SIZE(pw%pw_grid%gsq))
8447 IF (.NOT.
PRESENT(scale)) c = 0.0_dp
8449 IF (
PRESENT(scale))
THEN
8452 c_1d(g_hatmap(gpt, 1) + 1) = scale* cmplx(pw%array(gpt), 0.0_dp, kind=dp)
8458 c_1d(g_hatmap(gpt, 1) + 1) = cmplx(pw%array(gpt), 0.0_dp, kind=dp)
8463 IF (pw%pw_grid%grid_span == halfspace)
THEN
8465 IF (
PRESENT(scale))
THEN
8468 c_1d(g_hatmap(gpt, 2) + 1) = scale*( cmplx(pw%array(gpt), 0.0_dp, kind=dp))
8474 c_1d(g_hatmap(gpt, 2) + 1) = ( cmplx(pw%array(gpt), 0.0_dp, kind=dp))
8483 CALL timestop(handle)
8485 END SUBROUTINE pw_scatter_s_r1d_c3d
8511 SUBROUTINE fft_wrap_pw1pw2_r3d_c1d_rs_gs (pw1, pw2, debug)
8513 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
8514 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
8515 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
8517 CHARACTER(len=*),
PARAMETER :: routineN =
'fft_wrap_pw1pw2'
8519 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
8520 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
8521 INTEGER :: handle, handle2, my_pos, nrays, &
8523 INTEGER,
DIMENSION(3) :: nloc
8524#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
8525 LOGICAL :: use_pw_gpu
8527 INTEGER,
DIMENSION(:),
POINTER :: n
8530 CALL timeset(routinen, handle2)
8531 out_unit = cp_logger_get_default_io_unit()
8532 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
8533 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
8538 IF (
PRESENT(debug))
THEN
8545 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
8546 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
8547 cpabort(
"PW grids not compatible")
8549 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
8550 cpabort(
"PW grids have not compatible MPI groups")
8554 n => pw1%pw_grid%npts
8556 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
8562 IF (test .AND. out_unit > 0)
THEN
8563 WRITE (out_unit,
'(A)')
" FFT Protocol "
8564 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
8565 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
8566 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
8569#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
8570 CALL pw_gpu_r3dc1d_3d(pw1, pw2)
8571#elif defined (__PW_FPGA)
8572 ALLOCATE (c_out(n(1), n(2), n(3)))
8575 IF (pw_fpga_init_bitstream(n) == 1)
THEN
8577#if (__PW_FPGA_SP && __PW_FPGA)
8578 CALL pw_fpga_r3dc1d_3d_sp(n, c_out)
8580 CALL pw_fpga_r3dc1d_3d_dp(n, c_out)
8582 CALL zdscal(n(1)*n(2)*n(3), 1.0_dp/pw1%pw_grid%ngpts, c_out, 1)
8583 CALL pw_gather_s_c1d_c3d(pw2, c_out)
8586 CALL fft3d(fwfft, n, c_out, debug=test)
8587 CALL pw_gather_s_c1d_c3d(pw2, c_out)
8591 ALLOCATE (c_out(n(1), n(2), n(3)))
8594 CALL fft3d(fwfft, n, c_out, debug=test)
8595 CALL pw_gather_s_c1d_c3d(pw2, c_out)
8599 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
8607 IF (test .AND. out_unit > 0)
THEN
8608 WRITE (out_unit,
'(A)')
" FFT Protocol "
8609 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
8610 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
8611 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
8614 my_pos = pw1%pw_grid%para%group%mepos
8615 nrays = pw1%pw_grid%para%nyzray(my_pos)
8616 grays => pw1%pw_grid%grays
8618#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
8620 use_pw_gpu = pw1%pw_grid%para%ray_distribution
8621 IF (use_pw_gpu)
THEN
8622 CALL pw_gpu_r3dc1d_3d_ps(pw1, pw2)
8626 nloc = pw1%pw_grid%npts_local
8627 ALLOCATE (c_in(nloc(1), nloc(2), nloc(3)))
8631 IF (pw1%pw_grid%para%ray_distribution)
THEN
8632 CALL fft3d(fwfft, n, c_in, grays, pw1%pw_grid%para%group, &
8633 pw1%pw_grid%para%yzp, pw1%pw_grid%para%nyzray, &
8634 pw1%pw_grid%para%bo, debug=test)
8636 CALL fft3d(fwfft, n, c_in, grays, pw1%pw_grid%para%group, &
8637 pw1%pw_grid%para%bo, debug=test)
8640 IF (test .AND. out_unit > 0) &
8641 WRITE (out_unit,
'(A)')
" PW_GATHER : 2d -> 1d "
8642 CALL pw_gather_p_c1d (pw2, grays)
8645#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
8650 IF (test .AND. out_unit > 0)
THEN
8651 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
8654 CALL timestop(handle)
8655 CALL timestop(handle2)
8657 END SUBROUTINE fft_wrap_pw1pw2_r3d_c1d_rs_gs
8676 SUBROUTINE fft_wrap_pw1pw2_r3d_c3d_rs_gs (pw1, pw2, debug)
8678 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
8679 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
8680 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
8682 CHARACTER(len=*),
PARAMETER :: routineN =
'fft_wrap_pw1pw2'
8684 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
8685 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
8686 INTEGER :: handle, handle2, my_pos, nrays, &
8688 INTEGER,
DIMENSION(:),
POINTER :: n
8691 CALL timeset(routinen, handle2)
8692 out_unit = cp_logger_get_default_io_unit()
8693 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
8694 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
8699 IF (
PRESENT(debug))
THEN
8706 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
8707 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
8708 cpabort(
"PW grids not compatible")
8710 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
8711 cpabort(
"PW grids have not compatible MPI groups")
8715 n => pw1%pw_grid%npts
8717 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
8723 IF (test .AND. out_unit > 0)
THEN
8724 WRITE (out_unit,
'(A)')
" FFT Protocol "
8725 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
8726 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
8727 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
8730 pw2%array = cmplx(pw1%array, 0.0_dp, kind=dp)
8732 CALL fft3d(fwfft, n, c_out, debug=test)
8734 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
8742 IF (test .AND. out_unit > 0)
THEN
8743 WRITE (out_unit,
'(A)')
" FFT Protocol "
8744 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
8745 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
8746 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
8749 my_pos = pw1%pw_grid%para%group%mepos
8750 nrays = pw1%pw_grid%para%nyzray(my_pos)
8751 grays => pw1%pw_grid%grays
8755 IF (test .AND. out_unit > 0)
THEN
8756 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
8759 CALL timestop(handle)
8760 CALL timestop(handle2)
8762 END SUBROUTINE fft_wrap_pw1pw2_r3d_c3d_rs_gs
8787 SUBROUTINE fft_wrap_pw1pw2_c1d_r3d_gs_rs (pw1, pw2, debug)
8789 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
8790 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw2
8791 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
8793 CHARACTER(len=*),
PARAMETER :: routineN =
'fft_wrap_pw1pw2'
8795 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
8796 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
8797 INTEGER :: handle, handle2, my_pos, nrays, &
8799 INTEGER,
DIMENSION(3) :: nloc
8800#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
8801 LOGICAL :: use_pw_gpu
8803 INTEGER,
DIMENSION(:),
POINTER :: n
8806 CALL timeset(routinen, handle2)
8807 out_unit = cp_logger_get_default_io_unit()
8808 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
8809 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
8814 IF (
PRESENT(debug))
THEN
8821 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
8822 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
8823 cpabort(
"PW grids not compatible")
8825 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
8826 cpabort(
"PW grids have not compatible MPI groups")
8830 n => pw1%pw_grid%npts
8832 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
8838 IF (test .AND. out_unit > 0)
THEN
8839 WRITE (out_unit,
'(A)')
" FFT Protocol "
8840 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
8841 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
8842 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
8845#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
8846 CALL pw_gpu_c1dr3d_3d(pw1, pw2)
8847#elif defined (__PW_FPGA)
8848 ALLOCATE (c_out(n(1), n(2), n(3)))
8851 IF (pw_fpga_init_bitstream(n) == 1)
THEN
8852 CALL pw_scatter_s_c1d_c3d(pw1, c_out)
8854#if (__PW_FPGA_SP && __PW_FPGA)
8855 CALL pw_fpga_c1dr3d_3d_sp(n, c_out)
8857 CALL pw_fpga_c1dr3d_3d_dp(n, c_out)
8859 CALL zdscal(n(1)*n(2)*n(3), 1.0_dp, c_out, 1)
8863 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" PW_SCATTER : 3d -> 1d "
8864 CALL pw_scatter_s_c1d_c3d(pw1, c_out)
8866 CALL fft3d(bwfft, n, c_out, debug=test)
8868 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" REAL part "
8873 ALLOCATE (c_out(n(1), n(2), n(3)))
8874 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" PW_SCATTER : 3d -> 1d "
8875 CALL pw_scatter_s_c1d_c3d(pw1, c_out)
8877 CALL fft3d(bwfft, n, c_out, debug=test)
8879 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" REAL part "
8884 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
8892 IF (test .AND. out_unit > 0)
THEN
8893 WRITE (out_unit,
'(A)')
" FFT Protocol "
8894 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
8895 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
8896 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
8899 my_pos = pw1%pw_grid%para%group%mepos
8900 nrays = pw1%pw_grid%para%nyzray(my_pos)
8901 grays => pw1%pw_grid%grays
8903#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
8905 use_pw_gpu = pw1%pw_grid%para%ray_distribution
8906 IF (use_pw_gpu)
THEN
8907 CALL pw_gpu_c1dr3d_3d_ps(pw1, pw2)
8911 IF (test .AND. out_unit > 0) &
8912 WRITE (out_unit,
'(A)')
" PW_SCATTER : 2d -> 1d "
8914 CALL pw_scatter_p_c1d (pw1, grays)
8915 nloc = pw2%pw_grid%npts_local
8916 ALLOCATE (c_in(nloc(1), nloc(2), nloc(3)))
8918 IF (pw1%pw_grid%para%ray_distribution)
THEN
8919 CALL fft3d(bwfft, n, c_in, grays, pw1%pw_grid%para%group, &
8920 pw1%pw_grid%para%yzp, pw1%pw_grid%para%nyzray, &
8921 pw1%pw_grid%para%bo, debug=test)
8923 CALL fft3d(bwfft, n, c_in, grays, pw1%pw_grid%para%group, &
8924 pw1%pw_grid%para%bo, debug=test)
8927 IF (test .AND. out_unit > 0) &
8928 WRITE (out_unit,
'(A)')
" Real part "
8931#if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
8936 IF (test .AND. out_unit > 0)
THEN
8937 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
8940 CALL timestop(handle)
8941 CALL timestop(handle2)
8943 END SUBROUTINE fft_wrap_pw1pw2_c1d_r3d_gs_rs
8956 SUBROUTINE pw_gather_s_c1d_r3d_2(pw1, pw2, scale)
8958 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
8959 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
8960 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
8962 CALL pw_gather_s_c1d_r3d (pw2, pw1%array, scale)
8964 END SUBROUTINE pw_gather_s_c1d_r3d_2
8975 SUBROUTINE pw_gather_s_c1d_r3d (pw, c, scale)
8977 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
8978 REAL(kind=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(IN),
TARGET :: c
8979 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
8981 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_gather_s'
8983 INTEGER :: gpt, handle
8984 REAL(kind=dp),
CONTIGUOUS,
DIMENSION(:),
POINTER :: c_1d
8986 CALL timeset(routinen, handle)
8989 c_1d(1:
SIZE(c)) => c
8990 associate(ngpts =>
SIZE(pw%pw_grid%gsq), g_hatmap => pw%pw_grid%g_hatmap)
8992 IF (
PRESENT(scale))
THEN
8995 pw%array(gpt) = scale* cmplx(c_1d(g_hatmap(gpt, 1) + 1), 0.0_dp, kind=dp)
9001 pw%array(gpt) = cmplx(c_1d(g_hatmap(gpt, 1) + 1), 0.0_dp, kind=dp)
9008 CALL timestop(handle)
9010 END SUBROUTINE pw_gather_s_c1d_r3d
9021 SUBROUTINE pw_scatter_s_c1d_r3d_2(pw1, pw2, scale)
9023 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
9024 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw2
9025 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
9027 CALL pw_scatter_s_c1d_r3d (pw1, pw2%array, scale)
9029 END SUBROUTINE pw_scatter_s_c1d_r3d_2
9040 SUBROUTINE pw_scatter_s_c1d_r3d (pw, c, scale)
9042 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
9043 REAL(kind=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(INOUT),
TARGET :: c
9044 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
9046 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_scatter_s'
9048 INTEGER :: gpt, handle
9049 REAL(kind=dp),
CONTIGUOUS,
DIMENSION(:),
POINTER :: c_1d
9051 CALL timeset(routinen, handle)
9054 c_1d(1:
SIZE(c)) => c
9055 associate(g_hatmap => pw%pw_grid%g_hatmap, ngpts =>
SIZE(pw%pw_grid%gsq))
9058 IF (.NOT.
PRESENT(scale)) c = 0.0_dp
9060 IF (
PRESENT(scale))
THEN
9063 c_1d(g_hatmap(gpt, 1) + 1) = scale* real(pw%array(gpt), kind=dp)
9069 c_1d(g_hatmap(gpt, 1) + 1) = real(pw%array(gpt), kind=dp)
9074 IF (pw%pw_grid%grid_span == halfspace)
THEN
9076 IF (
PRESENT(scale))
THEN
9079 c_1d(g_hatmap(gpt, 2) + 1) = scale*( real(pw%array(gpt), kind=dp))
9085 c_1d(g_hatmap(gpt, 2) + 1) = ( real(pw%array(gpt), kind=dp))
9094 CALL timestop(handle)
9096 END SUBROUTINE pw_scatter_s_c1d_r3d
9118 SUBROUTINE fft_wrap_pw1pw2_c1d_c3d_gs_rs (pw1, pw2, debug)
9120 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
9121 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw2
9122 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
9124 CHARACTER(len=*),
PARAMETER :: routineN =
'fft_wrap_pw1pw2'
9126 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
9127 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
9128 INTEGER :: handle, handle2, my_pos, nrays, &
9130 INTEGER,
DIMENSION(:),
POINTER :: n
9133 CALL timeset(routinen, handle2)
9134 out_unit = cp_logger_get_default_io_unit()
9135 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
9136 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
9141 IF (
PRESENT(debug))
THEN
9148 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
9149 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
9150 cpabort(
"PW grids not compatible")
9152 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
9153 cpabort(
"PW grids have not compatible MPI groups")
9157 n => pw1%pw_grid%npts
9159 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
9165 IF (test .AND. out_unit > 0)
THEN
9166 WRITE (out_unit,
'(A)')
" FFT Protocol "
9167 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
9168 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
9169 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
9173 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" PW_SCATTER : 3d -> 1d "
9174 CALL pw_scatter_s_c1d_c3d(pw1, c_out)
9175 CALL fft3d(bwfft, n, c_out, debug=test)
9177 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9185 IF (test .AND. out_unit > 0)
THEN
9186 WRITE (out_unit,
'(A)')
" FFT Protocol "
9187 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
9188 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
9189 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
9192 my_pos = pw1%pw_grid%para%group%mepos
9193 nrays = pw1%pw_grid%para%nyzray(my_pos)
9194 grays => pw1%pw_grid%grays
9197 IF (test .AND. out_unit > 0) &
9198 WRITE (out_unit,
'(A)')
" PW_SCATTER : 2d -> 1d "
9200 CALL pw_scatter_p_c1d (pw1, grays)
9203 IF (pw1%pw_grid%para%ray_distribution)
THEN
9204 CALL fft3d(bwfft, n, c_in, grays, pw1%pw_grid%para%group, &
9205 pw1%pw_grid%para%yzp, pw1%pw_grid%para%nyzray, &
9206 pw1%pw_grid%para%bo, debug=test)
9208 CALL fft3d(bwfft, n, c_in, grays, pw1%pw_grid%para%group, &
9209 pw1%pw_grid%para%bo, debug=test)
9214 IF (test .AND. out_unit > 0)
THEN
9215 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9218 CALL timestop(handle)
9219 CALL timestop(handle2)
9221 END SUBROUTINE fft_wrap_pw1pw2_c1d_c3d_gs_rs
9234 SUBROUTINE pw_gather_s_c1d_c3d_2(pw1, pw2, scale)
9236 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
9237 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
9238 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
9240 CALL pw_gather_s_c1d_c3d (pw2, pw1%array, scale)
9242 END SUBROUTINE pw_gather_s_c1d_c3d_2
9253 SUBROUTINE pw_gather_s_c1d_c3d (pw, c, scale)
9255 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
9256 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(IN),
TARGET :: c
9257 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
9259 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_gather_s'
9261 INTEGER :: gpt, handle
9262 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:),
POINTER :: c_1d
9264 CALL timeset(routinen, handle)
9267 c_1d(1:
SIZE(c)) => c
9268 associate(ngpts =>
SIZE(pw%pw_grid%gsq), g_hatmap => pw%pw_grid%g_hatmap)
9270 IF (
PRESENT(scale))
THEN
9273 pw%array(gpt) = scale* c_1d(g_hatmap(gpt, 1) + 1)
9279 pw%array(gpt) = c_1d(g_hatmap(gpt, 1) + 1)
9286 CALL timestop(handle)
9288 END SUBROUTINE pw_gather_s_c1d_c3d
9299 SUBROUTINE pw_scatter_s_c1d_c3d_2(pw1, pw2, scale)
9301 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
9302 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
9303 REAL(KIND=dp),
INTENT(IN),
OPTIONAL :: scale
9305 CALL pw_scatter_s_c1d_c3d (pw1, pw2%array, scale)
9307 END SUBROUTINE pw_scatter_s_c1d_c3d_2
9318 SUBROUTINE pw_scatter_s_c1d_c3d (pw, c, scale)
9320 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
9321 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(INOUT),
TARGET :: c
9322 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
9324 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_scatter_s'
9326 INTEGER :: gpt, handle
9327 COMPLEX(KIND=dp),
CONTIGUOUS,
DIMENSION(:),
POINTER :: c_1d
9329 CALL timeset(routinen, handle)
9332 c_1d(1:
SIZE(c)) => c
9333 associate(g_hatmap => pw%pw_grid%g_hatmap, ngpts =>
SIZE(pw%pw_grid%gsq))
9336 IF (.NOT.
PRESENT(scale)) c = 0.0_dp
9338 IF (
PRESENT(scale))
THEN
9341 c_1d(g_hatmap(gpt, 1) + 1) = scale* pw%array(gpt)
9347 c_1d(g_hatmap(gpt, 1) + 1) = pw%array(gpt)
9352 IF (pw%pw_grid%grid_span == halfspace)
THEN
9354 IF (
PRESENT(scale))
THEN
9357 c_1d(g_hatmap(gpt, 2) + 1) = scale*conjg( pw%array(gpt))
9363 c_1d(g_hatmap(gpt, 2) + 1) = conjg( pw%array(gpt))
9372 CALL timestop(handle)
9374 END SUBROUTINE pw_scatter_s_c1d_c3d
9396 SUBROUTINE fft_wrap_pw1pw2_c3d_r3d_gs_rs (pw1, pw2, debug)
9398 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
9399 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw2
9400 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
9402 CHARACTER(len=*),
PARAMETER :: routineN =
'fft_wrap_pw1pw2'
9404 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
9405 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
9406 INTEGER :: handle, handle2, my_pos, nrays, &
9408 INTEGER,
DIMENSION(:),
POINTER :: n
9411 CALL timeset(routinen, handle2)
9412 out_unit = cp_logger_get_default_io_unit()
9413 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
9414 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
9419 IF (
PRESENT(debug))
THEN
9426 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
9427 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
9428 cpabort(
"PW grids not compatible")
9430 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
9431 cpabort(
"PW grids have not compatible MPI groups")
9435 n => pw1%pw_grid%npts
9437 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
9443 IF (test .AND. out_unit > 0)
THEN
9444 WRITE (out_unit,
'(A)')
" FFT Protocol "
9445 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
9446 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
9447 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
9451 ALLOCATE (c_out(n(1), n(2), n(3)))
9452 CALL fft3d(bwfft, n, c_in, c_out, debug=test)
9454 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" REAL part "
9455 pw2%array = real(c_out, kind=dp)
9458 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9466 IF (test .AND. out_unit > 0)
THEN
9467 WRITE (out_unit,
'(A)')
" FFT Protocol "
9468 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
9469 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
9470 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
9473 my_pos = pw1%pw_grid%para%group%mepos
9474 nrays = pw1%pw_grid%para%nyzray(my_pos)
9475 grays => pw1%pw_grid%grays
9479 IF (test .AND. out_unit > 0)
THEN
9480 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9483 CALL timestop(handle)
9484 CALL timestop(handle2)
9486 END SUBROUTINE fft_wrap_pw1pw2_c3d_r3d_gs_rs
9504 SUBROUTINE fft_wrap_pw1pw2_c3d_c1d_rs_gs (pw1, pw2, debug)
9506 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
9507 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
9508 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
9510 CHARACTER(len=*),
PARAMETER :: routineN =
'fft_wrap_pw1pw2'
9512 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
9513 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
9514 INTEGER :: handle, handle2, my_pos, nrays, &
9516 INTEGER,
DIMENSION(:),
POINTER :: n
9519 CALL timeset(routinen, handle2)
9520 out_unit = cp_logger_get_default_io_unit()
9521 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
9522 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
9527 IF (
PRESENT(debug))
THEN
9534 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
9535 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
9536 cpabort(
"PW grids not compatible")
9538 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
9539 cpabort(
"PW grids have not compatible MPI groups")
9543 n => pw1%pw_grid%npts
9545 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
9551 IF (test .AND. out_unit > 0)
THEN
9552 WRITE (out_unit,
'(A)')
" FFT Protocol "
9553 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
9554 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
9555 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
9559 ALLOCATE (c_out(n(1), n(2), n(3)))
9561 CALL fft3d(fwfft, n, c_in, c_out, debug=test)
9563 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" PW_GATHER : 3d -> 1d "
9564 CALL pw_gather_s_c1d_c3d(pw2, c_out)
9567 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9575 IF (test .AND. out_unit > 0)
THEN
9576 WRITE (out_unit,
'(A)')
" FFT Protocol "
9577 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
9578 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
9579 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
9582 my_pos = pw1%pw_grid%para%group%mepos
9583 nrays = pw1%pw_grid%para%nyzray(my_pos)
9584 grays => pw1%pw_grid%grays
9590 IF (pw1%pw_grid%para%ray_distribution)
THEN
9591 CALL fft3d(fwfft, n, c_in, grays, pw1%pw_grid%para%group, &
9592 pw1%pw_grid%para%yzp, pw1%pw_grid%para%nyzray, &
9593 pw1%pw_grid%para%bo, debug=test)
9595 CALL fft3d(fwfft, n, c_in, grays, pw1%pw_grid%para%group, &
9596 pw1%pw_grid%para%bo, debug=test)
9599 IF (test .AND. out_unit > 0) &
9600 WRITE (out_unit,
'(A)')
" PW_GATHER : 2d -> 1d "
9601 CALL pw_gather_p_c1d (pw2, grays)
9604 IF (test .AND. out_unit > 0)
THEN
9605 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9608 CALL timestop(handle)
9609 CALL timestop(handle2)
9611 END SUBROUTINE fft_wrap_pw1pw2_c3d_c1d_rs_gs
9630 SUBROUTINE fft_wrap_pw1pw2_c3d_c3d_rs_gs (pw1, pw2, debug)
9632 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
9633 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
9634 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
9636 CHARACTER(len=*),
PARAMETER :: routineN =
'fft_wrap_pw1pw2'
9638 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
9639 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
9640 INTEGER :: handle, handle2, my_pos, nrays, &
9642 INTEGER,
DIMENSION(:),
POINTER :: n
9645 CALL timeset(routinen, handle2)
9646 out_unit = cp_logger_get_default_io_unit()
9647 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
9648 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
9653 IF (
PRESENT(debug))
THEN
9660 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
9661 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
9662 cpabort(
"PW grids not compatible")
9664 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
9665 cpabort(
"PW grids have not compatible MPI groups")
9669 n => pw1%pw_grid%npts
9671 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
9677 IF (test .AND. out_unit > 0)
THEN
9678 WRITE (out_unit,
'(A)')
" FFT Protocol "
9679 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
9680 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
9681 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
9686 CALL fft3d(fwfft, n, c_in, c_out, debug=test)
9688 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9696 IF (test .AND. out_unit > 0)
THEN
9697 WRITE (out_unit,
'(A)')
" FFT Protocol "
9698 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
9699 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
9700 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
9703 my_pos = pw1%pw_grid%para%group%mepos
9704 nrays = pw1%pw_grid%para%nyzray(my_pos)
9705 grays => pw1%pw_grid%grays
9709 IF (test .AND. out_unit > 0)
THEN
9710 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9713 CALL timestop(handle)
9714 CALL timestop(handle2)
9716 END SUBROUTINE fft_wrap_pw1pw2_c3d_c3d_rs_gs
9731 SUBROUTINE fft_wrap_pw1pw2_c3d_c3d_gs_rs (pw1, pw2, debug)
9733 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
9734 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw2
9735 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
9737 CHARACTER(len=*),
PARAMETER :: routineN =
'fft_wrap_pw1pw2'
9739 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
9740 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
9741 INTEGER :: handle, handle2, my_pos, nrays, &
9743 INTEGER,
DIMENSION(:),
POINTER :: n
9746 CALL timeset(routinen, handle2)
9747 out_unit = cp_logger_get_default_io_unit()
9748 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
9749 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
9754 IF (
PRESENT(debug))
THEN
9761 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
9762 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
9763 cpabort(
"PW grids not compatible")
9765 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
9766 cpabort(
"PW grids have not compatible MPI groups")
9770 n => pw1%pw_grid%npts
9772 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
9778 IF (test .AND. out_unit > 0)
THEN
9779 WRITE (out_unit,
'(A)')
" FFT Protocol "
9780 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
9781 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
9782 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
9787 CALL fft3d(bwfft, n, c_in, c_out, debug=test)
9789 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9797 IF (test .AND. out_unit > 0)
THEN
9798 WRITE (out_unit,
'(A)')
" FFT Protocol "
9799 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
9800 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
9801 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
9804 my_pos = pw1%pw_grid%para%group%mepos
9805 nrays = pw1%pw_grid%para%nyzray(my_pos)
9806 grays => pw1%pw_grid%grays
9810 IF (test .AND. out_unit > 0)
THEN
9811 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9814 CALL timestop(handle)
9815 CALL timestop(handle2)
9817 END SUBROUTINE fft_wrap_pw1pw2_c3d_c3d_gs_rs
9836 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
9837 REAL(kind=dp),
INTENT(IN) :: omega
9839 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_gauss_damp'
9841 INTEGER :: handle, i
9842 REAL(kind=dp) :: omega_2, tmp
9844 CALL timeset(routinen, handle)
9845 cpassert(omega >= 0)
9847 omega_2 = omega*omega
9848 omega_2 = 0.25_dp/omega_2
9851 DO i = 1,
SIZE(pw%array)
9852 tmp = exp(-pw%pw_grid%gsq(i)*omega_2)
9853 pw%array(i) = pw%array(i)*tmp
9857 CALL timestop(handle)
9870 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
9871 REAL(kind=dp),
INTENT(IN) :: omega
9873 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_log_deriv_gauss'
9875 INTEGER :: handle, i
9876 REAL(kind=dp) :: omega_2, tmp
9878 CALL timeset(routinen, handle)
9879 cpassert(omega >= 0)
9881 omega_2 = omega*omega
9882 omega_2 = 0.25_dp/omega_2
9885 DO i = 1,
SIZE(pw%array)
9886 tmp = 1.0_dp + omega_2*pw%pw_grid%gsq(i)
9887 pw%array(i) = pw%array(i)*tmp
9891 CALL timestop(handle)
9909 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
9910 REAL(kind=dp),
INTENT(IN) :: omega
9912 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_compl_gauss_damp'
9914 INTEGER :: cnt, handle, i
9915 REAL(kind=dp) :: omega_2, tmp, tmp2
9917 CALL timeset(routinen, handle)
9919 omega_2 = omega*omega
9920 omega_2 = 0.25_dp/omega_2
9922 cnt =
SIZE(pw%array)
9926 tmp = -omega_2*pw%pw_grid%gsq(i)
9927 IF (abs(tmp) > 1.0e-5_dp)
THEN
9928 tmp2 = 1.0_dp - exp(tmp)
9930 tmp2 = tmp + 0.5_dp*tmp*(tmp + (1.0_dp/3.0_dp)*tmp**2)
9932 pw%array(i) = pw%array(i)*tmp2
9936 CALL timestop(handle)
9948 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
9949 REAL(kind=dp),
INTENT(IN) :: omega
9951 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_log_deriv_compl_gauss'
9953 INTEGER :: handle, i
9954 REAL(kind=dp) :: omega_2, tmp, tmp2
9956 CALL timeset(routinen, handle)
9958 omega_2 = omega*omega
9959 omega_2 = 0.25_dp/omega_2
9963 DO i = 1,
SIZE(pw%array)
9964 tmp = omega_2*pw%pw_grid%gsq(i)
9966 IF (abs(tmp) >= 0.003_dp)
THEN
9967 tmp2 = 1.0_dp - tmp*exp(-tmp)/(1.0_dp - exp(-tmp))
9969 tmp2 = 0.5_dp*tmp - tmp**2/12.0_dp
9971 pw%array(i) = pw%array(i)*tmp2
9975 CALL timestop(handle)
9996 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
9997 REAL(kind=dp),
INTENT(IN) :: omega, scale_coul, scale_long
9999 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_gauss_damp_mix'
10001 INTEGER :: handle, i
10002 REAL(kind=dp) :: omega_2, tmp
10004 CALL timeset(routinen, handle)
10006 omega_2 = omega*omega
10007 omega_2 = 0.25_dp/omega_2
10010 DO i = 1,
SIZE(pw%array)
10011 tmp = scale_coul + scale_long*exp(-pw%pw_grid%gsq(i)*omega_2)
10012 pw%array(i) = pw%array(i)*tmp
10016 CALL timestop(handle)
10031 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
10032 REAL(kind=dp),
INTENT(IN) :: omega, scale_coul, scale_long
10034 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_log_deriv_mix_cl'
10036 INTEGER :: handle, i
10037 REAL(kind=dp) :: omega_2, tmp, tmp2
10039 CALL timeset(routinen, handle)
10041 omega_2 = omega*omega
10042 omega_2 = 0.25_dp/omega_2
10046 DO i = 1,
SIZE(pw%array)
10047 tmp = omega_2*pw%pw_grid%gsq(i)
10048 tmp2 = 1.0_dp + scale_long*tmp*exp(-tmp)/(scale_coul + scale_long*exp(-tmp))
10049 pw%array(i) = pw%array(i)*tmp2
10053 CALL timestop(handle)
10072 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
10073 REAL(kind=dp),
INTENT(IN) :: rcutoff
10075 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_truncated'
10077 INTEGER :: handle, i
10078 REAL(kind=dp) :: tmp, tmp2
10080 CALL timeset(routinen, handle)
10081 cpassert(rcutoff >= 0)
10084 DO i = 1,
SIZE(pw%array)
10085 tmp = sqrt(pw%pw_grid%gsq(i))*rcutoff
10086 IF (tmp >= 0.005_dp)
THEN
10087 tmp2 = 1.0_dp - cos(tmp)
10089 tmp2 = tmp**2/2.0_dp*(1.0_dp - tmp**2/12.0_dp)
10091 pw%array(i) = pw%array(i)*tmp2
10095 CALL timestop(handle)
10108 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
10109 REAL(kind=dp),
INTENT(IN) :: rcutoff
10111 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_log_deriv_trunc'
10113 INTEGER :: handle, i
10114 REAL(kind=dp) :: rchalf, tmp, tmp2
10116 CALL timeset(routinen, handle)
10117 cpassert(rcutoff >= 0)
10119 rchalf = 0.5_dp*rcutoff
10122 DO i = 1,
SIZE(pw%array)
10123 tmp = rchalf*sqrt(pw%pw_grid%gsq(i))
10125 IF (abs(tmp) >= 0.0001_dp)
THEN
10126 tmp2 = 1.0_dp - tmp/tan(tmp)
10128 tmp2 = tmp**2*(1.0_dp + tmp**2/15.0_dp)/3.0_dp
10130 pw%array(i) = pw%array(i)*tmp2
10134 CALL timestop(handle)
10150 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
10151 INTEGER,
DIMENSION(3),
INTENT(IN) :: n
10153 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_derive'
10155 COMPLEX(KIND=dp) :: im
10156 INTEGER :: handle, m,
idx, idir
10158 CALL timeset(routinen, handle)
10161 cpabort(
"Nonnegative exponents are not supported!")
10167 IF (n(idir) == 1)
THEN
10169 DO idx = 1,
SIZE(pw%array)
10170 pw%array(
idx) = pw%array(
idx)*pw%pw_grid%g(idir,
idx)
10173 ELSE IF (n(idir) > 1)
THEN
10175 DO idx = 1,
SIZE(pw%array)
10176 pw%array(
idx) = pw%array(
idx)*pw%pw_grid%g(idir,
idx)**n(idir)
10184 IF (abs(real(im, kind=dp) - 1.0_dp) > 1.0e-10_dp)
THEN
10186 pw%array(:) = im*pw%array(:)
10190 CALL timestop(handle)
10205 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
10207 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_laplace'
10211 CALL timeset(routinen, handle)
10214 pw%array(:) = -pw%array(:)*pw%pw_grid%gsq(:)
10217 CALL timestop(handle)
10234 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw, pwdr2
10235 INTEGER,
INTENT(IN) :: i, j
10237 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_dr2'
10239 INTEGER :: cnt, handle, ig
10240 REAL(kind=dp) :: gg, o3
10242 CALL timeset(routinen, handle)
10246 cnt =
SIZE(pw%array)
10251 gg = pw%pw_grid%g(i, ig)**2 - o3*pw%pw_grid%gsq(ig)
10252 pwdr2%array(ig) = gg*pw%array(ig)
10258 pwdr2%array(ig) = pw%array(ig)*(pw%pw_grid%g(i, ig)*pw%pw_grid%g(j, ig))
10263 CALL timestop(handle)
10282 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw, pwdr2_gg
10283 INTEGER,
INTENT(IN) :: i, j
10285 INTEGER :: cnt, handle, ig
10286 REAL(kind=dp) :: gg, o3
10287 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_dr2_gg'
10289 CALL timeset(routinen, handle)
10293 cnt =
SIZE(pw%array)
10297 DO ig = pw%pw_grid%first_gne0, cnt
10298 gg = pw%pw_grid%g(i, ig)**2 - o3*pw%pw_grid%gsq(ig)
10299 pwdr2_gg%array(ig) = gg/pw%pw_grid%gsq(ig)*pw%array(ig)
10304 DO ig = pw%pw_grid%first_gne0, cnt
10305 pwdr2_gg%array(ig) = pw%array(ig)*((pw%pw_grid%g(i, ig)*pw%pw_grid%g(j, ig)) &
10306 /pw%pw_grid%gsq(ig))
10311 IF (pw%pw_grid%have_g0) pwdr2_gg%array(1) = 0.0_dp
10313 CALL timestop(handle)
10330 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
10331 REAL(kind=dp),
INTENT(IN) :: ecut, sigma
10333 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_smoothing'
10335 INTEGER :: cnt, handle, ig
10336 REAL(kind=dp) :: arg, f
10338 CALL timeset(routinen, handle)
10340 cnt =
SIZE(pw%array)
10344 arg = (ecut - pw%pw_grid%gsq(ig))/sigma
10345 f = exp(arg)/(1 + exp(arg))
10346 pw%array(ig) = f*pw%array(ig)
10350 CALL timestop(handle)
10360 ELEMENTAL FUNCTION pw_compatible(grida, gridb)
RESULT(compat)
10362 TYPE(pw_grid_type),
INTENT(IN) :: grida, gridb
10367 IF (grida%id_nr == gridb%id_nr)
THEN
10369 ELSE IF (grida%reference == gridb%id_nr)
THEN
10371 ELSE IF (gridb%reference == grida%id_nr)
THEN
10375 END FUNCTION pw_compatible
10388 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: sf
10389 REAL(kind=dp),
DIMENSION(:),
INTENT(IN) :: r
10391 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_structure_factor'
10393 INTEGER :: cnt, handle, ig
10394 REAL(kind=dp) :: arg
10396 CALL timeset(routinen, handle)
10398 cnt =
SIZE(sf%array)
10402 arg = dot_product(sf%pw_grid%g(:, ig), r)
10403 sf%array(ig) = cmplx(cos(arg), -sin(arg), kind=dp)
10407 CALL timestop(handle)
static GRID_HOST_DEVICE int idx(const orbital a)
Return coset index of given orbital angular momentum.
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Defines the basic variable types.
integer, parameter, public dp
Machine interface based on Fortran 2003 and POSIX.
subroutine, public m_memory(mem)
Returns the total amount of memory [bytes] in use, if known, zero otherwise.
Definition of mathematical constants and functions.
complex(kind=dp), parameter, public gaussi
complex(kind=dp), parameter, public z_zero
subroutine, public pw_copy_match(pw1, pw2)
copy a pw type variable
integer function, public pw_fpga_init_bitstream(n)
Invoke the pw_fpga_check_bitstream C function passing the path to the data dir.
subroutine, public pw_fpga_r3dc1d_3d_dp(n, c_out)
perform an in-place double precision fft3d on the FPGA
subroutine, public pw_fpga_r3dc1d_3d_sp(n, c_out)
perform an in-place single precision fft3d on the FPGA
subroutine, public pw_fpga_c1dr3d_3d_dp(n, c_out)
perform an in-place double precision inverse fft3d on the FPGA
subroutine, public pw_fpga_c1dr3d_3d_sp(n, c_out)
perform an in-place single precision inverse fft3d on the FPGA
subroutine, public pw_gpu_c1dr3d_3d_ps(pw1, pw2)
perform an parallel scatter followed by a fft on the gpu
subroutine, public pw_gpu_c1dr3d_3d(pw1, pw2)
perform an scatter followed by a fft on the gpu
subroutine, public pw_gpu_r3dc1d_3d(pw1, pw2)
perform an fft followed by a gather on the gpu
subroutine, public pw_gpu_r3dc1d_3d_ps(pw1, pw2)
perform an parallel fft followed by a gather on the gpu
integer, parameter, public halfspace
integer, parameter, public pw_mode_local
integer, parameter, public pw_mode_distributed
subroutine, public pw_smoothing(pw, ecut, sigma)
Multiplies a G-space function with a smoothing factor of the form f(|G|) = exp((ecut - G^2)/sigma)/(1...
subroutine, public pw_log_deriv_trunc(pw, rcutoff)
Multiply all data points with the logarithmic derivative of the complementary cosine This function is...
subroutine, public pw_gauss_damp(pw, omega)
Multiply all data points with a Gaussian damping factor Needed for longrange Coulomb potential V(\vec...
subroutine, public pw_log_deriv_compl_gauss(pw, omega)
Multiply all data points with the logarithmic derivative of the complementary Gaussian damping factor...
subroutine, public pw_dr2_gg(pw, pwdr2_gg, i, j)
Calculate the tensorial 2nd derivative of a plane wave vector and divides by |G|^2....
subroutine, public pw_laplace(pw)
Calculate the Laplacian of a plane wave vector.
subroutine, public pw_log_deriv_gauss(pw, omega)
Multiply all data points with the logarithmic derivative of a Gaussian.
subroutine, public pw_gauss_damp_mix(pw, omega, scale_coul, scale_long)
Multiply all data points with a Gaussian damping factor and mixes it with the original function Neede...
subroutine, public pw_truncated(pw, rcutoff)
Multiply all data points with a complementary cosine Needed for truncated Coulomb potential V(\vec r)...
integer, parameter, public do_accurate_sum
integer, parameter, public do_standard_sum
subroutine, public pw_derive(pw, n)
Calculate the derivative of a plane wave vector.
subroutine, public pw_structure_factor(sf, r)
Calculate the structure factor for point r.
subroutine, public pw_compl_gauss_damp(pw, omega)
Multiply all data points with a Gaussian damping factor Needed for longrange Coulomb potential V(\vec...
subroutine, public pw_log_deriv_mix_cl(pw, omega, scale_coul, scale_long)
Multiply all data points with the logarithmic derivative of the mixed longrange/Coulomb potential Nee...
subroutine, public pw_dr2(pw, pwdr2, i, j)
Calculate the tensorial 2nd derivative of a plane wave vector.