33 USE kahan_sum,
ONLY: accurate_dot_product, &
60 #include "../base/base_uses.f90"
67 PUBLIC :: pw_copy, pw_axpy, pw_transfer, pw_scale
71 PUBLIC :: pw_integral_ab, pw_integral_a2b
72 PUBLIC ::
pw_dr2_gg, pw_integrate_function
74 PUBLIC :: pw_scatter, pw_gather
75 PUBLIC :: pw_copy_to_array, pw_copy_from_array
77 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'pw_methods'
78 LOGICAL,
PARAMETER,
PRIVATE :: debug_this_module = .false.
83 MODULE PROCEDURE pw_zero_r1d_rs
84 MODULE PROCEDURE pw_zero_r3d_rs
85 MODULE PROCEDURE pw_zero_c1d_rs
86 MODULE PROCEDURE pw_zero_c3d_rs
87 MODULE PROCEDURE pw_zero_r1d_gs
88 MODULE PROCEDURE pw_zero_r3d_gs
89 MODULE PROCEDURE pw_zero_c1d_gs
90 MODULE PROCEDURE pw_zero_c3d_gs
94 MODULE PROCEDURE pw_scale_r1d_rs
95 MODULE PROCEDURE pw_scale_r3d_rs
96 MODULE PROCEDURE pw_scale_c1d_rs
97 MODULE PROCEDURE pw_scale_c3d_rs
98 MODULE PROCEDURE pw_scale_r1d_gs
99 MODULE PROCEDURE pw_scale_r3d_gs
100 MODULE PROCEDURE pw_scale_c1d_gs
101 MODULE PROCEDURE pw_scale_c3d_gs
105 MODULE PROCEDURE pw_write_r1d_rs
106 MODULE PROCEDURE pw_write_r3d_rs
107 MODULE PROCEDURE pw_write_c1d_rs
108 MODULE PROCEDURE pw_write_c3d_rs
109 MODULE PROCEDURE pw_write_r1d_gs
110 MODULE PROCEDURE pw_write_r3d_gs
111 MODULE PROCEDURE pw_write_c1d_gs
112 MODULE PROCEDURE pw_write_c3d_gs
115 INTERFACE pw_integrate_function
116 MODULE PROCEDURE pw_integrate_function_r1d_rs
117 MODULE PROCEDURE pw_integrate_function_r3d_rs
118 MODULE PROCEDURE pw_integrate_function_c1d_rs
119 MODULE PROCEDURE pw_integrate_function_c3d_rs
120 MODULE PROCEDURE pw_integrate_function_r1d_gs
121 MODULE PROCEDURE pw_integrate_function_r3d_gs
122 MODULE PROCEDURE pw_integrate_function_c1d_gs
123 MODULE PROCEDURE pw_integrate_function_c3d_gs
127 MODULE PROCEDURE pw_set_value_r1d_rs
128 MODULE PROCEDURE pw_zero_r1d_rs
129 MODULE PROCEDURE pw_set_value_r3d_rs
130 MODULE PROCEDURE pw_zero_r3d_rs
131 MODULE PROCEDURE pw_set_value_c1d_rs
132 MODULE PROCEDURE pw_zero_c1d_rs
133 MODULE PROCEDURE pw_set_value_c3d_rs
134 MODULE PROCEDURE pw_zero_c3d_rs
135 MODULE PROCEDURE pw_set_value_r1d_gs
136 MODULE PROCEDURE pw_zero_r1d_gs
137 MODULE PROCEDURE pw_set_value_r3d_gs
138 MODULE PROCEDURE pw_zero_r3d_gs
139 MODULE PROCEDURE pw_set_value_c1d_gs
140 MODULE PROCEDURE pw_zero_c1d_gs
141 MODULE PROCEDURE pw_set_value_c3d_gs
142 MODULE PROCEDURE pw_zero_c3d_gs
146 MODULE PROCEDURE pw_copy_r1d_r1d_rs
147 MODULE PROCEDURE pw_copy_r1d_c1d_rs
148 MODULE PROCEDURE pw_copy_r3d_r3d_rs
149 MODULE PROCEDURE pw_copy_r3d_c3d_rs
150 MODULE PROCEDURE pw_copy_c1d_r1d_rs
151 MODULE PROCEDURE pw_copy_c1d_c1d_rs
152 MODULE PROCEDURE pw_copy_c3d_r3d_rs
153 MODULE PROCEDURE pw_copy_c3d_c3d_rs
154 MODULE PROCEDURE pw_copy_r1d_r1d_gs
155 MODULE PROCEDURE pw_copy_r1d_c1d_gs
156 MODULE PROCEDURE pw_copy_r3d_r3d_gs
157 MODULE PROCEDURE pw_copy_r3d_c3d_gs
158 MODULE PROCEDURE pw_copy_c1d_r1d_gs
159 MODULE PROCEDURE pw_copy_c1d_c1d_gs
160 MODULE PROCEDURE pw_copy_c3d_r3d_gs
161 MODULE PROCEDURE pw_copy_c3d_c3d_gs
165 MODULE PROCEDURE pw_axpy_r1d_r1d_rs
166 MODULE PROCEDURE pw_axpy_r1d_c1d_rs
167 MODULE PROCEDURE pw_axpy_r3d_r3d_rs
168 MODULE PROCEDURE pw_axpy_r3d_c3d_rs
169 MODULE PROCEDURE pw_axpy_c1d_r1d_rs
170 MODULE PROCEDURE pw_axpy_c1d_c1d_rs
171 MODULE PROCEDURE pw_axpy_c3d_r3d_rs
172 MODULE PROCEDURE pw_axpy_c3d_c3d_rs
173 MODULE PROCEDURE pw_axpy_r1d_r1d_gs
174 MODULE PROCEDURE pw_axpy_r1d_c1d_gs
175 MODULE PROCEDURE pw_axpy_r3d_r3d_gs
176 MODULE PROCEDURE pw_axpy_r3d_c3d_gs
177 MODULE PROCEDURE pw_axpy_c1d_r1d_gs
178 MODULE PROCEDURE pw_axpy_c1d_c1d_gs
179 MODULE PROCEDURE pw_axpy_c3d_r3d_gs
180 MODULE PROCEDURE pw_axpy_c3d_c3d_gs
183 INTERFACE pw_multiply
184 MODULE PROCEDURE pw_multiply_r1d_r1d_rs
185 MODULE PROCEDURE pw_multiply_r1d_c1d_rs
186 MODULE PROCEDURE pw_multiply_r3d_r3d_rs
187 MODULE PROCEDURE pw_multiply_r3d_c3d_rs
188 MODULE PROCEDURE pw_multiply_c1d_r1d_rs
189 MODULE PROCEDURE pw_multiply_c1d_c1d_rs
190 MODULE PROCEDURE pw_multiply_c3d_r3d_rs
191 MODULE PROCEDURE pw_multiply_c3d_c3d_rs
192 MODULE PROCEDURE pw_multiply_r1d_r1d_gs
193 MODULE PROCEDURE pw_multiply_r1d_c1d_gs
194 MODULE PROCEDURE pw_multiply_r3d_r3d_gs
195 MODULE PROCEDURE pw_multiply_r3d_c3d_gs
196 MODULE PROCEDURE pw_multiply_c1d_r1d_gs
197 MODULE PROCEDURE pw_multiply_c1d_c1d_gs
198 MODULE PROCEDURE pw_multiply_c3d_r3d_gs
199 MODULE PROCEDURE pw_multiply_c3d_c3d_gs
202 INTERFACE pw_multiply_with
203 MODULE PROCEDURE pw_multiply_with_r1d_r1d_rs
204 MODULE PROCEDURE pw_multiply_with_r1d_c1d_rs
205 MODULE PROCEDURE pw_multiply_with_r3d_r3d_rs
206 MODULE PROCEDURE pw_multiply_with_r3d_c3d_rs
207 MODULE PROCEDURE pw_multiply_with_c1d_r1d_rs
208 MODULE PROCEDURE pw_multiply_with_c1d_c1d_rs
209 MODULE PROCEDURE pw_multiply_with_c3d_r3d_rs
210 MODULE PROCEDURE pw_multiply_with_c3d_c3d_rs
211 MODULE PROCEDURE pw_multiply_with_r1d_r1d_gs
212 MODULE PROCEDURE pw_multiply_with_r1d_c1d_gs
213 MODULE PROCEDURE pw_multiply_with_r3d_r3d_gs
214 MODULE PROCEDURE pw_multiply_with_r3d_c3d_gs
215 MODULE PROCEDURE pw_multiply_with_c1d_r1d_gs
216 MODULE PROCEDURE pw_multiply_with_c1d_c1d_gs
217 MODULE PROCEDURE pw_multiply_with_c3d_r3d_gs
218 MODULE PROCEDURE pw_multiply_with_c3d_c3d_gs
221 INTERFACE pw_integral_ab
222 MODULE PROCEDURE pw_integral_ab_r1d_r1d_rs
223 MODULE PROCEDURE pw_integral_ab_r1d_c1d_rs
224 MODULE PROCEDURE pw_integral_ab_r3d_r3d_rs
225 MODULE PROCEDURE pw_integral_ab_r3d_c3d_rs
226 MODULE PROCEDURE pw_integral_ab_c1d_r1d_rs
227 MODULE PROCEDURE pw_integral_ab_c1d_c1d_rs
228 MODULE PROCEDURE pw_integral_ab_c3d_r3d_rs
229 MODULE PROCEDURE pw_integral_ab_c3d_c3d_rs
230 MODULE PROCEDURE pw_integral_ab_r1d_r1d_gs
231 MODULE PROCEDURE pw_integral_ab_r1d_c1d_gs
232 MODULE PROCEDURE pw_integral_ab_r3d_r3d_gs
233 MODULE PROCEDURE pw_integral_ab_r3d_c3d_gs
234 MODULE PROCEDURE pw_integral_ab_c1d_r1d_gs
235 MODULE PROCEDURE pw_integral_ab_c1d_c1d_gs
236 MODULE PROCEDURE pw_integral_ab_c3d_r3d_gs
237 MODULE PROCEDURE pw_integral_ab_c3d_c3d_gs
240 INTERFACE pw_integral_a2b
241 MODULE PROCEDURE pw_integral_a2b_r1d_r1d
242 MODULE PROCEDURE pw_integral_a2b_r1d_c1d
243 MODULE PROCEDURE pw_integral_a2b_c1d_r1d
244 MODULE PROCEDURE pw_integral_a2b_c1d_c1d
248 MODULE PROCEDURE pw_gather_p_r1d
249 MODULE PROCEDURE pw_gather_p_c1d
250 MODULE PROCEDURE pw_gather_s_r1d_r3d
251 MODULE PROCEDURE pw_gather_s_r1d_c3d
252 MODULE PROCEDURE pw_gather_s_c1d_r3d
253 MODULE PROCEDURE pw_gather_s_c1d_c3d
257 MODULE PROCEDURE pw_scatter_p_r1d
258 MODULE PROCEDURE pw_scatter_p_c1d
259 MODULE PROCEDURE pw_scatter_s_r1d_r3d
260 MODULE PROCEDURE pw_scatter_s_r1d_c3d
261 MODULE PROCEDURE pw_scatter_s_c1d_r3d
262 MODULE PROCEDURE pw_scatter_s_c1d_c3d
265 INTERFACE pw_copy_to_array
266 MODULE PROCEDURE pw_copy_to_array_r1d_r1d_rs
267 MODULE PROCEDURE pw_copy_to_array_r1d_c1d_rs
268 MODULE PROCEDURE pw_copy_to_array_r3d_r3d_rs
269 MODULE PROCEDURE pw_copy_to_array_r3d_c3d_rs
270 MODULE PROCEDURE pw_copy_to_array_c1d_r1d_rs
271 MODULE PROCEDURE pw_copy_to_array_c1d_c1d_rs
272 MODULE PROCEDURE pw_copy_to_array_c3d_r3d_rs
273 MODULE PROCEDURE pw_copy_to_array_c3d_c3d_rs
274 MODULE PROCEDURE pw_copy_to_array_r1d_r1d_gs
275 MODULE PROCEDURE pw_copy_to_array_r1d_c1d_gs
276 MODULE PROCEDURE pw_copy_to_array_r3d_r3d_gs
277 MODULE PROCEDURE pw_copy_to_array_r3d_c3d_gs
278 MODULE PROCEDURE pw_copy_to_array_c1d_r1d_gs
279 MODULE PROCEDURE pw_copy_to_array_c1d_c1d_gs
280 MODULE PROCEDURE pw_copy_to_array_c3d_r3d_gs
281 MODULE PROCEDURE pw_copy_to_array_c3d_c3d_gs
284 INTERFACE pw_copy_from_array
285 MODULE PROCEDURE pw_copy_from_array_r1d_r1d_rs
286 MODULE PROCEDURE pw_copy_from_array_r1d_c1d_rs
287 MODULE PROCEDURE pw_copy_from_array_r3d_r3d_rs
288 MODULE PROCEDURE pw_copy_from_array_r3d_c3d_rs
289 MODULE PROCEDURE pw_copy_from_array_c1d_r1d_rs
290 MODULE PROCEDURE pw_copy_from_array_c1d_c1d_rs
291 MODULE PROCEDURE pw_copy_from_array_c3d_r3d_rs
292 MODULE PROCEDURE pw_copy_from_array_c3d_c3d_rs
293 MODULE PROCEDURE pw_copy_from_array_r1d_r1d_gs
294 MODULE PROCEDURE pw_copy_from_array_r1d_c1d_gs
295 MODULE PROCEDURE pw_copy_from_array_r3d_r3d_gs
296 MODULE PROCEDURE pw_copy_from_array_r3d_c3d_gs
297 MODULE PROCEDURE pw_copy_from_array_c1d_r1d_gs
298 MODULE PROCEDURE pw_copy_from_array_c1d_c1d_gs
299 MODULE PROCEDURE pw_copy_from_array_c3d_r3d_gs
300 MODULE PROCEDURE pw_copy_from_array_c3d_c3d_gs
303 INTERFACE pw_transfer
304 MODULE PROCEDURE pw_copy_r1d_r1d_rs
305 MODULE PROCEDURE pw_copy_r1d_r1d_gs
306 MODULE PROCEDURE pw_gather_s_r1d_r3d_2
307 MODULE PROCEDURE pw_scatter_s_r1d_r3d_2
308 MODULE PROCEDURE pw_copy_r1d_c1d_rs
309 MODULE PROCEDURE pw_copy_r1d_c1d_gs
310 MODULE PROCEDURE pw_gather_s_r1d_c3d_2
311 MODULE PROCEDURE pw_scatter_s_r1d_c3d_2
312 MODULE PROCEDURE pw_copy_r3d_r3d_rs
313 MODULE PROCEDURE pw_copy_r3d_r3d_gs
314 MODULE PROCEDURE fft_wrap_pw1pw2_r3d_c1d_rs_gs
315 MODULE PROCEDURE pw_copy_r3d_c3d_rs
316 MODULE PROCEDURE pw_copy_r3d_c3d_gs
317 MODULE PROCEDURE fft_wrap_pw1pw2_r3d_c3d_rs_gs
318 MODULE PROCEDURE pw_copy_c1d_r1d_rs
319 MODULE PROCEDURE pw_copy_c1d_r1d_gs
320 MODULE PROCEDURE pw_gather_s_c1d_r3d_2
321 MODULE PROCEDURE pw_scatter_s_c1d_r3d_2
322 MODULE PROCEDURE fft_wrap_pw1pw2_c1d_r3d_gs_rs
323 MODULE PROCEDURE pw_copy_c1d_c1d_rs
324 MODULE PROCEDURE pw_copy_c1d_c1d_gs
325 MODULE PROCEDURE pw_gather_s_c1d_c3d_2
326 MODULE PROCEDURE pw_scatter_s_c1d_c3d_2
327 MODULE PROCEDURE fft_wrap_pw1pw2_c1d_c3d_gs_rs
328 MODULE PROCEDURE pw_copy_c3d_r3d_rs
329 MODULE PROCEDURE pw_copy_c3d_r3d_gs
330 MODULE PROCEDURE fft_wrap_pw1pw2_c3d_r3d_gs_rs
331 MODULE PROCEDURE fft_wrap_pw1pw2_c3d_c1d_rs_gs
332 MODULE PROCEDURE pw_copy_c3d_c3d_rs
333 MODULE PROCEDURE pw_copy_c3d_c3d_gs
334 MODULE PROCEDURE fft_wrap_pw1pw2_c3d_c3d_rs_gs
335 MODULE PROCEDURE fft_wrap_pw1pw2_c3d_c3d_gs_rs
346 SUBROUTINE pw_zero_r1d_rs (pw)
348 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw
350 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_zero'
354 CALL timeset(routinen, handle)
360 CALL timestop(handle)
362 END SUBROUTINE pw_zero_r1d_rs
371 SUBROUTINE pw_scale_r1d_rs (pw, a)
373 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw
374 REAL(kind=
dp),
INTENT(IN) :: a
376 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_scale'
380 CALL timeset(routinen, handle)
383 pw%array = a*pw%array
386 CALL timestop(handle)
388 END SUBROUTINE pw_scale_r1d_rs
400 SUBROUTINE pw_write_r1d_rs (pw, unit_nr)
402 TYPE(pw_r1d_rs_type),
INTENT(in) :: pw
403 INTEGER,
INTENT(in) :: unit_nr
407 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
409 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=r1d"
410 IF (
ASSOCIATED(pw%array))
THEN
411 WRITE (unit=unit_nr, fmt=
"(' array=<r(',i8,':',i8,')>')") &
412 lbound(pw%array, 1), ubound(pw%array, 1)
414 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
417 END SUBROUTINE pw_write_r1d_rs
426 FUNCTION pw_integrate_function_r1d_rs (fun, isign, oprt)
RESULT(total_fun)
428 TYPE(pw_r1d_rs_type),
INTENT(IN) :: fun
429 INTEGER,
INTENT(IN),
OPTIONAL :: isign
430 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
431 REAL(kind=
dp) :: total_fun
437 IF (
PRESENT(oprt))
THEN
442 cpabort(
"Unknown operator")
450 total_fun = fun%pw_grid%dvol*accurate_sum(abs(fun%array))
452 total_fun = fun%pw_grid%dvol*accurate_sum( fun%array)
456 CALL fun%pw_grid%para%group%sum(total_fun)
459 IF (
PRESENT(isign))
THEN
460 total_fun = total_fun*sign(1._dp, real(isign,
dp))
463 END FUNCTION pw_integrate_function_r1d_rs
470 SUBROUTINE pw_set_value_r1d_rs (pw, value)
471 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw
472 REAL(kind=
dp),
INTENT(IN) ::
value
474 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_set_value'
478 CALL timeset(routinen, handle)
484 CALL timestop(handle)
486 END SUBROUTINE pw_set_value_r1d_rs
494 SUBROUTINE pw_zero_r1d_gs (pw)
496 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw
498 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_zero'
502 CALL timeset(routinen, handle)
508 CALL timestop(handle)
510 END SUBROUTINE pw_zero_r1d_gs
519 SUBROUTINE pw_scale_r1d_gs (pw, a)
521 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw
522 REAL(kind=
dp),
INTENT(IN) :: a
524 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_scale'
528 CALL timeset(routinen, handle)
531 pw%array = a*pw%array
534 CALL timestop(handle)
536 END SUBROUTINE pw_scale_r1d_gs
548 SUBROUTINE pw_write_r1d_gs (pw, unit_nr)
550 TYPE(pw_r1d_gs_type),
INTENT(in) :: pw
551 INTEGER,
INTENT(in) :: unit_nr
555 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
557 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=r1d"
558 IF (
ASSOCIATED(pw%array))
THEN
559 WRITE (unit=unit_nr, fmt=
"(' array=<r(',i8,':',i8,')>')") &
560 lbound(pw%array, 1), ubound(pw%array, 1)
562 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
565 END SUBROUTINE pw_write_r1d_gs
574 FUNCTION pw_integrate_function_r1d_gs (fun, isign, oprt)
RESULT(total_fun)
576 TYPE(pw_r1d_gs_type),
INTENT(IN) :: fun
577 INTEGER,
INTENT(IN),
OPTIONAL :: isign
578 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
579 REAL(kind=
dp) :: total_fun
585 IF (
PRESENT(oprt))
THEN
590 cpabort(
"Unknown operator")
597 cpabort(
"Operator ABS not implemented")
598 IF (fun%pw_grid%have_g0) total_fun = fun%pw_grid%vol* fun%array(1)
601 CALL fun%pw_grid%para%group%sum(total_fun)
604 IF (
PRESENT(isign))
THEN
605 total_fun = total_fun*sign(1._dp, real(isign,
dp))
608 END FUNCTION pw_integrate_function_r1d_gs
615 SUBROUTINE pw_set_value_r1d_gs (pw, value)
616 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw
617 REAL(kind=
dp),
INTENT(IN) ::
value
619 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_set_value'
623 CALL timeset(routinen, handle)
629 CALL timestop(handle)
631 END SUBROUTINE pw_set_value_r1d_gs
639 SUBROUTINE pw_gather_p_r1d (pw, c, scale)
641 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw
642 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
INTENT(IN) :: c
643 REAL(kind=
dp),
INTENT(IN),
OPTIONAL :: scale
645 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_gather_p'
647 INTEGER :: gpt, handle, l, m, mn, n
649 CALL timeset(routinen, handle)
652 cpabort(
"This grid type is not distributed")
655 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
656 ngpts =>
SIZE(pw%pw_grid%gsq), ghat => pw%pw_grid%g_hat, yzq => pw%pw_grid%para%yzq)
658 IF (
PRESENT(scale))
THEN
663 l = mapl(ghat(1, gpt)) + 1
664 m = mapm(ghat(2, gpt)) + 1
665 n = mapn(ghat(3, gpt)) + 1
667 pw%array(gpt) = scale* real(c(l, mn), kind=
dp)
675 l = mapl(ghat(1, gpt)) + 1
676 m = mapm(ghat(2, gpt)) + 1
677 n = mapn(ghat(3, gpt)) + 1
679 pw%array(gpt) = real(c(l, mn), kind=
dp)
686 CALL timestop(handle)
688 END SUBROUTINE pw_gather_p_r1d
696 SUBROUTINE pw_scatter_p_r1d (pw, c, scale)
697 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw
698 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
INTENT(INOUT) :: c
699 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
701 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_scatter_p'
703 INTEGER :: gpt, handle, l, m, mn, n
705 CALL timeset(routinen, handle)
707 IF (pw%pw_grid%para%mode /= pw_mode_distributed)
THEN
708 cpabort(
"This grid type is not distributed")
711 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
712 ghat => pw%pw_grid%g_hat, yzq => pw%pw_grid%para%yzq, ngpts =>
SIZE(pw%pw_grid%gsq))
714 IF (.NOT.
PRESENT(scale)) c = z_zero
716 IF (
PRESENT(scale))
THEN
721 l = mapl(ghat(1, gpt)) + 1
722 m = mapm(ghat(2, gpt)) + 1
723 n = mapn(ghat(3, gpt)) + 1
725 c(l, mn) = cmplx(scale*pw%array(gpt), 0.0_dp, kind=dp)
733 l = mapl(ghat(1, gpt)) + 1
734 m = mapm(ghat(2, gpt)) + 1
735 n = mapn(ghat(3, gpt)) + 1
737 c(l, mn) = cmplx(pw%array(gpt), 0.0_dp, kind=dp)
744 IF (pw%pw_grid%grid_span == halfspace)
THEN
746 associate(mapm => pw%pw_grid%mapm%neg, mapn => pw%pw_grid%mapn%neg, mapl => pw%pw_grid%mapl%neg, &
747 ghat => pw%pw_grid%g_hat, ngpts =>
SIZE(pw%pw_grid%gsq), yzq => pw%pw_grid%para%yzq)
749 IF (
PRESENT(scale))
THEN
754 l = mapl(ghat(1, gpt)) + 1
755 m = mapm(ghat(2, gpt)) + 1
756 n = mapn(ghat(3, gpt)) + 1
758 c(l, mn) = scale*( cmplx(pw%array(gpt), 0.0_dp, kind=dp))
766 l = mapl(ghat(1, gpt)) + 1
767 m = mapm(ghat(2, gpt)) + 1
768 n = mapn(ghat(3, gpt)) + 1
770 c(l, mn) = ( cmplx(pw%array(gpt), 0.0_dp, kind=dp))
777 CALL timestop(handle)
779 END SUBROUTINE pw_scatter_p_r1d
787 SUBROUTINE pw_zero_r3d_rs (pw)
789 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw
791 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_zero'
795 CALL timeset(routinen, handle)
801 CALL timestop(handle)
803 END SUBROUTINE pw_zero_r3d_rs
812 SUBROUTINE pw_scale_r3d_rs (pw, a)
814 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw
815 REAL(kind=dp),
INTENT(IN) :: a
817 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_scale'
821 CALL timeset(routinen, handle)
824 pw%array = a*pw%array
827 CALL timestop(handle)
829 END SUBROUTINE pw_scale_r3d_rs
841 SUBROUTINE pw_write_r3d_rs (pw, unit_nr)
843 TYPE(pw_r3d_rs_type),
INTENT(in) :: pw
844 INTEGER,
INTENT(in) :: unit_nr
848 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
850 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=r3d"
851 IF (
ASSOCIATED(pw%array))
THEN
852 WRITE (unit=unit_nr, fmt=
"(' array=<r(',i8,':',i8,',',i8,':',i8,',',i8,':',i8,')>')") &
853 lbound(pw%array, 1), ubound(pw%array, 1), lbound(pw%array, 2), ubound(pw%array, 2), &
854 lbound(pw%array, 3), ubound(pw%array, 3)
856 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
859 END SUBROUTINE pw_write_r3d_rs
868 FUNCTION pw_integrate_function_r3d_rs (fun, isign, oprt)
RESULT(total_fun)
870 TYPE(pw_r3d_rs_type),
INTENT(IN) :: fun
871 INTEGER,
INTENT(IN),
OPTIONAL :: isign
872 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
873 REAL(kind=dp) :: total_fun
879 IF (
PRESENT(oprt))
THEN
884 cpabort(
"Unknown operator")
892 total_fun = fun%pw_grid%dvol*accurate_sum(abs(fun%array))
894 total_fun = fun%pw_grid%dvol*accurate_sum( fun%array)
897 IF (fun%pw_grid%para%mode /= pw_mode_local)
THEN
898 CALL fun%pw_grid%para%group%sum(total_fun)
901 IF (
PRESENT(isign))
THEN
902 total_fun = total_fun*sign(1._dp, real(isign, dp))
905 END FUNCTION pw_integrate_function_r3d_rs
912 SUBROUTINE pw_set_value_r3d_rs (pw, value)
913 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw
914 REAL(kind=dp),
INTENT(IN) ::
value
916 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_set_value'
920 CALL timeset(routinen, handle)
926 CALL timestop(handle)
928 END SUBROUTINE pw_set_value_r3d_rs
936 SUBROUTINE pw_zero_r3d_gs (pw)
938 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw
940 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_zero'
944 CALL timeset(routinen, handle)
950 CALL timestop(handle)
952 END SUBROUTINE pw_zero_r3d_gs
961 SUBROUTINE pw_scale_r3d_gs (pw, a)
963 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw
964 REAL(kind=dp),
INTENT(IN) :: a
966 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_scale'
970 CALL timeset(routinen, handle)
973 pw%array = a*pw%array
976 CALL timestop(handle)
978 END SUBROUTINE pw_scale_r3d_gs
990 SUBROUTINE pw_write_r3d_gs (pw, unit_nr)
992 TYPE(pw_r3d_gs_type),
INTENT(in) :: pw
993 INTEGER,
INTENT(in) :: unit_nr
997 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
999 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=r3d"
1000 IF (
ASSOCIATED(pw%array))
THEN
1001 WRITE (unit=unit_nr, fmt=
"(' array=<r(',i8,':',i8,',',i8,':',i8,',',i8,':',i8,')>')") &
1002 lbound(pw%array, 1), ubound(pw%array, 1), lbound(pw%array, 2), ubound(pw%array, 2), &
1003 lbound(pw%array, 3), ubound(pw%array, 3)
1005 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
1008 END SUBROUTINE pw_write_r3d_gs
1017 FUNCTION pw_integrate_function_r3d_gs (fun, isign, oprt)
RESULT(total_fun)
1019 TYPE(pw_r3d_gs_type),
INTENT(IN) :: fun
1020 INTEGER,
INTENT(IN),
OPTIONAL :: isign
1021 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
1022 REAL(kind=dp) :: total_fun
1028 IF (
PRESENT(oprt))
THEN
1033 cpabort(
"Unknown operator")
1040 cpabort(
"Operator ABS not implemented")
1041 cpabort(
"Reciprocal space integration for 3D grids not implemented!")
1043 IF (fun%pw_grid%para%mode /= pw_mode_local)
THEN
1044 CALL fun%pw_grid%para%group%sum(total_fun)
1047 IF (
PRESENT(isign))
THEN
1048 total_fun = total_fun*sign(1._dp, real(isign, dp))
1051 END FUNCTION pw_integrate_function_r3d_gs
1058 SUBROUTINE pw_set_value_r3d_gs (pw, value)
1059 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw
1060 REAL(kind=dp),
INTENT(IN) ::
value
1062 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_set_value'
1066 CALL timeset(routinen, handle)
1072 CALL timestop(handle)
1074 END SUBROUTINE pw_set_value_r3d_gs
1083 SUBROUTINE pw_zero_c1d_rs (pw)
1085 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw
1087 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_zero'
1091 CALL timeset(routinen, handle)
1094 pw%array = cmplx(0.0_dp, 0.0_dp, kind=dp)
1097 CALL timestop(handle)
1099 END SUBROUTINE pw_zero_c1d_rs
1108 SUBROUTINE pw_scale_c1d_rs (pw, a)
1110 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw
1111 REAL(kind=dp),
INTENT(IN) :: a
1113 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_scale'
1117 CALL timeset(routinen, handle)
1120 pw%array = a*pw%array
1123 CALL timestop(handle)
1125 END SUBROUTINE pw_scale_c1d_rs
1137 SUBROUTINE pw_write_c1d_rs (pw, unit_nr)
1139 TYPE(pw_c1d_rs_type),
INTENT(in) :: pw
1140 INTEGER,
INTENT(in) :: unit_nr
1144 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
1146 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=c1d"
1147 IF (
ASSOCIATED(pw%array))
THEN
1148 WRITE (unit=unit_nr, fmt=
"(' array=<c(',i8,':',i8,')>')") &
1149 lbound(pw%array, 1), ubound(pw%array, 1)
1151 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
1154 END SUBROUTINE pw_write_c1d_rs
1163 FUNCTION pw_integrate_function_c1d_rs (fun, isign, oprt)
RESULT(total_fun)
1165 TYPE(pw_c1d_rs_type),
INTENT(IN) :: fun
1166 INTEGER,
INTENT(IN),
OPTIONAL :: isign
1167 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
1168 REAL(kind=dp) :: total_fun
1174 IF (
PRESENT(oprt))
THEN
1179 cpabort(
"Unknown operator")
1187 total_fun = fun%pw_grid%dvol*accurate_sum(abs(fun%array))
1189 total_fun = fun%pw_grid%dvol*accurate_sum( real(fun%array, kind=dp))
1192 IF (fun%pw_grid%para%mode /= pw_mode_local)
THEN
1193 CALL fun%pw_grid%para%group%sum(total_fun)
1196 IF (
PRESENT(isign))
THEN
1197 total_fun = total_fun*sign(1._dp, real(isign, dp))
1200 END FUNCTION pw_integrate_function_c1d_rs
1207 SUBROUTINE pw_set_value_c1d_rs (pw, value)
1208 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw
1209 REAL(kind=dp),
INTENT(IN) ::
value
1211 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_set_value'
1215 CALL timeset(routinen, handle)
1218 pw%array = cmplx(
value, 0.0_dp, kind=dp)
1221 CALL timestop(handle)
1223 END SUBROUTINE pw_set_value_c1d_rs
1231 SUBROUTINE pw_zero_c1d_gs (pw)
1233 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
1235 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_zero'
1239 CALL timeset(routinen, handle)
1242 pw%array = cmplx(0.0_dp, 0.0_dp, kind=dp)
1245 CALL timestop(handle)
1247 END SUBROUTINE pw_zero_c1d_gs
1256 SUBROUTINE pw_scale_c1d_gs (pw, a)
1258 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
1259 REAL(kind=dp),
INTENT(IN) :: a
1261 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_scale'
1265 CALL timeset(routinen, handle)
1268 pw%array = a*pw%array
1271 CALL timestop(handle)
1273 END SUBROUTINE pw_scale_c1d_gs
1285 SUBROUTINE pw_write_c1d_gs (pw, unit_nr)
1287 TYPE(pw_c1d_gs_type),
INTENT(in) :: pw
1288 INTEGER,
INTENT(in) :: unit_nr
1292 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
1294 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=c1d"
1295 IF (
ASSOCIATED(pw%array))
THEN
1296 WRITE (unit=unit_nr, fmt=
"(' array=<c(',i8,':',i8,')>')") &
1297 lbound(pw%array, 1), ubound(pw%array, 1)
1299 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
1302 END SUBROUTINE pw_write_c1d_gs
1311 FUNCTION pw_integrate_function_c1d_gs (fun, isign, oprt)
RESULT(total_fun)
1313 TYPE(pw_c1d_gs_type),
INTENT(IN) :: fun
1314 INTEGER,
INTENT(IN),
OPTIONAL :: isign
1315 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
1316 REAL(kind=dp) :: total_fun
1322 IF (
PRESENT(oprt))
THEN
1327 cpabort(
"Unknown operator")
1334 cpabort(
"Operator ABS not implemented")
1335 IF (fun%pw_grid%have_g0) total_fun = fun%pw_grid%vol* real(fun%array(1), kind=dp)
1337 IF (fun%pw_grid%para%mode /= pw_mode_local)
THEN
1338 CALL fun%pw_grid%para%group%sum(total_fun)
1341 IF (
PRESENT(isign))
THEN
1342 total_fun = total_fun*sign(1._dp, real(isign, dp))
1345 END FUNCTION pw_integrate_function_c1d_gs
1352 SUBROUTINE pw_set_value_c1d_gs (pw, value)
1353 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
1354 REAL(kind=dp),
INTENT(IN) ::
value
1356 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_set_value'
1360 CALL timeset(routinen, handle)
1363 pw%array = cmplx(
value, 0.0_dp, kind=dp)
1366 CALL timestop(handle)
1368 END SUBROUTINE pw_set_value_c1d_gs
1376 SUBROUTINE pw_gather_p_c1d (pw, c, scale)
1378 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
1379 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
INTENT(IN) :: c
1380 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
1382 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_gather_p'
1384 INTEGER :: gpt, handle, l, m, mn, n
1386 CALL timeset(routinen, handle)
1388 IF (pw%pw_grid%para%mode /= pw_mode_distributed)
THEN
1389 cpabort(
"This grid type is not distributed")
1392 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
1393 ngpts =>
SIZE(pw%pw_grid%gsq), ghat => pw%pw_grid%g_hat, yzq => pw%pw_grid%para%yzq)
1395 IF (
PRESENT(scale))
THEN
1400 l = mapl(ghat(1, gpt)) + 1
1401 m = mapm(ghat(2, gpt)) + 1
1402 n = mapn(ghat(3, gpt)) + 1
1404 pw%array(gpt) = scale* c(l, mn)
1412 l = mapl(ghat(1, gpt)) + 1
1413 m = mapm(ghat(2, gpt)) + 1
1414 n = mapn(ghat(3, gpt)) + 1
1416 pw%array(gpt) = c(l, mn)
1423 CALL timestop(handle)
1425 END SUBROUTINE pw_gather_p_c1d
1433 SUBROUTINE pw_scatter_p_c1d (pw, c, scale)
1434 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
1435 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
INTENT(INOUT) :: c
1436 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
1438 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_scatter_p'
1440 INTEGER :: gpt, handle, l, m, mn, n
1442 CALL timeset(routinen, handle)
1444 IF (pw%pw_grid%para%mode /= pw_mode_distributed)
THEN
1445 cpabort(
"This grid type is not distributed")
1448 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
1449 ghat => pw%pw_grid%g_hat, yzq => pw%pw_grid%para%yzq, ngpts =>
SIZE(pw%pw_grid%gsq))
1451 IF (.NOT.
PRESENT(scale)) c = z_zero
1453 IF (
PRESENT(scale))
THEN
1458 l = mapl(ghat(1, gpt)) + 1
1459 m = mapm(ghat(2, gpt)) + 1
1460 n = mapn(ghat(3, gpt)) + 1
1462 c(l, mn) = scale*pw%array(gpt)
1470 l = mapl(ghat(1, gpt)) + 1
1471 m = mapm(ghat(2, gpt)) + 1
1472 n = mapn(ghat(3, gpt)) + 1
1474 c(l, mn) = pw%array(gpt)
1481 IF (pw%pw_grid%grid_span == halfspace)
THEN
1483 associate(mapm => pw%pw_grid%mapm%neg, mapn => pw%pw_grid%mapn%neg, mapl => pw%pw_grid%mapl%neg, &
1484 ghat => pw%pw_grid%g_hat, ngpts =>
SIZE(pw%pw_grid%gsq), yzq => pw%pw_grid%para%yzq)
1486 IF (
PRESENT(scale))
THEN
1491 l = mapl(ghat(1, gpt)) + 1
1492 m = mapm(ghat(2, gpt)) + 1
1493 n = mapn(ghat(3, gpt)) + 1
1495 c(l, mn) = scale*conjg( pw%array(gpt))
1503 l = mapl(ghat(1, gpt)) + 1
1504 m = mapm(ghat(2, gpt)) + 1
1505 n = mapn(ghat(3, gpt)) + 1
1507 c(l, mn) = conjg( pw%array(gpt))
1514 CALL timestop(handle)
1516 END SUBROUTINE pw_scatter_p_c1d
1524 SUBROUTINE pw_zero_c3d_rs (pw)
1526 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw
1528 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_zero'
1532 CALL timeset(routinen, handle)
1535 pw%array = cmplx(0.0_dp, 0.0_dp, kind=dp)
1538 CALL timestop(handle)
1540 END SUBROUTINE pw_zero_c3d_rs
1549 SUBROUTINE pw_scale_c3d_rs (pw, a)
1551 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw
1552 REAL(kind=dp),
INTENT(IN) :: a
1554 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_scale'
1558 CALL timeset(routinen, handle)
1561 pw%array = a*pw%array
1564 CALL timestop(handle)
1566 END SUBROUTINE pw_scale_c3d_rs
1578 SUBROUTINE pw_write_c3d_rs (pw, unit_nr)
1580 TYPE(pw_c3d_rs_type),
INTENT(in) :: pw
1581 INTEGER,
INTENT(in) :: unit_nr
1585 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
1587 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=c3d"
1588 IF (
ASSOCIATED(pw%array))
THEN
1589 WRITE (unit=unit_nr, fmt=
"(' array=<c(',i8,':',i8,',',i8,':',i8,',',i8,':',i8,')>')") &
1590 lbound(pw%array, 1), ubound(pw%array, 1), lbound(pw%array, 2), ubound(pw%array, 2), &
1591 lbound(pw%array, 3), ubound(pw%array, 3)
1593 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
1596 END SUBROUTINE pw_write_c3d_rs
1605 FUNCTION pw_integrate_function_c3d_rs (fun, isign, oprt)
RESULT(total_fun)
1607 TYPE(pw_c3d_rs_type),
INTENT(IN) :: fun
1608 INTEGER,
INTENT(IN),
OPTIONAL :: isign
1609 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
1610 REAL(kind=dp) :: total_fun
1616 IF (
PRESENT(oprt))
THEN
1621 cpabort(
"Unknown operator")
1629 total_fun = fun%pw_grid%dvol*accurate_sum(abs(fun%array))
1631 total_fun = fun%pw_grid%dvol*accurate_sum( real(fun%array, kind=dp))
1634 IF (fun%pw_grid%para%mode /= pw_mode_local)
THEN
1635 CALL fun%pw_grid%para%group%sum(total_fun)
1638 IF (
PRESENT(isign))
THEN
1639 total_fun = total_fun*sign(1._dp, real(isign, dp))
1642 END FUNCTION pw_integrate_function_c3d_rs
1649 SUBROUTINE pw_set_value_c3d_rs (pw, value)
1650 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw
1651 REAL(kind=dp),
INTENT(IN) ::
value
1653 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_set_value'
1657 CALL timeset(routinen, handle)
1660 pw%array = cmplx(
value, 0.0_dp, kind=dp)
1663 CALL timestop(handle)
1665 END SUBROUTINE pw_set_value_c3d_rs
1673 SUBROUTINE pw_zero_c3d_gs (pw)
1675 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw
1677 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_zero'
1681 CALL timeset(routinen, handle)
1684 pw%array = cmplx(0.0_dp, 0.0_dp, kind=dp)
1687 CALL timestop(handle)
1689 END SUBROUTINE pw_zero_c3d_gs
1698 SUBROUTINE pw_scale_c3d_gs (pw, a)
1700 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw
1701 REAL(kind=dp),
INTENT(IN) :: a
1703 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_scale'
1707 CALL timeset(routinen, handle)
1710 pw%array = a*pw%array
1713 CALL timestop(handle)
1715 END SUBROUTINE pw_scale_c3d_gs
1727 SUBROUTINE pw_write_c3d_gs (pw, unit_nr)
1729 TYPE(pw_c3d_gs_type),
INTENT(in) :: pw
1730 INTEGER,
INTENT(in) :: unit_nr
1734 WRITE (unit=unit_nr, fmt=
"('<pw>:{ ')", iostat=iostatus)
1736 WRITE (unit=unit_nr, fmt=
"(a)", iostat=iostatus)
" in_use=c3d"
1737 IF (
ASSOCIATED(pw%array))
THEN
1738 WRITE (unit=unit_nr, fmt=
"(' array=<c(',i8,':',i8,',',i8,':',i8,',',i8,':',i8,')>')") &
1739 lbound(pw%array, 1), ubound(pw%array, 1), lbound(pw%array, 2), ubound(pw%array, 2), &
1740 lbound(pw%array, 3), ubound(pw%array, 3)
1742 WRITE (unit=unit_nr, fmt=
"(' array=*null*')")
1745 END SUBROUTINE pw_write_c3d_gs
1754 FUNCTION pw_integrate_function_c3d_gs (fun, isign, oprt)
RESULT(total_fun)
1756 TYPE(pw_c3d_gs_type),
INTENT(IN) :: fun
1757 INTEGER,
INTENT(IN),
OPTIONAL :: isign
1758 CHARACTER(len=*),
INTENT(IN),
OPTIONAL :: oprt
1759 REAL(kind=dp) :: total_fun
1765 IF (
PRESENT(oprt))
THEN
1770 cpabort(
"Unknown operator")
1777 cpabort(
"Operator ABS not implemented")
1778 cpabort(
"Reciprocal space integration for 3D grids not implemented!")
1780 IF (fun%pw_grid%para%mode /= pw_mode_local)
THEN
1781 CALL fun%pw_grid%para%group%sum(total_fun)
1784 IF (
PRESENT(isign))
THEN
1785 total_fun = total_fun*sign(1._dp, real(isign, dp))
1788 END FUNCTION pw_integrate_function_c3d_gs
1795 SUBROUTINE pw_set_value_c3d_gs (pw, value)
1796 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw
1797 REAL(kind=dp),
INTENT(IN) ::
value
1799 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_set_value'
1803 CALL timeset(routinen, handle)
1806 pw%array = cmplx(
value, 0.0_dp, kind=dp)
1809 CALL timestop(handle)
1811 END SUBROUTINE pw_set_value_c3d_gs
1827 SUBROUTINE pw_copy_r1d_r1d_rs (pw1, pw2)
1829 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
1830 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw2
1832 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy'
1835 INTEGER :: i, j, ng, ng1, ng2, ns
1837 CALL timeset(routinen, handle)
1839 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
1840 cpabort(
"Both grids must be either spherical or non-spherical!")
1841 IF (pw1%pw_grid%spherical) &
1842 cpabort(
"Spherical grids only exist in reciprocal space!")
1844 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
1845 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
1846 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
1847 ng1 =
SIZE(pw1%array)
1848 ng2 =
SIZE(pw2%array)
1851 pw2%array(1:ng) = pw1%array(1:ng)
1855 pw2%array(ng + 1:ng2) = 0.0_dp
1859 cpabort(
"Copies between spherical grids require compatible grids!")
1862 ng1 =
SIZE(pw1%array)
1863 ng2 =
SIZE(pw2%array)
1864 ns = 2*max(ng1, ng2)
1866 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
1867 IF (ng1 >= ng2)
THEN
1870 j = pw2%pw_grid%gidx(i)
1871 pw2%array(i) = pw1%array(j)
1878 j = pw2%pw_grid%gidx(i)
1879 pw2%array(j) = pw1%array(i)
1883 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
1884 IF (ng1 >= ng2)
THEN
1887 j = pw1%pw_grid%gidx(i)
1888 pw2%array(i) = pw1%array(j)
1895 j = pw1%pw_grid%gidx(i)
1896 pw2%array(j) = pw1%array(i)
1901 cpabort(
"Copy not implemented!")
1908 pw2%array = pw1%array
1912 CALL timestop(handle)
1914 END SUBROUTINE pw_copy_r1d_r1d_rs
1921 SUBROUTINE pw_copy_to_array_r1d_r1d_rs (pw, array)
1922 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw
1923 REAL(kind=dp),
DIMENSION(:),
INTENT(INOUT) :: array
1925 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_to_array'
1929 CALL timeset(routinen, handle)
1932 array(:) = pw%array(:)
1935 CALL timestop(handle)
1936 END SUBROUTINE pw_copy_to_array_r1d_r1d_rs
1943 SUBROUTINE pw_copy_from_array_r1d_r1d_rs (pw, array)
1944 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw
1945 REAL(kind=dp),
DIMENSION(:),
INTENT(IN) :: array
1947 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_from_array'
1951 CALL timeset(routinen, handle)
1957 CALL timestop(handle)
1958 END SUBROUTINE pw_copy_from_array_r1d_r1d_rs
1976 SUBROUTINE pw_axpy_r1d_r1d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
1978 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
1979 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw2
1980 REAL(kind=dp),
INTENT(in),
OPTIONAL :: alpha, beta
1981 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
1983 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_axpy'
1986 LOGICAL :: my_allow_noncompatible_grids
1987 REAL(kind=dp) :: my_alpha, my_beta
1988 INTEGER :: i, j, ng, ng1, ng2
1990 CALL timeset(routinen, handle)
1993 IF (
PRESENT(alpha)) my_alpha = alpha
1996 IF (
PRESENT(beta)) my_beta = beta
1998 my_allow_noncompatible_grids = .false.
1999 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
2001 IF (my_beta /= 1.0_dp)
THEN
2002 IF (my_beta == 0.0_dp)
THEN
2006 pw2%array = pw2%array*my_beta
2011 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
2013 IF (my_alpha == 1.0_dp)
THEN
2015 pw2%array = pw2%array + pw1%array
2017 ELSE IF (my_alpha /= 0.0_dp)
THEN
2019 pw2%array = pw2%array + my_alpha* pw1%array
2023 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
2025 ng1 =
SIZE(pw1%array)
2026 ng2 =
SIZE(pw2%array)
2029 IF (pw1%pw_grid%spherical)
THEN
2030 IF (my_alpha == 1.0_dp)
THEN
2033 pw2%array(i) = pw2%array(i) + pw1%array(i)
2036 ELSE IF (my_alpha /= 0.0_dp)
THEN
2039 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(i)
2043 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
2044 IF (ng1 >= ng2)
THEN
2045 IF (my_alpha == 1.0_dp)
THEN
2048 j = pw2%pw_grid%gidx(i)
2049 pw2%array(i) = pw2%array(i) + pw1%array(j)
2052 ELSE IF (my_alpha /= 0.0_dp)
THEN
2055 j = pw2%pw_grid%gidx(i)
2056 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
2061 IF (my_alpha == 1.0_dp)
THEN
2064 j = pw2%pw_grid%gidx(i)
2065 pw2%array(j) = pw2%array(j) + pw1%array(i)
2068 ELSE IF (my_alpha /= 0.0_dp)
THEN
2071 j = pw2%pw_grid%gidx(i)
2072 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
2077 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
2078 IF (ng1 >= ng2)
THEN
2079 IF (my_alpha == 1.0_dp)
THEN
2082 j = pw1%pw_grid%gidx(i)
2083 pw2%array(i) = pw2%array(i) + pw1%array(j)
2086 ELSE IF (my_alpha /= 0.0_dp)
THEN
2089 j = pw1%pw_grid%gidx(i)
2090 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
2095 IF (my_alpha == 1.0_dp)
THEN
2098 j = pw1%pw_grid%gidx(i)
2099 pw2%array(j) = pw2%array(j) + pw1%array(i)
2102 ELSE IF (my_alpha /= 0.0_dp)
THEN
2105 j = pw1%pw_grid%gidx(i)
2106 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
2112 cpabort(
"Grids not compatible")
2117 cpabort(
"Grids not compatible")
2121 CALL timestop(handle)
2123 END SUBROUTINE pw_axpy_r1d_r1d_rs
2134 SUBROUTINE pw_multiply_r1d_r1d_rs (pw_out, pw1, pw2, alpha)
2136 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw_out
2137 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
2138 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw2
2139 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: alpha
2141 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_multiply'
2144 REAL(kind=dp) :: my_alpha
2146 CALL timeset(routinen, handle)
2149 IF (
PRESENT(alpha)) my_alpha = alpha
2151 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
2152 cpabort(
"pw_multiply not implemented for non-identical grids!")
2154 IF (my_alpha == 1.0_dp)
THEN
2156 pw_out%array = pw_out%array + pw1%array* pw2%array
2160 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
2164 CALL timestop(handle)
2166 END SUBROUTINE pw_multiply_r1d_r1d_rs
2173 SUBROUTINE pw_multiply_with_r1d_r1d_rs (pw1, pw2)
2174 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw1
2175 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw2
2177 CHARACTER(LEN=*),
PARAMETER :: routinen =
'pw_multiply_with'
2181 CALL timeset(routinen, handle)
2183 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
2184 cpabort(
"Incompatible grids!")
2187 pw1%array = pw1%array* pw2%array
2190 CALL timestop(handle)
2192 END SUBROUTINE pw_multiply_with_r1d_r1d_rs
2207 FUNCTION pw_integral_ab_r1d_r1d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
2209 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
2210 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw2
2211 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
2212 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
2213 REAL(kind=dp) :: integral_value
2215 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab'
2217 INTEGER :: handle, loc_sumtype
2218 LOGICAL :: my_just_sum, my_local_only
2220 CALL timeset(routinen, handle)
2223 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
2225 my_local_only = .false.
2226 IF (
PRESENT(local_only)) my_local_only = local_only
2228 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
2229 cpabort(
"Grids incompatible")
2232 my_just_sum = .false.
2233 IF (
PRESENT(just_sum)) my_just_sum = just_sum
2240 integral_value = dot_product(pw1%array, pw2%array)
2245 integral_value = accurate_dot_product(pw1%array, pw2%array)
2249 IF (.NOT. my_just_sum)
THEN
2250 integral_value = integral_value*pw1%pw_grid%dvol
2253 IF (pw1%pw_grid%grid_span == halfspace)
THEN
2254 integral_value = 2.0_dp*integral_value
2255 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
2256 pw1%array(1)*pw2%array(1)
2259 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
2260 CALL pw1%pw_grid%para%group%sum(integral_value)
2262 CALL timestop(handle)
2264 END FUNCTION pw_integral_ab_r1d_r1d_rs
2278 SUBROUTINE pw_copy_r1d_r1d_gs (pw1, pw2)
2280 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
2281 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw2
2283 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy'
2286 INTEGER :: i, j, ng, ng1, ng2, ns
2288 CALL timeset(routinen, handle)
2290 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
2291 cpabort(
"Both grids must be either spherical or non-spherical!")
2293 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
2294 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
2295 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
2296 ng1 =
SIZE(pw1%array)
2297 ng2 =
SIZE(pw2%array)
2300 pw2%array(1:ng) = pw1%array(1:ng)
2304 pw2%array(ng + 1:ng2) = 0.0_dp
2308 cpabort(
"Copies between spherical grids require compatible grids!")
2311 ng1 =
SIZE(pw1%array)
2312 ng2 =
SIZE(pw2%array)
2313 ns = 2*max(ng1, ng2)
2315 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
2316 IF (ng1 >= ng2)
THEN
2319 j = pw2%pw_grid%gidx(i)
2320 pw2%array(i) = pw1%array(j)
2327 j = pw2%pw_grid%gidx(i)
2328 pw2%array(j) = pw1%array(i)
2332 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
2333 IF (ng1 >= ng2)
THEN
2336 j = pw1%pw_grid%gidx(i)
2337 pw2%array(i) = pw1%array(j)
2344 j = pw1%pw_grid%gidx(i)
2345 pw2%array(j) = pw1%array(i)
2350 cpabort(
"Copy not implemented!")
2357 pw2%array = pw1%array
2361 CALL timestop(handle)
2363 END SUBROUTINE pw_copy_r1d_r1d_gs
2370 SUBROUTINE pw_copy_to_array_r1d_r1d_gs (pw, array)
2371 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw
2372 REAL(kind=dp),
DIMENSION(:),
INTENT(INOUT) :: array
2374 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_to_array'
2378 CALL timeset(routinen, handle)
2381 array(:) = pw%array(:)
2384 CALL timestop(handle)
2385 END SUBROUTINE pw_copy_to_array_r1d_r1d_gs
2392 SUBROUTINE pw_copy_from_array_r1d_r1d_gs (pw, array)
2393 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw
2394 REAL(kind=dp),
DIMENSION(:),
INTENT(IN) :: array
2396 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_from_array'
2400 CALL timeset(routinen, handle)
2406 CALL timestop(handle)
2407 END SUBROUTINE pw_copy_from_array_r1d_r1d_gs
2425 SUBROUTINE pw_axpy_r1d_r1d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
2427 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
2428 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw2
2429 REAL(kind=dp),
INTENT(in),
OPTIONAL :: alpha, beta
2430 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
2432 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_axpy'
2435 LOGICAL :: my_allow_noncompatible_grids
2436 REAL(kind=dp) :: my_alpha, my_beta
2437 INTEGER :: i, j, ng, ng1, ng2
2439 CALL timeset(routinen, handle)
2442 IF (
PRESENT(alpha)) my_alpha = alpha
2445 IF (
PRESENT(beta)) my_beta = beta
2447 my_allow_noncompatible_grids = .false.
2448 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
2450 IF (my_beta /= 1.0_dp)
THEN
2451 IF (my_beta == 0.0_dp)
THEN
2455 pw2%array = pw2%array*my_beta
2460 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
2462 IF (my_alpha == 1.0_dp)
THEN
2464 pw2%array = pw2%array + pw1%array
2466 ELSE IF (my_alpha /= 0.0_dp)
THEN
2468 pw2%array = pw2%array + my_alpha* pw1%array
2472 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
2474 ng1 =
SIZE(pw1%array)
2475 ng2 =
SIZE(pw2%array)
2478 IF (pw1%pw_grid%spherical)
THEN
2479 IF (my_alpha == 1.0_dp)
THEN
2482 pw2%array(i) = pw2%array(i) + pw1%array(i)
2485 ELSE IF (my_alpha /= 0.0_dp)
THEN
2488 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(i)
2492 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
2493 IF (ng1 >= ng2)
THEN
2494 IF (my_alpha == 1.0_dp)
THEN
2497 j = pw2%pw_grid%gidx(i)
2498 pw2%array(i) = pw2%array(i) + pw1%array(j)
2501 ELSE IF (my_alpha /= 0.0_dp)
THEN
2504 j = pw2%pw_grid%gidx(i)
2505 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
2510 IF (my_alpha == 1.0_dp)
THEN
2513 j = pw2%pw_grid%gidx(i)
2514 pw2%array(j) = pw2%array(j) + pw1%array(i)
2517 ELSE IF (my_alpha /= 0.0_dp)
THEN
2520 j = pw2%pw_grid%gidx(i)
2521 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
2526 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
2527 IF (ng1 >= ng2)
THEN
2528 IF (my_alpha == 1.0_dp)
THEN
2531 j = pw1%pw_grid%gidx(i)
2532 pw2%array(i) = pw2%array(i) + pw1%array(j)
2535 ELSE IF (my_alpha /= 0.0_dp)
THEN
2538 j = pw1%pw_grid%gidx(i)
2539 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
2544 IF (my_alpha == 1.0_dp)
THEN
2547 j = pw1%pw_grid%gidx(i)
2548 pw2%array(j) = pw2%array(j) + pw1%array(i)
2551 ELSE IF (my_alpha /= 0.0_dp)
THEN
2554 j = pw1%pw_grid%gidx(i)
2555 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
2561 cpabort(
"Grids not compatible")
2566 cpabort(
"Grids not compatible")
2570 CALL timestop(handle)
2572 END SUBROUTINE pw_axpy_r1d_r1d_gs
2583 SUBROUTINE pw_multiply_r1d_r1d_gs (pw_out, pw1, pw2, alpha)
2585 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw_out
2586 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
2587 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
2588 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: alpha
2590 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_multiply'
2593 REAL(kind=dp) :: my_alpha
2595 CALL timeset(routinen, handle)
2598 IF (
PRESENT(alpha)) my_alpha = alpha
2600 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
2601 cpabort(
"pw_multiply not implemented for non-identical grids!")
2603 IF (my_alpha == 1.0_dp)
THEN
2605 pw_out%array = pw_out%array + pw1%array* pw2%array
2609 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
2613 CALL timestop(handle)
2615 END SUBROUTINE pw_multiply_r1d_r1d_gs
2622 SUBROUTINE pw_multiply_with_r1d_r1d_gs (pw1, pw2)
2623 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw1
2624 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
2626 CHARACTER(LEN=*),
PARAMETER :: routinen =
'pw_multiply_with'
2630 CALL timeset(routinen, handle)
2632 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
2633 cpabort(
"Incompatible grids!")
2636 pw1%array = pw1%array* pw2%array
2639 CALL timestop(handle)
2641 END SUBROUTINE pw_multiply_with_r1d_r1d_gs
2656 FUNCTION pw_integral_ab_r1d_r1d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
2658 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
2659 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
2660 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
2661 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
2662 REAL(kind=dp) :: integral_value
2664 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab'
2666 INTEGER :: handle, loc_sumtype
2667 LOGICAL :: my_just_sum, my_local_only
2669 CALL timeset(routinen, handle)
2672 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
2674 my_local_only = .false.
2675 IF (
PRESENT(local_only)) my_local_only = local_only
2677 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
2678 cpabort(
"Grids incompatible")
2681 my_just_sum = .false.
2682 IF (
PRESENT(just_sum)) my_just_sum = just_sum
2689 integral_value = dot_product(pw1%array, pw2%array)
2694 integral_value = accurate_dot_product(pw1%array, pw2%array)
2698 IF (.NOT. my_just_sum)
THEN
2699 integral_value = integral_value*pw1%pw_grid%vol
2702 IF (pw1%pw_grid%grid_span == halfspace)
THEN
2703 integral_value = 2.0_dp*integral_value
2704 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
2705 pw1%array(1)*pw2%array(1)
2708 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
2709 CALL pw1%pw_grid%para%group%sum(integral_value)
2711 CALL timestop(handle)
2713 END FUNCTION pw_integral_ab_r1d_r1d_gs
2721 FUNCTION pw_integral_a2b_r1d_r1d (pw1, pw2)
RESULT(integral_value)
2723 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
2724 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
2725 REAL(kind=dp) :: integral_value
2727 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_a2b'
2731 CALL timeset(routinen, handle)
2733 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
2734 cpabort(
"Grids incompatible")
2737 integral_value = accurate_sum(pw1%array*pw2%array*pw1%pw_grid%gsq)
2738 IF (pw1%pw_grid%grid_span == halfspace)
THEN
2739 integral_value = 2.0_dp*integral_value
2742 integral_value = integral_value*pw1%pw_grid%vol
2744 IF (pw1%pw_grid%para%mode == pw_mode_distributed) &
2745 CALL pw1%pw_grid%para%group%sum(integral_value)
2746 CALL timestop(handle)
2748 END FUNCTION pw_integral_a2b_r1d_r1d
2762 SUBROUTINE pw_copy_r1d_c1d_rs (pw1, pw2)
2764 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
2765 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw2
2767 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy'
2770 INTEGER :: i, j, ng, ng1, ng2, ns
2772 CALL timeset(routinen, handle)
2774 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
2775 cpabort(
"Both grids must be either spherical or non-spherical!")
2776 IF (pw1%pw_grid%spherical) &
2777 cpabort(
"Spherical grids only exist in reciprocal space!")
2779 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
2780 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
2781 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
2782 ng1 =
SIZE(pw1%array)
2783 ng2 =
SIZE(pw2%array)
2786 pw2%array(1:ng) = cmplx(pw1%array(1:ng), 0.0_dp, kind=dp)
2790 pw2%array(ng + 1:ng2) = cmplx(0.0_dp, 0.0_dp, kind=dp)
2794 cpabort(
"Copies between spherical grids require compatible grids!")
2797 ng1 =
SIZE(pw1%array)
2798 ng2 =
SIZE(pw2%array)
2799 ns = 2*max(ng1, ng2)
2801 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
2802 IF (ng1 >= ng2)
THEN
2805 j = pw2%pw_grid%gidx(i)
2806 pw2%array(i) = cmplx(pw1%array(j), 0.0_dp, kind=dp)
2813 j = pw2%pw_grid%gidx(i)
2814 pw2%array(j) = cmplx(pw1%array(i), 0.0_dp, kind=dp)
2818 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
2819 IF (ng1 >= ng2)
THEN
2822 j = pw1%pw_grid%gidx(i)
2823 pw2%array(i) = cmplx(pw1%array(j), 0.0_dp, kind=dp)
2830 j = pw1%pw_grid%gidx(i)
2831 pw2%array(j) = cmplx(pw1%array(i), 0.0_dp, kind=dp)
2836 cpabort(
"Copy not implemented!")
2843 pw2%array = cmplx(pw1%array, 0.0_dp, kind=dp)
2847 CALL timestop(handle)
2849 END SUBROUTINE pw_copy_r1d_c1d_rs
2856 SUBROUTINE pw_copy_to_array_r1d_c1d_rs (pw, array)
2857 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw
2858 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(INOUT) :: array
2860 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_to_array'
2864 CALL timeset(routinen, handle)
2867 array(:) = cmplx(pw%array(:), 0.0_dp, kind=dp)
2870 CALL timestop(handle)
2871 END SUBROUTINE pw_copy_to_array_r1d_c1d_rs
2878 SUBROUTINE pw_copy_from_array_r1d_c1d_rs (pw, array)
2879 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw
2880 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(IN) :: array
2882 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_from_array'
2886 CALL timeset(routinen, handle)
2889 pw%array = real(array, kind=dp)
2892 CALL timestop(handle)
2893 END SUBROUTINE pw_copy_from_array_r1d_c1d_rs
2911 SUBROUTINE pw_axpy_r1d_c1d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
2913 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
2914 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw2
2915 REAL(kind=dp),
INTENT(in),
OPTIONAL :: alpha, beta
2916 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
2918 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_axpy'
2921 LOGICAL :: my_allow_noncompatible_grids
2922 REAL(kind=dp) :: my_alpha, my_beta
2923 INTEGER :: i, j, ng, ng1, ng2
2925 CALL timeset(routinen, handle)
2928 IF (
PRESENT(alpha)) my_alpha = alpha
2931 IF (
PRESENT(beta)) my_beta = beta
2933 my_allow_noncompatible_grids = .false.
2934 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
2936 IF (my_beta /= 1.0_dp)
THEN
2937 IF (my_beta == 0.0_dp)
THEN
2941 pw2%array = pw2%array*my_beta
2946 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
2948 IF (my_alpha == 1.0_dp)
THEN
2950 pw2%array = pw2%array + cmplx(pw1%array, 0.0_dp, kind=dp)
2952 ELSE IF (my_alpha /= 0.0_dp)
THEN
2954 pw2%array = pw2%array + my_alpha* cmplx(pw1%array, 0.0_dp, kind=dp)
2958 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
2960 ng1 =
SIZE(pw1%array)
2961 ng2 =
SIZE(pw2%array)
2964 IF (pw1%pw_grid%spherical)
THEN
2965 IF (my_alpha == 1.0_dp)
THEN
2968 pw2%array(i) = pw2%array(i) + cmplx(pw1%array(i), 0.0_dp, kind=dp)
2971 ELSE IF (my_alpha /= 0.0_dp)
THEN
2974 pw2%array(i) = pw2%array(i) + my_alpha* cmplx(pw1%array(i), 0.0_dp, kind=dp)
2978 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
2979 IF (ng1 >= ng2)
THEN
2980 IF (my_alpha == 1.0_dp)
THEN
2983 j = pw2%pw_grid%gidx(i)
2984 pw2%array(i) = pw2%array(i) + cmplx(pw1%array(j), 0.0_dp, kind=dp)
2987 ELSE IF (my_alpha /= 0.0_dp)
THEN
2990 j = pw2%pw_grid%gidx(i)
2991 pw2%array(i) = pw2%array(i) + my_alpha* cmplx(pw1%array(j), 0.0_dp, kind=dp)
2996 IF (my_alpha == 1.0_dp)
THEN
2999 j = pw2%pw_grid%gidx(i)
3000 pw2%array(j) = pw2%array(j) + cmplx(pw1%array(i), 0.0_dp, kind=dp)
3003 ELSE IF (my_alpha /= 0.0_dp)
THEN
3006 j = pw2%pw_grid%gidx(i)
3007 pw2%array(j) = pw2%array(j) + my_alpha* cmplx(pw1%array(i), 0.0_dp, kind=dp)
3012 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
3013 IF (ng1 >= ng2)
THEN
3014 IF (my_alpha == 1.0_dp)
THEN
3017 j = pw1%pw_grid%gidx(i)
3018 pw2%array(i) = pw2%array(i) + cmplx(pw1%array(j), 0.0_dp, kind=dp)
3021 ELSE IF (my_alpha /= 0.0_dp)
THEN
3024 j = pw1%pw_grid%gidx(i)
3025 pw2%array(i) = pw2%array(i) + my_alpha* cmplx(pw1%array(j), 0.0_dp, kind=dp)
3030 IF (my_alpha == 1.0_dp)
THEN
3033 j = pw1%pw_grid%gidx(i)
3034 pw2%array(j) = pw2%array(j) + cmplx(pw1%array(i), 0.0_dp, kind=dp)
3037 ELSE IF (my_alpha /= 0.0_dp)
THEN
3040 j = pw1%pw_grid%gidx(i)
3041 pw2%array(j) = pw2%array(j) + my_alpha* cmplx(pw1%array(i), 0.0_dp, kind=dp)
3047 cpabort(
"Grids not compatible")
3052 cpabort(
"Grids not compatible")
3056 CALL timestop(handle)
3058 END SUBROUTINE pw_axpy_r1d_c1d_rs
3069 SUBROUTINE pw_multiply_r1d_c1d_rs (pw_out, pw1, pw2, alpha)
3071 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw_out
3072 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
3073 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw2
3074 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: alpha
3076 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_multiply'
3079 REAL(kind=dp) :: my_alpha
3081 CALL timeset(routinen, handle)
3084 IF (
PRESENT(alpha)) my_alpha = alpha
3086 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
3087 cpabort(
"pw_multiply not implemented for non-identical grids!")
3089 IF (my_alpha == 1.0_dp)
THEN
3091 pw_out%array = pw_out%array + pw1%array* real(pw2%array, kind=dp)
3095 pw_out%array = pw_out%array + my_alpha*pw1%array* real(pw2%array, kind=dp)
3099 CALL timestop(handle)
3101 END SUBROUTINE pw_multiply_r1d_c1d_rs
3108 SUBROUTINE pw_multiply_with_r1d_c1d_rs (pw1, pw2)
3109 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw1
3110 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw2
3112 CHARACTER(LEN=*),
PARAMETER :: routinen =
'pw_multiply_with'
3116 CALL timeset(routinen, handle)
3118 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
3119 cpabort(
"Incompatible grids!")
3122 pw1%array = pw1%array* real(pw2%array, kind=dp)
3125 CALL timestop(handle)
3127 END SUBROUTINE pw_multiply_with_r1d_c1d_rs
3142 FUNCTION pw_integral_ab_r1d_c1d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
3144 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw1
3145 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw2
3146 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
3147 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
3148 REAL(kind=dp) :: integral_value
3150 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab'
3152 INTEGER :: handle, loc_sumtype
3153 LOGICAL :: my_just_sum, my_local_only
3155 CALL timeset(routinen, handle)
3158 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
3160 my_local_only = .false.
3161 IF (
PRESENT(local_only)) my_local_only = local_only
3163 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
3164 cpabort(
"Grids incompatible")
3167 my_just_sum = .false.
3168 IF (
PRESENT(just_sum)) my_just_sum = just_sum
3175 integral_value = sum(pw1%array*real(pw2%array, kind=dp))
3180 integral_value = accurate_sum(pw1%array*real(pw2%array, kind=dp))
3184 IF (.NOT. my_just_sum)
THEN
3185 integral_value = integral_value*pw1%pw_grid%dvol
3188 IF (pw1%pw_grid%grid_span == halfspace)
THEN
3189 integral_value = 2.0_dp*integral_value
3190 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
3191 pw1%array(1)*real(pw2%array(1), kind=dp)
3194 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
3195 CALL pw1%pw_grid%para%group%sum(integral_value)
3197 CALL timestop(handle)
3199 END FUNCTION pw_integral_ab_r1d_c1d_rs
3213 SUBROUTINE pw_copy_r1d_c1d_gs (pw1, pw2)
3215 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
3216 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
3218 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy'
3221 INTEGER :: i, j, ng, ng1, ng2, ns
3223 CALL timeset(routinen, handle)
3225 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
3226 cpabort(
"Both grids must be either spherical or non-spherical!")
3228 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
3229 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
3230 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
3231 ng1 =
SIZE(pw1%array)
3232 ng2 =
SIZE(pw2%array)
3235 pw2%array(1:ng) = cmplx(pw1%array(1:ng), 0.0_dp, kind=dp)
3239 pw2%array(ng + 1:ng2) = cmplx(0.0_dp, 0.0_dp, kind=dp)
3243 cpabort(
"Copies between spherical grids require compatible grids!")
3246 ng1 =
SIZE(pw1%array)
3247 ng2 =
SIZE(pw2%array)
3248 ns = 2*max(ng1, ng2)
3250 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
3251 IF (ng1 >= ng2)
THEN
3254 j = pw2%pw_grid%gidx(i)
3255 pw2%array(i) = cmplx(pw1%array(j), 0.0_dp, kind=dp)
3262 j = pw2%pw_grid%gidx(i)
3263 pw2%array(j) = cmplx(pw1%array(i), 0.0_dp, kind=dp)
3267 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
3268 IF (ng1 >= ng2)
THEN
3271 j = pw1%pw_grid%gidx(i)
3272 pw2%array(i) = cmplx(pw1%array(j), 0.0_dp, kind=dp)
3279 j = pw1%pw_grid%gidx(i)
3280 pw2%array(j) = cmplx(pw1%array(i), 0.0_dp, kind=dp)
3285 cpabort(
"Copy not implemented!")
3292 pw2%array = cmplx(pw1%array, 0.0_dp, kind=dp)
3296 CALL timestop(handle)
3298 END SUBROUTINE pw_copy_r1d_c1d_gs
3305 SUBROUTINE pw_copy_to_array_r1d_c1d_gs (pw, array)
3306 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw
3307 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(INOUT) :: array
3309 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_to_array'
3313 CALL timeset(routinen, handle)
3316 array(:) = cmplx(pw%array(:), 0.0_dp, kind=dp)
3319 CALL timestop(handle)
3320 END SUBROUTINE pw_copy_to_array_r1d_c1d_gs
3327 SUBROUTINE pw_copy_from_array_r1d_c1d_gs (pw, array)
3328 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw
3329 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(IN) :: array
3331 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_from_array'
3335 CALL timeset(routinen, handle)
3338 pw%array = real(array, kind=dp)
3341 CALL timestop(handle)
3342 END SUBROUTINE pw_copy_from_array_r1d_c1d_gs
3360 SUBROUTINE pw_axpy_r1d_c1d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
3362 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
3363 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
3364 REAL(kind=dp),
INTENT(in),
OPTIONAL :: alpha, beta
3365 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
3367 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_axpy'
3370 LOGICAL :: my_allow_noncompatible_grids
3371 REAL(kind=dp) :: my_alpha, my_beta
3372 INTEGER :: i, j, ng, ng1, ng2
3374 CALL timeset(routinen, handle)
3377 IF (
PRESENT(alpha)) my_alpha = alpha
3380 IF (
PRESENT(beta)) my_beta = beta
3382 my_allow_noncompatible_grids = .false.
3383 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
3385 IF (my_beta /= 1.0_dp)
THEN
3386 IF (my_beta == 0.0_dp)
THEN
3390 pw2%array = pw2%array*my_beta
3395 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
3397 IF (my_alpha == 1.0_dp)
THEN
3399 pw2%array = pw2%array + cmplx(pw1%array, 0.0_dp, kind=dp)
3401 ELSE IF (my_alpha /= 0.0_dp)
THEN
3403 pw2%array = pw2%array + my_alpha* cmplx(pw1%array, 0.0_dp, kind=dp)
3407 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
3409 ng1 =
SIZE(pw1%array)
3410 ng2 =
SIZE(pw2%array)
3413 IF (pw1%pw_grid%spherical)
THEN
3414 IF (my_alpha == 1.0_dp)
THEN
3417 pw2%array(i) = pw2%array(i) + cmplx(pw1%array(i), 0.0_dp, kind=dp)
3420 ELSE IF (my_alpha /= 0.0_dp)
THEN
3423 pw2%array(i) = pw2%array(i) + my_alpha* cmplx(pw1%array(i), 0.0_dp, kind=dp)
3427 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
3428 IF (ng1 >= ng2)
THEN
3429 IF (my_alpha == 1.0_dp)
THEN
3432 j = pw2%pw_grid%gidx(i)
3433 pw2%array(i) = pw2%array(i) + cmplx(pw1%array(j), 0.0_dp, kind=dp)
3436 ELSE IF (my_alpha /= 0.0_dp)
THEN
3439 j = pw2%pw_grid%gidx(i)
3440 pw2%array(i) = pw2%array(i) + my_alpha* cmplx(pw1%array(j), 0.0_dp, kind=dp)
3445 IF (my_alpha == 1.0_dp)
THEN
3448 j = pw2%pw_grid%gidx(i)
3449 pw2%array(j) = pw2%array(j) + cmplx(pw1%array(i), 0.0_dp, kind=dp)
3452 ELSE IF (my_alpha /= 0.0_dp)
THEN
3455 j = pw2%pw_grid%gidx(i)
3456 pw2%array(j) = pw2%array(j) + my_alpha* cmplx(pw1%array(i), 0.0_dp, kind=dp)
3461 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
3462 IF (ng1 >= ng2)
THEN
3463 IF (my_alpha == 1.0_dp)
THEN
3466 j = pw1%pw_grid%gidx(i)
3467 pw2%array(i) = pw2%array(i) + cmplx(pw1%array(j), 0.0_dp, kind=dp)
3470 ELSE IF (my_alpha /= 0.0_dp)
THEN
3473 j = pw1%pw_grid%gidx(i)
3474 pw2%array(i) = pw2%array(i) + my_alpha* cmplx(pw1%array(j), 0.0_dp, kind=dp)
3479 IF (my_alpha == 1.0_dp)
THEN
3482 j = pw1%pw_grid%gidx(i)
3483 pw2%array(j) = pw2%array(j) + cmplx(pw1%array(i), 0.0_dp, kind=dp)
3486 ELSE IF (my_alpha /= 0.0_dp)
THEN
3489 j = pw1%pw_grid%gidx(i)
3490 pw2%array(j) = pw2%array(j) + my_alpha* cmplx(pw1%array(i), 0.0_dp, kind=dp)
3496 cpabort(
"Grids not compatible")
3501 cpabort(
"Grids not compatible")
3505 CALL timestop(handle)
3507 END SUBROUTINE pw_axpy_r1d_c1d_gs
3518 SUBROUTINE pw_multiply_r1d_c1d_gs (pw_out, pw1, pw2, alpha)
3520 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw_out
3521 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
3522 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
3523 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: alpha
3525 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_multiply'
3528 REAL(kind=dp) :: my_alpha
3530 CALL timeset(routinen, handle)
3533 IF (
PRESENT(alpha)) my_alpha = alpha
3535 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
3536 cpabort(
"pw_multiply not implemented for non-identical grids!")
3538 IF (my_alpha == 1.0_dp)
THEN
3540 pw_out%array = pw_out%array + pw1%array* real(pw2%array, kind=dp)
3544 pw_out%array = pw_out%array + my_alpha*pw1%array* real(pw2%array, kind=dp)
3548 CALL timestop(handle)
3550 END SUBROUTINE pw_multiply_r1d_c1d_gs
3557 SUBROUTINE pw_multiply_with_r1d_c1d_gs (pw1, pw2)
3558 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw1
3559 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
3561 CHARACTER(LEN=*),
PARAMETER :: routinen =
'pw_multiply_with'
3565 CALL timeset(routinen, handle)
3567 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
3568 cpabort(
"Incompatible grids!")
3571 pw1%array = pw1%array* real(pw2%array, kind=dp)
3574 CALL timestop(handle)
3576 END SUBROUTINE pw_multiply_with_r1d_c1d_gs
3591 FUNCTION pw_integral_ab_r1d_c1d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
3593 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
3594 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
3595 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
3596 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
3597 REAL(kind=dp) :: integral_value
3599 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab'
3601 INTEGER :: handle, loc_sumtype
3602 LOGICAL :: my_just_sum, my_local_only
3604 CALL timeset(routinen, handle)
3607 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
3609 my_local_only = .false.
3610 IF (
PRESENT(local_only)) my_local_only = local_only
3612 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
3613 cpabort(
"Grids incompatible")
3616 my_just_sum = .false.
3617 IF (
PRESENT(just_sum)) my_just_sum = just_sum
3624 integral_value = sum(pw1%array*real(pw2%array, kind=dp))
3629 integral_value = accurate_sum(pw1%array*real(pw2%array, kind=dp))
3633 IF (.NOT. my_just_sum)
THEN
3634 integral_value = integral_value*pw1%pw_grid%vol
3637 IF (pw1%pw_grid%grid_span == halfspace)
THEN
3638 integral_value = 2.0_dp*integral_value
3639 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
3640 pw1%array(1)*real(pw2%array(1), kind=dp)
3643 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
3644 CALL pw1%pw_grid%para%group%sum(integral_value)
3646 CALL timestop(handle)
3648 END FUNCTION pw_integral_ab_r1d_c1d_gs
3656 FUNCTION pw_integral_a2b_r1d_c1d (pw1, pw2)
RESULT(integral_value)
3658 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
3659 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
3660 REAL(kind=dp) :: integral_value
3662 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_a2b'
3666 CALL timeset(routinen, handle)
3668 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
3669 cpabort(
"Grids incompatible")
3672 integral_value = accurate_sum(pw1%array*real(pw2%array, kind=dp)*pw1%pw_grid%gsq)
3673 IF (pw1%pw_grid%grid_span == halfspace)
THEN
3674 integral_value = 2.0_dp*integral_value
3677 integral_value = integral_value*pw1%pw_grid%vol
3679 IF (pw1%pw_grid%para%mode == pw_mode_distributed) &
3680 CALL pw1%pw_grid%para%group%sum(integral_value)
3681 CALL timestop(handle)
3683 END FUNCTION pw_integral_a2b_r1d_c1d
3697 SUBROUTINE pw_copy_r3d_r3d_rs (pw1, pw2)
3699 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
3700 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw2
3702 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy'
3706 CALL timeset(routinen, handle)
3708 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
3709 cpabort(
"Both grids must be either spherical or non-spherical!")
3710 IF (pw1%pw_grid%spherical) &
3711 cpabort(
"Spherical grids only exist in reciprocal space!")
3713 IF (any(shape(pw2%array) /= shape(pw1%array))) &
3714 cpabort(
"3D grids must be compatible!")
3715 IF (pw1%pw_grid%spherical) &
3716 cpabort(
"3D grids must not be spherical!")
3718 pw2%array(:, :, :) = pw1%array(:, :, :)
3721 CALL timestop(handle)
3723 END SUBROUTINE pw_copy_r3d_r3d_rs
3730 SUBROUTINE pw_copy_to_array_r3d_r3d_rs (pw, array)
3731 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw
3732 REAL(kind=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
3734 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_to_array'
3738 CALL timeset(routinen, handle)
3741 array(:, :, :) = pw%array(:, :, :)
3744 CALL timestop(handle)
3745 END SUBROUTINE pw_copy_to_array_r3d_r3d_rs
3752 SUBROUTINE pw_copy_from_array_r3d_r3d_rs (pw, array)
3753 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw
3754 REAL(kind=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
3756 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_from_array'
3760 CALL timeset(routinen, handle)
3766 CALL timestop(handle)
3767 END SUBROUTINE pw_copy_from_array_r3d_r3d_rs
3785 SUBROUTINE pw_axpy_r3d_r3d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
3787 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
3788 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw2
3789 REAL(kind=dp),
INTENT(in),
OPTIONAL :: alpha, beta
3790 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
3792 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_axpy'
3795 LOGICAL :: my_allow_noncompatible_grids
3796 REAL(kind=dp) :: my_alpha, my_beta
3798 CALL timeset(routinen, handle)
3801 IF (
PRESENT(alpha)) my_alpha = alpha
3804 IF (
PRESENT(beta)) my_beta = beta
3806 my_allow_noncompatible_grids = .false.
3807 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
3809 IF (my_beta /= 1.0_dp)
THEN
3810 IF (my_beta == 0.0_dp)
THEN
3814 pw2%array = pw2%array*my_beta
3819 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
3820 IF (my_alpha == 1.0_dp)
THEN
3822 pw2%array = pw2%array + pw1%array
3824 ELSE IF (my_alpha /= 0.0_dp)
THEN
3826 pw2%array = pw2%array + my_alpha* pw1%array
3830 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
3832 IF (any(shape(pw1%array) /= shape(pw2%array))) &
3833 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
3835 IF (my_alpha == 1.0_dp)
THEN
3837 pw2%array = pw2%array + pw1%array
3841 pw2%array = pw2%array + my_alpha* pw1%array
3847 cpabort(
"Grids not compatible")
3851 CALL timestop(handle)
3853 END SUBROUTINE pw_axpy_r3d_r3d_rs
3864 SUBROUTINE pw_multiply_r3d_r3d_rs (pw_out, pw1, pw2, alpha)
3866 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw_out
3867 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
3868 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw2
3869 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: alpha
3871 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_multiply'
3874 REAL(kind=dp) :: my_alpha
3876 CALL timeset(routinen, handle)
3879 IF (
PRESENT(alpha)) my_alpha = alpha
3881 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
3882 cpabort(
"pw_multiply not implemented for non-identical grids!")
3884 IF (my_alpha == 1.0_dp)
THEN
3886 pw_out%array = pw_out%array + pw1%array* pw2%array
3890 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
3894 CALL timestop(handle)
3896 END SUBROUTINE pw_multiply_r3d_r3d_rs
3903 SUBROUTINE pw_multiply_with_r3d_r3d_rs (pw1, pw2)
3904 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw1
3905 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw2
3907 CHARACTER(LEN=*),
PARAMETER :: routinen =
'pw_multiply_with'
3911 CALL timeset(routinen, handle)
3913 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
3914 cpabort(
"Incompatible grids!")
3917 pw1%array = pw1%array* pw2%array
3920 CALL timestop(handle)
3922 END SUBROUTINE pw_multiply_with_r3d_r3d_rs
3937 FUNCTION pw_integral_ab_r3d_r3d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
3939 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
3940 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw2
3941 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
3942 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
3943 REAL(kind=dp) :: integral_value
3945 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab'
3947 INTEGER :: handle, loc_sumtype
3948 LOGICAL :: my_just_sum, my_local_only
3950 CALL timeset(routinen, handle)
3953 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
3955 my_local_only = .false.
3956 IF (
PRESENT(local_only)) my_local_only = local_only
3958 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
3959 cpabort(
"Grids incompatible")
3962 my_just_sum = .false.
3963 IF (
PRESENT(just_sum)) my_just_sum = just_sum
3970 integral_value = sum(pw1%array*pw2%array)
3975 integral_value = accurate_dot_product(pw1%array, pw2%array)
3979 IF (.NOT. my_just_sum)
THEN
3980 integral_value = integral_value*pw1%pw_grid%dvol
3984 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
3985 CALL pw1%pw_grid%para%group%sum(integral_value)
3987 CALL timestop(handle)
3989 END FUNCTION pw_integral_ab_r3d_r3d_rs
4003 SUBROUTINE pw_copy_r3d_r3d_gs (pw1, pw2)
4005 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4006 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw2
4008 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy'
4012 CALL timeset(routinen, handle)
4014 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
4015 cpabort(
"Both grids must be either spherical or non-spherical!")
4017 IF (any(shape(pw2%array) /= shape(pw1%array))) &
4018 cpabort(
"3D grids must be compatible!")
4019 IF (pw1%pw_grid%spherical) &
4020 cpabort(
"3D grids must not be spherical!")
4022 pw2%array(:, :, :) = pw1%array(:, :, :)
4025 CALL timestop(handle)
4027 END SUBROUTINE pw_copy_r3d_r3d_gs
4034 SUBROUTINE pw_copy_to_array_r3d_r3d_gs (pw, array)
4035 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw
4036 REAL(kind=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
4038 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_to_array'
4042 CALL timeset(routinen, handle)
4045 array(:, :, :) = pw%array(:, :, :)
4048 CALL timestop(handle)
4049 END SUBROUTINE pw_copy_to_array_r3d_r3d_gs
4056 SUBROUTINE pw_copy_from_array_r3d_r3d_gs (pw, array)
4057 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw
4058 REAL(kind=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
4060 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_from_array'
4064 CALL timeset(routinen, handle)
4070 CALL timestop(handle)
4071 END SUBROUTINE pw_copy_from_array_r3d_r3d_gs
4089 SUBROUTINE pw_axpy_r3d_r3d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
4091 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4092 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw2
4093 REAL(kind=dp),
INTENT(in),
OPTIONAL :: alpha, beta
4094 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
4096 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_axpy'
4099 LOGICAL :: my_allow_noncompatible_grids
4100 REAL(kind=dp) :: my_alpha, my_beta
4102 CALL timeset(routinen, handle)
4105 IF (
PRESENT(alpha)) my_alpha = alpha
4108 IF (
PRESENT(beta)) my_beta = beta
4110 my_allow_noncompatible_grids = .false.
4111 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
4113 IF (my_beta /= 1.0_dp)
THEN
4114 IF (my_beta == 0.0_dp)
THEN
4118 pw2%array = pw2%array*my_beta
4123 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
4124 IF (my_alpha == 1.0_dp)
THEN
4126 pw2%array = pw2%array + pw1%array
4128 ELSE IF (my_alpha /= 0.0_dp)
THEN
4130 pw2%array = pw2%array + my_alpha* pw1%array
4134 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
4136 IF (any(shape(pw1%array) /= shape(pw2%array))) &
4137 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
4139 IF (my_alpha == 1.0_dp)
THEN
4141 pw2%array = pw2%array + pw1%array
4145 pw2%array = pw2%array + my_alpha* pw1%array
4151 cpabort(
"Grids not compatible")
4155 CALL timestop(handle)
4157 END SUBROUTINE pw_axpy_r3d_r3d_gs
4168 SUBROUTINE pw_multiply_r3d_r3d_gs (pw_out, pw1, pw2, alpha)
4170 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw_out
4171 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4172 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw2
4173 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: alpha
4175 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_multiply'
4178 REAL(kind=dp) :: my_alpha
4180 CALL timeset(routinen, handle)
4183 IF (
PRESENT(alpha)) my_alpha = alpha
4185 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
4186 cpabort(
"pw_multiply not implemented for non-identical grids!")
4188 IF (my_alpha == 1.0_dp)
THEN
4190 pw_out%array = pw_out%array + pw1%array* pw2%array
4194 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
4198 CALL timestop(handle)
4200 END SUBROUTINE pw_multiply_r3d_r3d_gs
4207 SUBROUTINE pw_multiply_with_r3d_r3d_gs (pw1, pw2)
4208 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw1
4209 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw2
4211 CHARACTER(LEN=*),
PARAMETER :: routinen =
'pw_multiply_with'
4215 CALL timeset(routinen, handle)
4217 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
4218 cpabort(
"Incompatible grids!")
4221 pw1%array = pw1%array* pw2%array
4224 CALL timestop(handle)
4226 END SUBROUTINE pw_multiply_with_r3d_r3d_gs
4241 FUNCTION pw_integral_ab_r3d_r3d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
4243 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4244 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw2
4245 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
4246 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
4247 REAL(kind=dp) :: integral_value
4249 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab'
4251 INTEGER :: handle, loc_sumtype
4252 LOGICAL :: my_just_sum, my_local_only
4254 CALL timeset(routinen, handle)
4257 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
4259 my_local_only = .false.
4260 IF (
PRESENT(local_only)) my_local_only = local_only
4262 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
4263 cpabort(
"Grids incompatible")
4266 my_just_sum = .false.
4267 IF (
PRESENT(just_sum)) my_just_sum = just_sum
4274 integral_value = sum(pw1%array*pw2%array)
4279 integral_value = accurate_dot_product(pw1%array, pw2%array)
4283 IF (.NOT. my_just_sum)
THEN
4284 integral_value = integral_value*pw1%pw_grid%vol
4288 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
4289 CALL pw1%pw_grid%para%group%sum(integral_value)
4291 CALL timestop(handle)
4293 END FUNCTION pw_integral_ab_r3d_r3d_gs
4308 SUBROUTINE pw_copy_r3d_c3d_rs (pw1, pw2)
4310 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
4311 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw2
4313 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy'
4317 CALL timeset(routinen, handle)
4319 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
4320 cpabort(
"Both grids must be either spherical or non-spherical!")
4321 IF (pw1%pw_grid%spherical) &
4322 cpabort(
"Spherical grids only exist in reciprocal space!")
4324 IF (any(shape(pw2%array) /= shape(pw1%array))) &
4325 cpabort(
"3D grids must be compatible!")
4326 IF (pw1%pw_grid%spherical) &
4327 cpabort(
"3D grids must not be spherical!")
4329 pw2%array(:, :, :) = cmplx(pw1%array(:, :, :), 0.0_dp, kind=dp)
4332 CALL timestop(handle)
4334 END SUBROUTINE pw_copy_r3d_c3d_rs
4341 SUBROUTINE pw_copy_to_array_r3d_c3d_rs (pw, array)
4342 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw
4343 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
4345 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_to_array'
4349 CALL timeset(routinen, handle)
4352 array(:, :, :) = cmplx(pw%array(:, :, :), 0.0_dp, kind=dp)
4355 CALL timestop(handle)
4356 END SUBROUTINE pw_copy_to_array_r3d_c3d_rs
4363 SUBROUTINE pw_copy_from_array_r3d_c3d_rs (pw, array)
4364 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw
4365 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
4367 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_from_array'
4371 CALL timeset(routinen, handle)
4374 pw%array = real(array, kind=dp)
4377 CALL timestop(handle)
4378 END SUBROUTINE pw_copy_from_array_r3d_c3d_rs
4396 SUBROUTINE pw_axpy_r3d_c3d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
4398 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
4399 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw2
4400 REAL(kind=dp),
INTENT(in),
OPTIONAL :: alpha, beta
4401 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
4403 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_axpy'
4406 LOGICAL :: my_allow_noncompatible_grids
4407 REAL(kind=dp) :: my_alpha, my_beta
4409 CALL timeset(routinen, handle)
4412 IF (
PRESENT(alpha)) my_alpha = alpha
4415 IF (
PRESENT(beta)) my_beta = beta
4417 my_allow_noncompatible_grids = .false.
4418 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
4420 IF (my_beta /= 1.0_dp)
THEN
4421 IF (my_beta == 0.0_dp)
THEN
4425 pw2%array = pw2%array*my_beta
4430 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
4431 IF (my_alpha == 1.0_dp)
THEN
4433 pw2%array = pw2%array + cmplx(pw1%array, 0.0_dp, kind=dp)
4435 ELSE IF (my_alpha /= 0.0_dp)
THEN
4437 pw2%array = pw2%array + my_alpha* cmplx(pw1%array, 0.0_dp, kind=dp)
4441 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
4443 IF (any(shape(pw1%array) /= shape(pw2%array))) &
4444 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
4446 IF (my_alpha == 1.0_dp)
THEN
4448 pw2%array = pw2%array + cmplx(pw1%array, 0.0_dp, kind=dp)
4452 pw2%array = pw2%array + my_alpha* cmplx(pw1%array, 0.0_dp, kind=dp)
4458 cpabort(
"Grids not compatible")
4462 CALL timestop(handle)
4464 END SUBROUTINE pw_axpy_r3d_c3d_rs
4475 SUBROUTINE pw_multiply_r3d_c3d_rs (pw_out, pw1, pw2, alpha)
4477 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw_out
4478 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
4479 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw2
4480 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: alpha
4482 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_multiply'
4485 REAL(kind=dp) :: my_alpha
4487 CALL timeset(routinen, handle)
4490 IF (
PRESENT(alpha)) my_alpha = alpha
4492 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
4493 cpabort(
"pw_multiply not implemented for non-identical grids!")
4495 IF (my_alpha == 1.0_dp)
THEN
4497 pw_out%array = pw_out%array + pw1%array* real(pw2%array, kind=dp)
4501 pw_out%array = pw_out%array + my_alpha*pw1%array* real(pw2%array, kind=dp)
4505 CALL timestop(handle)
4507 END SUBROUTINE pw_multiply_r3d_c3d_rs
4514 SUBROUTINE pw_multiply_with_r3d_c3d_rs (pw1, pw2)
4515 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw1
4516 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw2
4518 CHARACTER(LEN=*),
PARAMETER :: routinen =
'pw_multiply_with'
4522 CALL timeset(routinen, handle)
4524 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
4525 cpabort(
"Incompatible grids!")
4528 pw1%array = pw1%array* real(pw2%array, kind=dp)
4531 CALL timestop(handle)
4533 END SUBROUTINE pw_multiply_with_r3d_c3d_rs
4548 FUNCTION pw_integral_ab_r3d_c3d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
4550 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
4551 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw2
4552 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
4553 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
4554 REAL(kind=dp) :: integral_value
4556 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab'
4558 INTEGER :: handle, loc_sumtype
4559 LOGICAL :: my_just_sum, my_local_only
4561 CALL timeset(routinen, handle)
4564 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
4566 my_local_only = .false.
4567 IF (
PRESENT(local_only)) my_local_only = local_only
4569 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
4570 cpabort(
"Grids incompatible")
4573 my_just_sum = .false.
4574 IF (
PRESENT(just_sum)) my_just_sum = just_sum
4581 integral_value = sum(pw1%array*real(pw2%array, kind=dp))
4586 integral_value = accurate_sum(pw1%array*real(pw2%array, kind=dp))
4590 IF (.NOT. my_just_sum)
THEN
4591 integral_value = integral_value*pw1%pw_grid%dvol
4595 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
4596 CALL pw1%pw_grid%para%group%sum(integral_value)
4598 CALL timestop(handle)
4600 END FUNCTION pw_integral_ab_r3d_c3d_rs
4614 SUBROUTINE pw_copy_r3d_c3d_gs (pw1, pw2)
4616 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4617 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
4619 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy'
4623 CALL timeset(routinen, handle)
4625 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
4626 cpabort(
"Both grids must be either spherical or non-spherical!")
4628 IF (any(shape(pw2%array) /= shape(pw1%array))) &
4629 cpabort(
"3D grids must be compatible!")
4630 IF (pw1%pw_grid%spherical) &
4631 cpabort(
"3D grids must not be spherical!")
4633 pw2%array(:, :, :) = cmplx(pw1%array(:, :, :), 0.0_dp, kind=dp)
4636 CALL timestop(handle)
4638 END SUBROUTINE pw_copy_r3d_c3d_gs
4645 SUBROUTINE pw_copy_to_array_r3d_c3d_gs (pw, array)
4646 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw
4647 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
4649 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_to_array'
4653 CALL timeset(routinen, handle)
4656 array(:, :, :) = cmplx(pw%array(:, :, :), 0.0_dp, kind=dp)
4659 CALL timestop(handle)
4660 END SUBROUTINE pw_copy_to_array_r3d_c3d_gs
4667 SUBROUTINE pw_copy_from_array_r3d_c3d_gs (pw, array)
4668 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw
4669 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
4671 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_from_array'
4675 CALL timeset(routinen, handle)
4678 pw%array = real(array, kind=dp)
4681 CALL timestop(handle)
4682 END SUBROUTINE pw_copy_from_array_r3d_c3d_gs
4700 SUBROUTINE pw_axpy_r3d_c3d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
4702 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4703 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
4704 REAL(kind=dp),
INTENT(in),
OPTIONAL :: alpha, beta
4705 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
4707 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_axpy'
4710 LOGICAL :: my_allow_noncompatible_grids
4711 REAL(kind=dp) :: my_alpha, my_beta
4713 CALL timeset(routinen, handle)
4716 IF (
PRESENT(alpha)) my_alpha = alpha
4719 IF (
PRESENT(beta)) my_beta = beta
4721 my_allow_noncompatible_grids = .false.
4722 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
4724 IF (my_beta /= 1.0_dp)
THEN
4725 IF (my_beta == 0.0_dp)
THEN
4729 pw2%array = pw2%array*my_beta
4734 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
4735 IF (my_alpha == 1.0_dp)
THEN
4737 pw2%array = pw2%array + cmplx(pw1%array, 0.0_dp, kind=dp)
4739 ELSE IF (my_alpha /= 0.0_dp)
THEN
4741 pw2%array = pw2%array + my_alpha* cmplx(pw1%array, 0.0_dp, kind=dp)
4745 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
4747 IF (any(shape(pw1%array) /= shape(pw2%array))) &
4748 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
4750 IF (my_alpha == 1.0_dp)
THEN
4752 pw2%array = pw2%array + cmplx(pw1%array, 0.0_dp, kind=dp)
4756 pw2%array = pw2%array + my_alpha* cmplx(pw1%array, 0.0_dp, kind=dp)
4762 cpabort(
"Grids not compatible")
4766 CALL timestop(handle)
4768 END SUBROUTINE pw_axpy_r3d_c3d_gs
4779 SUBROUTINE pw_multiply_r3d_c3d_gs (pw_out, pw1, pw2, alpha)
4781 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw_out
4782 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4783 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw2
4784 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: alpha
4786 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_multiply'
4789 REAL(kind=dp) :: my_alpha
4791 CALL timeset(routinen, handle)
4794 IF (
PRESENT(alpha)) my_alpha = alpha
4796 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
4797 cpabort(
"pw_multiply not implemented for non-identical grids!")
4799 IF (my_alpha == 1.0_dp)
THEN
4801 pw_out%array = pw_out%array + pw1%array* real(pw2%array, kind=dp)
4805 pw_out%array = pw_out%array + my_alpha*pw1%array* real(pw2%array, kind=dp)
4809 CALL timestop(handle)
4811 END SUBROUTINE pw_multiply_r3d_c3d_gs
4818 SUBROUTINE pw_multiply_with_r3d_c3d_gs (pw1, pw2)
4819 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw1
4820 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw2
4822 CHARACTER(LEN=*),
PARAMETER :: routinen =
'pw_multiply_with'
4826 CALL timeset(routinen, handle)
4828 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
4829 cpabort(
"Incompatible grids!")
4832 pw1%array = pw1%array* real(pw2%array, kind=dp)
4835 CALL timestop(handle)
4837 END SUBROUTINE pw_multiply_with_r3d_c3d_gs
4852 FUNCTION pw_integral_ab_r3d_c3d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
4854 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
4855 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw2
4856 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
4857 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
4858 REAL(kind=dp) :: integral_value
4860 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab'
4862 INTEGER :: handle, loc_sumtype
4863 LOGICAL :: my_just_sum, my_local_only
4865 CALL timeset(routinen, handle)
4868 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
4870 my_local_only = .false.
4871 IF (
PRESENT(local_only)) my_local_only = local_only
4873 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
4874 cpabort(
"Grids incompatible")
4877 my_just_sum = .false.
4878 IF (
PRESENT(just_sum)) my_just_sum = just_sum
4885 integral_value = sum(pw1%array*real(pw2%array, kind=dp))
4890 integral_value = accurate_sum(pw1%array*real(pw2%array, kind=dp))
4894 IF (.NOT. my_just_sum)
THEN
4895 integral_value = integral_value*pw1%pw_grid%vol
4899 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
4900 CALL pw1%pw_grid%para%group%sum(integral_value)
4902 CALL timestop(handle)
4904 END FUNCTION pw_integral_ab_r3d_c3d_gs
4919 SUBROUTINE pw_copy_c1d_r1d_rs (pw1, pw2)
4921 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
4922 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw2
4924 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy'
4927 INTEGER :: i, j, ng, ng1, ng2, ns
4929 CALL timeset(routinen, handle)
4931 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
4932 cpabort(
"Both grids must be either spherical or non-spherical!")
4933 IF (pw1%pw_grid%spherical) &
4934 cpabort(
"Spherical grids only exist in reciprocal space!")
4936 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
4937 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
4938 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
4939 ng1 =
SIZE(pw1%array)
4940 ng2 =
SIZE(pw2%array)
4943 pw2%array(1:ng) = real(pw1%array(1:ng), kind=dp)
4947 pw2%array(ng + 1:ng2) = 0.0_dp
4951 cpabort(
"Copies between spherical grids require compatible grids!")
4954 ng1 =
SIZE(pw1%array)
4955 ng2 =
SIZE(pw2%array)
4956 ns = 2*max(ng1, ng2)
4958 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
4959 IF (ng1 >= ng2)
THEN
4962 j = pw2%pw_grid%gidx(i)
4963 pw2%array(i) = real(pw1%array(j), kind=dp)
4970 j = pw2%pw_grid%gidx(i)
4971 pw2%array(j) = real(pw1%array(i), kind=dp)
4975 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
4976 IF (ng1 >= ng2)
THEN
4979 j = pw1%pw_grid%gidx(i)
4980 pw2%array(i) = real(pw1%array(j), kind=dp)
4987 j = pw1%pw_grid%gidx(i)
4988 pw2%array(j) = real(pw1%array(i), kind=dp)
4993 cpabort(
"Copy not implemented!")
5000 pw2%array = real(pw1%array, kind=dp)
5004 CALL timestop(handle)
5006 END SUBROUTINE pw_copy_c1d_r1d_rs
5013 SUBROUTINE pw_copy_to_array_c1d_r1d_rs (pw, array)
5014 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw
5015 REAL(kind=dp),
DIMENSION(:),
INTENT(INOUT) :: array
5017 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_to_array'
5021 CALL timeset(routinen, handle)
5024 array(:) = real(pw%array(:), kind=dp)
5027 CALL timestop(handle)
5028 END SUBROUTINE pw_copy_to_array_c1d_r1d_rs
5035 SUBROUTINE pw_copy_from_array_c1d_r1d_rs (pw, array)
5036 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw
5037 REAL(kind=dp),
DIMENSION(:),
INTENT(IN) :: array
5039 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_from_array'
5043 CALL timeset(routinen, handle)
5046 pw%array = cmplx(array, 0.0_dp, kind=dp)
5049 CALL timestop(handle)
5050 END SUBROUTINE pw_copy_from_array_c1d_r1d_rs
5068 SUBROUTINE pw_axpy_c1d_r1d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
5070 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
5071 TYPE(pw_r1d_rs_type),
INTENT(INOUT) :: pw2
5072 REAL(kind=dp),
INTENT(in),
OPTIONAL :: alpha, beta
5073 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
5075 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_axpy'
5078 LOGICAL :: my_allow_noncompatible_grids
5079 REAL(kind=dp) :: my_alpha, my_beta
5080 INTEGER :: i, j, ng, ng1, ng2
5082 CALL timeset(routinen, handle)
5085 IF (
PRESENT(alpha)) my_alpha = alpha
5088 IF (
PRESENT(beta)) my_beta = beta
5090 my_allow_noncompatible_grids = .false.
5091 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
5093 IF (my_beta /= 1.0_dp)
THEN
5094 IF (my_beta == 0.0_dp)
THEN
5098 pw2%array = pw2%array*my_beta
5103 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
5105 IF (my_alpha == 1.0_dp)
THEN
5107 pw2%array = pw2%array + real(pw1%array, kind=dp)
5109 ELSE IF (my_alpha /= 0.0_dp)
THEN
5111 pw2%array = pw2%array + my_alpha* real(pw1%array, kind=dp)
5115 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
5117 ng1 =
SIZE(pw1%array)
5118 ng2 =
SIZE(pw2%array)
5121 IF (pw1%pw_grid%spherical)
THEN
5122 IF (my_alpha == 1.0_dp)
THEN
5125 pw2%array(i) = pw2%array(i) + real(pw1%array(i), kind=dp)
5128 ELSE IF (my_alpha /= 0.0_dp)
THEN
5131 pw2%array(i) = pw2%array(i) + my_alpha* real(pw1%array(i), kind=dp)
5135 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
5136 IF (ng1 >= ng2)
THEN
5137 IF (my_alpha == 1.0_dp)
THEN
5140 j = pw2%pw_grid%gidx(i)
5141 pw2%array(i) = pw2%array(i) + real(pw1%array(j), kind=dp)
5144 ELSE IF (my_alpha /= 0.0_dp)
THEN
5147 j = pw2%pw_grid%gidx(i)
5148 pw2%array(i) = pw2%array(i) + my_alpha* real(pw1%array(j), kind=dp)
5153 IF (my_alpha == 1.0_dp)
THEN
5156 j = pw2%pw_grid%gidx(i)
5157 pw2%array(j) = pw2%array(j) + real(pw1%array(i), kind=dp)
5160 ELSE IF (my_alpha /= 0.0_dp)
THEN
5163 j = pw2%pw_grid%gidx(i)
5164 pw2%array(j) = pw2%array(j) + my_alpha* real(pw1%array(i), kind=dp)
5169 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
5170 IF (ng1 >= ng2)
THEN
5171 IF (my_alpha == 1.0_dp)
THEN
5174 j = pw1%pw_grid%gidx(i)
5175 pw2%array(i) = pw2%array(i) + real(pw1%array(j), kind=dp)
5178 ELSE IF (my_alpha /= 0.0_dp)
THEN
5181 j = pw1%pw_grid%gidx(i)
5182 pw2%array(i) = pw2%array(i) + my_alpha* real(pw1%array(j), kind=dp)
5187 IF (my_alpha == 1.0_dp)
THEN
5190 j = pw1%pw_grid%gidx(i)
5191 pw2%array(j) = pw2%array(j) + real(pw1%array(i), kind=dp)
5194 ELSE IF (my_alpha /= 0.0_dp)
THEN
5197 j = pw1%pw_grid%gidx(i)
5198 pw2%array(j) = pw2%array(j) + my_alpha* real(pw1%array(i), kind=dp)
5204 cpabort(
"Grids not compatible")
5209 cpabort(
"Grids not compatible")
5213 CALL timestop(handle)
5215 END SUBROUTINE pw_axpy_c1d_r1d_rs
5226 SUBROUTINE pw_multiply_c1d_r1d_rs (pw_out, pw1, pw2, alpha)
5228 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw_out
5229 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
5230 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw2
5231 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: alpha
5233 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_multiply'
5236 REAL(kind=dp) :: my_alpha
5238 CALL timeset(routinen, handle)
5241 IF (
PRESENT(alpha)) my_alpha = alpha
5243 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
5244 cpabort(
"pw_multiply not implemented for non-identical grids!")
5246 IF (my_alpha == 1.0_dp)
THEN
5248 pw_out%array = pw_out%array + pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5252 pw_out%array = pw_out%array + my_alpha*pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5256 CALL timestop(handle)
5258 END SUBROUTINE pw_multiply_c1d_r1d_rs
5265 SUBROUTINE pw_multiply_with_c1d_r1d_rs (pw1, pw2)
5266 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw1
5267 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw2
5269 CHARACTER(LEN=*),
PARAMETER :: routinen =
'pw_multiply_with'
5273 CALL timeset(routinen, handle)
5275 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
5276 cpabort(
"Incompatible grids!")
5279 pw1%array = pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5282 CALL timestop(handle)
5284 END SUBROUTINE pw_multiply_with_c1d_r1d_rs
5299 FUNCTION pw_integral_ab_c1d_r1d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
5301 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
5302 TYPE(pw_r1d_rs_type),
INTENT(IN) :: pw2
5303 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
5304 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
5305 REAL(kind=dp) :: integral_value
5307 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab'
5309 INTEGER :: handle, loc_sumtype
5310 LOGICAL :: my_just_sum, my_local_only
5312 CALL timeset(routinen, handle)
5315 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
5317 my_local_only = .false.
5318 IF (
PRESENT(local_only)) my_local_only = local_only
5320 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
5321 cpabort(
"Grids incompatible")
5324 my_just_sum = .false.
5325 IF (
PRESENT(just_sum)) my_just_sum = just_sum
5332 integral_value = sum(real(pw1%array, kind=dp)*pw2%array)
5337 integral_value = accurate_sum(real(pw1%array, kind=dp)*pw2%array)
5341 IF (.NOT. my_just_sum)
THEN
5342 integral_value = integral_value*pw1%pw_grid%dvol
5345 IF (pw1%pw_grid%grid_span == halfspace)
THEN
5346 integral_value = 2.0_dp*integral_value
5347 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
5348 REAL(conjg(pw1%array(1))*pw2%array(1), kind=dp)
5351 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
5352 CALL pw1%pw_grid%para%group%sum(integral_value)
5354 CALL timestop(handle)
5356 END FUNCTION pw_integral_ab_c1d_r1d_rs
5370 SUBROUTINE pw_copy_c1d_r1d_gs (pw1, pw2)
5372 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
5373 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw2
5375 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy'
5378 INTEGER :: i, j, ng, ng1, ng2, ns
5380 CALL timeset(routinen, handle)
5382 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
5383 cpabort(
"Both grids must be either spherical or non-spherical!")
5385 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
5386 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
5387 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
5388 ng1 =
SIZE(pw1%array)
5389 ng2 =
SIZE(pw2%array)
5392 pw2%array(1:ng) = real(pw1%array(1:ng), kind=dp)
5396 pw2%array(ng + 1:ng2) = 0.0_dp
5400 cpabort(
"Copies between spherical grids require compatible grids!")
5403 ng1 =
SIZE(pw1%array)
5404 ng2 =
SIZE(pw2%array)
5405 ns = 2*max(ng1, ng2)
5407 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
5408 IF (ng1 >= ng2)
THEN
5411 j = pw2%pw_grid%gidx(i)
5412 pw2%array(i) = real(pw1%array(j), kind=dp)
5419 j = pw2%pw_grid%gidx(i)
5420 pw2%array(j) = real(pw1%array(i), kind=dp)
5424 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
5425 IF (ng1 >= ng2)
THEN
5428 j = pw1%pw_grid%gidx(i)
5429 pw2%array(i) = real(pw1%array(j), kind=dp)
5436 j = pw1%pw_grid%gidx(i)
5437 pw2%array(j) = real(pw1%array(i), kind=dp)
5442 cpabort(
"Copy not implemented!")
5449 pw2%array = real(pw1%array, kind=dp)
5453 CALL timestop(handle)
5455 END SUBROUTINE pw_copy_c1d_r1d_gs
5462 SUBROUTINE pw_copy_to_array_c1d_r1d_gs (pw, array)
5463 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
5464 REAL(kind=dp),
DIMENSION(:),
INTENT(INOUT) :: array
5466 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_to_array'
5470 CALL timeset(routinen, handle)
5473 array(:) = real(pw%array(:), kind=dp)
5476 CALL timestop(handle)
5477 END SUBROUTINE pw_copy_to_array_c1d_r1d_gs
5484 SUBROUTINE pw_copy_from_array_c1d_r1d_gs (pw, array)
5485 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
5486 REAL(kind=dp),
DIMENSION(:),
INTENT(IN) :: array
5488 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_from_array'
5492 CALL timeset(routinen, handle)
5495 pw%array = cmplx(array, 0.0_dp, kind=dp)
5498 CALL timestop(handle)
5499 END SUBROUTINE pw_copy_from_array_c1d_r1d_gs
5517 SUBROUTINE pw_axpy_c1d_r1d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
5519 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
5520 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw2
5521 REAL(kind=dp),
INTENT(in),
OPTIONAL :: alpha, beta
5522 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
5524 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_axpy'
5527 LOGICAL :: my_allow_noncompatible_grids
5528 REAL(kind=dp) :: my_alpha, my_beta
5529 INTEGER :: i, j, ng, ng1, ng2
5531 CALL timeset(routinen, handle)
5534 IF (
PRESENT(alpha)) my_alpha = alpha
5537 IF (
PRESENT(beta)) my_beta = beta
5539 my_allow_noncompatible_grids = .false.
5540 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
5542 IF (my_beta /= 1.0_dp)
THEN
5543 IF (my_beta == 0.0_dp)
THEN
5547 pw2%array = pw2%array*my_beta
5552 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
5554 IF (my_alpha == 1.0_dp)
THEN
5556 pw2%array = pw2%array + real(pw1%array, kind=dp)
5558 ELSE IF (my_alpha /= 0.0_dp)
THEN
5560 pw2%array = pw2%array + my_alpha* real(pw1%array, kind=dp)
5564 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
5566 ng1 =
SIZE(pw1%array)
5567 ng2 =
SIZE(pw2%array)
5570 IF (pw1%pw_grid%spherical)
THEN
5571 IF (my_alpha == 1.0_dp)
THEN
5574 pw2%array(i) = pw2%array(i) + real(pw1%array(i), kind=dp)
5577 ELSE IF (my_alpha /= 0.0_dp)
THEN
5580 pw2%array(i) = pw2%array(i) + my_alpha* real(pw1%array(i), kind=dp)
5584 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
5585 IF (ng1 >= ng2)
THEN
5586 IF (my_alpha == 1.0_dp)
THEN
5589 j = pw2%pw_grid%gidx(i)
5590 pw2%array(i) = pw2%array(i) + real(pw1%array(j), kind=dp)
5593 ELSE IF (my_alpha /= 0.0_dp)
THEN
5596 j = pw2%pw_grid%gidx(i)
5597 pw2%array(i) = pw2%array(i) + my_alpha* real(pw1%array(j), kind=dp)
5602 IF (my_alpha == 1.0_dp)
THEN
5605 j = pw2%pw_grid%gidx(i)
5606 pw2%array(j) = pw2%array(j) + real(pw1%array(i), kind=dp)
5609 ELSE IF (my_alpha /= 0.0_dp)
THEN
5612 j = pw2%pw_grid%gidx(i)
5613 pw2%array(j) = pw2%array(j) + my_alpha* real(pw1%array(i), kind=dp)
5618 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
5619 IF (ng1 >= ng2)
THEN
5620 IF (my_alpha == 1.0_dp)
THEN
5623 j = pw1%pw_grid%gidx(i)
5624 pw2%array(i) = pw2%array(i) + real(pw1%array(j), kind=dp)
5627 ELSE IF (my_alpha /= 0.0_dp)
THEN
5630 j = pw1%pw_grid%gidx(i)
5631 pw2%array(i) = pw2%array(i) + my_alpha* real(pw1%array(j), kind=dp)
5636 IF (my_alpha == 1.0_dp)
THEN
5639 j = pw1%pw_grid%gidx(i)
5640 pw2%array(j) = pw2%array(j) + real(pw1%array(i), kind=dp)
5643 ELSE IF (my_alpha /= 0.0_dp)
THEN
5646 j = pw1%pw_grid%gidx(i)
5647 pw2%array(j) = pw2%array(j) + my_alpha* real(pw1%array(i), kind=dp)
5653 cpabort(
"Grids not compatible")
5658 cpabort(
"Grids not compatible")
5662 CALL timestop(handle)
5664 END SUBROUTINE pw_axpy_c1d_r1d_gs
5675 SUBROUTINE pw_multiply_c1d_r1d_gs (pw_out, pw1, pw2, alpha)
5677 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw_out
5678 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
5679 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
5680 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: alpha
5682 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_multiply'
5685 REAL(kind=dp) :: my_alpha
5687 CALL timeset(routinen, handle)
5690 IF (
PRESENT(alpha)) my_alpha = alpha
5692 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
5693 cpabort(
"pw_multiply not implemented for non-identical grids!")
5695 IF (my_alpha == 1.0_dp)
THEN
5697 pw_out%array = pw_out%array + pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5701 pw_out%array = pw_out%array + my_alpha*pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5705 CALL timestop(handle)
5707 END SUBROUTINE pw_multiply_c1d_r1d_gs
5714 SUBROUTINE pw_multiply_with_c1d_r1d_gs (pw1, pw2)
5715 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw1
5716 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
5718 CHARACTER(LEN=*),
PARAMETER :: routinen =
'pw_multiply_with'
5722 CALL timeset(routinen, handle)
5724 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
5725 cpabort(
"Incompatible grids!")
5728 pw1%array = pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
5731 CALL timestop(handle)
5733 END SUBROUTINE pw_multiply_with_c1d_r1d_gs
5748 FUNCTION pw_integral_ab_c1d_r1d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
5750 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
5751 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
5752 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
5753 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
5754 REAL(kind=dp) :: integral_value
5756 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab'
5758 INTEGER :: handle, loc_sumtype
5759 LOGICAL :: my_just_sum, my_local_only
5761 CALL timeset(routinen, handle)
5764 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
5766 my_local_only = .false.
5767 IF (
PRESENT(local_only)) my_local_only = local_only
5769 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
5770 cpabort(
"Grids incompatible")
5773 my_just_sum = .false.
5774 IF (
PRESENT(just_sum)) my_just_sum = just_sum
5781 integral_value = sum(real(pw1%array, kind=dp)*pw2%array)
5786 integral_value = accurate_sum(real(pw1%array, kind=dp)*pw2%array)
5790 IF (.NOT. my_just_sum)
THEN
5791 integral_value = integral_value*pw1%pw_grid%vol
5794 IF (pw1%pw_grid%grid_span == halfspace)
THEN
5795 integral_value = 2.0_dp*integral_value
5796 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
5797 REAL(conjg(pw1%array(1))*pw2%array(1), kind=dp)
5800 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
5801 CALL pw1%pw_grid%para%group%sum(integral_value)
5803 CALL timestop(handle)
5805 END FUNCTION pw_integral_ab_c1d_r1d_gs
5813 FUNCTION pw_integral_a2b_c1d_r1d (pw1, pw2)
RESULT(integral_value)
5815 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
5816 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw2
5817 REAL(kind=dp) :: integral_value
5819 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_a2b'
5823 CALL timeset(routinen, handle)
5825 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
5826 cpabort(
"Grids incompatible")
5829 integral_value = accurate_sum(real(conjg(pw1%array), kind=dp)*pw2%array*pw1%pw_grid%gsq)
5830 IF (pw1%pw_grid%grid_span == halfspace)
THEN
5831 integral_value = 2.0_dp*integral_value
5834 integral_value = integral_value*pw1%pw_grid%vol
5836 IF (pw1%pw_grid%para%mode == pw_mode_distributed) &
5837 CALL pw1%pw_grid%para%group%sum(integral_value)
5838 CALL timestop(handle)
5840 END FUNCTION pw_integral_a2b_c1d_r1d
5854 SUBROUTINE pw_copy_c1d_c1d_rs (pw1, pw2)
5856 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
5857 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw2
5859 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy'
5862 INTEGER :: i, j, ng, ng1, ng2, ns
5864 CALL timeset(routinen, handle)
5866 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
5867 cpabort(
"Both grids must be either spherical or non-spherical!")
5868 IF (pw1%pw_grid%spherical) &
5869 cpabort(
"Spherical grids only exist in reciprocal space!")
5871 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
5872 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
5873 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
5874 ng1 =
SIZE(pw1%array)
5875 ng2 =
SIZE(pw2%array)
5878 pw2%array(1:ng) = pw1%array(1:ng)
5882 pw2%array(ng + 1:ng2) = cmplx(0.0_dp, 0.0_dp, kind=dp)
5886 cpabort(
"Copies between spherical grids require compatible grids!")
5889 ng1 =
SIZE(pw1%array)
5890 ng2 =
SIZE(pw2%array)
5891 ns = 2*max(ng1, ng2)
5893 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
5894 IF (ng1 >= ng2)
THEN
5897 j = pw2%pw_grid%gidx(i)
5898 pw2%array(i) = pw1%array(j)
5905 j = pw2%pw_grid%gidx(i)
5906 pw2%array(j) = pw1%array(i)
5910 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
5911 IF (ng1 >= ng2)
THEN
5914 j = pw1%pw_grid%gidx(i)
5915 pw2%array(i) = pw1%array(j)
5922 j = pw1%pw_grid%gidx(i)
5923 pw2%array(j) = pw1%array(i)
5928 cpabort(
"Copy not implemented!")
5935 pw2%array = pw1%array
5939 CALL timestop(handle)
5941 END SUBROUTINE pw_copy_c1d_c1d_rs
5948 SUBROUTINE pw_copy_to_array_c1d_c1d_rs (pw, array)
5949 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw
5950 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(INOUT) :: array
5952 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_to_array'
5956 CALL timeset(routinen, handle)
5959 array(:) = pw%array(:)
5962 CALL timestop(handle)
5963 END SUBROUTINE pw_copy_to_array_c1d_c1d_rs
5970 SUBROUTINE pw_copy_from_array_c1d_c1d_rs (pw, array)
5971 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw
5972 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(IN) :: array
5974 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_from_array'
5978 CALL timeset(routinen, handle)
5984 CALL timestop(handle)
5985 END SUBROUTINE pw_copy_from_array_c1d_c1d_rs
6003 SUBROUTINE pw_axpy_c1d_c1d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
6005 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
6006 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw2
6007 REAL(kind=dp),
INTENT(in),
OPTIONAL :: alpha, beta
6008 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
6010 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_axpy'
6013 LOGICAL :: my_allow_noncompatible_grids
6014 REAL(kind=dp) :: my_alpha, my_beta
6015 INTEGER :: i, j, ng, ng1, ng2
6017 CALL timeset(routinen, handle)
6020 IF (
PRESENT(alpha)) my_alpha = alpha
6023 IF (
PRESENT(beta)) my_beta = beta
6025 my_allow_noncompatible_grids = .false.
6026 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
6028 IF (my_beta /= 1.0_dp)
THEN
6029 IF (my_beta == 0.0_dp)
THEN
6033 pw2%array = pw2%array*my_beta
6038 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
6040 IF (my_alpha == 1.0_dp)
THEN
6042 pw2%array = pw2%array + pw1%array
6044 ELSE IF (my_alpha /= 0.0_dp)
THEN
6046 pw2%array = pw2%array + my_alpha* pw1%array
6050 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
6052 ng1 =
SIZE(pw1%array)
6053 ng2 =
SIZE(pw2%array)
6056 IF (pw1%pw_grid%spherical)
THEN
6057 IF (my_alpha == 1.0_dp)
THEN
6060 pw2%array(i) = pw2%array(i) + pw1%array(i)
6063 ELSE IF (my_alpha /= 0.0_dp)
THEN
6066 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(i)
6070 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
6071 IF (ng1 >= ng2)
THEN
6072 IF (my_alpha == 1.0_dp)
THEN
6075 j = pw2%pw_grid%gidx(i)
6076 pw2%array(i) = pw2%array(i) + pw1%array(j)
6079 ELSE IF (my_alpha /= 0.0_dp)
THEN
6082 j = pw2%pw_grid%gidx(i)
6083 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
6088 IF (my_alpha == 1.0_dp)
THEN
6091 j = pw2%pw_grid%gidx(i)
6092 pw2%array(j) = pw2%array(j) + pw1%array(i)
6095 ELSE IF (my_alpha /= 0.0_dp)
THEN
6098 j = pw2%pw_grid%gidx(i)
6099 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
6104 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
6105 IF (ng1 >= ng2)
THEN
6106 IF (my_alpha == 1.0_dp)
THEN
6109 j = pw1%pw_grid%gidx(i)
6110 pw2%array(i) = pw2%array(i) + pw1%array(j)
6113 ELSE IF (my_alpha /= 0.0_dp)
THEN
6116 j = pw1%pw_grid%gidx(i)
6117 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
6122 IF (my_alpha == 1.0_dp)
THEN
6125 j = pw1%pw_grid%gidx(i)
6126 pw2%array(j) = pw2%array(j) + pw1%array(i)
6129 ELSE IF (my_alpha /= 0.0_dp)
THEN
6132 j = pw1%pw_grid%gidx(i)
6133 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
6139 cpabort(
"Grids not compatible")
6144 cpabort(
"Grids not compatible")
6148 CALL timestop(handle)
6150 END SUBROUTINE pw_axpy_c1d_c1d_rs
6161 SUBROUTINE pw_multiply_c1d_c1d_rs (pw_out, pw1, pw2, alpha)
6163 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw_out
6164 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
6165 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw2
6166 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: alpha
6168 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_multiply'
6171 REAL(kind=dp) :: my_alpha
6173 CALL timeset(routinen, handle)
6176 IF (
PRESENT(alpha)) my_alpha = alpha
6178 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
6179 cpabort(
"pw_multiply not implemented for non-identical grids!")
6181 IF (my_alpha == 1.0_dp)
THEN
6183 pw_out%array = pw_out%array + pw1%array* pw2%array
6187 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
6191 CALL timestop(handle)
6193 END SUBROUTINE pw_multiply_c1d_c1d_rs
6200 SUBROUTINE pw_multiply_with_c1d_c1d_rs (pw1, pw2)
6201 TYPE(pw_c1d_rs_type),
INTENT(INOUT) :: pw1
6202 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw2
6204 CHARACTER(LEN=*),
PARAMETER :: routinen =
'pw_multiply_with'
6208 CALL timeset(routinen, handle)
6210 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
6211 cpabort(
"Incompatible grids!")
6214 pw1%array = pw1%array* pw2%array
6217 CALL timestop(handle)
6219 END SUBROUTINE pw_multiply_with_c1d_c1d_rs
6234 FUNCTION pw_integral_ab_c1d_c1d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
6236 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw1
6237 TYPE(pw_c1d_rs_type),
INTENT(IN) :: pw2
6238 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
6239 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
6240 REAL(kind=dp) :: integral_value
6242 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab'
6244 INTEGER :: handle, loc_sumtype
6245 LOGICAL :: my_just_sum, my_local_only
6247 CALL timeset(routinen, handle)
6250 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
6252 my_local_only = .false.
6253 IF (
PRESENT(local_only)) my_local_only = local_only
6255 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
6256 cpabort(
"Grids incompatible")
6259 my_just_sum = .false.
6260 IF (
PRESENT(just_sum)) my_just_sum = just_sum
6267 integral_value = sum(real(conjg(pw1%array) &
6268 *pw2%array, kind=dp))
6273 integral_value = accurate_sum(real(conjg(pw1%array) &
6274 *pw2%array, kind=dp))
6278 IF (.NOT. my_just_sum)
THEN
6279 integral_value = integral_value*pw1%pw_grid%dvol
6282 IF (pw1%pw_grid%grid_span == halfspace)
THEN
6283 integral_value = 2.0_dp*integral_value
6284 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
6285 REAL(conjg(pw1%array(1))*pw2%array(1), kind=dp)
6288 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
6289 CALL pw1%pw_grid%para%group%sum(integral_value)
6291 CALL timestop(handle)
6293 END FUNCTION pw_integral_ab_c1d_c1d_rs
6307 SUBROUTINE pw_copy_c1d_c1d_gs (pw1, pw2)
6309 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
6310 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
6312 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy'
6315 INTEGER :: i, j, ng, ng1, ng2, ns
6317 CALL timeset(routinen, handle)
6319 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
6320 cpabort(
"Both grids must be either spherical or non-spherical!")
6322 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
6323 IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical)
THEN
6324 IF (pw_compatible(pw1%pw_grid, pw2%pw_grid))
THEN
6325 ng1 =
SIZE(pw1%array)
6326 ng2 =
SIZE(pw2%array)
6329 pw2%array(1:ng) = pw1%array(1:ng)
6333 pw2%array(ng + 1:ng2) = cmplx(0.0_dp, 0.0_dp, kind=dp)
6337 cpabort(
"Copies between spherical grids require compatible grids!")
6340 ng1 =
SIZE(pw1%array)
6341 ng2 =
SIZE(pw2%array)
6342 ns = 2*max(ng1, ng2)
6344 IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
6345 IF (ng1 >= ng2)
THEN
6348 j = pw2%pw_grid%gidx(i)
6349 pw2%array(i) = pw1%array(j)
6356 j = pw2%pw_grid%gidx(i)
6357 pw2%array(j) = pw1%array(i)
6361 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
6362 IF (ng1 >= ng2)
THEN
6365 j = pw1%pw_grid%gidx(i)
6366 pw2%array(i) = pw1%array(j)
6373 j = pw1%pw_grid%gidx(i)
6374 pw2%array(j) = pw1%array(i)
6379 CALL pw_copy_match(pw1, pw2)
6386 pw2%array = pw1%array
6390 CALL timestop(handle)
6392 END SUBROUTINE pw_copy_c1d_c1d_gs
6399 SUBROUTINE pw_copy_to_array_c1d_c1d_gs (pw, array)
6400 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
6401 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(INOUT) :: array
6403 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_to_array'
6407 CALL timeset(routinen, handle)
6410 array(:) = pw%array(:)
6413 CALL timestop(handle)
6414 END SUBROUTINE pw_copy_to_array_c1d_c1d_gs
6421 SUBROUTINE pw_copy_from_array_c1d_c1d_gs (pw, array)
6422 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
6423 COMPLEX(KIND=dp),
DIMENSION(:),
INTENT(IN) :: array
6425 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_from_array'
6429 CALL timeset(routinen, handle)
6435 CALL timestop(handle)
6436 END SUBROUTINE pw_copy_from_array_c1d_c1d_gs
6454 SUBROUTINE pw_axpy_c1d_c1d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
6456 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
6457 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
6458 REAL(kind=dp),
INTENT(in),
OPTIONAL :: alpha, beta
6459 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
6461 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_axpy'
6464 LOGICAL :: my_allow_noncompatible_grids
6465 REAL(kind=dp) :: my_alpha, my_beta
6466 INTEGER :: i, j, ng, ng1, ng2
6468 CALL timeset(routinen, handle)
6471 IF (
PRESENT(alpha)) my_alpha = alpha
6474 IF (
PRESENT(beta)) my_beta = beta
6476 my_allow_noncompatible_grids = .false.
6477 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
6479 IF (my_beta /= 1.0_dp)
THEN
6480 IF (my_beta == 0.0_dp)
THEN
6484 pw2%array = pw2%array*my_beta
6489 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
6491 IF (my_alpha == 1.0_dp)
THEN
6493 pw2%array = pw2%array + pw1%array
6495 ELSE IF (my_alpha /= 0.0_dp)
THEN
6497 pw2%array = pw2%array + my_alpha* pw1%array
6501 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
6503 ng1 =
SIZE(pw1%array)
6504 ng2 =
SIZE(pw2%array)
6507 IF (pw1%pw_grid%spherical)
THEN
6508 IF (my_alpha == 1.0_dp)
THEN
6511 pw2%array(i) = pw2%array(i) + pw1%array(i)
6514 ELSE IF (my_alpha /= 0.0_dp)
THEN
6517 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(i)
6521 ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference))
THEN
6522 IF (ng1 >= ng2)
THEN
6523 IF (my_alpha == 1.0_dp)
THEN
6526 j = pw2%pw_grid%gidx(i)
6527 pw2%array(i) = pw2%array(i) + pw1%array(j)
6530 ELSE IF (my_alpha /= 0.0_dp)
THEN
6533 j = pw2%pw_grid%gidx(i)
6534 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
6539 IF (my_alpha == 1.0_dp)
THEN
6542 j = pw2%pw_grid%gidx(i)
6543 pw2%array(j) = pw2%array(j) + pw1%array(i)
6546 ELSE IF (my_alpha /= 0.0_dp)
THEN
6549 j = pw2%pw_grid%gidx(i)
6550 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
6555 ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference))
THEN
6556 IF (ng1 >= ng2)
THEN
6557 IF (my_alpha == 1.0_dp)
THEN
6560 j = pw1%pw_grid%gidx(i)
6561 pw2%array(i) = pw2%array(i) + pw1%array(j)
6564 ELSE IF (my_alpha /= 0.0_dp)
THEN
6567 j = pw1%pw_grid%gidx(i)
6568 pw2%array(i) = pw2%array(i) + my_alpha* pw1%array(j)
6573 IF (my_alpha == 1.0_dp)
THEN
6576 j = pw1%pw_grid%gidx(i)
6577 pw2%array(j) = pw2%array(j) + pw1%array(i)
6580 ELSE IF (my_alpha /= 0.0_dp)
THEN
6583 j = pw1%pw_grid%gidx(i)
6584 pw2%array(j) = pw2%array(j) + my_alpha* pw1%array(i)
6590 cpabort(
"Grids not compatible")
6595 cpabort(
"Grids not compatible")
6599 CALL timestop(handle)
6601 END SUBROUTINE pw_axpy_c1d_c1d_gs
6612 SUBROUTINE pw_multiply_c1d_c1d_gs (pw_out, pw1, pw2, alpha)
6614 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw_out
6615 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
6616 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
6617 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: alpha
6619 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_multiply'
6622 REAL(kind=dp) :: my_alpha
6624 CALL timeset(routinen, handle)
6627 IF (
PRESENT(alpha)) my_alpha = alpha
6629 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
6630 cpabort(
"pw_multiply not implemented for non-identical grids!")
6632 IF (my_alpha == 1.0_dp)
THEN
6634 pw_out%array = pw_out%array + pw1%array* pw2%array
6638 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
6642 CALL timestop(handle)
6644 END SUBROUTINE pw_multiply_c1d_c1d_gs
6651 SUBROUTINE pw_multiply_with_c1d_c1d_gs (pw1, pw2)
6652 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw1
6653 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
6655 CHARACTER(LEN=*),
PARAMETER :: routinen =
'pw_multiply_with'
6659 CALL timeset(routinen, handle)
6661 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
6662 cpabort(
"Incompatible grids!")
6665 pw1%array = pw1%array* pw2%array
6668 CALL timestop(handle)
6670 END SUBROUTINE pw_multiply_with_c1d_c1d_gs
6685 FUNCTION pw_integral_ab_c1d_c1d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
6687 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
6688 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
6689 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
6690 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
6691 REAL(kind=dp) :: integral_value
6693 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab'
6695 INTEGER :: handle, loc_sumtype
6696 LOGICAL :: my_just_sum, my_local_only
6698 CALL timeset(routinen, handle)
6701 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
6703 my_local_only = .false.
6704 IF (
PRESENT(local_only)) my_local_only = local_only
6706 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
6707 cpabort(
"Grids incompatible")
6710 my_just_sum = .false.
6711 IF (
PRESENT(just_sum)) my_just_sum = just_sum
6718 integral_value = sum(real(conjg(pw1%array) &
6719 *pw2%array, kind=dp))
6724 integral_value = accurate_sum(real(conjg(pw1%array) &
6725 *pw2%array, kind=dp))
6729 IF (.NOT. my_just_sum)
THEN
6730 integral_value = integral_value*pw1%pw_grid%vol
6733 IF (pw1%pw_grid%grid_span == halfspace)
THEN
6734 integral_value = 2.0_dp*integral_value
6735 IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
6736 REAL(conjg(pw1%array(1))*pw2%array(1), kind=dp)
6739 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
6740 CALL pw1%pw_grid%para%group%sum(integral_value)
6742 CALL timestop(handle)
6744 END FUNCTION pw_integral_ab_c1d_c1d_gs
6752 FUNCTION pw_integral_a2b_c1d_c1d (pw1, pw2)
RESULT(integral_value)
6754 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
6755 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw2
6756 REAL(kind=dp) :: integral_value
6758 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_a2b'
6762 CALL timeset(routinen, handle)
6764 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
6765 cpabort(
"Grids incompatible")
6768 integral_value = accurate_sum(real(conjg(pw1%array)*pw2%array, kind=dp)*pw1%pw_grid%gsq)
6769 IF (pw1%pw_grid%grid_span == halfspace)
THEN
6770 integral_value = 2.0_dp*integral_value
6773 integral_value = integral_value*pw1%pw_grid%vol
6775 IF (pw1%pw_grid%para%mode == pw_mode_distributed) &
6776 CALL pw1%pw_grid%para%group%sum(integral_value)
6777 CALL timestop(handle)
6779 END FUNCTION pw_integral_a2b_c1d_c1d
6793 SUBROUTINE pw_copy_c3d_r3d_rs (pw1, pw2)
6795 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
6796 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw2
6798 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy'
6802 CALL timeset(routinen, handle)
6804 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
6805 cpabort(
"Both grids must be either spherical or non-spherical!")
6806 IF (pw1%pw_grid%spherical) &
6807 cpabort(
"Spherical grids only exist in reciprocal space!")
6809 IF (any(shape(pw2%array) /= shape(pw1%array))) &
6810 cpabort(
"3D grids must be compatible!")
6811 IF (pw1%pw_grid%spherical) &
6812 cpabort(
"3D grids must not be spherical!")
6814 pw2%array(:, :, :) = real(pw1%array(:, :, :), kind=dp)
6817 CALL timestop(handle)
6819 END SUBROUTINE pw_copy_c3d_r3d_rs
6826 SUBROUTINE pw_copy_to_array_c3d_r3d_rs (pw, array)
6827 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw
6828 REAL(kind=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
6830 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_to_array'
6834 CALL timeset(routinen, handle)
6837 array(:, :, :) = real(pw%array(:, :, :), kind=dp)
6840 CALL timestop(handle)
6841 END SUBROUTINE pw_copy_to_array_c3d_r3d_rs
6848 SUBROUTINE pw_copy_from_array_c3d_r3d_rs (pw, array)
6849 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw
6850 REAL(kind=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
6852 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_from_array'
6856 CALL timeset(routinen, handle)
6859 pw%array = cmplx(array, 0.0_dp, kind=dp)
6862 CALL timestop(handle)
6863 END SUBROUTINE pw_copy_from_array_c3d_r3d_rs
6881 SUBROUTINE pw_axpy_c3d_r3d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
6883 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
6884 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw2
6885 REAL(kind=dp),
INTENT(in),
OPTIONAL :: alpha, beta
6886 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
6888 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_axpy'
6891 LOGICAL :: my_allow_noncompatible_grids
6892 REAL(kind=dp) :: my_alpha, my_beta
6894 CALL timeset(routinen, handle)
6897 IF (
PRESENT(alpha)) my_alpha = alpha
6900 IF (
PRESENT(beta)) my_beta = beta
6902 my_allow_noncompatible_grids = .false.
6903 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
6905 IF (my_beta /= 1.0_dp)
THEN
6906 IF (my_beta == 0.0_dp)
THEN
6910 pw2%array = pw2%array*my_beta
6915 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
6916 IF (my_alpha == 1.0_dp)
THEN
6918 pw2%array = pw2%array + real(pw1%array, kind=dp)
6920 ELSE IF (my_alpha /= 0.0_dp)
THEN
6922 pw2%array = pw2%array + my_alpha* real(pw1%array, kind=dp)
6926 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
6928 IF (any(shape(pw1%array) /= shape(pw2%array))) &
6929 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
6931 IF (my_alpha == 1.0_dp)
THEN
6933 pw2%array = pw2%array + real(pw1%array, kind=dp)
6937 pw2%array = pw2%array + my_alpha* real(pw1%array, kind=dp)
6943 cpabort(
"Grids not compatible")
6947 CALL timestop(handle)
6949 END SUBROUTINE pw_axpy_c3d_r3d_rs
6960 SUBROUTINE pw_multiply_c3d_r3d_rs (pw_out, pw1, pw2, alpha)
6962 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw_out
6963 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
6964 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw2
6965 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: alpha
6967 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_multiply'
6970 REAL(kind=dp) :: my_alpha
6972 CALL timeset(routinen, handle)
6975 IF (
PRESENT(alpha)) my_alpha = alpha
6977 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
6978 cpabort(
"pw_multiply not implemented for non-identical grids!")
6980 IF (my_alpha == 1.0_dp)
THEN
6982 pw_out%array = pw_out%array + pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
6986 pw_out%array = pw_out%array + my_alpha*pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
6990 CALL timestop(handle)
6992 END SUBROUTINE pw_multiply_c3d_r3d_rs
6999 SUBROUTINE pw_multiply_with_c3d_r3d_rs (pw1, pw2)
7000 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw1
7001 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw2
7003 CHARACTER(LEN=*),
PARAMETER :: routinen =
'pw_multiply_with'
7007 CALL timeset(routinen, handle)
7009 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
7010 cpabort(
"Incompatible grids!")
7013 pw1%array = pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7016 CALL timestop(handle)
7018 END SUBROUTINE pw_multiply_with_c3d_r3d_rs
7033 FUNCTION pw_integral_ab_c3d_r3d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
7035 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
7036 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw2
7037 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
7038 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
7039 REAL(kind=dp) :: integral_value
7041 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab'
7043 INTEGER :: handle, loc_sumtype
7044 LOGICAL :: my_just_sum, my_local_only
7046 CALL timeset(routinen, handle)
7049 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
7051 my_local_only = .false.
7052 IF (
PRESENT(local_only)) my_local_only = local_only
7054 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
7055 cpabort(
"Grids incompatible")
7058 my_just_sum = .false.
7059 IF (
PRESENT(just_sum)) my_just_sum = just_sum
7066 integral_value = sum(real(pw1%array, kind=dp)*pw2%array)
7071 integral_value = accurate_sum(real(pw1%array, kind=dp)*pw2%array)
7075 IF (.NOT. my_just_sum)
THEN
7076 integral_value = integral_value*pw1%pw_grid%dvol
7080 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
7081 CALL pw1%pw_grid%para%group%sum(integral_value)
7083 CALL timestop(handle)
7085 END FUNCTION pw_integral_ab_c3d_r3d_rs
7099 SUBROUTINE pw_copy_c3d_r3d_gs (pw1, pw2)
7101 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
7102 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw2
7104 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy'
7108 CALL timeset(routinen, handle)
7110 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
7111 cpabort(
"Both grids must be either spherical or non-spherical!")
7113 IF (any(shape(pw2%array) /= shape(pw1%array))) &
7114 cpabort(
"3D grids must be compatible!")
7115 IF (pw1%pw_grid%spherical) &
7116 cpabort(
"3D grids must not be spherical!")
7118 pw2%array(:, :, :) = real(pw1%array(:, :, :), kind=dp)
7121 CALL timestop(handle)
7123 END SUBROUTINE pw_copy_c3d_r3d_gs
7130 SUBROUTINE pw_copy_to_array_c3d_r3d_gs (pw, array)
7131 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw
7132 REAL(kind=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
7134 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_to_array'
7138 CALL timeset(routinen, handle)
7141 array(:, :, :) = real(pw%array(:, :, :), kind=dp)
7144 CALL timestop(handle)
7145 END SUBROUTINE pw_copy_to_array_c3d_r3d_gs
7152 SUBROUTINE pw_copy_from_array_c3d_r3d_gs (pw, array)
7153 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw
7154 REAL(kind=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
7156 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_from_array'
7160 CALL timeset(routinen, handle)
7163 pw%array = cmplx(array, 0.0_dp, kind=dp)
7166 CALL timestop(handle)
7167 END SUBROUTINE pw_copy_from_array_c3d_r3d_gs
7185 SUBROUTINE pw_axpy_c3d_r3d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
7187 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
7188 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw2
7189 REAL(kind=dp),
INTENT(in),
OPTIONAL :: alpha, beta
7190 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
7192 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_axpy'
7195 LOGICAL :: my_allow_noncompatible_grids
7196 REAL(kind=dp) :: my_alpha, my_beta
7198 CALL timeset(routinen, handle)
7201 IF (
PRESENT(alpha)) my_alpha = alpha
7204 IF (
PRESENT(beta)) my_beta = beta
7206 my_allow_noncompatible_grids = .false.
7207 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
7209 IF (my_beta /= 1.0_dp)
THEN
7210 IF (my_beta == 0.0_dp)
THEN
7214 pw2%array = pw2%array*my_beta
7219 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
7220 IF (my_alpha == 1.0_dp)
THEN
7222 pw2%array = pw2%array + real(pw1%array, kind=dp)
7224 ELSE IF (my_alpha /= 0.0_dp)
THEN
7226 pw2%array = pw2%array + my_alpha* real(pw1%array, kind=dp)
7230 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
7232 IF (any(shape(pw1%array) /= shape(pw2%array))) &
7233 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
7235 IF (my_alpha == 1.0_dp)
THEN
7237 pw2%array = pw2%array + real(pw1%array, kind=dp)
7241 pw2%array = pw2%array + my_alpha* real(pw1%array, kind=dp)
7247 cpabort(
"Grids not compatible")
7251 CALL timestop(handle)
7253 END SUBROUTINE pw_axpy_c3d_r3d_gs
7264 SUBROUTINE pw_multiply_c3d_r3d_gs (pw_out, pw1, pw2, alpha)
7266 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw_out
7267 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
7268 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw2
7269 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: alpha
7271 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_multiply'
7274 REAL(kind=dp) :: my_alpha
7276 CALL timeset(routinen, handle)
7279 IF (
PRESENT(alpha)) my_alpha = alpha
7281 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
7282 cpabort(
"pw_multiply not implemented for non-identical grids!")
7284 IF (my_alpha == 1.0_dp)
THEN
7286 pw_out%array = pw_out%array + pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7290 pw_out%array = pw_out%array + my_alpha*pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7294 CALL timestop(handle)
7296 END SUBROUTINE pw_multiply_c3d_r3d_gs
7303 SUBROUTINE pw_multiply_with_c3d_r3d_gs (pw1, pw2)
7304 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw1
7305 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw2
7307 CHARACTER(LEN=*),
PARAMETER :: routinen =
'pw_multiply_with'
7311 CALL timeset(routinen, handle)
7313 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
7314 cpabort(
"Incompatible grids!")
7317 pw1%array = pw1%array* cmplx(pw2%array, 0.0_dp, kind=dp)
7320 CALL timestop(handle)
7322 END SUBROUTINE pw_multiply_with_c3d_r3d_gs
7337 FUNCTION pw_integral_ab_c3d_r3d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
7339 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
7340 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw2
7341 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
7342 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
7343 REAL(kind=dp) :: integral_value
7345 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab'
7347 INTEGER :: handle, loc_sumtype
7348 LOGICAL :: my_just_sum, my_local_only
7350 CALL timeset(routinen, handle)
7353 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
7355 my_local_only = .false.
7356 IF (
PRESENT(local_only)) my_local_only = local_only
7358 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
7359 cpabort(
"Grids incompatible")
7362 my_just_sum = .false.
7363 IF (
PRESENT(just_sum)) my_just_sum = just_sum
7370 integral_value = sum(real(pw1%array, kind=dp)*pw2%array)
7375 integral_value = accurate_sum(real(pw1%array, kind=dp)*pw2%array)
7379 IF (.NOT. my_just_sum)
THEN
7380 integral_value = integral_value*pw1%pw_grid%vol
7384 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
7385 CALL pw1%pw_grid%para%group%sum(integral_value)
7387 CALL timestop(handle)
7389 END FUNCTION pw_integral_ab_c3d_r3d_gs
7404 SUBROUTINE pw_copy_c3d_c3d_rs (pw1, pw2)
7406 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
7407 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw2
7409 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy'
7413 CALL timeset(routinen, handle)
7415 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
7416 cpabort(
"Both grids must be either spherical or non-spherical!")
7417 IF (pw1%pw_grid%spherical) &
7418 cpabort(
"Spherical grids only exist in reciprocal space!")
7420 IF (any(shape(pw2%array) /= shape(pw1%array))) &
7421 cpabort(
"3D grids must be compatible!")
7422 IF (pw1%pw_grid%spherical) &
7423 cpabort(
"3D grids must not be spherical!")
7425 pw2%array(:, :, :) = pw1%array(:, :, :)
7428 CALL timestop(handle)
7430 END SUBROUTINE pw_copy_c3d_c3d_rs
7437 SUBROUTINE pw_copy_to_array_c3d_c3d_rs (pw, array)
7438 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw
7439 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
7441 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_to_array'
7445 CALL timeset(routinen, handle)
7448 array(:, :, :) = pw%array(:, :, :)
7451 CALL timestop(handle)
7452 END SUBROUTINE pw_copy_to_array_c3d_c3d_rs
7459 SUBROUTINE pw_copy_from_array_c3d_c3d_rs (pw, array)
7460 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw
7461 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
7463 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_from_array'
7467 CALL timeset(routinen, handle)
7473 CALL timestop(handle)
7474 END SUBROUTINE pw_copy_from_array_c3d_c3d_rs
7492 SUBROUTINE pw_axpy_c3d_c3d_rs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
7494 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
7495 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw2
7496 REAL(kind=dp),
INTENT(in),
OPTIONAL :: alpha, beta
7497 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
7499 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_axpy'
7502 LOGICAL :: my_allow_noncompatible_grids
7503 REAL(kind=dp) :: my_alpha, my_beta
7505 CALL timeset(routinen, handle)
7508 IF (
PRESENT(alpha)) my_alpha = alpha
7511 IF (
PRESENT(beta)) my_beta = beta
7513 my_allow_noncompatible_grids = .false.
7514 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
7516 IF (my_beta /= 1.0_dp)
THEN
7517 IF (my_beta == 0.0_dp)
THEN
7521 pw2%array = pw2%array*my_beta
7526 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
7527 IF (my_alpha == 1.0_dp)
THEN
7529 pw2%array = pw2%array + pw1%array
7531 ELSE IF (my_alpha /= 0.0_dp)
THEN
7533 pw2%array = pw2%array + my_alpha* pw1%array
7537 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
7539 IF (any(shape(pw1%array) /= shape(pw2%array))) &
7540 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
7542 IF (my_alpha == 1.0_dp)
THEN
7544 pw2%array = pw2%array + pw1%array
7548 pw2%array = pw2%array + my_alpha* pw1%array
7554 cpabort(
"Grids not compatible")
7558 CALL timestop(handle)
7560 END SUBROUTINE pw_axpy_c3d_c3d_rs
7571 SUBROUTINE pw_multiply_c3d_c3d_rs (pw_out, pw1, pw2, alpha)
7573 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw_out
7574 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
7575 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw2
7576 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: alpha
7578 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_multiply'
7581 REAL(kind=dp) :: my_alpha
7583 CALL timeset(routinen, handle)
7586 IF (
PRESENT(alpha)) my_alpha = alpha
7588 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
7589 cpabort(
"pw_multiply not implemented for non-identical grids!")
7591 IF (my_alpha == 1.0_dp)
THEN
7593 pw_out%array = pw_out%array + pw1%array* pw2%array
7597 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
7601 CALL timestop(handle)
7603 END SUBROUTINE pw_multiply_c3d_c3d_rs
7610 SUBROUTINE pw_multiply_with_c3d_c3d_rs (pw1, pw2)
7611 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw1
7612 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw2
7614 CHARACTER(LEN=*),
PARAMETER :: routinen =
'pw_multiply_with'
7618 CALL timeset(routinen, handle)
7620 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
7621 cpabort(
"Incompatible grids!")
7624 pw1%array = pw1%array* pw2%array
7627 CALL timestop(handle)
7629 END SUBROUTINE pw_multiply_with_c3d_c3d_rs
7644 FUNCTION pw_integral_ab_c3d_c3d_rs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
7646 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
7647 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw2
7648 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
7649 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
7650 REAL(kind=dp) :: integral_value
7652 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab'
7654 INTEGER :: handle, loc_sumtype
7655 LOGICAL :: my_just_sum, my_local_only
7657 CALL timeset(routinen, handle)
7660 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
7662 my_local_only = .false.
7663 IF (
PRESENT(local_only)) my_local_only = local_only
7665 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
7666 cpabort(
"Grids incompatible")
7669 my_just_sum = .false.
7670 IF (
PRESENT(just_sum)) my_just_sum = just_sum
7677 integral_value = sum(real(conjg(pw1%array) &
7678 *pw2%array, kind=dp))
7683 integral_value = accurate_sum(real(conjg(pw1%array) &
7684 *pw2%array, kind=dp))
7688 IF (.NOT. my_just_sum)
THEN
7689 integral_value = integral_value*pw1%pw_grid%dvol
7693 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
7694 CALL pw1%pw_grid%para%group%sum(integral_value)
7696 CALL timestop(handle)
7698 END FUNCTION pw_integral_ab_c3d_c3d_rs
7712 SUBROUTINE pw_copy_c3d_c3d_gs (pw1, pw2)
7714 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
7715 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
7717 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy'
7721 CALL timeset(routinen, handle)
7723 IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
7724 cpabort(
"Both grids must be either spherical or non-spherical!")
7726 IF (any(shape(pw2%array) /= shape(pw1%array))) &
7727 cpabort(
"3D grids must be compatible!")
7728 IF (pw1%pw_grid%spherical) &
7729 cpabort(
"3D grids must not be spherical!")
7731 pw2%array(:, :, :) = pw1%array(:, :, :)
7734 CALL timestop(handle)
7736 END SUBROUTINE pw_copy_c3d_c3d_gs
7743 SUBROUTINE pw_copy_to_array_c3d_c3d_gs (pw, array)
7744 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw
7745 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: array
7747 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_to_array'
7751 CALL timeset(routinen, handle)
7754 array(:, :, :) = pw%array(:, :, :)
7757 CALL timestop(handle)
7758 END SUBROUTINE pw_copy_to_array_c3d_c3d_gs
7765 SUBROUTINE pw_copy_from_array_c3d_c3d_gs (pw, array)
7766 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw
7767 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
INTENT(IN) :: array
7769 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_copy_from_array'
7773 CALL timeset(routinen, handle)
7779 CALL timestop(handle)
7780 END SUBROUTINE pw_copy_from_array_c3d_c3d_gs
7798 SUBROUTINE pw_axpy_c3d_c3d_gs (pw1, pw2, alpha, beta, allow_noncompatible_grids)
7800 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
7801 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
7802 REAL(kind=dp),
INTENT(in),
OPTIONAL :: alpha, beta
7803 LOGICAL,
INTENT(IN),
OPTIONAL :: allow_noncompatible_grids
7805 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_axpy'
7808 LOGICAL :: my_allow_noncompatible_grids
7809 REAL(kind=dp) :: my_alpha, my_beta
7811 CALL timeset(routinen, handle)
7814 IF (
PRESENT(alpha)) my_alpha = alpha
7817 IF (
PRESENT(beta)) my_beta = beta
7819 my_allow_noncompatible_grids = .false.
7820 IF (
PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
7822 IF (my_beta /= 1.0_dp)
THEN
7823 IF (my_beta == 0.0_dp)
THEN
7827 pw2%array = pw2%array*my_beta
7832 IF (
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
7833 IF (my_alpha == 1.0_dp)
THEN
7835 pw2%array = pw2%array + pw1%array
7837 ELSE IF (my_alpha /= 0.0_dp)
THEN
7839 pw2%array = pw2%array + my_alpha* pw1%array
7843 ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids)
THEN
7845 IF (any(shape(pw1%array) /= shape(pw2%array))) &
7846 cpabort(
"Noncommensurate grids not implemented for 3D grids!")
7848 IF (my_alpha == 1.0_dp)
THEN
7850 pw2%array = pw2%array + pw1%array
7854 pw2%array = pw2%array + my_alpha* pw1%array
7860 cpabort(
"Grids not compatible")
7864 CALL timestop(handle)
7866 END SUBROUTINE pw_axpy_c3d_c3d_gs
7877 SUBROUTINE pw_multiply_c3d_c3d_gs (pw_out, pw1, pw2, alpha)
7879 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw_out
7880 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
7881 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw2
7882 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: alpha
7884 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_multiply'
7887 REAL(kind=dp) :: my_alpha
7889 CALL timeset(routinen, handle)
7892 IF (
PRESENT(alpha)) my_alpha = alpha
7894 IF (.NOT.
ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT.
ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
7895 cpabort(
"pw_multiply not implemented for non-identical grids!")
7897 IF (my_alpha == 1.0_dp)
THEN
7899 pw_out%array = pw_out%array + pw1%array* pw2%array
7903 pw_out%array = pw_out%array + my_alpha*pw1%array* pw2%array
7907 CALL timestop(handle)
7909 END SUBROUTINE pw_multiply_c3d_c3d_gs
7916 SUBROUTINE pw_multiply_with_c3d_c3d_gs (pw1, pw2)
7917 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw1
7918 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw2
7920 CHARACTER(LEN=*),
PARAMETER :: routinen =
'pw_multiply_with'
7924 CALL timeset(routinen, handle)
7926 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
7927 cpabort(
"Incompatible grids!")
7930 pw1%array = pw1%array* pw2%array
7933 CALL timestop(handle)
7935 END SUBROUTINE pw_multiply_with_c3d_c3d_gs
7950 FUNCTION pw_integral_ab_c3d_c3d_gs (pw1, pw2, sumtype, just_sum, local_only)
RESULT(integral_value)
7952 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
7953 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw2
7954 INTEGER,
INTENT(IN),
OPTIONAL :: sumtype
7955 LOGICAL,
INTENT(IN),
OPTIONAL :: just_sum, local_only
7956 REAL(kind=dp) :: integral_value
7958 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_integral_ab'
7960 INTEGER :: handle, loc_sumtype
7961 LOGICAL :: my_just_sum, my_local_only
7963 CALL timeset(routinen, handle)
7966 IF (
PRESENT(sumtype)) loc_sumtype = sumtype
7968 my_local_only = .false.
7969 IF (
PRESENT(local_only)) my_local_only = local_only
7971 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
7972 cpabort(
"Grids incompatible")
7975 my_just_sum = .false.
7976 IF (
PRESENT(just_sum)) my_just_sum = just_sum
7983 integral_value = sum(real(conjg(pw1%array) &
7984 *pw2%array, kind=dp))
7989 integral_value = accurate_sum(real(conjg(pw1%array) &
7990 *pw2%array, kind=dp))
7994 IF (.NOT. my_just_sum)
THEN
7995 integral_value = integral_value*pw1%pw_grid%vol
7999 IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == pw_mode_distributed) &
8000 CALL pw1%pw_grid%para%group%sum(integral_value)
8002 CALL timestop(handle)
8004 END FUNCTION pw_integral_ab_c3d_c3d_gs
8027 SUBROUTINE pw_gather_s_r1d_r3d_2(pw1, pw2, scale)
8029 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
8030 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw2
8031 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
8033 CALL pw_gather_s_r1d_r3d (pw2, pw1%array, scale)
8035 END SUBROUTINE pw_gather_s_r1d_r3d_2
8046 SUBROUTINE pw_gather_s_r1d_r3d (pw, c, scale)
8048 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw
8049 REAL(kind=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(IN) :: c
8050 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
8052 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_gather_s'
8054 INTEGER :: gpt, handle, l, m, n
8056 CALL timeset(routinen, handle)
8058 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
8059 ngpts =>
SIZE(pw%pw_grid%gsq), ghat => pw%pw_grid%g_hat)
8061 IF (
PRESENT(scale))
THEN
8064 l = mapl(ghat(1, gpt)) + 1
8065 m = mapm(ghat(2, gpt)) + 1
8066 n = mapn(ghat(3, gpt)) + 1
8067 pw%array(gpt) = scale* c(l, m, n)
8073 l = mapl(ghat(1, gpt)) + 1
8074 m = mapm(ghat(2, gpt)) + 1
8075 n = mapn(ghat(3, gpt)) + 1
8076 pw%array(gpt) = c(l, m, n)
8083 CALL timestop(handle)
8085 END SUBROUTINE pw_gather_s_r1d_r3d
8096 SUBROUTINE pw_scatter_s_r1d_r3d_2(pw1, pw2, scale)
8098 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
8099 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw2
8100 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
8102 CALL pw_scatter_s_r1d_r3d (pw1, pw2%array, scale)
8104 END SUBROUTINE pw_scatter_s_r1d_r3d_2
8115 SUBROUTINE pw_scatter_s_r1d_r3d (pw, c, scale)
8117 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw
8118 REAL(kind=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(INOUT) :: c
8119 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
8121 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_scatter_s'
8123 INTEGER :: gpt, handle, l, m, n
8125 CALL timeset(routinen, handle)
8127 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
8128 ghat => pw%pw_grid%g_hat, ngpts =>
SIZE(pw%pw_grid%gsq))
8131 IF (.NOT.
PRESENT(scale)) c = 0.0_dp
8133 IF (
PRESENT(scale))
THEN
8136 l = mapl(ghat(1, gpt)) + 1
8137 m = mapm(ghat(2, gpt)) + 1
8138 n = mapn(ghat(3, gpt)) + 1
8139 c(l, m, n) = scale* pw%array(gpt)
8145 l = mapl(ghat(1, gpt)) + 1
8146 m = mapm(ghat(2, gpt)) + 1
8147 n = mapn(ghat(3, gpt)) + 1
8148 c(l, m, n) = pw%array(gpt)
8155 IF (pw%pw_grid%grid_span == halfspace)
THEN
8157 associate(mapl => pw%pw_grid%mapl%neg, mapm => pw%pw_grid%mapm%neg, mapn => pw%pw_grid%mapn%neg, &
8158 ghat => pw%pw_grid%g_hat, ngpts =>
SIZE(pw%pw_grid%gsq))
8160 IF (
PRESENT(scale))
THEN
8163 l = mapl(ghat(1, gpt)) + 1
8164 m = mapm(ghat(2, gpt)) + 1
8165 n = mapn(ghat(3, gpt)) + 1
8166 c(l, m, n) = scale*( pw%array(gpt))
8172 l = mapl(ghat(1, gpt)) + 1
8173 m = mapm(ghat(2, gpt)) + 1
8174 n = mapn(ghat(3, gpt)) + 1
8175 c(l, m, n) = ( pw%array(gpt))
8184 CALL timestop(handle)
8186 END SUBROUTINE pw_scatter_s_r1d_r3d
8207 SUBROUTINE pw_gather_s_r1d_c3d_2(pw1, pw2, scale)
8209 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
8210 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw2
8211 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
8213 CALL pw_gather_s_r1d_c3d (pw2, pw1%array, scale)
8215 END SUBROUTINE pw_gather_s_r1d_c3d_2
8226 SUBROUTINE pw_gather_s_r1d_c3d (pw, c, scale)
8228 TYPE(pw_r1d_gs_type),
INTENT(INOUT) :: pw
8229 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(IN) :: c
8230 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
8232 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_gather_s'
8234 INTEGER :: gpt, handle, l, m, n
8236 CALL timeset(routinen, handle)
8238 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
8239 ngpts =>
SIZE(pw%pw_grid%gsq), ghat => pw%pw_grid%g_hat)
8241 IF (
PRESENT(scale))
THEN
8244 l = mapl(ghat(1, gpt)) + 1
8245 m = mapm(ghat(2, gpt)) + 1
8246 n = mapn(ghat(3, gpt)) + 1
8247 pw%array(gpt) = scale* real(c(l, m, n), kind=dp)
8253 l = mapl(ghat(1, gpt)) + 1
8254 m = mapm(ghat(2, gpt)) + 1
8255 n = mapn(ghat(3, gpt)) + 1
8256 pw%array(gpt) = real(c(l, m, n), kind=dp)
8263 CALL timestop(handle)
8265 END SUBROUTINE pw_gather_s_r1d_c3d
8276 SUBROUTINE pw_scatter_s_r1d_c3d_2(pw1, pw2, scale)
8278 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw1
8279 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
8280 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
8282 CALL pw_scatter_s_r1d_c3d (pw1, pw2%array, scale)
8284 END SUBROUTINE pw_scatter_s_r1d_c3d_2
8295 SUBROUTINE pw_scatter_s_r1d_c3d (pw, c, scale)
8297 TYPE(pw_r1d_gs_type),
INTENT(IN) :: pw
8298 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(INOUT) :: c
8299 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
8301 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_scatter_s'
8303 INTEGER :: gpt, handle, l, m, n
8305 CALL timeset(routinen, handle)
8307 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
8308 ghat => pw%pw_grid%g_hat, ngpts =>
SIZE(pw%pw_grid%gsq))
8311 IF (.NOT.
PRESENT(scale)) c = 0.0_dp
8313 IF (
PRESENT(scale))
THEN
8316 l = mapl(ghat(1, gpt)) + 1
8317 m = mapm(ghat(2, gpt)) + 1
8318 n = mapn(ghat(3, gpt)) + 1
8319 c(l, m, n) = scale* cmplx(pw%array(gpt), 0.0_dp, kind=dp)
8325 l = mapl(ghat(1, gpt)) + 1
8326 m = mapm(ghat(2, gpt)) + 1
8327 n = mapn(ghat(3, gpt)) + 1
8328 c(l, m, n) = cmplx(pw%array(gpt), 0.0_dp, kind=dp)
8335 IF (pw%pw_grid%grid_span == halfspace)
THEN
8337 associate(mapl => pw%pw_grid%mapl%neg, mapm => pw%pw_grid%mapm%neg, mapn => pw%pw_grid%mapn%neg, &
8338 ghat => pw%pw_grid%g_hat, ngpts =>
SIZE(pw%pw_grid%gsq))
8340 IF (
PRESENT(scale))
THEN
8343 l = mapl(ghat(1, gpt)) + 1
8344 m = mapm(ghat(2, gpt)) + 1
8345 n = mapn(ghat(3, gpt)) + 1
8346 c(l, m, n) = scale*( cmplx(pw%array(gpt), 0.0_dp, kind=dp))
8352 l = mapl(ghat(1, gpt)) + 1
8353 m = mapm(ghat(2, gpt)) + 1
8354 n = mapn(ghat(3, gpt)) + 1
8355 c(l, m, n) = ( cmplx(pw%array(gpt), 0.0_dp, kind=dp))
8364 CALL timestop(handle)
8366 END SUBROUTINE pw_scatter_s_r1d_c3d
8392 SUBROUTINE fft_wrap_pw1pw2_r3d_c1d_rs_gs (pw1, pw2, debug)
8394 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
8395 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
8396 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
8398 CHARACTER(len=*),
PARAMETER :: routinen =
'fft_wrap_pw1pw2'
8400 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
8401 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
8402 INTEGER :: handle, handle2, my_pos, nrays, &
8404 INTEGER,
DIMENSION(3) :: nloc
8405 #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
8406 LOGICAL :: use_pw_gpu
8408 INTEGER,
DIMENSION(:),
POINTER :: n
8410 REAL(kind=dp) :: norm
8412 CALL timeset(routinen, handle2)
8413 out_unit = cp_logger_get_default_io_unit()
8414 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
8415 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
8420 IF (
PRESENT(debug))
THEN
8427 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
8428 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
8429 cpabort(
"PW grids not compatible")
8431 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
8432 cpabort(
"PW grids have not compatible MPI groups")
8437 norm = 1.0_dp/pw1%pw_grid%ngpts
8439 n => pw1%pw_grid%npts
8441 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
8447 IF (test .AND. out_unit > 0)
THEN
8448 WRITE (out_unit,
'(A)')
" FFT Protocol "
8449 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
8450 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
8451 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
8452 WRITE (out_unit,
'(A,T66,E15.6)')
" scale factor", norm
8455 #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
8456 CALL pw_gpu_r3dc1d_3d(pw1, pw2, scale=norm)
8457 #elif defined (__PW_FPGA)
8458 ALLOCATE (c_out(n(1), n(2), n(3)))
8461 IF (pw_fpga_init_bitstream(n) == 1)
THEN
8462 CALL pw_copy_to_array(pw1, c_out)
8463 #if (__PW_FPGA_SP && __PW_FPGA)
8464 CALL pw_fpga_r3dc1d_3d_sp(n, c_out)
8466 CALL pw_fpga_r3dc1d_3d_dp(n, c_out)
8468 CALL zdscal(n(1)*n(2)*n(3), norm, c_out, 1)
8469 CALL pw_gather_s_c1d_c3d(pw2, c_out)
8471 CALL pw_copy_to_array(pw1, c_out)
8472 CALL fft3d(fwfft, n, c_out, scale=norm, debug=test)
8473 CALL pw_gather_s_c1d_c3d(pw2, c_out)
8477 ALLOCATE (c_out(n(1), n(2), n(3)))
8479 CALL pw_copy_to_array(pw1, c_out)
8480 CALL fft3d(fwfft, n, c_out, scale=norm, debug=test)
8481 CALL pw_gather_s_c1d_c3d(pw2, c_out)
8485 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
8493 IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0)
THEN
8494 WRITE (out_unit,
'(A)')
" FFT Protocol "
8495 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
8496 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
8497 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
8498 WRITE (out_unit,
'(A,T66,E15.6)')
" scale factor", norm
8501 my_pos = pw1%pw_grid%para%my_pos
8502 nrays = pw1%pw_grid%para%nyzray(my_pos)
8503 grays => pw1%pw_grid%grays
8505 #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
8507 use_pw_gpu = pw1%pw_grid%para%ray_distribution
8508 IF (use_pw_gpu)
THEN
8509 CALL pw_gpu_r3dc1d_3d_ps(pw1, pw2, scale=norm)
8513 nloc = pw1%pw_grid%npts_local
8514 ALLOCATE (c_in(nloc(1), nloc(2), nloc(3)))
8515 CALL pw_copy_to_array(pw1, c_in)
8518 IF (pw1%pw_grid%para%ray_distribution)
THEN
8519 CALL fft3d(fwfft, n, c_in, grays, pw1%pw_grid%para%group, &
8520 pw1%pw_grid%para%rs_group, &
8521 pw1%pw_grid%para%yzp, pw1%pw_grid%para%nyzray, &
8522 pw1%pw_grid%para%bo, scale=norm, debug=test)
8524 CALL fft3d(fwfft, n, c_in, grays, pw1%pw_grid%para%rs_group, &
8525 pw1%pw_grid%para%bo, scale=norm, debug=test)
8528 IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0) &
8529 WRITE (out_unit,
'(A)')
" PW_GATHER : 2d -> 1d "
8530 CALL pw_gather_p_c1d (pw2, grays)
8533 #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
8538 IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0)
THEN
8539 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
8542 CALL timestop(handle)
8543 CALL timestop(handle2)
8545 END SUBROUTINE fft_wrap_pw1pw2_r3d_c1d_rs_gs
8564 SUBROUTINE fft_wrap_pw1pw2_r3d_c3d_rs_gs (pw1, pw2, debug)
8566 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw1
8567 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
8568 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
8570 CHARACTER(len=*),
PARAMETER :: routinen =
'fft_wrap_pw1pw2'
8572 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
8573 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
8574 INTEGER :: handle, handle2, my_pos, nrays, &
8576 INTEGER,
DIMENSION(:),
POINTER :: n
8578 REAL(kind=dp) :: norm
8580 CALL timeset(routinen, handle2)
8581 out_unit = cp_logger_get_default_io_unit()
8582 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
8583 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
8588 IF (
PRESENT(debug))
THEN
8595 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
8596 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
8597 cpabort(
"PW grids not compatible")
8599 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
8600 cpabort(
"PW grids have not compatible MPI groups")
8605 norm = 1.0_dp/pw1%pw_grid%ngpts
8607 n => pw1%pw_grid%npts
8609 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
8615 IF (test .AND. out_unit > 0)
THEN
8616 WRITE (out_unit,
'(A)')
" FFT Protocol "
8617 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
8618 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
8619 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
8620 WRITE (out_unit,
'(A,T66,E15.6)')
" scale factor", norm
8623 pw2%array = cmplx(pw1%array, 0.0_dp, kind=dp)
8625 CALL fft3d(fwfft, n, c_out, scale=norm, debug=test)
8627 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
8635 IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0)
THEN
8636 WRITE (out_unit,
'(A)')
" FFT Protocol "
8637 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
8638 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
8639 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
8640 WRITE (out_unit,
'(A,T66,E15.6)')
" scale factor", norm
8643 my_pos = pw1%pw_grid%para%my_pos
8644 nrays = pw1%pw_grid%para%nyzray(my_pos)
8645 grays => pw1%pw_grid%grays
8649 IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0)
THEN
8650 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
8653 CALL timestop(handle)
8654 CALL timestop(handle2)
8656 END SUBROUTINE fft_wrap_pw1pw2_r3d_c3d_rs_gs
8681 SUBROUTINE fft_wrap_pw1pw2_c1d_r3d_gs_rs (pw1, pw2, debug)
8683 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
8684 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw2
8685 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
8687 CHARACTER(len=*),
PARAMETER :: routinen =
'fft_wrap_pw1pw2'
8689 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
8690 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
8691 INTEGER :: handle, handle2, my_pos, nrays, &
8693 INTEGER,
DIMENSION(3) :: nloc
8694 #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
8695 LOGICAL :: use_pw_gpu
8697 INTEGER,
DIMENSION(:),
POINTER :: n
8699 REAL(kind=dp) :: norm
8701 CALL timeset(routinen, handle2)
8702 out_unit = cp_logger_get_default_io_unit()
8703 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
8704 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
8709 IF (
PRESENT(debug))
THEN
8716 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
8717 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
8718 cpabort(
"PW grids not compatible")
8720 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
8721 cpabort(
"PW grids have not compatible MPI groups")
8728 n => pw1%pw_grid%npts
8730 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
8736 IF (test .AND. out_unit > 0)
THEN
8737 WRITE (out_unit,
'(A)')
" FFT Protocol "
8738 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
8739 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
8740 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
8741 WRITE (out_unit,
'(A,T66,E15.6)')
" scale factor", norm
8744 #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
8745 CALL pw_gpu_c1dr3d_3d(pw1, pw2, scale=norm)
8746 #elif defined (__PW_FPGA)
8747 ALLOCATE (c_out(n(1), n(2), n(3)))
8750 IF (pw_fpga_init_bitstream(n) == 1)
THEN
8751 CALL pw_scatter_s_c1d_c3d(pw1, c_out)
8753 #if (__PW_FPGA_SP && __PW_FPGA)
8754 CALL pw_fpga_c1dr3d_3d_sp(n, c_out)
8756 CALL pw_fpga_c1dr3d_3d_dp(n, c_out)
8758 CALL zdscal(n(1)*n(2)*n(3), norm, c_out, 1)
8760 CALL pw_copy_from_array(pw2, c_out)
8762 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" PW_SCATTER : 3d -> 1d "
8763 CALL pw_scatter_s_c1d_c3d(pw1, c_out)
8765 CALL fft3d(bwfft, n, c_out, scale=norm, debug=test)
8767 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" REAL part "
8768 CALL pw_copy_from_array(pw2, c_out)
8772 ALLOCATE (c_out(n(1), n(2), n(3)))
8773 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" PW_SCATTER : 3d -> 1d "
8774 CALL pw_scatter_s_c1d_c3d(pw1, c_out)
8776 CALL fft3d(bwfft, n, c_out, scale=norm, debug=test)
8778 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" REAL part "
8779 CALL pw_copy_from_array(pw2, c_out)
8783 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
8791 IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0)
THEN
8792 WRITE (out_unit,
'(A)')
" FFT Protocol "
8793 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
8794 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
8795 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
8796 WRITE (out_unit,
'(A,T66,E15.6)')
" scale factor", norm
8799 my_pos = pw1%pw_grid%para%my_pos
8800 nrays = pw1%pw_grid%para%nyzray(my_pos)
8801 grays => pw1%pw_grid%grays
8803 #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
8805 use_pw_gpu = pw1%pw_grid%para%ray_distribution
8806 IF (use_pw_gpu)
THEN
8807 CALL pw_gpu_c1dr3d_3d_ps(pw1, pw2, scale=norm)
8811 IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0) &
8812 WRITE (out_unit,
'(A)')
" PW_SCATTER : 2d -> 1d "
8814 CALL pw_scatter_p_c1d (pw1, grays)
8815 nloc = pw2%pw_grid%npts_local
8816 ALLOCATE (c_in(nloc(1), nloc(2), nloc(3)))
8818 IF (pw1%pw_grid%para%ray_distribution)
THEN
8819 CALL fft3d(bwfft, n, c_in, grays, pw1%pw_grid%para%group, &
8820 pw1%pw_grid%para%rs_group, &
8821 pw1%pw_grid%para%yzp, pw1%pw_grid%para%nyzray, &
8822 pw1%pw_grid%para%bo, scale=norm, debug=test)
8824 CALL fft3d(bwfft, n, c_in, grays, pw1%pw_grid%para%rs_group, &
8825 pw1%pw_grid%para%bo, scale=norm, debug=test)
8828 IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0) &
8829 WRITE (out_unit,
'(A)')
" Real part "
8830 CALL pw_copy_from_array(pw2, c_in)
8832 #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
8837 IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0)
THEN
8838 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
8841 CALL timestop(handle)
8842 CALL timestop(handle2)
8844 END SUBROUTINE fft_wrap_pw1pw2_c1d_r3d_gs_rs
8857 SUBROUTINE pw_gather_s_c1d_r3d_2(pw1, pw2, scale)
8859 TYPE(pw_r3d_gs_type),
INTENT(IN) :: pw1
8860 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
8861 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
8863 CALL pw_gather_s_c1d_r3d (pw2, pw1%array, scale)
8865 END SUBROUTINE pw_gather_s_c1d_r3d_2
8876 SUBROUTINE pw_gather_s_c1d_r3d (pw, c, scale)
8878 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
8879 REAL(kind=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(IN) :: c
8880 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
8882 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_gather_s'
8884 INTEGER :: gpt, handle, l, m, n
8886 CALL timeset(routinen, handle)
8888 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
8889 ngpts =>
SIZE(pw%pw_grid%gsq), ghat => pw%pw_grid%g_hat)
8891 IF (
PRESENT(scale))
THEN
8894 l = mapl(ghat(1, gpt)) + 1
8895 m = mapm(ghat(2, gpt)) + 1
8896 n = mapn(ghat(3, gpt)) + 1
8897 pw%array(gpt) = scale* cmplx(c(l, m, n), 0.0_dp, kind=dp)
8903 l = mapl(ghat(1, gpt)) + 1
8904 m = mapm(ghat(2, gpt)) + 1
8905 n = mapn(ghat(3, gpt)) + 1
8906 pw%array(gpt) = cmplx(c(l, m, n), 0.0_dp, kind=dp)
8913 CALL timestop(handle)
8915 END SUBROUTINE pw_gather_s_c1d_r3d
8926 SUBROUTINE pw_scatter_s_c1d_r3d_2(pw1, pw2, scale)
8928 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
8929 TYPE(pw_r3d_gs_type),
INTENT(INOUT) :: pw2
8930 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
8932 CALL pw_scatter_s_c1d_r3d (pw1, pw2%array, scale)
8934 END SUBROUTINE pw_scatter_s_c1d_r3d_2
8945 SUBROUTINE pw_scatter_s_c1d_r3d (pw, c, scale)
8947 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
8948 REAL(kind=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(INOUT) :: c
8949 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
8951 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_scatter_s'
8953 INTEGER :: gpt, handle, l, m, n
8955 CALL timeset(routinen, handle)
8957 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
8958 ghat => pw%pw_grid%g_hat, ngpts =>
SIZE(pw%pw_grid%gsq))
8961 IF (.NOT.
PRESENT(scale)) c = 0.0_dp
8963 IF (
PRESENT(scale))
THEN
8966 l = mapl(ghat(1, gpt)) + 1
8967 m = mapm(ghat(2, gpt)) + 1
8968 n = mapn(ghat(3, gpt)) + 1
8969 c(l, m, n) = scale* real(pw%array(gpt), kind=dp)
8975 l = mapl(ghat(1, gpt)) + 1
8976 m = mapm(ghat(2, gpt)) + 1
8977 n = mapn(ghat(3, gpt)) + 1
8978 c(l, m, n) = real(pw%array(gpt), kind=dp)
8985 IF (pw%pw_grid%grid_span == halfspace)
THEN
8987 associate(mapl => pw%pw_grid%mapl%neg, mapm => pw%pw_grid%mapm%neg, mapn => pw%pw_grid%mapn%neg, &
8988 ghat => pw%pw_grid%g_hat, ngpts =>
SIZE(pw%pw_grid%gsq))
8990 IF (
PRESENT(scale))
THEN
8993 l = mapl(ghat(1, gpt)) + 1
8994 m = mapm(ghat(2, gpt)) + 1
8995 n = mapn(ghat(3, gpt)) + 1
8996 c(l, m, n) = scale*( real(pw%array(gpt), kind=dp))
9002 l = mapl(ghat(1, gpt)) + 1
9003 m = mapm(ghat(2, gpt)) + 1
9004 n = mapn(ghat(3, gpt)) + 1
9005 c(l, m, n) = ( real(pw%array(gpt), kind=dp))
9014 CALL timestop(handle)
9016 END SUBROUTINE pw_scatter_s_c1d_r3d
9038 SUBROUTINE fft_wrap_pw1pw2_c1d_c3d_gs_rs (pw1, pw2, debug)
9040 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
9041 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw2
9042 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
9044 CHARACTER(len=*),
PARAMETER :: routinen =
'fft_wrap_pw1pw2'
9046 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
9047 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
9048 INTEGER :: handle, handle2, my_pos, nrays, &
9050 INTEGER,
DIMENSION(:),
POINTER :: n
9052 REAL(kind=dp) :: norm
9054 CALL timeset(routinen, handle2)
9055 out_unit = cp_logger_get_default_io_unit()
9056 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
9057 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
9062 IF (
PRESENT(debug))
THEN
9069 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
9070 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
9071 cpabort(
"PW grids not compatible")
9073 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
9074 cpabort(
"PW grids have not compatible MPI groups")
9081 n => pw1%pw_grid%npts
9083 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
9089 IF (test .AND. out_unit > 0)
THEN
9090 WRITE (out_unit,
'(A)')
" FFT Protocol "
9091 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
9092 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
9093 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
9094 WRITE (out_unit,
'(A,T66,E15.6)')
" scale factor", norm
9098 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" PW_SCATTER : 3d -> 1d "
9099 CALL pw_scatter_s_c1d_c3d(pw1, c_out)
9100 CALL fft3d(bwfft, n, c_out, scale=norm, debug=test)
9102 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9110 IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0)
THEN
9111 WRITE (out_unit,
'(A)')
" FFT Protocol "
9112 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
9113 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
9114 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
9115 WRITE (out_unit,
'(A,T66,E15.6)')
" scale factor", norm
9118 my_pos = pw1%pw_grid%para%my_pos
9119 nrays = pw1%pw_grid%para%nyzray(my_pos)
9120 grays => pw1%pw_grid%grays
9123 IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0) &
9124 WRITE (out_unit,
'(A)')
" PW_SCATTER : 2d -> 1d "
9126 CALL pw_scatter_p_c1d (pw1, grays)
9129 IF (pw1%pw_grid%para%ray_distribution)
THEN
9130 CALL fft3d(bwfft, n, c_in, grays, pw1%pw_grid%para%group, &
9131 pw1%pw_grid%para%rs_group, &
9132 pw1%pw_grid%para%yzp, pw1%pw_grid%para%nyzray, &
9133 pw1%pw_grid%para%bo, scale=norm, debug=test)
9135 CALL fft3d(bwfft, n, c_in, grays, pw1%pw_grid%para%rs_group, &
9136 pw1%pw_grid%para%bo, scale=norm, debug=test)
9141 IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0)
THEN
9142 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9145 CALL timestop(handle)
9146 CALL timestop(handle2)
9148 END SUBROUTINE fft_wrap_pw1pw2_c1d_c3d_gs_rs
9161 SUBROUTINE pw_gather_s_c1d_c3d_2(pw1, pw2, scale)
9163 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
9164 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
9165 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
9167 CALL pw_gather_s_c1d_c3d (pw2, pw1%array, scale)
9169 END SUBROUTINE pw_gather_s_c1d_c3d_2
9180 SUBROUTINE pw_gather_s_c1d_c3d (pw, c, scale)
9182 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
9183 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(IN) :: c
9184 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
9186 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_gather_s'
9188 INTEGER :: gpt, handle, l, m, n
9190 CALL timeset(routinen, handle)
9192 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
9193 ngpts =>
SIZE(pw%pw_grid%gsq), ghat => pw%pw_grid%g_hat)
9195 IF (
PRESENT(scale))
THEN
9198 l = mapl(ghat(1, gpt)) + 1
9199 m = mapm(ghat(2, gpt)) + 1
9200 n = mapn(ghat(3, gpt)) + 1
9201 pw%array(gpt) = scale* c(l, m, n)
9207 l = mapl(ghat(1, gpt)) + 1
9208 m = mapm(ghat(2, gpt)) + 1
9209 n = mapn(ghat(3, gpt)) + 1
9210 pw%array(gpt) = c(l, m, n)
9217 CALL timestop(handle)
9219 END SUBROUTINE pw_gather_s_c1d_c3d
9230 SUBROUTINE pw_scatter_s_c1d_c3d_2(pw1, pw2, scale)
9232 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw1
9233 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
9234 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
9236 CALL pw_scatter_s_c1d_c3d (pw1, pw2%array, scale)
9238 END SUBROUTINE pw_scatter_s_c1d_c3d_2
9249 SUBROUTINE pw_scatter_s_c1d_c3d (pw, c, scale)
9251 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
9252 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
INTENT(INOUT) :: c
9253 REAL(kind=dp),
INTENT(IN),
OPTIONAL :: scale
9255 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_scatter_s'
9257 INTEGER :: gpt, handle, l, m, n
9259 CALL timeset(routinen, handle)
9261 associate(mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
9262 ghat => pw%pw_grid%g_hat, ngpts =>
SIZE(pw%pw_grid%gsq))
9265 IF (.NOT.
PRESENT(scale)) c = 0.0_dp
9267 IF (
PRESENT(scale))
THEN
9270 l = mapl(ghat(1, gpt)) + 1
9271 m = mapm(ghat(2, gpt)) + 1
9272 n = mapn(ghat(3, gpt)) + 1
9273 c(l, m, n) = scale* pw%array(gpt)
9279 l = mapl(ghat(1, gpt)) + 1
9280 m = mapm(ghat(2, gpt)) + 1
9281 n = mapn(ghat(3, gpt)) + 1
9282 c(l, m, n) = pw%array(gpt)
9289 IF (pw%pw_grid%grid_span == halfspace)
THEN
9291 associate(mapl => pw%pw_grid%mapl%neg, mapm => pw%pw_grid%mapm%neg, mapn => pw%pw_grid%mapn%neg, &
9292 ghat => pw%pw_grid%g_hat, ngpts =>
SIZE(pw%pw_grid%gsq))
9294 IF (
PRESENT(scale))
THEN
9297 l = mapl(ghat(1, gpt)) + 1
9298 m = mapm(ghat(2, gpt)) + 1
9299 n = mapn(ghat(3, gpt)) + 1
9300 c(l, m, n) = scale*conjg( pw%array(gpt))
9306 l = mapl(ghat(1, gpt)) + 1
9307 m = mapm(ghat(2, gpt)) + 1
9308 n = mapn(ghat(3, gpt)) + 1
9309 c(l, m, n) = conjg( pw%array(gpt))
9318 CALL timestop(handle)
9320 END SUBROUTINE pw_scatter_s_c1d_c3d
9342 SUBROUTINE fft_wrap_pw1pw2_c3d_r3d_gs_rs (pw1, pw2, debug)
9344 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
9345 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw2
9346 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
9348 CHARACTER(len=*),
PARAMETER :: routinen =
'fft_wrap_pw1pw2'
9350 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
9351 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
9352 INTEGER :: handle, handle2, my_pos, nrays, &
9354 INTEGER,
DIMENSION(:),
POINTER :: n
9356 REAL(kind=dp) :: norm
9358 CALL timeset(routinen, handle2)
9359 out_unit = cp_logger_get_default_io_unit()
9360 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
9361 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
9366 IF (
PRESENT(debug))
THEN
9373 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
9374 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
9375 cpabort(
"PW grids not compatible")
9377 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
9378 cpabort(
"PW grids have not compatible MPI groups")
9385 n => pw1%pw_grid%npts
9387 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
9393 IF (test .AND. out_unit > 0)
THEN
9394 WRITE (out_unit,
'(A)')
" FFT Protocol "
9395 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
9396 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
9397 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
9398 WRITE (out_unit,
'(A,T66,E15.6)')
" scale factor", norm
9402 ALLOCATE (c_out(n(1), n(2), n(3)))
9403 CALL fft3d(bwfft, n, c_in, c_out, scale=norm, debug=test)
9405 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" REAL part "
9406 pw2%array = real(c_out, kind=dp)
9409 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9417 IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0)
THEN
9418 WRITE (out_unit,
'(A)')
" FFT Protocol "
9419 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
9420 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
9421 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
9422 WRITE (out_unit,
'(A,T66,E15.6)')
" scale factor", norm
9425 my_pos = pw1%pw_grid%para%my_pos
9426 nrays = pw1%pw_grid%para%nyzray(my_pos)
9427 grays => pw1%pw_grid%grays
9431 IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0)
THEN
9432 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9435 CALL timestop(handle)
9436 CALL timestop(handle2)
9438 END SUBROUTINE fft_wrap_pw1pw2_c3d_r3d_gs_rs
9456 SUBROUTINE fft_wrap_pw1pw2_c3d_c1d_rs_gs (pw1, pw2, debug)
9458 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
9459 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw2
9460 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
9462 CHARACTER(len=*),
PARAMETER :: routinen =
'fft_wrap_pw1pw2'
9464 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
9465 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
9466 INTEGER :: handle, handle2, my_pos, nrays, &
9468 INTEGER,
DIMENSION(:),
POINTER :: n
9470 REAL(kind=dp) :: norm
9472 CALL timeset(routinen, handle2)
9473 out_unit = cp_logger_get_default_io_unit()
9474 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
9475 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
9480 IF (
PRESENT(debug))
THEN
9487 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
9488 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
9489 cpabort(
"PW grids not compatible")
9491 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
9492 cpabort(
"PW grids have not compatible MPI groups")
9497 norm = 1.0_dp/pw1%pw_grid%ngpts
9499 n => pw1%pw_grid%npts
9501 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
9507 IF (test .AND. out_unit > 0)
THEN
9508 WRITE (out_unit,
'(A)')
" FFT Protocol "
9509 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
9510 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
9511 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
9512 WRITE (out_unit,
'(A,T66,E15.6)')
" scale factor", norm
9516 ALLOCATE (c_out(n(1), n(2), n(3)))
9518 CALL fft3d(fwfft, n, c_in, c_out, scale=norm, debug=test)
9520 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" PW_GATHER : 3d -> 1d "
9521 CALL pw_gather_s_c1d_c3d(pw2, c_out)
9524 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9532 IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0)
THEN
9533 WRITE (out_unit,
'(A)')
" FFT Protocol "
9534 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
9535 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
9536 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
9537 WRITE (out_unit,
'(A,T66,E15.6)')
" scale factor", norm
9540 my_pos = pw1%pw_grid%para%my_pos
9541 nrays = pw1%pw_grid%para%nyzray(my_pos)
9542 grays => pw1%pw_grid%grays
9548 IF (pw1%pw_grid%para%ray_distribution)
THEN
9549 CALL fft3d(fwfft, n, c_in, grays, pw1%pw_grid%para%group, &
9550 pw1%pw_grid%para%rs_group, &
9551 pw1%pw_grid%para%yzp, pw1%pw_grid%para%nyzray, &
9552 pw1%pw_grid%para%bo, scale=norm, debug=test)
9554 CALL fft3d(fwfft, n, c_in, grays, pw1%pw_grid%para%rs_group, &
9555 pw1%pw_grid%para%bo, scale=norm, debug=test)
9558 IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0) &
9559 WRITE (out_unit,
'(A)')
" PW_GATHER : 2d -> 1d "
9560 CALL pw_gather_p_c1d (pw2, grays)
9563 IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0)
THEN
9564 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9567 CALL timestop(handle)
9568 CALL timestop(handle2)
9570 END SUBROUTINE fft_wrap_pw1pw2_c3d_c1d_rs_gs
9589 SUBROUTINE fft_wrap_pw1pw2_c3d_c3d_rs_gs (pw1, pw2, debug)
9591 TYPE(pw_c3d_rs_type),
INTENT(IN) :: pw1
9592 TYPE(pw_c3d_gs_type),
INTENT(INOUT) :: pw2
9593 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
9595 CHARACTER(len=*),
PARAMETER :: routinen =
'fft_wrap_pw1pw2'
9597 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
9598 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
9599 INTEGER :: handle, handle2, my_pos, nrays, &
9601 INTEGER,
DIMENSION(:),
POINTER :: n
9603 REAL(kind=dp) :: norm
9605 CALL timeset(routinen, handle2)
9606 out_unit = cp_logger_get_default_io_unit()
9607 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
9608 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
9613 IF (
PRESENT(debug))
THEN
9620 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
9621 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
9622 cpabort(
"PW grids not compatible")
9624 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
9625 cpabort(
"PW grids have not compatible MPI groups")
9630 norm = 1.0_dp/pw1%pw_grid%ngpts
9632 n => pw1%pw_grid%npts
9634 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
9640 IF (test .AND. out_unit > 0)
THEN
9641 WRITE (out_unit,
'(A)')
" FFT Protocol "
9642 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
9643 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
9644 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
9645 WRITE (out_unit,
'(A,T66,E15.6)')
" scale factor", norm
9650 CALL fft3d(fwfft, n, c_in, c_out, scale=norm, debug=test)
9652 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9660 IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0)
THEN
9661 WRITE (out_unit,
'(A)')
" FFT Protocol "
9662 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"FWFFT"
9663 WRITE (out_unit,
'(A,T72,A)')
" in space ",
"REALSPACE"
9664 WRITE (out_unit,
'(A,T72,A)')
" out space ",
"REALSPACE"
9665 WRITE (out_unit,
'(A,T66,E15.6)')
" scale factor", norm
9668 my_pos = pw1%pw_grid%para%my_pos
9669 nrays = pw1%pw_grid%para%nyzray(my_pos)
9670 grays => pw1%pw_grid%grays
9674 IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0)
THEN
9675 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9678 CALL timestop(handle)
9679 CALL timestop(handle2)
9681 END SUBROUTINE fft_wrap_pw1pw2_c3d_c3d_rs_gs
9696 SUBROUTINE fft_wrap_pw1pw2_c3d_c3d_gs_rs (pw1, pw2, debug)
9698 TYPE(pw_c3d_gs_type),
INTENT(IN) :: pw1
9699 TYPE(pw_c3d_rs_type),
INTENT(INOUT) :: pw2
9700 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
9702 CHARACTER(len=*),
PARAMETER :: routinen =
'fft_wrap_pw1pw2'
9704 COMPLEX(KIND=dp),
DIMENSION(:, :),
CONTIGUOUS,
POINTER :: grays
9705 COMPLEX(KIND=dp),
DIMENSION(:, :, :),
CONTIGUOUS,
POINTER :: c_in, c_out
9706 INTEGER :: handle, handle2, my_pos, nrays, &
9708 INTEGER,
DIMENSION(:),
POINTER :: n
9710 REAL(kind=dp) :: norm
9712 CALL timeset(routinen, handle2)
9713 out_unit = cp_logger_get_default_io_unit()
9714 CALL timeset(routinen//
"_"//trim(adjustl(cp_to_string( &
9715 ceiling(pw1%pw_grid%cutoff/10)*10))), handle)
9720 IF (
PRESENT(debug))
THEN
9727 IF (.NOT.
ASSOCIATED(pw1%pw_grid, pw2%pw_grid))
THEN
9728 IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol)
THEN
9729 cpabort(
"PW grids not compatible")
9731 IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group)
THEN
9732 cpabort(
"PW grids have not compatible MPI groups")
9739 n => pw1%pw_grid%npts
9741 IF (pw1%pw_grid%para%mode == pw_mode_local)
THEN
9747 IF (test .AND. out_unit > 0)
THEN
9748 WRITE (out_unit,
'(A)')
" FFT Protocol "
9749 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
9750 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
9751 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
9752 WRITE (out_unit,
'(A,T66,E15.6)')
" scale factor", norm
9757 CALL fft3d(bwfft, n, c_in, c_out, scale=norm, debug=test)
9759 IF (test .AND. out_unit > 0)
WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9767 IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0)
THEN
9768 WRITE (out_unit,
'(A)')
" FFT Protocol "
9769 WRITE (out_unit,
'(A,T76,A)')
" Transform direction ",
"BWFFT"
9770 WRITE (out_unit,
'(A,T66,A)')
" in space ",
"RECIPROCALSPACE"
9771 WRITE (out_unit,
'(A,T66,A)')
" out space ",
"RECIPROCALSPACE"
9772 WRITE (out_unit,
'(A,T66,E15.6)')
" scale factor", norm
9775 my_pos = pw1%pw_grid%para%my_pos
9776 nrays = pw1%pw_grid%para%nyzray(my_pos)
9777 grays => pw1%pw_grid%grays
9781 IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0)
THEN
9782 WRITE (out_unit,
'(A)')
" End of FFT Protocol "
9785 CALL timestop(handle)
9786 CALL timestop(handle2)
9788 END SUBROUTINE fft_wrap_pw1pw2_c3d_c3d_gs_rs
9807 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
9808 REAL(kind=dp),
INTENT(IN) :: omega
9810 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_gauss_damp'
9813 REAL(kind=dp) :: omega_2
9815 CALL timeset(routinen, handle)
9816 cpassert(omega >= 0)
9818 omega_2 = omega*omega
9819 omega_2 = 0.25_dp/omega_2
9822 pw%array = pw%array*exp(-pw%pw_grid%gsq*omega_2)
9825 CALL timestop(handle)
9838 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
9839 REAL(kind=dp),
INTENT(IN) :: omega
9841 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_log_deriv_gauss'
9844 REAL(kind=dp) :: omega_2
9846 CALL timeset(routinen, handle)
9847 cpassert(omega >= 0)
9849 omega_2 = omega*omega
9850 omega_2 = 0.25_dp/omega_2
9853 pw%array = pw%array*(1.0_dp + omega_2*pw%pw_grid%gsq)
9856 CALL timestop(handle)
9874 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
9875 REAL(kind=dp),
INTENT(IN) :: omega
9877 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_compl_gauss_damp'
9879 INTEGER :: cnt, handle, i
9880 REAL(kind=dp) :: omega_2, tmp
9882 CALL timeset(routinen, handle)
9884 omega_2 = omega*omega
9885 omega_2 = 0.25_dp/omega_2
9887 cnt =
SIZE(pw%array)
9891 tmp = -omega_2*pw%pw_grid%gsq(i)
9892 IF (abs(tmp) > 1.0e-5_dp)
THEN
9893 pw%array(i) = pw%array(i)*(1.0_dp - exp(tmp))
9895 pw%array(i) = pw%array(i)*(tmp + 0.5_dp*tmp*(tmp + (1.0_dp/3.0_dp)*tmp**2))
9900 CALL timestop(handle)
9912 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
9913 REAL(kind=dp),
INTENT(IN) :: omega
9915 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_log_deriv_compl_gauss'
9917 INTEGER :: handle, i
9918 REAL(kind=dp) :: omega_2, tmp
9920 CALL timeset(routinen, handle)
9922 omega_2 = omega*omega
9923 omega_2 = 0.25_dp/omega_2
9927 DO i = 1,
SIZE(pw%array)
9928 tmp = omega_2*pw%pw_grid%gsq(i)
9930 IF (abs(tmp) >= 0.003_dp)
THEN
9931 pw%array(i) = pw%array(i)*(1.0_dp - tmp*exp(-tmp)/(1.0_dp - exp(-tmp)))
9933 pw%array(i) = pw%array(i)*(0.5_dp*tmp - tmp**2/12.0_dp)
9938 CALL timestop(handle)
9959 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
9960 REAL(kind=dp),
INTENT(IN) :: omega, scale_coul, scale_long
9962 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_gauss_damp_mix'
9965 REAL(kind=dp) :: omega_2
9967 CALL timeset(routinen, handle)
9969 omega_2 = omega*omega
9970 omega_2 = 0.25_dp/omega_2
9973 pw%array = pw%array*(scale_coul + scale_long*exp(-pw%pw_grid%gsq*omega_2))
9976 CALL timestop(handle)
9991 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
9992 REAL(kind=dp),
INTENT(IN) :: omega, scale_coul, scale_long
9994 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_log_deriv_mix_cl'
9996 INTEGER :: handle, i
9997 REAL(kind=dp) :: omega_2, tmp
9999 CALL timeset(routinen, handle)
10001 omega_2 = omega*omega
10002 omega_2 = 0.25_dp/omega_2
10006 DO i = 1,
SIZE(pw%array)
10007 tmp = omega_2*pw%pw_grid%gsq(i)
10008 pw%array(i) = pw%array(i)*(1.0_dp + scale_long*tmp*exp(-tmp)/(scale_coul + scale_long*exp(-tmp)))
10012 CALL timestop(handle)
10031 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
10032 REAL(kind=dp),
INTENT(IN) :: rcutoff
10034 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_truncated'
10036 INTEGER :: handle, i
10037 REAL(kind=dp) :: tmp
10039 CALL timeset(routinen, handle)
10040 cpassert(rcutoff >= 0)
10043 DO i = 1,
SIZE(pw%array)
10044 tmp = sqrt(pw%pw_grid%gsq(i))*rcutoff
10045 IF (tmp >= 0.005_dp)
THEN
10046 pw%array(i) = pw%array(i)*(1.0_dp - cos(tmp))
10048 pw%array(i) = pw%array(i)*tmp**2/2.0_dp*(1.0 - tmp**2/12.0_dp)
10053 CALL timestop(handle)
10066 TYPE(pw_c1d_gs_type),
INTENT(IN) :: pw
10067 REAL(kind=dp),
INTENT(IN) :: rcutoff
10069 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_log_deriv_trunc'
10071 INTEGER :: handle, i
10072 REAL(kind=dp) :: rchalf, tmp
10074 CALL timeset(routinen, handle)
10075 cpassert(rcutoff >= 0)
10077 rchalf = 0.5_dp*rcutoff
10080 DO i = 1,
SIZE(pw%array)
10081 tmp = rchalf*sqrt(pw%pw_grid%gsq(i))
10083 IF (abs(tmp) >= 0.0001_dp)
THEN
10084 pw%array(i) = pw%array(i)*(1.0_dp - tmp/tan(tmp))
10086 pw%array(i) = pw%array(i)*tmp**2*(1.0_dp + tmp**2/15.0_dp)/3.0_dp
10091 CALL timestop(handle)
10107 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
10108 INTEGER,
DIMENSION(3),
INTENT(IN) :: n
10110 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_derive'
10112 COMPLEX(KIND=dp) :: im
10113 INTEGER :: handle, m
10115 CALL timeset(routinen, handle)
10118 cpabort(
"Nonnegative exponents are not supported!")
10121 im = cmplx(0.0_dp, 1.0_dp, kind=dp)**m
10123 IF (n(1) == 1)
THEN
10125 pw%array(:) = pw%array(:)*pw%pw_grid%g(1, :)
10127 ELSE IF (n(1) > 1)
THEN
10129 pw%array(:) = pw%array(:)*(pw%pw_grid%g(1, :)**n(1))
10133 IF (n(2) == 1)
THEN
10135 pw%array(:) = pw%array(:)*pw%pw_grid%g(2, :)
10137 ELSE IF (n(2) > 1)
THEN
10139 pw%array(:) = pw%array(:)*(pw%pw_grid%g(2, :)**n(2))
10143 IF (n(3) == 1)
THEN
10145 pw%array(:) = pw%array(:)*pw%pw_grid%g(3, :)
10147 ELSE IF (n(3) > 1)
THEN
10149 pw%array(:) = pw%array(:)*(pw%pw_grid%g(3, :)**n(3))
10155 IF (abs(real(im, kind=dp) - 1.0_dp) > 1.0e-10_dp)
THEN
10157 pw%array(:) = im*pw%array(:)
10161 CALL timestop(handle)
10176 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
10178 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_laplace'
10182 CALL timeset(routinen, handle)
10185 pw%array(:) = -pw%array(:)*pw%pw_grid%gsq(:)
10188 CALL timestop(handle)
10205 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw, pwdr2
10206 INTEGER,
INTENT(IN) :: i, j
10208 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_dr2'
10210 INTEGER :: cnt, handle, ig
10211 REAL(kind=dp) :: gg, o3
10213 CALL timeset(routinen, handle)
10217 cnt =
SIZE(pw%array)
10222 gg = pw%pw_grid%g(i, ig)**2 - o3*pw%pw_grid%gsq(ig)
10223 pwdr2%array(ig) = gg*pw%array(ig)
10229 pwdr2%array(ig) = pw%array(ig)*(pw%pw_grid%g(i, ig)*pw%pw_grid%g(j, ig))
10234 CALL timestop(handle)
10253 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw, pwdr2_gg
10254 INTEGER,
INTENT(IN) :: i, j
10256 INTEGER :: cnt, handle, ig
10257 REAL(kind=dp) :: gg, o3
10258 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_dr2_gg'
10260 CALL timeset(routinen, handle)
10264 cnt =
SIZE(pw%array)
10268 DO ig = pw%pw_grid%first_gne0, cnt
10269 gg = pw%pw_grid%g(i, ig)**2 - o3*pw%pw_grid%gsq(ig)
10270 pwdr2_gg%array(ig) = gg*pw%array(ig)/pw%pw_grid%gsq(ig)
10275 DO ig = pw%pw_grid%first_gne0, cnt
10276 pwdr2_gg%array(ig) = pw%array(ig)*(pw%pw_grid%g(i, ig)*pw%pw_grid%g(j, ig)) &
10277 /pw%pw_grid%gsq(ig)
10282 IF (pw%pw_grid%have_g0) pwdr2_gg%array(1) = 0.0_dp
10284 CALL timestop(handle)
10301 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: pw
10302 REAL(kind=dp),
INTENT(IN) :: ecut, sigma
10304 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_smoothing'
10306 INTEGER :: cnt, handle, ig
10307 REAL(kind=dp) :: arg, f
10309 CALL timeset(routinen, handle)
10311 cnt =
SIZE(pw%array)
10315 arg = (ecut - pw%pw_grid%gsq(ig))/sigma
10316 f = exp(arg)/(1 + exp(arg))
10317 pw%array(ig) = f*pw%array(ig)
10321 CALL timestop(handle)
10331 ELEMENTAL FUNCTION pw_compatible(grida, gridb)
RESULT(compat)
10333 TYPE(pw_grid_type),
INTENT(IN) :: grida, gridb
10338 IF (grida%id_nr == gridb%id_nr)
THEN
10340 ELSE IF (grida%reference == gridb%id_nr)
THEN
10342 ELSE IF (gridb%reference == grida%id_nr)
THEN
10346 END FUNCTION pw_compatible
10359 TYPE(pw_c1d_gs_type),
INTENT(INOUT) :: sf
10360 REAL(kind=dp),
DIMENSION(:),
INTENT(IN) :: r
10362 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_structure_factor'
10364 INTEGER :: cnt, handle, ig
10365 REAL(kind=dp) :: arg
10367 CALL timeset(routinen, handle)
10369 cnt =
SIZE(sf%array)
10373 arg = dot_product(sf%pw_grid%g(:, ig), r)
10374 sf%array(ig) = cmplx(cos(arg), -sin(arg), kind=dp)
10378 CALL timestop(handle)
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Defines the basic variable types.
integer, parameter, public dp
Machine interface based on Fortran 2003 and POSIX.
subroutine, public m_memory(mem)
Returns the total amount of memory [bytes] in use, if known, zero otherwise.
Definition of mathematical constants and functions.
complex(kind=dp), parameter, public z_zero
subroutine, public pw_copy_match(pw1, pw2)
copy a pw type variable
integer function, public pw_fpga_init_bitstream(n)
Invoke the pw_fpga_check_bitstream C function passing the path to the data dir.
subroutine, public pw_fpga_r3dc1d_3d_dp(n, c_out)
perform an in-place double precision fft3d on the FPGA
subroutine, public pw_fpga_r3dc1d_3d_sp(n, c_out)
perform an in-place single precision fft3d on the FPGA
subroutine, public pw_fpga_c1dr3d_3d_dp(n, c_out)
perform an in-place double precision inverse fft3d on the FPGA
subroutine, public pw_fpga_c1dr3d_3d_sp(n, c_out)
perform an in-place single precision inverse fft3d on the FPGA
subroutine, public pw_gpu_c1dr3d_3d(pw1, pw2, scale)
perform an scatter followed by a fft on the gpu
subroutine, public pw_gpu_c1dr3d_3d_ps(pw1, pw2, scale)
perform an parallel scatter followed by a fft on the gpu
subroutine, public pw_gpu_r3dc1d_3d(pw1, pw2, scale)
perform an fft followed by a gather on the gpu
subroutine, public pw_gpu_r3dc1d_3d_ps(pw1, pw2, scale)
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.