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