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