Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \note
10 : !> If parallel mode is distributed certain combination of
11 : !> "in_use" and "in_space" can not be used.
12 : !> For performance reasons it would be better to have the loops
13 : !> over g-vectors in the gather/scatter routines in new subprograms
14 : !> with the actual arrays (also the addressing) in the parameter list
15 : !> \par History
16 : !> JGH (29-Dec-2000) : Changes for parallel use
17 : !> JGH (13-Mar-2001) : added timing calls
18 : !> JGH (26-Feb-2003) : OpenMP enabled
19 : !> JGH (17-Nov-2007) : Removed mass arrays
20 : !> JGH (01-Dec-2007) : Removed and renamed routines
21 : !> JGH (04-Jul-2019) : added pw_multiply routine
22 : !> 03.2008 [tlaino] : Splitting pw_types into pw_types and pw_methods
23 : !> \author apsi
24 : ! **************************************************************************************************
25 : MODULE pw_methods
26 : #:include "pw_types.fypp"
27 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit, &
28 : cp_to_string
29 : USE fft_tools, ONLY: BWFFT, &
30 : FWFFT, &
31 : fft3d
32 : USE kahan_sum, ONLY: accurate_dot_product, &
33 : accurate_sum
34 : USE kinds, ONLY: dp
35 : USE machine, ONLY: m_memory
36 : USE mathconstants, ONLY: z_zero
37 : USE pw_copy_all, ONLY: pw_copy_match
38 : USE pw_fpga, ONLY: pw_fpga_c1dr3d_3d_dp, &
39 : pw_fpga_c1dr3d_3d_sp, &
40 : pw_fpga_init_bitstream, &
41 : pw_fpga_r3dc1d_3d_dp, &
42 : pw_fpga_r3dc1d_3d_sp
43 : USE pw_gpu, ONLY: pw_gpu_c1dr3d_3d, &
44 : pw_gpu_c1dr3d_3d_ps, &
45 : pw_gpu_r3dc1d_3d, &
46 : pw_gpu_r3dc1d_3d_ps
47 : USE pw_grid_types, ONLY: HALFSPACE, &
48 : PW_MODE_DISTRIBUTED, &
49 : PW_MODE_LOCAL, &
50 : pw_grid_type
51 : #:for space in pw_spaces
52 : #:for kind in pw_kinds
53 : USE pw_types, ONLY: pw_${kind}$_${space}$_type
54 : #:endfor
55 : #:endfor
56 : #include "../base/base_uses.f90"
57 :
58 : IMPLICIT NONE
59 :
60 : PRIVATE
61 :
62 : PUBLIC :: pw_zero, pw_structure_factor, pw_smoothing
63 : PUBLIC :: pw_copy, pw_axpy, pw_transfer, pw_scale
64 : PUBLIC :: pw_gauss_damp, pw_compl_gauss_damp, pw_derive, pw_laplace, pw_dr2, pw_write, pw_multiply
65 : PUBLIC :: pw_log_deriv_gauss, pw_log_deriv_compl_gauss, pw_log_deriv_mix_cl, pw_log_deriv_trunc
66 : PUBLIC :: pw_gauss_damp_mix, pw_multiply_with
67 : PUBLIC :: pw_integral_ab, pw_integral_a2b
68 : PUBLIC :: pw_dr2_gg, pw_integrate_function
69 : PUBLIC :: pw_set, pw_truncated
70 : PUBLIC :: pw_scatter, pw_gather
71 : PUBLIC :: pw_copy_to_array, pw_copy_from_array
72 :
73 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pw_methods'
74 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
75 : INTEGER, PARAMETER, PUBLIC :: do_accurate_sum = 0, &
76 : do_standard_sum = 1
77 :
78 : INTERFACE pw_zero
79 : #:for space in pw_spaces
80 : #:for kind in pw_kinds
81 : MODULE PROCEDURE pw_zero_${kind}$_${space}$
82 : #:endfor
83 : #:endfor
84 : END INTERFACE
85 :
86 : INTERFACE pw_scale
87 : #:for space in pw_spaces
88 : #:for kind in pw_kinds
89 : MODULE PROCEDURE pw_scale_${kind}$_${space}$
90 : #:endfor
91 : #:endfor
92 : END INTERFACE
93 :
94 : INTERFACE pw_write
95 : #:for space in pw_spaces
96 : #:for kind in pw_kinds
97 : MODULE PROCEDURE pw_write_${kind}$_${space}$
98 : #:endfor
99 : #:endfor
100 : END INTERFACE
101 :
102 : INTERFACE pw_integrate_function
103 : #:for space in pw_spaces
104 : #:for kind in pw_kinds
105 : MODULE PROCEDURE pw_integrate_function_${kind}$_${space}$
106 : #:endfor
107 : #:endfor
108 : END INTERFACE
109 :
110 : INTERFACE pw_set
111 : #:for space in pw_spaces
112 : #:for kind in pw_kinds
113 : MODULE PROCEDURE pw_set_value_${kind}$_${space}$
114 : MODULE PROCEDURE pw_zero_${kind}$_${space}$
115 : #:endfor
116 : #:endfor
117 : END INTERFACE
118 :
119 : INTERFACE pw_copy
120 : #:for space in pw_spaces
121 : #:for kind, kind2 in pw_kinds2_sameD
122 : MODULE PROCEDURE pw_copy_${kind}$_${kind2}$_${space}$
123 : #:endfor
124 : #:endfor
125 : END INTERFACE
126 :
127 : INTERFACE pw_axpy
128 : #:for space in pw_spaces
129 : #:for kind, kind2 in pw_kinds2_sameD
130 : MODULE PROCEDURE pw_axpy_${kind}$_${kind2}$_${space}$
131 : #:endfor
132 : #:endfor
133 : END INTERFACE
134 :
135 : INTERFACE pw_multiply
136 : #:for space in pw_spaces
137 : #:for kind, kind2 in pw_kinds2_sameD
138 : MODULE PROCEDURE pw_multiply_${kind}$_${kind2}$_${space}$
139 : #:endfor
140 : #:endfor
141 : END INTERFACE
142 :
143 : INTERFACE pw_multiply_with
144 : #:for space in pw_spaces
145 : #:for kind, kind2 in pw_kinds2_sameD
146 : MODULE PROCEDURE pw_multiply_with_${kind}$_${kind2}$_${space}$
147 : #:endfor
148 : #:endfor
149 : END INTERFACE
150 :
151 : INTERFACE pw_integral_ab
152 : #:for space in pw_spaces
153 : #:for kind, kind2 in pw_kinds2_sameD
154 : MODULE PROCEDURE pw_integral_ab_${kind}$_${kind2}$_${space}$
155 : #:endfor
156 : #:endfor
157 : END INTERFACE
158 :
159 : INTERFACE pw_integral_a2b
160 : #:for kind, kind2 in pw_kinds2_sameD
161 : #:if kind[1]=="1"
162 : MODULE PROCEDURE pw_integral_a2b_${kind}$_${kind2}$
163 : #:endif
164 : #:endfor
165 : END INTERFACE
166 :
167 : INTERFACE pw_gather
168 : #:for kind in pw_kinds
169 : #:if kind[1]=="1"
170 : MODULE PROCEDURE pw_gather_p_${kind}$
171 : #:endif
172 : #:endfor
173 : #:for kind, kind2 in pw_kinds2
174 : #:if kind[1]=="1" and kind2[1]=="3"
175 : MODULE PROCEDURE pw_gather_s_${kind}$_${kind2}$
176 : #:endif
177 : #:endfor
178 : END INTERFACE
179 :
180 : INTERFACE pw_scatter
181 : #:for kind in pw_kinds
182 : #:if kind[1]=="1"
183 : MODULE PROCEDURE pw_scatter_p_${kind}$
184 : #:endif
185 : #:endfor
186 : #:for kind, kind2 in pw_kinds2
187 : #:if kind[1]=="1" and kind2[1]=="3"
188 : MODULE PROCEDURE pw_scatter_s_${kind}$_${kind2}$
189 : #:endif
190 : #:endfor
191 : END INTERFACE
192 :
193 : INTERFACE pw_copy_to_array
194 : #:for space in pw_spaces
195 : #:for kind, kind2 in pw_kinds2_sameD
196 : MODULE PROCEDURE pw_copy_to_array_${kind}$_${kind2}$_${space}$
197 : #:endfor
198 : #:endfor
199 : END INTERFACE
200 :
201 : INTERFACE pw_copy_from_array
202 : #:for space in pw_spaces
203 : #:for kind, kind2 in pw_kinds2_sameD
204 : MODULE PROCEDURE pw_copy_from_array_${kind}$_${kind2}$_${space}$
205 : #:endfor
206 : #:endfor
207 : END INTERFACE
208 :
209 : INTERFACE pw_transfer
210 : #:for kind, kind2 in pw_kinds2
211 : #:if kind[1]=="1" and kind2[1]=="3"
212 : MODULE PROCEDURE pw_gather_s_${kind}$_${kind2}$_2
213 : MODULE PROCEDURE pw_scatter_s_${kind}$_${kind2}$_2
214 : #:endif
215 : #:for space in pw_spaces
216 : #:if kind[1]==kind2[1]
217 : MODULE PROCEDURE pw_copy_${kind}$_${kind2}$_${space}$
218 : #:endif
219 : #:endfor
220 : #:if kind2[0]=="c" and kind[1]=="3"
221 : MODULE PROCEDURE fft_wrap_pw1pw2_${kind}$_${kind2}$_rs_gs
222 : #:endif
223 : #:if kind[0]=="c" and kind2[1]=="3"
224 : MODULE PROCEDURE fft_wrap_pw1pw2_${kind}$_${kind2}$_gs_rs
225 : #:endif
226 : #:endfor
227 : END INTERFACE
228 :
229 : CONTAINS
230 : #:for kind, type in pw_list
231 : #:for space in pw_spaces
232 : ! **************************************************************************************************
233 : !> \brief Set values of a pw type to zero
234 : !> \param pw ...
235 : !> \par History
236 : !> none
237 : !> \author apsi
238 : ! **************************************************************************************************
239 1108706 : SUBROUTINE pw_zero_${kind}$_${space}$ (pw)
240 :
241 : TYPE(pw_${kind}$_${space}$_type), INTENT(INOUT) :: pw
242 :
243 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_zero'
244 :
245 : INTEGER :: handle
246 :
247 1108706 : CALL timeset(routineN, handle)
248 :
249 : #if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
250 1108706 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw)
251 : pw%array = ${type2type("0.0_dp", "r3d", kind)}$
252 : !$OMP END PARALLEL WORKSHARE
253 : #else
254 : pw%array = ${type2type("0.0_dp", "r3d", kind)}$
255 : #endif
256 :
257 1108706 : CALL timestop(handle)
258 :
259 1108706 : END SUBROUTINE pw_zero_${kind}$_${space}$
260 :
261 : ! **************************************************************************************************
262 : !> \brief multiplies pw coeffs with a number
263 : !> \param pw ...
264 : !> \param a ...
265 : !> \par History
266 : !> 11.2004 created [Joost VandeVondele]
267 : ! **************************************************************************************************
268 458802 : SUBROUTINE pw_scale_${kind}$_${space}$ (pw, a)
269 :
270 : TYPE(pw_${kind}$_${space}$_type), INTENT(INOUT) :: pw
271 : REAL(KIND=dp), INTENT(IN) :: a
272 :
273 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_scale'
274 :
275 : INTEGER :: handle
276 :
277 458802 : CALL timeset(routineN, handle)
278 :
279 458802 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(A, pw)
280 : pw%array = a*pw%array
281 : !$OMP END PARALLEL WORKSHARE
282 :
283 458802 : CALL timestop(handle)
284 :
285 458802 : END SUBROUTINE pw_scale_${kind}$_${space}$
286 :
287 : ! **************************************************************************************************
288 : !> \brief writes a small description of the actual grid
289 : !> (change to output the data as cube file, maybe with an
290 : !> optional long_description arg?)
291 : !> \param pw the pw data to output
292 : !> \param unit_nr the unit to output to
293 : !> \par History
294 : !> 08.2002 created [fawzi]
295 : !> \author Fawzi Mohamed
296 : ! **************************************************************************************************
297 0 : SUBROUTINE pw_write_${kind}$_${space}$ (pw, unit_nr)
298 :
299 : TYPE(pw_${kind}$_${space}$_type), INTENT(in) :: pw
300 : INTEGER, INTENT(in) :: unit_nr
301 :
302 : INTEGER :: iostatus
303 :
304 0 : WRITE (unit=unit_nr, fmt="('<pw>:{ ')", iostat=iostatus)
305 :
306 0 : WRITE (unit=unit_nr, fmt="(a)", iostat=iostatus) " in_use=${kind}$"
307 0 : IF (ASSOCIATED(pw%array)) THEN
308 : #:if kind[1]=='1'
309 : WRITE (unit=unit_nr, fmt="(' array=<${kind[0]}$(',i8,':',i8,')>')") &
310 0 : LBOUND(pw%array, 1), UBOUND(pw%array, 1)
311 : #:else
312 : WRITE (unit=unit_nr, fmt="(' array=<${kind[0]}$(',i8,':',i8,',',i8,':',i8,',',i8,':',i8,')>')") &
313 0 : LBOUND(pw%array, 1), UBOUND(pw%array, 1), LBOUND(pw%array, 2), UBOUND(pw%array, 2), &
314 0 : LBOUND(pw%array, 3), UBOUND(pw%array, 3)
315 : #:endif
316 : ELSE
317 0 : WRITE (unit=unit_nr, fmt="(' array=*null*')")
318 : END IF
319 :
320 0 : END SUBROUTINE pw_write_${kind}$_${space}$
321 :
322 : ! **************************************************************************************************
323 : !> \brief ...
324 : !> \param fun ...
325 : !> \param isign ...
326 : !> \param oprt ...
327 : !> \return ...
328 : ! **************************************************************************************************
329 426460 : FUNCTION pw_integrate_function_${kind}$_${space}$ (fun, isign, oprt) RESULT(total_fun)
330 :
331 : TYPE(pw_${kind}$_${space}$_type), INTENT(IN) :: fun
332 : INTEGER, INTENT(IN), OPTIONAL :: isign
333 : CHARACTER(len=*), INTENT(IN), OPTIONAL :: oprt
334 : REAL(KIND=dp) :: total_fun
335 :
336 : INTEGER :: iop
337 :
338 426460 : iop = 0
339 :
340 426460 : IF (PRESENT(oprt)) THEN
341 : SELECT CASE (oprt)
342 : CASE ("ABS", "abs")
343 0 : iop = 1
344 : CASE DEFAULT
345 2028 : CPABORT("Unknown operator")
346 : END SELECT
347 : END IF
348 :
349 102823 : total_fun = 0.0_dp
350 :
351 : #:if space == "rs"
352 : ! do reduction using maximum accuracy
353 : IF (iop == 1) THEN
354 106902165 : total_fun = fun%pw_grid%dvol*accurate_sum(ABS(fun%array))
355 : ELSE
356 321609 : total_fun = fun%pw_grid%dvol*accurate_sum(${type2type("fun%array", kind, "r1d")}$)
357 : END IF
358 : #:elif space == "gs"
359 : IF (iop == 1) &
360 0 : CPABORT("Operator ABS not implemented")
361 : #:if kind[1:]=="1d"
362 102823 : IF (fun%pw_grid%have_g0) total_fun = fun%pw_grid%vol*${type2type("fun%array(1)", kind, "r1d")}$
363 : #:else
364 0 : CPABORT("Reciprocal space integration for 3D grids not implemented!")
365 : #:endif
366 : #:else
367 : CPABORT("No space defined")
368 : #:endif
369 :
370 426460 : IF (fun%pw_grid%para%mode /= PW_MODE_LOCAL) THEN
371 421020 : CALL fun%pw_grid%para%group%sum(total_fun)
372 : END IF
373 :
374 426460 : IF (PRESENT(isign)) THEN
375 309036 : total_fun = total_fun*SIGN(1._dp, REAL(isign, dp))
376 : END IF
377 :
378 426460 : END FUNCTION pw_integrate_function_${kind}$_${space}$
379 :
380 : ! **************************************************************************************************
381 : !> \brief ...
382 : !> \param pw ...
383 : !> \param value ...
384 : ! **************************************************************************************************
385 118 : SUBROUTINE pw_set_value_${kind}$_${space}$ (pw, value)
386 : TYPE(pw_${kind}$_${space}$_type), INTENT(IN) :: pw
387 : REAL(KIND=dp), INTENT(IN) :: value
388 :
389 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_set_value'
390 :
391 : INTEGER :: handle
392 :
393 118 : CALL timeset(routineN, handle)
394 :
395 118 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw,value)
396 : pw%array = ${type2type("value", "r3d", kind)}$
397 : !$OMP END PARALLEL WORKSHARE
398 :
399 118 : CALL timestop(handle)
400 :
401 118 : END SUBROUTINE pw_set_value_${kind}$_${space}$
402 : #:endfor
403 :
404 : #:if kind[1]=="1"
405 : ! **************************************************************************************************
406 : !> \brief ...
407 : !> \param pw ...
408 : !> \param c ...
409 : !> \param scale ...
410 : ! **************************************************************************************************
411 1307046 : SUBROUTINE pw_gather_p_${kind}$ (pw, c, scale)
412 :
413 : TYPE(pw_${kind}$_gs_type), INTENT(INOUT) :: pw
414 : COMPLEX(KIND=dp), DIMENSION(:, :), CONTIGUOUS, INTENT(IN) :: c
415 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: scale
416 :
417 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_gather_p'
418 :
419 : INTEGER :: gpt, handle, l, m, mn, n
420 :
421 1307046 : CALL timeset(routineN, handle)
422 :
423 1307046 : IF (pw%pw_grid%para%mode /= PW_MODE_DISTRIBUTED) THEN
424 0 : CPABORT("This grid type is not distributed")
425 : END IF
426 :
427 : ASSOCIATE (mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
428 : ngpts => SIZE(pw%pw_grid%gsq), ghat => pw%pw_grid%g_hat, yzq => pw%pw_grid%para%yzq)
429 :
430 1307046 : IF (PRESENT(scale)) THEN
431 : !$OMP PARALLEL DO DEFAULT(NONE) &
432 : !$OMP PRIVATE(l, m, mn, n) &
433 0 : !$OMP SHARED(c, pw, scale)
434 : DO gpt = 1, ngpts
435 : l = mapl(ghat(1, gpt)) + 1
436 : m = mapm(ghat(2, gpt)) + 1
437 : n = mapn(ghat(3, gpt)) + 1
438 : mn = yzq(m, n)
439 : pw%array(gpt) = scale*${type2type("c(l, mn)", "c3d", kind)}$
440 : END DO
441 : !$OMP END PARALLEL DO
442 : ELSE
443 : !$OMP PARALLEL DO DEFAULT(NONE) &
444 : !$OMP PRIVATE(l, m, mn, n) &
445 1307046 : !$OMP SHARED(c, pw)
446 : DO gpt = 1, ngpts
447 : l = mapl(ghat(1, gpt)) + 1
448 : m = mapm(ghat(2, gpt)) + 1
449 : n = mapn(ghat(3, gpt)) + 1
450 : mn = yzq(m, n)
451 : pw%array(gpt) = ${type2type("c(l, mn)", "c3d", kind)}$
452 : END DO
453 : !$OMP END PARALLEL DO
454 : END IF
455 :
456 : END ASSOCIATE
457 :
458 1307046 : CALL timestop(handle)
459 :
460 1307046 : END SUBROUTINE pw_gather_p_${kind}$
461 :
462 : ! **************************************************************************************************
463 : !> \brief ...
464 : !> \param pw ...
465 : !> \param c ...
466 : !> \param scale ...
467 : ! **************************************************************************************************
468 1336308 : SUBROUTINE pw_scatter_p_${kind}$ (pw, c, scale)
469 : TYPE(pw_${kind}$_gs_type), INTENT(IN) :: pw
470 : COMPLEX(KIND=dp), DIMENSION(:, :), CONTIGUOUS, INTENT(INOUT) :: c
471 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: scale
472 :
473 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_scatter_p'
474 :
475 : INTEGER :: gpt, handle, l, m, mn, n
476 :
477 1336308 : CALL timeset(routineN, handle)
478 :
479 1336308 : IF (pw%pw_grid%para%mode /= PW_MODE_DISTRIBUTED) THEN
480 0 : CPABORT("This grid type is not distributed")
481 : END IF
482 :
483 : ASSOCIATE (mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
484 : ghat => pw%pw_grid%g_hat, yzq => pw%pw_grid%para%yzq, ngpts => SIZE(pw%pw_grid%gsq))
485 :
486 41333396938 : IF (.NOT. PRESENT(scale)) c = z_zero
487 :
488 1336308 : IF (PRESENT(scale)) THEN
489 : !$OMP PARALLEL DO DEFAULT(NONE) &
490 : !$OMP PRIVATE(l, m, mn, n) &
491 0 : !$OMP SHARED(c, pw, scale)
492 : DO gpt = 1, ngpts
493 : l = mapl(ghat(1, gpt)) + 1
494 : m = mapm(ghat(2, gpt)) + 1
495 : n = mapn(ghat(3, gpt)) + 1
496 : mn = yzq(m, n)
497 : c(l, mn) = ${type2type("scale*pw%array(gpt)", kind, "c3d")}$
498 : END DO
499 : !$OMP END PARALLEL DO
500 : ELSE
501 : !$OMP PARALLEL DO DEFAULT(NONE) &
502 : !$OMP PRIVATE(l, m, mn, n) &
503 1336308 : !$OMP SHARED(c, pw)
504 : DO gpt = 1, ngpts
505 : l = mapl(ghat(1, gpt)) + 1
506 : m = mapm(ghat(2, gpt)) + 1
507 : n = mapn(ghat(3, gpt)) + 1
508 : mn = yzq(m, n)
509 : c(l, mn) = ${type2type("pw%array(gpt)", kind, "c3d")}$
510 : END DO
511 : !$OMP END PARALLEL DO
512 : END IF
513 :
514 : END ASSOCIATE
515 :
516 1336308 : IF (pw%pw_grid%grid_span == HALFSPACE) THEN
517 :
518 : ASSOCIATE (mapm => pw%pw_grid%mapm%neg, mapn => pw%pw_grid%mapn%neg, mapl => pw%pw_grid%mapl%neg, &
519 : ghat => pw%pw_grid%g_hat, ngpts => SIZE(pw%pw_grid%gsq), yzq => pw%pw_grid%para%yzq)
520 :
521 119554 : IF (PRESENT(scale)) THEN
522 : !$OMP PARALLEL DO DEFAULT(NONE) &
523 : !$OMP PRIVATE(l, m, mn, n) &
524 0 : !$OMP SHARED(c, pw, scale)
525 : DO gpt = 1, ngpts
526 : l = mapl(ghat(1, gpt)) + 1
527 : m = mapm(ghat(2, gpt)) + 1
528 : n = mapn(ghat(3, gpt)) + 1
529 : mn = yzq(m, n)
530 : c(l, mn) = scale*#{if kind[0]=="c"}#CONJG#{endif}#(${type2type("pw%array(gpt)", kind, "c3d")}$)
531 : END DO
532 : !$OMP END PARALLEL DO
533 : ELSE
534 : !$OMP PARALLEL DO DEFAULT(NONE) &
535 : !$OMP PRIVATE(l, m, mn, n) &
536 119554 : !$OMP SHARED(c, pw)
537 : DO gpt = 1, ngpts
538 : l = mapl(ghat(1, gpt)) + 1
539 : m = mapm(ghat(2, gpt)) + 1
540 : n = mapn(ghat(3, gpt)) + 1
541 : mn = yzq(m, n)
542 : c(l, mn) = #{if kind[0]=="c"}#CONJG#{endif}#(${type2type("pw%array(gpt)", kind, "c3d")}$)
543 : END DO
544 : !$OMP END PARALLEL DO
545 : END IF
546 : END ASSOCIATE
547 : END IF
548 :
549 1336308 : CALL timestop(handle)
550 :
551 1336308 : END SUBROUTINE pw_scatter_p_${kind}$
552 : #:endif
553 : #:endfor
554 :
555 : #:for kind, type, kind2, type2 in pw_list2_sameD
556 : #:for space in pw_spaces
557 : ! **************************************************************************************************
558 : !> \brief copy a pw type variable
559 : !> \param pw1 ...
560 : !> \param pw2 ...
561 : !> \par History
562 : !> JGH (7-Mar-2001) : check for pw_grid %id_nr, allow copy if
563 : !> in_use == COMPLEXDATA1D and in_space == RECIPROCALSPACE
564 : !> JGH (21-Feb-2003) : Code for generalized reference grids
565 : !> \author apsi
566 : !> \note
567 : !> Currently only copying of respective types allowed,
568 : !> in order to avoid errors
569 : ! **************************************************************************************************
570 3032124 : SUBROUTINE pw_copy_${kind}$_${kind2}$_${space}$ (pw1, pw2)
571 :
572 : TYPE(pw_${kind}$_${space}$_type), INTENT(IN) :: pw1
573 : TYPE(pw_${kind2}$_${space}$_type), INTENT(INOUT) :: pw2
574 :
575 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_copy'
576 :
577 : INTEGER :: handle
578 : #:if kind[1:]=='1d'
579 : INTEGER :: i, j, ng, ng1, ng2, ns
580 : #:endif
581 :
582 3032124 : CALL timeset(routineN, handle)
583 :
584 3032124 : IF (pw1%pw_grid%spherical .neqv. pw2%pw_grid%spherical) &
585 0 : CPABORT("Both grids must be either spherical or non-spherical!")
586 : #:if space != "gs"
587 615401 : IF (pw1%pw_grid%spherical) &
588 0 : CPABORT("Spherical grids only exist in reciprocal space!")
589 : #:endif
590 :
591 : #:if kind[1]=='1'
592 2416723 : IF (.NOT. ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) THEN
593 656880 : IF (pw1%pw_grid%spherical .AND. pw2%pw_grid%spherical) THEN
594 0 : IF (pw_compatible(pw1%pw_grid, pw2%pw_grid)) THEN
595 0 : ng1 = SIZE(pw1%array)
596 0 : ng2 = SIZE(pw2%array)
597 0 : ng = MIN(ng1, ng2)
598 0 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(ng, pw1, pw2)
599 : pw2%array(1:ng) = ${type2type("pw1%array(1:ng)", kind, kind2)}$
600 : !$OMP END PARALLEL WORKSHARE
601 0 : IF (ng2 > ng) THEN
602 0 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(ng, ng2, pw2)
603 : pw2%array(ng + 1:ng2) = ${type2type("0.0_dp", "r3d", kind2)}$
604 : !$OMP END PARALLEL WORKSHARE
605 : END IF
606 : ELSE
607 0 : CPABORT("Copies between spherical grids require compatible grids!")
608 : END IF
609 : ELSE
610 656880 : ng1 = SIZE(pw1%array)
611 656880 : ng2 = SIZE(pw2%array)
612 656880 : ns = 2*MAX(ng1, ng2)
613 :
614 656880 : IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference)) THEN
615 656512 : IF (ng1 >= ng2) THEN
616 : !$OMP PARALLEL DO DEFAULT(NONE) &
617 : !$OMP PRIVATE(i,j) &
618 656512 : !$OMP SHARED(ng2, pw1, pw2)
619 : DO i = 1, ng2
620 : j = pw2%pw_grid%gidx(i)
621 : pw2%array(i) = ${type2type("pw1%array(j)", kind, kind2)}$
622 : END DO
623 : !$OMP END PARALLEL DO
624 : ELSE
625 0 : CALL pw_zero(pw2)
626 : !$OMP PARALLEL DO DEFAULT(NONE) &
627 : !$OMP PRIVATE(i,j) &
628 0 : !$OMP SHARED(ng1, pw1, pw2)
629 : DO i = 1, ng1
630 : j = pw2%pw_grid%gidx(i)
631 : pw2%array(j) = ${type2type("pw1%array(i)", kind, kind2)}$
632 : END DO
633 : !$OMP END PARALLEL DO
634 : END IF
635 368 : ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference)) THEN
636 366 : IF (ng1 >= ng2) THEN
637 0 : !$OMP PARALLEL DO PRIVATE(i,j) DEFAULT(NONE) SHARED(ng2, pw1, pw2)
638 : DO i = 1, ng2
639 : j = pw1%pw_grid%gidx(i)
640 : pw2%array(i) = ${type2type("pw1%array(j)", kind, kind2)}$
641 : END DO
642 : !$OMP END PARALLEL DO
643 : ELSE
644 366 : CALL pw_zero(pw2)
645 366 : !$OMP PARALLEL DO PRIVATE(i,j) DEFAULT(NONE) SHARED(ng1, pw1, pw2)
646 : DO i = 1, ng1
647 : j = pw1%pw_grid%gidx(i)
648 : pw2%array(j) = ${type2type("pw1%array(i)", kind, kind2)}$
649 : END DO
650 : !$OMP END PARALLEL DO
651 : END IF
652 : ELSE
653 : #:if kind=='c1d' and kind2=='c1d' and space=="gs"
654 2 : CALL pw_copy_match(pw1, pw2)
655 : #:else
656 0 : CPABORT("Copy not implemented!")
657 : #:endif
658 : END IF
659 :
660 : END IF
661 :
662 : ELSE
663 1759843 : !$OMP PARALLEL WORKSHARE PRIVATE(i) DEFAULT(NONE) SHARED(pw1, pw2)
664 : pw2%array = ${type2type("pw1%array", kind, kind2)}$
665 : !$OMP END PARALLEL WORKSHARE
666 : END IF
667 : #:else
668 2461604 : IF (ANY(SHAPE(pw2%array) /= SHAPE(pw1%array))) &
669 0 : CPABORT("3D grids must be compatible!")
670 615401 : IF (pw1%pw_grid%spherical) &
671 0 : CPABORT("3D grids must not be spherical!")
672 615401 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw1,pw2)
673 : pw2%array(:, :, :) = ${type2type("pw1%array(:, :, :)", kind, kind2)}$
674 : !$OMP END PARALLEL WORKSHARE
675 : #:endif
676 :
677 3032124 : CALL timestop(handle)
678 :
679 3032124 : END SUBROUTINE pw_copy_${kind}$_${kind2}$_${space}$
680 :
681 : ! **************************************************************************************************
682 : !> \brief ...
683 : !> \param pw ...
684 : !> \param array ...
685 : ! **************************************************************************************************
686 1597931 : SUBROUTINE pw_copy_to_array_${kind}$_${kind2}$_${space}$ (pw, array)
687 : TYPE(pw_${kind}$_${space}$_type), INTENT(IN) :: pw
688 : ${type2}$, INTENT(INOUT) :: array
689 :
690 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_copy_to_array'
691 :
692 : INTEGER :: handle
693 :
694 1597931 : CALL timeset(routineN, handle)
695 :
696 : #:if kind[1]=="1"
697 0 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw, array)
698 : array(:) = ${type2type("pw%array(:)", kind, kind2)}$
699 : !$OMP END PARALLEL WORKSHARE
700 : #:else
701 1597931 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw, array)
702 : array(:, :, :) = ${type2type("pw%array(:, :, :)", kind, kind2)}$
703 : !$OMP END PARALLEL WORKSHARE
704 : #:endif
705 :
706 1597931 : CALL timestop(handle)
707 1597931 : END SUBROUTINE pw_copy_to_array_${kind}$_${kind2}$_${space}$
708 :
709 : ! **************************************************************************************************
710 : !> \brief ...
711 : !> \param pw ...
712 : !> \param array ...
713 : ! **************************************************************************************************
714 1642926 : SUBROUTINE pw_copy_from_array_${kind}$_${kind2}$_${space}$ (pw, array)
715 : TYPE(pw_${kind}$_${space}$_type), INTENT(IN) :: pw
716 : ${type2}$, INTENT(IN) :: array
717 :
718 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_copy_from_array'
719 :
720 : INTEGER :: handle
721 :
722 1642926 : CALL timeset(routineN, handle)
723 :
724 1642926 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw, array)
725 : pw%array = ${type2type("array", kind2, kind)}$
726 : !$OMP END PARALLEL WORKSHARE
727 :
728 1642926 : CALL timestop(handle)
729 1642926 : END SUBROUTINE pw_copy_from_array_${kind}$_${kind2}$_${space}$
730 :
731 : ! **************************************************************************************************
732 : !> \brief pw2 = alpha*pw1 + beta*pw2
733 : !> alpha defaults to 1, beta defaults to 1
734 : !> \param pw1 ...
735 : !> \param pw2 ...
736 : !> \param alpha ...
737 : !> \param beta ...
738 : !> \param allow_noncompatible_grids ...
739 : !> \par History
740 : !> JGH (21-Feb-2003) : added reference grid functionality
741 : !> JGH (01-Dec-2007) : rename and remove complex alpha
742 : !> \author apsi
743 : !> \note
744 : !> Currently only summing up of respective types allowed,
745 : !> in order to avoid errors
746 : ! **************************************************************************************************
747 1869142 : SUBROUTINE pw_axpy_${kind}$_${kind2}$_${space}$ (pw1, pw2, alpha, beta, allow_noncompatible_grids)
748 :
749 : TYPE(pw_${kind}$_${space}$_type), INTENT(IN) :: pw1
750 : TYPE(pw_${kind2}$_${space}$_type), INTENT(INOUT) :: pw2
751 : REAL(KIND=dp), INTENT(in), OPTIONAL :: alpha, beta
752 : LOGICAL, INTENT(IN), OPTIONAL :: allow_noncompatible_grids
753 :
754 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_axpy'
755 :
756 : INTEGER :: handle
757 : LOGICAL :: my_allow_noncompatible_grids
758 : REAL(KIND=dp) :: my_alpha, my_beta
759 : #:if kind[1]=='1'
760 : INTEGER :: i, j, ng, ng1, ng2
761 : #:endif
762 :
763 1869142 : CALL timeset(routineN, handle)
764 :
765 1869142 : my_alpha = 1.0_dp
766 1869142 : IF (PRESENT(alpha)) my_alpha = alpha
767 :
768 1869142 : my_beta = 1.0_dp
769 1869142 : IF (PRESENT(beta)) my_beta = beta
770 :
771 1869142 : my_allow_noncompatible_grids = .FALSE.
772 1869142 : IF (PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
773 :
774 1869142 : IF (my_beta /= 1.0_dp) THEN
775 96305 : IF (my_beta == 0.0_dp) THEN
776 1524 : CALL pw_zero(pw2)
777 : ELSE
778 94781 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw2,my_beta)
779 : pw2%array = pw2%array*my_beta
780 : !$OMP END PARALLEL WORKSHARE
781 : END IF
782 : END IF
783 :
784 : #:if kind[1]=='1'
785 1367212 : IF (ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) THEN
786 :
787 643584 : IF (my_alpha == 1.0_dp) THEN
788 625484 : !$OMP PARALLEL WORKSHARE PRIVATE(i) DEFAULT(NONE) SHARED(pw1, pw2)
789 : pw2%array = pw2%array + ${type2type("pw1%array", kind, kind2)}$
790 : !$OMP END PARALLEL WORKSHARE
791 18100 : ELSE IF (my_alpha /= 0.0_dp) THEN
792 18100 : !$OMP PARALLEL WORKSHARE PRIVATE(i) DEFAULT(NONE) SHARED(pw2,pw1,my_alpha)
793 : pw2%array = pw2%array + my_alpha*${type2type("pw1%array", kind, kind2)}$
794 : !$OMP END PARALLEL WORKSHARE
795 : END IF
796 :
797 723628 : ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids) THEN
798 :
799 723628 : ng1 = SIZE(pw1%array)
800 723628 : ng2 = SIZE(pw2%array)
801 723628 : ng = MIN(ng1, ng2)
802 :
803 723628 : IF (pw1%pw_grid%spherical) THEN
804 0 : IF (my_alpha == 1.0_dp) THEN
805 0 : !$OMP PARALLEL DO PRIVATE(i) DEFAULT(NONE) SHARED(ng, pw1, pw2)
806 : DO i = 1, ng
807 : pw2%array(i) = pw2%array(i) + ${type2type("pw1%array(i)", kind, kind2)}$
808 : END DO
809 : !$OMP END PARALLEL DO
810 0 : ELSE IF (my_alpha /= 0.0_dp) THEN
811 0 : !$OMP PARALLEL DO PRIVATE(i) DEFAULT(NONE) SHARED(pw2,pw1,my_alpha,ng)
812 : DO i = 1, ng
813 : pw2%array(i) = pw2%array(i) + my_alpha*${type2type("pw1%array(i)", kind, kind2)}$
814 : END DO
815 : !$OMP END PARALLEL DO
816 : END IF
817 723628 : ELSE IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference)) THEN
818 1106 : IF (ng1 >= ng2) THEN
819 1106 : IF (my_alpha == 1.0_dp) THEN
820 1106 : !$OMP PARALLEL DO PRIVATE(i,j) DEFAULT(NONE) SHARED(ng, pw1, pw2)
821 : DO i = 1, ng
822 : j = pw2%pw_grid%gidx(i)
823 : pw2%array(i) = pw2%array(i) + ${type2type("pw1%array(j)", kind, kind2)}$
824 : END DO
825 : !$OMP END PARALLEL DO
826 0 : ELSE IF (my_alpha /= 0.0_dp) THEN
827 0 : !$OMP PARALLEL DO PRIVATE(i,j) DEFAULT(NONE) SHARED(my_alpha, ng, pw1, pw2)
828 : DO i = 1, ng
829 : j = pw2%pw_grid%gidx(i)
830 : pw2%array(i) = pw2%array(i) + my_alpha*${type2type("pw1%array(j)", kind, kind2)}$
831 : END DO
832 : !$OMP END PARALLEL DO
833 : END IF
834 : ELSE
835 0 : IF (my_alpha == 1.0_dp) THEN
836 0 : !$OMP PARALLEL DO PRIVATE(i,j) DEFAULT(NONE) SHARED(ng, pw1, pw2)
837 : DO i = 1, ng
838 : j = pw2%pw_grid%gidx(i)
839 : pw2%array(j) = pw2%array(j) + ${type2type("pw1%array(i)", kind, kind2)}$
840 : END DO
841 : !$OMP END PARALLEL DO
842 0 : ELSE IF (my_alpha /= 0.0_dp) THEN
843 0 : !$OMP PARALLEL DO PRIVATE(i,j) DEFAULT(NONE) SHARED(my_alpha, ng, pw1, pw2)
844 : DO i = 1, ng
845 : j = pw2%pw_grid%gidx(i)
846 : pw2%array(j) = pw2%array(j) + my_alpha*${type2type("pw1%array(i)", kind, kind2)}$
847 : END DO
848 : !$OMP END PARALLEL DO
849 : END IF
850 : END IF
851 722522 : ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference)) THEN
852 722522 : IF (ng1 >= ng2) THEN
853 0 : IF (my_alpha == 1.0_dp) THEN
854 0 : !$OMP PARALLEL DO PRIVATE(i,j) DEFAULT(NONE) SHARED(ng, pw1, pw2)
855 : DO i = 1, ng
856 : j = pw1%pw_grid%gidx(i)
857 : pw2%array(i) = pw2%array(i) + ${type2type("pw1%array(j)", kind, kind2)}$
858 : END DO
859 : !$OMP END PARALLEL DO
860 0 : ELSE IF (my_alpha /= 0.0_dp) THEN
861 0 : !$OMP PARALLEL DO PRIVATE(i,j) DEFAULT(NONE) SHARED(my_alpha, ng, pw1, pw2)
862 : DO i = 1, ng
863 : j = pw1%pw_grid%gidx(i)
864 : pw2%array(i) = pw2%array(i) + my_alpha*${type2type("pw1%array(j)", kind, kind2)}$
865 : END DO
866 : !$OMP END PARALLEL DO
867 : END IF
868 : ELSE
869 722522 : IF (my_alpha == 1.0_dp) THEN
870 722522 : !$OMP PARALLEL DO PRIVATE(i,j) DEFAULT(NONE) SHARED(ng, pw1, pw2)
871 : DO i = 1, ng
872 : j = pw1%pw_grid%gidx(i)
873 : pw2%array(j) = pw2%array(j) + ${type2type("pw1%array(i)", kind, kind2)}$
874 : END DO
875 : !$OMP END PARALLEL DO
876 0 : ELSE IF (my_alpha /= 0.0_dp) THEN
877 0 : !$OMP PARALLEL DO PRIVATE(i,j) DEFAULT(NONE) SHARED(my_alpha, ng, pw1, pw2)
878 : DO i = 1, ng
879 : j = pw1%pw_grid%gidx(i)
880 : pw2%array(j) = pw2%array(j) + my_alpha*${type2type("pw1%array(i)", kind, kind2)}$
881 : END DO
882 : !$OMP END PARALLEL DO
883 : END IF
884 : END IF
885 : ELSE
886 0 : CPABORT("Grids not compatible")
887 : END IF
888 : #:else
889 501930 : IF (ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) THEN
890 501400 : IF (my_alpha == 1.0_dp) THEN
891 290438 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw1, pw2)
892 : pw2%array = pw2%array + ${type2type("pw1%array", kind, kind2)}$
893 : !$OMP END PARALLEL WORKSHARE
894 210962 : ELSE IF (my_alpha /= 0.0_dp) THEN
895 209534 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw1, pw2, my_alpha)
896 : pw2%array = pw2%array + my_alpha*${type2type("pw1%array", kind, kind2)}$
897 : !$OMP END PARALLEL WORKSHARE
898 : END IF
899 :
900 530 : ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids) THEN
901 :
902 2120 : IF (ANY(SHAPE(pw1%array) /= SHAPE(pw2%array))) &
903 0 : CPABORT("Noncommensurate grids not implemented for 3D grids!")
904 :
905 530 : IF (my_alpha == 1.0_dp) THEN
906 444 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw1,pw2)
907 : pw2%array = pw2%array + ${type2type("pw1%array", kind, kind2)}$
908 : !$OMP END PARALLEL WORKSHARE
909 : ELSE
910 86 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw1,pw2,my_alpha)
911 : pw2%array = pw2%array + my_alpha*${type2type("pw1%array", kind, kind2)}$
912 : !$OMP END PARALLEL WORKSHARE
913 : END IF
914 : #:endif
915 :
916 : ELSE
917 :
918 0 : CPABORT("Grids not compatible")
919 :
920 : END IF
921 :
922 1869142 : CALL timestop(handle)
923 :
924 1869142 : END SUBROUTINE pw_axpy_${kind}$_${kind2}$_${space}$
925 :
926 : ! **************************************************************************************************
927 : !> \brief pw_out = pw_out + alpha * pw1 * pw2
928 : !> alpha defaults to 1
929 : !> \param pw_out ...
930 : !> \param pw1 ...
931 : !> \param pw2 ...
932 : !> \param alpha ...
933 : !> \author JGH
934 : ! **************************************************************************************************
935 3849 : SUBROUTINE pw_multiply_${kind}$_${kind2}$_${space}$ (pw_out, pw1, pw2, alpha)
936 :
937 : TYPE(pw_${kind}$_${space}$_type), INTENT(INOUT) :: pw_out
938 : TYPE(pw_${kind}$_${space}$_type), INTENT(IN) :: pw1
939 : TYPE(pw_${kind2}$_${space}$_type), INTENT(IN) :: pw2
940 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: alpha
941 :
942 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_multiply'
943 :
944 : INTEGER :: handle
945 : REAL(KIND=dp) :: my_alpha
946 :
947 3849 : CALL timeset(routineN, handle)
948 :
949 3849 : my_alpha = 1.0_dp
950 3849 : IF (PRESENT(alpha)) my_alpha = alpha
951 :
952 7698 : IF (.NOT. ASSOCIATED(pw_out%pw_grid, pw1%pw_grid) .OR. .NOT. ASSOCIATED(pw_out%pw_grid, pw2%pw_grid)) &
953 0 : CPABORT("pw_multiply not implemented for non-identical grids!")
954 :
955 : #if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
956 3849 : IF (my_alpha == 1.0_dp) THEN
957 3785 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw_out, pw1, pw2)
958 : pw_out%array = pw_out%array + pw1%array*${type2type("pw2%array", kind2, kind)}$
959 : !$OMP END PARALLEL WORKSHARE
960 : ELSE
961 64 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(my_alpha, pw_out, pw1, pw2)
962 : pw_out%array = pw_out%array + my_alpha*pw1%array*${type2type("pw2%array", kind2, kind)}$
963 : !$OMP END PARALLEL WORKSHARE
964 : END IF
965 : #else
966 : IF (my_alpha == 1.0_dp) THEN
967 : pw_out%array = pw_out%array + pw1%array*${type2type("pw2%array", kind2, kind)}$
968 : ELSE
969 : pw_out%array = pw_out%array + my_alpha*pw1%array*${type2type("pw2%array", kind2, kind)}$
970 : END IF
971 : #endif
972 :
973 3849 : CALL timestop(handle)
974 :
975 3849 : END SUBROUTINE pw_multiply_${kind}$_${kind2}$_${space}$
976 :
977 : ! **************************************************************************************************
978 : !> \brief ...
979 : !> \param pw1 ...
980 : !> \param pw2 ...
981 : ! **************************************************************************************************
982 413206 : SUBROUTINE pw_multiply_with_${kind}$_${kind2}$_${space}$ (pw1, pw2)
983 : TYPE(pw_${kind}$_${space}$_type), INTENT(INOUT) :: pw1
984 : TYPE(pw_${kind2}$_${space}$_type), INTENT(IN) :: pw2
985 :
986 : CHARACTER(LEN=*), PARAMETER :: routineN = 'pw_multiply_with'
987 :
988 : INTEGER :: handle
989 :
990 413206 : CALL timeset(routineN, handle)
991 :
992 413206 : IF (.NOT. ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
993 0 : CPABORT("Incompatible grids!")
994 :
995 413206 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw1,pw2)
996 : pw1%array = pw1%array*${type2type("pw2%array", kind2, kind)}$
997 : !$OMP END PARALLEL WORKSHARE
998 :
999 413206 : CALL timestop(handle)
1000 :
1001 413206 : END SUBROUTINE pw_multiply_with_${kind}$_${kind2}$_${space}$
1002 :
1003 : ! **************************************************************************************************
1004 : !> \brief Calculate integral over unit cell for functions in plane wave basis
1005 : !> only returns the real part of it ......
1006 : !> \param pw1 ...
1007 : !> \param pw2 ...
1008 : !> \param sumtype ...
1009 : !> \param just_sum ...
1010 : !> \param local_only ...
1011 : !> \return ...
1012 : !> \par History
1013 : !> JGH (14-Mar-2001) : Parallel sum and some tests, HALFSPACE case
1014 : !> \author apsi
1015 : ! **************************************************************************************************
1016 702803 : FUNCTION pw_integral_ab_${kind}$_${kind2}$_${space}$ (pw1, pw2, sumtype, just_sum, local_only) RESULT(integral_value)
1017 :
1018 : TYPE(pw_${kind}$_${space}$_type), INTENT(IN) :: pw1
1019 : TYPE(pw_${kind2}$_${space}$_type), INTENT(IN) :: pw2
1020 : INTEGER, INTENT(IN), OPTIONAL :: sumtype
1021 : LOGICAL, INTENT(IN), OPTIONAL :: just_sum, local_only
1022 : REAL(KIND=dp) :: integral_value
1023 :
1024 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_integral_ab_${kind}$_${kind2}$_${space}$'
1025 :
1026 : INTEGER :: handle, loc_sumtype
1027 : LOGICAL :: my_just_sum, my_local_only
1028 :
1029 702803 : CALL timeset(routineN, handle)
1030 :
1031 702803 : loc_sumtype = do_accurate_sum
1032 702803 : IF (PRESENT(sumtype)) loc_sumtype = sumtype
1033 :
1034 702803 : my_local_only = .FALSE.
1035 702803 : IF (PRESENT(local_only)) my_local_only = local_only
1036 :
1037 702803 : IF (.NOT. ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) THEN
1038 0 : CPABORT("Grids incompatible")
1039 : END IF
1040 :
1041 702803 : my_just_sum = .FALSE.
1042 702803 : IF (PRESENT(just_sum)) my_just_sum = just_sum
1043 :
1044 : ! do standard sum
1045 702803 : IF (loc_sumtype == do_standard_sum) THEN
1046 :
1047 : ! Do standard sum
1048 :
1049 : #:if kind=="r1d" and kind2=="r1d"
1050 0 : integral_value = DOT_PRODUCT(pw1%array, pw2%array)
1051 : #:elif kind=="r3d" and kind2=="r3d"
1052 1060095956 : integral_value = SUM(pw1%array*pw2%array)
1053 : #:elif kind[0]=="r"
1054 0 : integral_value = SUM(pw1%array*REAL(pw2%array, KIND=dp)) !? complex bit
1055 : #:elif kind2[0]=="r"
1056 0 : integral_value = SUM(REAL(pw1%array, KIND=dp)*pw2%array) !? complex bit
1057 : #:else
1058 : integral_value = SUM(REAL(CONJG(pw1%array) &
1059 0 : *pw2%array, KIND=dp)) !? complex bit
1060 : #:endif
1061 :
1062 : ELSE
1063 :
1064 : ! Do accurate sum
1065 : #:if kind[0]=="r" and kind2[0]=="r"
1066 54840 : integral_value = accurate_dot_product(pw1%array, pw2%array)
1067 : #:elif kind[0]=="r"
1068 0 : integral_value = accurate_sum(pw1%array*REAL(pw2%array, KIND=dp)) !? complex bit
1069 : #:elif kind2[0]=="r"
1070 0 : integral_value = accurate_sum(REAL(pw1%array, KIND=dp)*pw2%array) !? complex bit
1071 : #:else
1072 8864704396 : integral_value = accurate_sum(REAL(CONJG(pw1%array)*pw2%array, KIND=dp))
1073 : #:endif
1074 :
1075 : END IF
1076 :
1077 702803 : IF (.NOT. my_just_sum) THEN
1078 : #:if space == "rs"
1079 306288 : integral_value = integral_value*pw1%pw_grid%dvol
1080 : #:elif space == "gs"
1081 396417 : integral_value = integral_value*pw1%pw_grid%vol
1082 : #:else
1083 : #:stop "Unknown space: "+space
1084 : #:endif
1085 : END IF
1086 :
1087 : #:if kind[1]=="1"
1088 396417 : IF (pw1%pw_grid%grid_span == HALFSPACE) THEN
1089 236846 : integral_value = 2.0_dp*integral_value
1090 236846 : IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
1091 : #:if kind[0]=="c"
1092 128832 : REAL(CONJG(pw1%array(1))*pw2%array(1), KIND=dp)
1093 : #:elif kind2[0]=="r"
1094 0 : pw1%array(1)*pw2%array(1)
1095 : #:else
1096 0 : pw1%array(1)*REAL(pw2%array(1), KIND=dp)
1097 : #:endif
1098 : END IF
1099 : #:endif
1100 :
1101 702803 : IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == PW_MODE_DISTRIBUTED) &
1102 383400 : CALL pw1%pw_grid%para%group%sum(integral_value)
1103 :
1104 702803 : CALL timestop(handle)
1105 :
1106 702803 : END FUNCTION pw_integral_ab_${kind}$_${kind2}$_${space}$
1107 : #:endfor
1108 :
1109 : #:if kind[1]=="1"
1110 : ! **************************************************************************************************
1111 : !> \brief ...
1112 : !> \param pw1 ...
1113 : !> \param pw2 ...
1114 : !> \return ...
1115 : ! **************************************************************************************************
1116 82980 : FUNCTION pw_integral_a2b_${kind}$_${kind2}$ (pw1, pw2) RESULT(integral_value)
1117 :
1118 : TYPE(pw_${kind}$_gs_type), INTENT(IN) :: pw1
1119 : TYPE(pw_${kind2}$_gs_type), INTENT(IN) :: pw2
1120 : REAL(KIND=dp) :: integral_value
1121 :
1122 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_integral_a2b'
1123 :
1124 : INTEGER :: handle
1125 :
1126 82980 : CALL timeset(routineN, handle)
1127 :
1128 82980 : IF (.NOT. ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) THEN
1129 0 : CPABORT("Grids incompatible")
1130 : END IF
1131 :
1132 : #:if kind[0]=="c"
1133 : #:if kind2[0]=="c"
1134 241320192 : integral_value = accurate_sum(REAL(CONJG(pw1%array)*pw2%array, KIND=dp)*pw1%pw_grid%gsq)
1135 : #:else
1136 0 : integral_value = accurate_sum(REAL(CONJG(pw1%array), KIND=dp)*pw2%array*pw1%pw_grid%gsq)
1137 : #:endif
1138 : #:elif kind2[0]=="c"
1139 0 : integral_value = accurate_sum(pw1%array*REAL(pw2%array, KIND=dp)*pw1%pw_grid%gsq)
1140 : #:else
1141 0 : integral_value = accurate_sum(pw1%array*pw2%array*pw1%pw_grid%gsq)
1142 : #:endif
1143 82980 : IF (pw1%pw_grid%grid_span == HALFSPACE) THEN
1144 82980 : integral_value = 2.0_dp*integral_value
1145 : END IF
1146 :
1147 82980 : integral_value = integral_value*pw1%pw_grid%vol
1148 :
1149 82980 : IF (pw1%pw_grid%para%mode == PW_MODE_DISTRIBUTED) &
1150 78252 : CALL pw1%pw_grid%para%group%sum(integral_value)
1151 82980 : CALL timestop(handle)
1152 :
1153 82980 : END FUNCTION pw_integral_a2b_${kind}$_${kind2}$
1154 : #:endif
1155 : #:endfor
1156 :
1157 : #:for kind, type, kind2, type2 in pw_list2
1158 : #:for space in pw_spaces
1159 : #:for space2 in pw_spaces
1160 :
1161 : #:if space != space2 and ((space=="rs" and kind[1]=="3" and kind2[0]=="c") or (space=="gs" and kind2[1]=="3" and kind[0]=="c"))
1162 : ! **************************************************************************************************
1163 : !> \brief Generic function for 3d FFT of a coefficient_type or pw_r3d_rs_type
1164 : !> \param pw1 ...
1165 : !> \param pw2 ...
1166 : !> \param debug ...
1167 : !> \par History
1168 : !> JGH (30-12-2000): New setup of functions and adaptation to parallelism
1169 : !> JGH (04-01-2001): Moved routine from pws to this module, only covers
1170 : !> pw_types, no more coefficient types
1171 : !> \author apsi
1172 : !> \note
1173 : !> fft_wrap_pw1pw2
1174 : ! **************************************************************************************************
1175 3241929 : SUBROUTINE fft_wrap_pw1pw2_${kind}$_${kind2}$_${space}$_${space2}$ (pw1, pw2, debug)
1176 :
1177 : TYPE(pw_${kind}$_${space}$_type), INTENT(IN) :: pw1
1178 : TYPE(pw_${kind2}$_${space2}$_type), INTENT(INOUT) :: pw2
1179 : LOGICAL, INTENT(IN), OPTIONAL :: debug
1180 :
1181 : CHARACTER(len=*), PARAMETER :: routineN = 'fft_wrap_pw1pw2'
1182 :
1183 3241929 : COMPLEX(KIND=dp), DIMENSION(:, :), CONTIGUOUS, POINTER :: grays
1184 3241929 : COMPLEX(KIND=dp), DIMENSION(:, :, :), CONTIGUOUS, POINTER :: c_in, c_out
1185 : INTEGER :: handle, handle2, my_pos, nrays, &
1186 : out_unit
1187 : #:if (space=="rs" and kind=="r3d" and kind2=="c1d") or (space=="gs" and kind=="c1d" and kind2=="r3d")
1188 : INTEGER, DIMENSION(3) :: nloc
1189 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
1190 : LOGICAL :: use_pw_gpu
1191 : #endif
1192 : #:endif
1193 3241929 : INTEGER, DIMENSION(:), POINTER :: n
1194 : LOGICAL :: test
1195 :
1196 3241929 : CALL timeset(routineN, handle2)
1197 3241929 : out_unit = cp_logger_get_default_io_unit()
1198 : CALL timeset(routineN//"_"//TRIM(ADJUSTL(cp_to_string( &
1199 3241929 : CEILING(pw1%pw_grid%cutoff/10)*10))), handle)
1200 :
1201 3241929 : NULLIFY (c_in)
1202 3241929 : NULLIFY (c_out)
1203 :
1204 3241929 : IF (PRESENT(debug)) THEN
1205 1072 : test = debug
1206 : ELSE
1207 3240857 : test = .FALSE.
1208 : END IF
1209 :
1210 : !..check if grids are compatible
1211 3241929 : IF (.NOT. ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) THEN
1212 0 : IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol) THEN
1213 0 : CPABORT("PW grids not compatible")
1214 : END IF
1215 0 : IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group) THEN
1216 0 : CPABORT("PW grids have not compatible MPI groups")
1217 : END IF
1218 : END IF
1219 :
1220 3241929 : n => pw1%pw_grid%npts
1221 :
1222 3241929 : IF (pw1%pw_grid%para%mode == PW_MODE_LOCAL) THEN
1223 :
1224 : !
1225 : !..replicated data, use local FFT
1226 : !
1227 :
1228 598575 : IF (test .AND. out_unit > 0) THEN
1229 0 : WRITE (out_unit, '(A)') " FFT Protocol "
1230 : #:if space=="rs"
1231 0 : WRITE (out_unit, '(A,T76,A)') " Transform direction ", "FWFFT"
1232 0 : WRITE (out_unit, '(A,T72,A)') " in space ", "REALSPACE"
1233 0 : WRITE (out_unit, '(A,T72,A)') " out space ", "REALSPACE"
1234 : #:else
1235 0 : WRITE (out_unit, '(A,T76,A)') " Transform direction ", "BWFFT"
1236 0 : WRITE (out_unit, '(A,T66,A)') " in space ", "RECIPROCALSPACE"
1237 0 : WRITE (out_unit, '(A,T66,A)') " out space ", "RECIPROCALSPACE"
1238 : #:endif
1239 : END IF
1240 :
1241 : #:if space=="rs"
1242 : #:if kind==kind2=="c3d"
1243 0 : c_in => pw1%array
1244 0 : c_out => pw2%array
1245 0 : CALL fft3d(FWFFT, n, c_in, c_out, debug=test)
1246 : #:elif kind=="r3d" and kind2=="c3d"
1247 0 : pw2%array = CMPLX(pw1%array, 0.0_dp, KIND=dp)
1248 0 : c_out => pw2%array
1249 0 : CALL fft3d(FWFFT, n, c_out, debug=test)
1250 : #:elif kind=="c3d" and kind2=="c1d"
1251 0 : c_in => pw1%array
1252 0 : ALLOCATE (c_out(n(1), n(2), n(3)))
1253 : ! transform
1254 0 : CALL fft3d(FWFFT, n, c_in, c_out, debug=test)
1255 : ! gather results
1256 0 : IF (test .AND. out_unit > 0) WRITE (out_unit, '(A)') " PW_GATHER : 3d -> 1d "
1257 0 : CALL pw_gather_s_${kind2}$_c3d(pw2, c_out)
1258 0 : DEALLOCATE (c_out)
1259 : #:elif kind=="r3d" and kind2=="c1d"
1260 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
1261 : CALL pw_gpu_r3dc1d_3d(pw1, pw2)
1262 : #elif defined (__PW_FPGA)
1263 : ALLOCATE (c_out(n(1), n(2), n(3)))
1264 : ! check if bitstream for the fft size is present
1265 : ! if not, perform fft3d in CPU
1266 : IF (pw_fpga_init_bitstream(n) == 1) THEN
1267 : CALL pw_copy_to_array(pw1, c_out)
1268 : #if (__PW_FPGA_SP && __PW_FPGA)
1269 : CALL pw_fpga_r3dc1d_3d_sp(n, c_out)
1270 : #else
1271 : CALL pw_fpga_r3dc1d_3d_dp(n, c_out)
1272 : #endif
1273 : CALL zdscal(n(1)*n(2)*n(3), 1.0_dp/pw1%pw_grid%ngpts, c_out, 1)
1274 : CALL pw_gather_s_${kind2}$_c3d(pw2, c_out)
1275 : ELSE
1276 : CALL pw_copy_to_array(pw1, c_out)
1277 : CALL fft3d(FWFFT, n, c_out, debug=test)
1278 : CALL pw_gather_s_${kind2}$_c3d(pw2, c_out)
1279 : END IF
1280 : DEALLOCATE (c_out)
1281 : #else
1282 1457105 : ALLOCATE (c_out(n(1), n(2), n(3)))
1283 2747198148 : c_out = 0.0_dp
1284 291421 : CALL pw_copy_to_array(pw1, c_out)
1285 291421 : CALL fft3d(FWFFT, n, c_out, debug=test)
1286 291421 : CALL pw_gather_s_${kind2}$_c3d(pw2, c_out)
1287 291421 : DEALLOCATE (c_out)
1288 : #endif
1289 : #:endif
1290 : #:else
1291 : #:if kind=="c3d" and kind2=="c3d"
1292 0 : c_in => pw1%array
1293 0 : c_out => pw2%array
1294 0 : CALL fft3d(BWFFT, n, c_in, c_out, debug=test)
1295 : #:elif kind=="c3d" and kind2=="r3d"
1296 0 : c_in => pw1%array
1297 0 : ALLOCATE (c_out(n(1), n(2), n(3)))
1298 0 : CALL fft3d(BWFFT, n, c_in, c_out, debug=test)
1299 : ! use real part only
1300 0 : IF (test .AND. out_unit > 0) WRITE (out_unit, '(A)') " REAL part "
1301 0 : pw2%array = REAL(c_out, KIND=dp)
1302 0 : DEALLOCATE (c_out)
1303 : #:elif kind=="c1d" and kind2=="c3d"
1304 0 : c_out => pw2%array
1305 0 : IF (test .AND. out_unit > 0) WRITE (out_unit, '(A)') " PW_SCATTER : 3d -> 1d "
1306 0 : CALL pw_scatter_s_${kind}$_c3d(pw1, c_out)
1307 0 : CALL fft3d(BWFFT, n, c_out, debug=test)
1308 : #:elif kind=="c1d" and kind2=="r3d"
1309 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
1310 : CALL pw_gpu_c1dr3d_3d(pw1, pw2)
1311 : #elif defined (__PW_FPGA)
1312 : ALLOCATE (c_out(n(1), n(2), n(3)))
1313 : ! check if bitstream for the fft size is present
1314 : ! if not, perform fft3d in CPU
1315 : IF (pw_fpga_init_bitstream(n) == 1) THEN
1316 : CALL pw_scatter_s_${kind}$_c3d(pw1, c_out)
1317 : ! transform using FPGA
1318 : #if (__PW_FPGA_SP && __PW_FPGA)
1319 : CALL pw_fpga_c1dr3d_3d_sp(n, c_out)
1320 : #else
1321 : CALL pw_fpga_c1dr3d_3d_dp(n, c_out)
1322 : #endif
1323 : CALL zdscal(n(1)*n(2)*n(3), 1.0_dp, c_out, 1)
1324 : ! use real part only
1325 : CALL pw_copy_from_array(pw2, c_out)
1326 : ELSE
1327 : IF (test .AND. out_unit > 0) WRITE (out_unit, '(A)') " PW_SCATTER : 3d -> 1d "
1328 : CALL pw_scatter_s_${kind}$_c3d(pw1, c_out)
1329 : ! transform
1330 : CALL fft3d(BWFFT, n, c_out, debug=test)
1331 : ! use real part only
1332 : IF (test .AND. out_unit > 0) WRITE (out_unit, '(A)') " REAL part "
1333 : CALL pw_copy_from_array(pw2, c_out)
1334 : END IF
1335 : DEALLOCATE (c_out)
1336 : #else
1337 1535770 : ALLOCATE (c_out(n(1), n(2), n(3)))
1338 307154 : IF (test .AND. out_unit > 0) WRITE (out_unit, '(A)') " PW_SCATTER : 3d -> 1d "
1339 307154 : CALL pw_scatter_s_${kind}$_c3d(pw1, c_out)
1340 : ! transform
1341 307154 : CALL fft3d(BWFFT, n, c_out, debug=test)
1342 : ! use real part only
1343 307154 : IF (test .AND. out_unit > 0) WRITE (out_unit, '(A)') " REAL part "
1344 307154 : CALL pw_copy_from_array(pw2, c_out)
1345 307154 : DEALLOCATE (c_out)
1346 : #endif
1347 : #:endif
1348 : #:endif
1349 :
1350 598575 : IF (test .AND. out_unit > 0) WRITE (out_unit, '(A)') " End of FFT Protocol "
1351 :
1352 : ELSE
1353 :
1354 : !
1355 : !..parallel FFT
1356 : !
1357 :
1358 2643354 : IF (test .AND. out_unit > 0) THEN
1359 8 : WRITE (out_unit, '(A)') " FFT Protocol "
1360 : #:if space=="rs"
1361 4 : WRITE (out_unit, '(A,T76,A)') " Transform direction ", "FWFFT"
1362 4 : WRITE (out_unit, '(A,T72,A)') " in space ", "REALSPACE"
1363 4 : WRITE (out_unit, '(A,T72,A)') " out space ", "REALSPACE"
1364 : #:else
1365 4 : WRITE (out_unit, '(A,T76,A)') " Transform direction ", "BWFFT"
1366 4 : WRITE (out_unit, '(A,T66,A)') " in space ", "RECIPROCALSPACE"
1367 4 : WRITE (out_unit, '(A,T66,A)') " out space ", "RECIPROCALSPACE"
1368 : #:endif
1369 : END IF
1370 :
1371 2643354 : my_pos = pw1%pw_grid%para%group%mepos
1372 2643354 : nrays = pw1%pw_grid%para%nyzray(my_pos)
1373 2643354 : grays => pw1%pw_grid%grays
1374 :
1375 : #:if space=="rs"
1376 : #:if kind=="c3d" and kind2=="c1d"
1377 : !..prepare input
1378 536 : c_in => pw1%array
1379 1918288 : grays = z_zero
1380 : !..transform
1381 536 : IF (pw1%pw_grid%para%ray_distribution) THEN
1382 : CALL fft3d(FWFFT, n, c_in, grays, pw1%pw_grid%para%group, &
1383 : pw1%pw_grid%para%yzp, pw1%pw_grid%para%nyzray, &
1384 432 : pw1%pw_grid%para%bo, debug=test)
1385 : ELSE
1386 : CALL fft3d(FWFFT, n, c_in, grays, pw1%pw_grid%para%group, &
1387 104 : pw1%pw_grid%para%bo, debug=test)
1388 : END IF
1389 : !..prepare output
1390 536 : IF (test .AND. out_unit > 0) &
1391 4 : WRITE (out_unit, '(A)') " PW_GATHER : 2d -> 1d "
1392 536 : CALL pw_gather_p_${kind2}$ (pw2, grays)
1393 : #:elif kind=="r3d" and kind2=="c1d"
1394 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
1395 : ! (no ray dist. is not efficient in CUDA)
1396 : use_pw_gpu = pw1%pw_grid%para%ray_distribution
1397 : IF (use_pw_gpu) THEN
1398 : CALL pw_gpu_r3dc1d_3d_ps(pw1, pw2)
1399 : ELSE
1400 : #endif
1401 : !.. prepare input
1402 5226040 : nloc = pw1%pw_grid%npts_local
1403 6532550 : ALLOCATE (c_in(nloc(1), nloc(2), nloc(3)))
1404 1306510 : CALL pw_copy_to_array(pw1, c_in)
1405 37395099527 : grays = z_zero
1406 : !..transform
1407 1306510 : IF (pw1%pw_grid%para%ray_distribution) THEN
1408 : CALL fft3d(FWFFT, n, c_in, grays, pw1%pw_grid%para%group, &
1409 : pw1%pw_grid%para%yzp, pw1%pw_grid%para%nyzray, &
1410 1306510 : pw1%pw_grid%para%bo, debug=test)
1411 : ELSE
1412 : CALL fft3d(FWFFT, n, c_in, grays, pw1%pw_grid%para%group, &
1413 0 : pw1%pw_grid%para%bo, debug=test)
1414 : END IF
1415 : !..prepare output
1416 1306510 : IF (test .AND. out_unit > 0) &
1417 0 : WRITE (out_unit, '(A)') " PW_GATHER : 2d -> 1d "
1418 1306510 : CALL pw_gather_p_${kind2}$ (pw2, grays)
1419 1306510 : DEALLOCATE (c_in)
1420 :
1421 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
1422 : END IF
1423 : #endif
1424 : #:endif
1425 : #:else
1426 : #:if kind=="c1d" and kind2=="c3d"
1427 : !..prepare input
1428 536 : IF (test .AND. out_unit > 0) &
1429 4 : WRITE (out_unit, '(A)') " PW_SCATTER : 2d -> 1d "
1430 1918288 : grays = z_zero
1431 536 : CALL pw_scatter_p_${kind}$ (pw1, grays)
1432 536 : c_in => pw2%array
1433 : !..transform
1434 536 : IF (pw1%pw_grid%para%ray_distribution) THEN
1435 : CALL fft3d(BWFFT, n, c_in, grays, pw1%pw_grid%para%group, &
1436 : pw1%pw_grid%para%yzp, pw1%pw_grid%para%nyzray, &
1437 432 : pw1%pw_grid%para%bo, debug=test)
1438 : ELSE
1439 : CALL fft3d(BWFFT, n, c_in, grays, pw1%pw_grid%para%group, &
1440 104 : pw1%pw_grid%para%bo, debug=test)
1441 : END IF
1442 : !..prepare output (nothing to do)
1443 : #:elif kind=="c1d" and kind2=="r3d"
1444 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
1445 : ! (no ray dist. is not efficient in CUDA)
1446 : use_pw_gpu = pw1%pw_grid%para%ray_distribution
1447 : IF (use_pw_gpu) THEN
1448 : CALL pw_gpu_c1dr3d_3d_ps(pw1, pw2)
1449 : ELSE
1450 : #endif
1451 : !.. prepare input
1452 1335772 : IF (test .AND. out_unit > 0) &
1453 0 : WRITE (out_unit, '(A)') " PW_SCATTER : 2d -> 1d "
1454 41331478650 : grays = z_zero
1455 1335772 : CALL pw_scatter_p_${kind}$ (pw1, grays)
1456 5343088 : nloc = pw2%pw_grid%npts_local
1457 6678860 : ALLOCATE (c_in(nloc(1), nloc(2), nloc(3)))
1458 : !..transform
1459 1335772 : IF (pw1%pw_grid%para%ray_distribution) THEN
1460 : CALL fft3d(BWFFT, n, c_in, grays, pw1%pw_grid%para%group, &
1461 : pw1%pw_grid%para%yzp, pw1%pw_grid%para%nyzray, &
1462 1335772 : pw1%pw_grid%para%bo, debug=test)
1463 : ELSE
1464 : CALL fft3d(BWFFT, n, c_in, grays, pw1%pw_grid%para%group, &
1465 0 : pw1%pw_grid%para%bo, debug=test)
1466 : END IF
1467 : !..prepare output
1468 1335772 : IF (test .AND. out_unit > 0) &
1469 0 : WRITE (out_unit, '(A)') " Real part "
1470 1335772 : CALL pw_copy_from_array(pw2, c_in)
1471 1335772 : DEALLOCATE (c_in)
1472 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
1473 : END IF
1474 : #endif
1475 : #:endif
1476 : #:endif
1477 : END IF
1478 :
1479 3241929 : IF (test .AND. out_unit > 0) THEN
1480 8 : WRITE (out_unit, '(A)') " End of FFT Protocol "
1481 : END IF
1482 :
1483 3241929 : CALL timestop(handle)
1484 3241929 : CALL timestop(handle2)
1485 :
1486 3241929 : END SUBROUTINE fft_wrap_pw1pw2_${kind}$_${kind2}$_${space}$_${space2}$
1487 : #:endif
1488 : #:endfor
1489 : #:endfor
1490 :
1491 : #:if kind[1]=='1' and kind2[1]=='3'
1492 :
1493 : ! **************************************************************************************************
1494 : !> \brief Gathers the pw vector from a 3d data field
1495 : !> \param pw ...
1496 : !> \param c ...
1497 : !> \param scale ...
1498 : !> \par History
1499 : !> none
1500 : !> \author JGH
1501 : ! **************************************************************************************************
1502 0 : SUBROUTINE pw_gather_s_${kind}$_${kind2}$_2(pw1, pw2, scale)
1503 :
1504 : TYPE(pw_${kind2}$_gs_type), INTENT(IN) :: pw1
1505 : TYPE(pw_${kind}$_gs_type), INTENT(INOUT) :: pw2
1506 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: scale
1507 :
1508 0 : CALL pw_gather_s_${kind}$_${kind2}$ (pw2, pw1%array, scale)
1509 :
1510 0 : END SUBROUTINE pw_gather_s_${kind}$_${kind2}$_2
1511 :
1512 : ! **************************************************************************************************
1513 : !> \brief Gathers the pw vector from a 3d data field
1514 : !> \param pw ...
1515 : !> \param c ...
1516 : !> \param scale ...
1517 : !> \par History
1518 : !> none
1519 : !> \author JGH
1520 : ! **************************************************************************************************
1521 291421 : SUBROUTINE pw_gather_s_${kind}$_${kind2}$ (pw, c, scale)
1522 :
1523 : TYPE(pw_${kind}$_gs_type), INTENT(INOUT) :: pw
1524 : ${type2}$, CONTIGUOUS, INTENT(IN) :: c
1525 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: scale
1526 :
1527 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_gather_s'
1528 :
1529 : INTEGER :: gpt, handle, l, m, n
1530 :
1531 291421 : CALL timeset(routineN, handle)
1532 :
1533 : ASSOCIATE (mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
1534 : ngpts => SIZE(pw%pw_grid%gsq), ghat => pw%pw_grid%g_hat)
1535 :
1536 291421 : IF (PRESENT(scale)) THEN
1537 0 : !$OMP PARALLEL DO PRIVATE(gpt, l, m, n) DEFAULT(NONE) SHARED(c, pw, scale)
1538 : DO gpt = 1, ngpts
1539 : l = mapl(ghat(1, gpt)) + 1
1540 : m = mapm(ghat(2, gpt)) + 1
1541 : n = mapn(ghat(3, gpt)) + 1
1542 : pw%array(gpt) = scale*${type2type("c(l, m, n)", kind2, kind)}$
1543 : END DO
1544 : !$OMP END PARALLEL DO
1545 : ELSE
1546 291421 : !$OMP PARALLEL DO PRIVATE(gpt, l, m, n) DEFAULT(NONE) SHARED(c, pw)
1547 : DO gpt = 1, ngpts
1548 : l = mapl(ghat(1, gpt)) + 1
1549 : m = mapm(ghat(2, gpt)) + 1
1550 : n = mapn(ghat(3, gpt)) + 1
1551 : pw%array(gpt) = ${type2type("c(l, m, n)", kind2, kind)}$
1552 : END DO
1553 : !$OMP END PARALLEL DO
1554 : END IF
1555 :
1556 : END ASSOCIATE
1557 :
1558 291421 : CALL timestop(handle)
1559 :
1560 291421 : END SUBROUTINE pw_gather_s_${kind}$_${kind2}$
1561 :
1562 : ! **************************************************************************************************
1563 : !> \brief Scatters a pw vector to a 3d data field
1564 : !> \param pw ...
1565 : !> \param c ...
1566 : !> \param scale ...
1567 : !> \par History
1568 : !> none
1569 : !> \author JGH
1570 : ! **************************************************************************************************
1571 0 : SUBROUTINE pw_scatter_s_${kind}$_${kind2}$_2(pw1, pw2, scale)
1572 :
1573 : TYPE(pw_${kind}$_gs_type), INTENT(IN) :: pw1
1574 : TYPE(pw_${kind2}$_gs_type), INTENT(INOUT) :: pw2
1575 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: scale
1576 :
1577 0 : CALL pw_scatter_s_${kind}$_${kind2}$ (pw1, pw2%array, scale)
1578 :
1579 0 : END SUBROUTINE pw_scatter_s_${kind}$_${kind2}$_2
1580 :
1581 : ! **************************************************************************************************
1582 : !> \brief Scatters a pw vector to a 3d data field
1583 : !> \param pw ...
1584 : !> \param c ...
1585 : !> \param scale ...
1586 : !> \par History
1587 : !> none
1588 : !> \author JGH
1589 : ! **************************************************************************************************
1590 307154 : SUBROUTINE pw_scatter_s_${kind}$_${kind2}$ (pw, c, scale)
1591 :
1592 : TYPE(pw_${kind}$_gs_type), INTENT(IN) :: pw
1593 : ${type2}$, CONTIGUOUS, INTENT(INOUT) :: c
1594 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: scale
1595 :
1596 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_scatter_s'
1597 :
1598 : INTEGER :: gpt, handle, l, m, n
1599 :
1600 307154 : CALL timeset(routineN, handle)
1601 :
1602 : ASSOCIATE (mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, &
1603 : ghat => pw%pw_grid%g_hat, ngpts => SIZE(pw%pw_grid%gsq))
1604 :
1605 : ! should only zero the unused bits (but the zero is needed)
1606 2656265070 : IF (.NOT. PRESENT(scale)) c = 0.0_dp
1607 :
1608 307154 : IF (PRESENT(scale)) THEN
1609 0 : !$OMP PARALLEL DO PRIVATE(gpt, l, m, n) DEFAULT(NONE) SHARED(c, pw, scale)
1610 : DO gpt = 1, ngpts
1611 : l = mapl(ghat(1, gpt)) + 1
1612 : m = mapm(ghat(2, gpt)) + 1
1613 : n = mapn(ghat(3, gpt)) + 1
1614 : c(l, m, n) = scale*${type2type("pw%array(gpt)", kind, kind2)}$
1615 : END DO
1616 : !$OMP END PARALLEL DO
1617 : ELSE
1618 307154 : !$OMP PARALLEL DO PRIVATE(gpt, l, m, n) DEFAULT(NONE) SHARED(c, pw)
1619 : DO gpt = 1, ngpts
1620 : l = mapl(ghat(1, gpt)) + 1
1621 : m = mapm(ghat(2, gpt)) + 1
1622 : n = mapn(ghat(3, gpt)) + 1
1623 : c(l, m, n) = ${type2type("pw%array(gpt)", kind, kind2)}$
1624 : END DO
1625 : !$OMP END PARALLEL DO
1626 : END IF
1627 :
1628 : END ASSOCIATE
1629 :
1630 307154 : IF (pw%pw_grid%grid_span == HALFSPACE) THEN
1631 :
1632 : ASSOCIATE (mapl => pw%pw_grid%mapl%neg, mapm => pw%pw_grid%mapm%neg, mapn => pw%pw_grid%mapn%neg, &
1633 : ghat => pw%pw_grid%g_hat, ngpts => SIZE(pw%pw_grid%gsq))
1634 :
1635 8922 : IF (PRESENT(scale)) THEN
1636 0 : !$OMP PARALLEL DO PRIVATE(gpt, l, m, n) DEFAULT(NONE) SHARED(c, pw, scale)
1637 : DO gpt = 1, ngpts
1638 : l = mapl(ghat(1, gpt)) + 1
1639 : m = mapm(ghat(2, gpt)) + 1
1640 : n = mapn(ghat(3, gpt)) + 1
1641 : c(l, m, n) = scale*#{if kind[0]=="c" and kind2[0]=="c"}#CONJG#{endif}#(${type2type("pw%array(gpt)", kind, kind2)}$)
1642 : END DO
1643 : !$OMP END PARALLEL DO
1644 : ELSE
1645 8922 : !$OMP PARALLEL DO PRIVATE(gpt, l, m, n) DEFAULT(NONE) SHARED(c, pw)
1646 : DO gpt = 1, ngpts
1647 : l = mapl(ghat(1, gpt)) + 1
1648 : m = mapm(ghat(2, gpt)) + 1
1649 : n = mapn(ghat(3, gpt)) + 1
1650 : c(l, m, n) = #{if kind[0]=="c" and kind2[0]=="c"}#CONJG#{endif}#(${type2type("pw%array(gpt)", kind, kind2)}$)
1651 : END DO
1652 : !$OMP END PARALLEL DO
1653 : END IF
1654 :
1655 : END ASSOCIATE
1656 :
1657 : END IF
1658 :
1659 307154 : CALL timestop(handle)
1660 :
1661 307154 : END SUBROUTINE pw_scatter_s_${kind}$_${kind2}$
1662 : #:endif
1663 : #:endfor
1664 :
1665 : ! **************************************************************************************************
1666 : !> \brief Multiply all data points with a Gaussian damping factor
1667 : !> Needed for longrange Coulomb potential
1668 : !> V(\vec r)=erf(omega*r)/r
1669 : !> V(\vec g)=\frac{4*\pi}{g**2}*exp(-g**2/omega**2)
1670 : !> \param pw ...
1671 : !> \param omega ...
1672 : !> \par History
1673 : !> Frederick Stein (12-04-2019) created
1674 : !> \author Frederick Stein (12-Apr-2019)
1675 : !> \note
1676 : !> Performs a Gaussian damping
1677 : ! **************************************************************************************************
1678 3841 : SUBROUTINE pw_gauss_damp(pw, omega)
1679 :
1680 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: pw
1681 : REAL(KIND=dp), INTENT(IN) :: omega
1682 :
1683 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_gauss_damp'
1684 :
1685 : INTEGER :: handle, i
1686 : REAL(KIND=dp) :: omega_2, tmp
1687 :
1688 3841 : CALL timeset(routineN, handle)
1689 3841 : CPASSERT(omega >= 0)
1690 :
1691 3841 : omega_2 = omega*omega
1692 3841 : omega_2 = 0.25_dp/omega_2
1693 :
1694 3841 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,tmp) SHARED(pw,omega_2)
1695 : DO i = 1, SIZE(pw%array)
1696 : tmp = EXP(-pw%pw_grid%gsq(i)*omega_2)
1697 : pw%array(i) = pw%array(i)*tmp
1698 : END DO
1699 : !$OMP END PARALLEL DO
1700 :
1701 3841 : CALL timestop(handle)
1702 :
1703 3841 : END SUBROUTINE pw_gauss_damp
1704 :
1705 : ! **************************************************************************************************
1706 : !> \brief Multiply all data points with the logarithmic derivative of a Gaussian
1707 : !> \param pw ...
1708 : !> \param omega ...
1709 : !> \note
1710 : !> Performs a Gaussian damping
1711 : ! **************************************************************************************************
1712 1329 : SUBROUTINE pw_log_deriv_gauss(pw, omega)
1713 :
1714 : TYPE(pw_c1d_gs_type), INTENT(IN) :: pw
1715 : REAL(KIND=dp), INTENT(IN) :: omega
1716 :
1717 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_log_deriv_gauss'
1718 :
1719 : INTEGER :: handle, i
1720 : REAL(KIND=dp) :: omega_2, tmp
1721 :
1722 1329 : CALL timeset(routineN, handle)
1723 1329 : CPASSERT(omega >= 0)
1724 :
1725 1329 : omega_2 = omega*omega
1726 1329 : omega_2 = 0.25_dp/omega_2
1727 :
1728 1329 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,tmp) SHARED(pw,omega_2)
1729 : DO i = 1, SIZE(pw%array)
1730 : tmp = 1.0_dp + omega_2*pw%pw_grid%gsq(i)
1731 : pw%array(i) = pw%array(i)*tmp
1732 : END DO
1733 : !$OMP END PARALLEL DO
1734 :
1735 1329 : CALL timestop(handle)
1736 1329 : END SUBROUTINE pw_log_deriv_gauss
1737 :
1738 : ! **************************************************************************************************
1739 : !> \brief Multiply all data points with a Gaussian damping factor
1740 : !> Needed for longrange Coulomb potential
1741 : !> V(\vec r)=erf(omega*r)/r
1742 : !> V(\vec g)=\frac{4*\pi}{g**2}*exp(-g**2/omega**2)
1743 : !> \param pw ...
1744 : !> \param omega ...
1745 : !> \par History
1746 : !> Frederick Stein (12-04-2019) created
1747 : !> \author Frederick Stein (12-Apr-2019)
1748 : !> \note
1749 : !> Performs a Gaussian damping
1750 : ! **************************************************************************************************
1751 0 : SUBROUTINE pw_compl_gauss_damp(pw, omega)
1752 :
1753 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: pw
1754 : REAL(KIND=dp), INTENT(IN) :: omega
1755 :
1756 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_compl_gauss_damp'
1757 :
1758 : INTEGER :: cnt, handle, i
1759 : REAL(KIND=dp) :: omega_2, tmp, tmp2
1760 :
1761 0 : CALL timeset(routineN, handle)
1762 :
1763 0 : omega_2 = omega*omega
1764 0 : omega_2 = 0.25_dp/omega_2
1765 :
1766 0 : cnt = SIZE(pw%array)
1767 :
1768 0 : !$OMP PARALLEL DO PRIVATE(i, tmp, tmp2) DEFAULT(NONE) SHARED(cnt, pw,omega_2)
1769 : DO i = 1, cnt
1770 : tmp = -omega_2*pw%pw_grid%gsq(i)
1771 : IF (ABS(tmp) > 1.0E-5_dp) THEN
1772 : tmp2 = 1.0_dp - EXP(tmp)
1773 : ELSE
1774 : tmp2 = tmp + 0.5_dp*tmp*(tmp + (1.0_dp/3.0_dp)*tmp**2)
1775 : END IF
1776 : pw%array(i) = pw%array(i)*tmp2
1777 : END DO
1778 : !$OMP END PARALLEL DO
1779 :
1780 0 : CALL timestop(handle)
1781 :
1782 0 : END SUBROUTINE pw_compl_gauss_damp
1783 :
1784 : ! **************************************************************************************************
1785 : !> \brief Multiply all data points with the logarithmic derivative of the complementary Gaussian damping factor
1786 : !> \param pw ...
1787 : !> \param omega ...
1788 : !> \note
1789 : ! **************************************************************************************************
1790 0 : SUBROUTINE pw_log_deriv_compl_gauss(pw, omega)
1791 :
1792 : TYPE(pw_c1d_gs_type), INTENT(IN) :: pw
1793 : REAL(KIND=dp), INTENT(IN) :: omega
1794 :
1795 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_log_deriv_compl_gauss'
1796 :
1797 : INTEGER :: handle, i
1798 : REAL(KIND=dp) :: omega_2, tmp, tmp2
1799 :
1800 0 : CALL timeset(routineN, handle)
1801 :
1802 0 : omega_2 = omega*omega
1803 0 : omega_2 = 0.25_dp/omega_2
1804 :
1805 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,tmp,tmp2) &
1806 0 : !$OMP SHARED(pw,omega_2)
1807 : DO i = 1, SIZE(pw%array)
1808 : tmp = omega_2*pw%pw_grid%gsq(i)
1809 : ! For too small arguments, use the Taylor polynomial to prevent division by zero
1810 : IF (ABS(tmp) >= 0.003_dp) THEN
1811 : tmp2 = 1.0_dp - tmp*EXP(-tmp)/(1.0_dp - EXP(-tmp))
1812 : ELSE
1813 : tmp2 = 0.5_dp*tmp - tmp**2/12.0_dp
1814 : END IF
1815 : pw%array(i) = pw%array(i)*tmp2
1816 : END DO
1817 : !$OMP END PARALLEL DO
1818 :
1819 0 : CALL timestop(handle)
1820 :
1821 0 : END SUBROUTINE pw_log_deriv_compl_gauss
1822 :
1823 : ! **************************************************************************************************
1824 : !> \brief Multiply all data points with a Gaussian damping factor and mixes it with the original function
1825 : !> Needed for mixed longrange/Coulomb potential
1826 : !> V(\vec r)=(a+b*erf(omega*r))/r
1827 : !> V(\vec g)=\frac{4*\pi}{g**2}*(a+b*exp(-g**2/omega**2))
1828 : !> \param pw ...
1829 : !> \param omega ...
1830 : !> \param scale_coul ...
1831 : !> \param scale_long ...
1832 : !> \par History
1833 : !> Frederick Stein (16-Dec-2021) created
1834 : !> \author Frederick Stein (16-Dec-2021)
1835 : !> \note
1836 : !> Performs a Gaussian damping
1837 : ! **************************************************************************************************
1838 332 : SUBROUTINE pw_gauss_damp_mix(pw, omega, scale_coul, scale_long)
1839 :
1840 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: pw
1841 : REAL(KIND=dp), INTENT(IN) :: omega, scale_coul, scale_long
1842 :
1843 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_gauss_damp_mix'
1844 :
1845 : INTEGER :: handle, i
1846 : REAL(KIND=dp) :: omega_2, tmp
1847 :
1848 332 : CALL timeset(routineN, handle)
1849 :
1850 332 : omega_2 = omega*omega
1851 332 : omega_2 = 0.25_dp/omega_2
1852 :
1853 332 : !$OMP PARALLEL DO DEFAULT(NONE) SHARED(pw, omega_2, scale_coul, scale_long) PRIVATE(i,tmp)
1854 : DO i = 1, SIZE(pw%array)
1855 : tmp = scale_coul + scale_long*EXP(-pw%pw_grid%gsq(i)*omega_2)
1856 : pw%array(i) = pw%array(i)*tmp
1857 : END DO
1858 : !$OMP END PARALLEL DO
1859 :
1860 332 : CALL timestop(handle)
1861 :
1862 332 : END SUBROUTINE pw_gauss_damp_mix
1863 :
1864 : ! **************************************************************************************************
1865 : !> \brief Multiply all data points with the logarithmic derivative of the mixed longrange/Coulomb potential
1866 : !> Needed for mixed longrange/Coulomb potential
1867 : !> \param pw ...
1868 : !> \param omega ...
1869 : !> \param scale_coul ...
1870 : !> \param scale_long ...
1871 : !> \note
1872 : ! **************************************************************************************************
1873 249 : SUBROUTINE pw_log_deriv_mix_cl(pw, omega, scale_coul, scale_long)
1874 :
1875 : TYPE(pw_c1d_gs_type), INTENT(IN) :: pw
1876 : REAL(KIND=dp), INTENT(IN) :: omega, scale_coul, scale_long
1877 :
1878 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_log_deriv_mix_cl'
1879 :
1880 : INTEGER :: handle, i
1881 : REAL(KIND=dp) :: omega_2, tmp, tmp2
1882 :
1883 249 : CALL timeset(routineN, handle)
1884 :
1885 249 : omega_2 = omega*omega
1886 249 : omega_2 = 0.25_dp/omega_2
1887 :
1888 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,tmp,tmp2) &
1889 249 : !$OMP SHARED(pw,omega_2,scale_long,scale_coul)
1890 : DO i = 1, SIZE(pw%array)
1891 : tmp = omega_2*pw%pw_grid%gsq(i)
1892 : tmp2 = 1.0_dp + scale_long*tmp*EXP(-tmp)/(scale_coul + scale_long*EXP(-tmp))
1893 : pw%array(i) = pw%array(i)*tmp2
1894 : END DO
1895 : !$OMP END PARALLEL DO
1896 :
1897 249 : CALL timestop(handle)
1898 :
1899 249 : END SUBROUTINE pw_log_deriv_mix_cl
1900 :
1901 : ! **************************************************************************************************
1902 : !> \brief Multiply all data points with a complementary cosine
1903 : !> Needed for truncated Coulomb potential
1904 : !> V(\vec r)=1/r if r<rc, 0 else
1905 : !> V(\vec g)=\frac{4*\pi}{g**2}*(1-cos(g*rc))
1906 : !> \param pw ...
1907 : !> \param rcutoff ...
1908 : !> \par History
1909 : !> Frederick Stein (07-06-2021) created
1910 : !> \author Frederick Stein (07-Jun-2021)
1911 : !> \note
1912 : !> Multiplies by complementary cosine
1913 : ! **************************************************************************************************
1914 0 : SUBROUTINE pw_truncated(pw, rcutoff)
1915 :
1916 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: pw
1917 : REAL(KIND=dp), INTENT(IN) :: rcutoff
1918 :
1919 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_truncated'
1920 :
1921 : INTEGER :: handle, i
1922 : REAL(KIND=dp) :: tmp, tmp2
1923 :
1924 0 : CALL timeset(routineN, handle)
1925 0 : CPASSERT(rcutoff >= 0)
1926 :
1927 0 : !$OMP PARALLEL DO PRIVATE(i,tmp,tmp2) DEFAULT(NONE) SHARED(pw, rcutoff)
1928 : DO i = 1, SIZE(pw%array)
1929 : tmp = SQRT(pw%pw_grid%gsq(i))*rcutoff
1930 : IF (tmp >= 0.005_dp) THEN
1931 : tmp2 = 1.0_dp - COS(tmp)
1932 : ELSE
1933 : tmp2 = tmp**2/2.0_dp*(1.0 - tmp**2/12.0_dp)
1934 : END IF
1935 : pw%array(i) = pw%array(i)*tmp2
1936 : END DO
1937 : !$OMP END PARALLEL DO
1938 :
1939 0 : CALL timestop(handle)
1940 :
1941 0 : END SUBROUTINE pw_truncated
1942 :
1943 : ! **************************************************************************************************
1944 : !> \brief Multiply all data points with the logarithmic derivative of the complementary cosine
1945 : !> This function is needed for virials using truncated Coulomb potentials
1946 : !> \param pw ...
1947 : !> \param rcutoff ...
1948 : !> \note
1949 : ! **************************************************************************************************
1950 0 : SUBROUTINE pw_log_deriv_trunc(pw, rcutoff)
1951 :
1952 : TYPE(pw_c1d_gs_type), INTENT(IN) :: pw
1953 : REAL(KIND=dp), INTENT(IN) :: rcutoff
1954 :
1955 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_log_deriv_trunc'
1956 :
1957 : INTEGER :: handle, i
1958 : REAL(KIND=dp) :: rchalf, tmp, tmp2
1959 :
1960 0 : CALL timeset(routineN, handle)
1961 0 : CPASSERT(rcutoff >= 0)
1962 :
1963 0 : rchalf = 0.5_dp*rcutoff
1964 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,tmp,tmp2) &
1965 0 : !$OMP SHARED(pw,rchalf)
1966 : DO i = 1, SIZE(pw%array)
1967 : tmp = rchalf*SQRT(pw%pw_grid%gsq(i))
1968 : ! For too small arguments, use the Taylor polynomial to prevent division by zero
1969 : IF (ABS(tmp) >= 0.0001_dp) THEN
1970 : tmp2 = 1.0_dp - tmp/TAN(tmp)
1971 : ELSE
1972 : tmp2 = tmp**2*(1.0_dp + tmp**2/15.0_dp)/3.0_dp
1973 : END IF
1974 : pw%array(i) = pw%array(i)*tmp2
1975 : END DO
1976 : !$OMP END PARALLEL DO
1977 :
1978 0 : CALL timestop(handle)
1979 :
1980 0 : END SUBROUTINE pw_log_deriv_trunc
1981 :
1982 : ! **************************************************************************************************
1983 : !> \brief Calculate the derivative of a plane wave vector
1984 : !> \param pw ...
1985 : !> \param n ...
1986 : !> \par History
1987 : !> JGH (06-10-2002) allow only for inplace derivatives
1988 : !> \author JGH (25-Feb-2001)
1989 : !> \note
1990 : !> Calculate the derivative dx^n(1) dy^n(2) dz^n(3) PW
1991 : ! **************************************************************************************************
1992 1263141 : SUBROUTINE pw_derive(pw, n)
1993 :
1994 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: pw
1995 : INTEGER, DIMENSION(3), INTENT(IN) :: n
1996 :
1997 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_derive'
1998 :
1999 : COMPLEX(KIND=dp) :: im
2000 : INTEGER :: handle, m, idx, idir
2001 :
2002 1263141 : CALL timeset(routineN, handle)
2003 :
2004 5052564 : IF (ANY(n < 0)) &
2005 0 : CPABORT("Nonnegative exponents are not supported!")
2006 :
2007 5052564 : m = SUM(n)
2008 1263141 : im = CMPLX(0.0_dp, 1.0_dp, KIND=dp)**m
2009 :
2010 5052564 : DO idir = 1, 3
2011 5052564 : IF (n(idir) == 1) THEN
2012 1263141 : !$OMP PARALLEL DO DEFAULT(NONE) SHARED(pw,idir) PRIVATE(idx)
2013 : DO idx = 1, SIZE(pw%array)
2014 : pw%array(idx) = pw%array(idx)*pw%pw_grid%g(idir, idx)
2015 : END DO
2016 : !$OMP END PARALLEL DO
2017 2526282 : ELSE IF (n(idir) > 1) THEN
2018 450 : !$OMP PARALLEL DO DEFAULT(NONE) SHARED(n, pw,idir) PRIVATE(idx)
2019 : DO idx = 1, SIZE(pw%array)
2020 : pw%array(idx) = pw%array(idx)*pw%pw_grid%g(idir, idx)**n(idir)
2021 : END DO
2022 : !$OMP END PARALLEL DO
2023 : END IF
2024 : END DO
2025 :
2026 : ! im can take the values 1, -1, i, -i
2027 : ! skip this if im == 1
2028 1263141 : IF (ABS(REAL(im, KIND=dp) - 1.0_dp) > 1.0E-10_dp) THEN
2029 1263141 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(im, pw)
2030 : pw%array(:) = im*pw%array(:)
2031 : !$OMP END PARALLEL WORKSHARE
2032 : END IF
2033 :
2034 1263141 : CALL timestop(handle)
2035 :
2036 1263141 : END SUBROUTINE pw_derive
2037 :
2038 : ! **************************************************************************************************
2039 : !> \brief Calculate the Laplacian of a plane wave vector
2040 : !> \param pw ...
2041 : !> \par History
2042 : !> Frederick Stein (01-02-2022) created
2043 : !> \author JGH (25-Feb-2001)
2044 : !> \note
2045 : !> Calculate the derivative DELTA PW
2046 : ! **************************************************************************************************
2047 2666 : SUBROUTINE pw_laplace(pw)
2048 :
2049 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: pw
2050 :
2051 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_laplace'
2052 :
2053 : INTEGER :: handle
2054 :
2055 2666 : CALL timeset(routineN, handle)
2056 :
2057 2666 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw)
2058 : pw%array(:) = -pw%array(:)*pw%pw_grid%gsq(:)
2059 : !$OMP END PARALLEL WORKSHARE
2060 :
2061 2666 : CALL timestop(handle)
2062 :
2063 2666 : END SUBROUTINE pw_laplace
2064 :
2065 : ! **************************************************************************************************
2066 : !> \brief Calculate the tensorial 2nd derivative of a plane wave vector
2067 : !> \param pw ...
2068 : !> \param pwdr2 ...
2069 : !> \param i ...
2070 : !> \param j ...
2071 : !> \par History
2072 : !> none
2073 : !> \author JGH (05-May-2006)
2074 : !> \note
2075 : ! **************************************************************************************************
2076 198 : SUBROUTINE pw_dr2(pw, pwdr2, i, j)
2077 :
2078 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: pw, pwdr2
2079 : INTEGER, INTENT(IN) :: i, j
2080 :
2081 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_dr2'
2082 :
2083 : INTEGER :: cnt, handle, ig
2084 : REAL(KIND=dp) :: gg, o3
2085 :
2086 198 : CALL timeset(routineN, handle)
2087 :
2088 198 : o3 = 1.0_dp/3.0_dp
2089 :
2090 198 : cnt = SIZE(pw%array)
2091 :
2092 198 : IF (i == j) THEN
2093 108 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(ig,gg) SHARED(cnt, i, o3, pw, pwdr2)
2094 : DO ig = 1, cnt
2095 : gg = pw%pw_grid%g(i, ig)**2 - o3*pw%pw_grid%gsq(ig)
2096 : pwdr2%array(ig) = gg*pw%array(ig)
2097 : END DO
2098 : !$OMP END PARALLEL DO
2099 : ELSE
2100 90 : !$OMP PARALLEL DO PRIVATE (ig) DEFAULT(NONE) SHARED(cnt, i, j, pw, pwdr2)
2101 : DO ig = 1, cnt
2102 : pwdr2%array(ig) = pw%array(ig)*(pw%pw_grid%g(i, ig)*pw%pw_grid%g(j, ig))
2103 : END DO
2104 : !$OMP END PARALLEL DO
2105 : END IF
2106 :
2107 198 : CALL timestop(handle)
2108 :
2109 198 : END SUBROUTINE pw_dr2
2110 :
2111 : ! **************************************************************************************************
2112 : !> \brief Calculate the tensorial 2nd derivative of a plane wave vector
2113 : !> and divides by |G|^2. pwdr2_gg(G=0) is put to zero.
2114 : !> \param pw ...
2115 : !> \param pwdr2_gg ...
2116 : !> \param i ...
2117 : !> \param j ...
2118 : !> \par History
2119 : !> none
2120 : !> \author RD (20-Nov-2006)
2121 : !> \note
2122 : !> Adapted from pw_dr2
2123 : ! **************************************************************************************************
2124 12 : SUBROUTINE pw_dr2_gg(pw, pwdr2_gg, i, j)
2125 :
2126 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: pw, pwdr2_gg
2127 : INTEGER, INTENT(IN) :: i, j
2128 :
2129 : INTEGER :: cnt, handle, ig
2130 : REAL(KIND=dp) :: gg, o3
2131 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_dr2_gg'
2132 :
2133 12 : CALL timeset(routineN, handle)
2134 :
2135 12 : o3 = 1.0_dp/3.0_dp
2136 :
2137 12 : cnt = SIZE(pw%array)
2138 :
2139 12 : IF (i == j) THEN
2140 6 : !$OMP PARALLEL DO PRIVATE (ig) DEFAULT(NONE) PRIVATE(gg) SHARED(cnt, i, o3, pw, pwdr2_gg)
2141 : DO ig = pw%pw_grid%first_gne0, cnt
2142 : gg = pw%pw_grid%g(i, ig)**2 - o3*pw%pw_grid%gsq(ig)
2143 : pwdr2_gg%array(ig) = gg/pw%pw_grid%gsq(ig)*pw%array(ig)
2144 : END DO
2145 : !$OMP END PARALLEL DO
2146 : ELSE
2147 6 : !$OMP PARALLEL DO PRIVATE (ig) DEFAULT(NONE) SHARED(cnt, i, j, pw, pwdr2_gg)
2148 : DO ig = pw%pw_grid%first_gne0, cnt
2149 : pwdr2_gg%array(ig) = pw%array(ig)*((pw%pw_grid%g(i, ig)*pw%pw_grid%g(j, ig)) &
2150 : /pw%pw_grid%gsq(ig))
2151 : END DO
2152 : !$OMP END PARALLEL DO
2153 : END IF
2154 :
2155 12 : IF (pw%pw_grid%have_g0) pwdr2_gg%array(1) = 0.0_dp
2156 :
2157 12 : CALL timestop(handle)
2158 :
2159 12 : END SUBROUTINE pw_dr2_gg
2160 :
2161 : ! **************************************************************************************************
2162 : !> \brief Multiplies a G-space function with a smoothing factor of the form
2163 : !> f(|G|) = exp((ecut - G^2)/sigma)/(1+exp((ecut - G^2)/sigma))
2164 : !> \param pw ...
2165 : !> \param ecut ...
2166 : !> \param sigma ...
2167 : !> \par History
2168 : !> none
2169 : !> \author JGH (09-June-2006)
2170 : !> \note
2171 : ! **************************************************************************************************
2172 30 : SUBROUTINE pw_smoothing(pw, ecut, sigma)
2173 :
2174 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: pw
2175 : REAL(KIND=dp), INTENT(IN) :: ecut, sigma
2176 :
2177 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_smoothing'
2178 :
2179 : INTEGER :: cnt, handle, ig
2180 : REAL(KIND=dp) :: arg, f
2181 :
2182 30 : CALL timeset(routineN, handle)
2183 :
2184 30 : cnt = SIZE(pw%array)
2185 :
2186 30 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(ig, arg, f) SHARED(cnt, ecut, pw, sigma)
2187 : DO ig = 1, cnt
2188 : arg = (ecut - pw%pw_grid%gsq(ig))/sigma
2189 : f = EXP(arg)/(1 + EXP(arg))
2190 : pw%array(ig) = f*pw%array(ig)
2191 : END DO
2192 : !$OMP END PARALLEL DO
2193 :
2194 30 : CALL timestop(handle)
2195 :
2196 30 : END SUBROUTINE pw_smoothing
2197 :
2198 : ! **************************************************************************************************
2199 : !> \brief ...
2200 : !> \param grida ...
2201 : !> \param gridb ...
2202 : !> \return ...
2203 : ! **************************************************************************************************
2204 724158 : ELEMENTAL FUNCTION pw_compatible(grida, gridb) RESULT(compat)
2205 :
2206 : TYPE(pw_grid_type), INTENT(IN) :: grida, gridb
2207 : LOGICAL :: compat
2208 :
2209 724158 : compat = .FALSE.
2210 :
2211 724158 : IF (grida%id_nr == gridb%id_nr) THEN
2212 : compat = .TRUE.
2213 724158 : ELSE IF (grida%reference == gridb%id_nr) THEN
2214 : compat = .TRUE.
2215 1636 : ELSE IF (gridb%reference == grida%id_nr) THEN
2216 1106 : compat = .TRUE.
2217 : END IF
2218 :
2219 724158 : END FUNCTION pw_compatible
2220 :
2221 : ! **************************************************************************************************
2222 : !> \brief Calculate the structure factor for point r
2223 : !> \param sf ...
2224 : !> \param r ...
2225 : !> \par History
2226 : !> none
2227 : !> \author JGH (05-May-2006)
2228 : !> \note
2229 : ! **************************************************************************************************
2230 76 : SUBROUTINE pw_structure_factor(sf, r)
2231 :
2232 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: sf
2233 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: r
2234 :
2235 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_structure_factor'
2236 :
2237 : INTEGER :: cnt, handle, ig
2238 : REAL(KIND=dp) :: arg
2239 :
2240 76 : CALL timeset(routineN, handle)
2241 :
2242 76 : cnt = SIZE(sf%array)
2243 :
2244 76 : !$OMP PARALLEL DO PRIVATE (ig, arg) DEFAULT(NONE) SHARED(cnt, r, sf)
2245 : DO ig = 1, cnt
2246 : arg = DOT_PRODUCT(sf%pw_grid%g(:, ig), r)
2247 : sf%array(ig) = CMPLX(COS(arg), -SIN(arg), KIND=dp)
2248 : END DO
2249 : !$OMP END PARALLEL DO
2250 :
2251 76 : CALL timestop(handle)
2252 :
2253 76 : END SUBROUTINE pw_structure_factor
2254 :
2255 : END MODULE pw_methods
|
