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 1111502 : 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 1111502 : CALL timeset(routineN, handle)
248 :
249 1111502 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw)
250 : pw%array = ${type2type("0.0_dp", "r3d", kind)}$
251 : !$OMP END PARALLEL WORKSHARE
252 :
253 1111502 : CALL timestop(handle)
254 :
255 1111502 : 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 452604 : 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 452604 : CALL timeset(routineN, handle)
274 :
275 452604 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(A, pw)
276 : pw%array = a*pw%array
277 : !$OMP END PARALLEL WORKSHARE
278 :
279 452604 : CALL timestop(handle)
280 :
281 452604 : 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 421678 : 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 421678 : iop = 0
335 :
336 421678 : IF (PRESENT(oprt)) THEN
337 : SELECT CASE (oprt)
338 : CASE ("ABS", "abs")
339 0 : iop = 1
340 : CASE DEFAULT
341 2016 : CPABORT("Unknown operator")
342 : END SELECT
343 : END IF
344 :
345 101161 : total_fun = 0.0_dp
346 :
347 : #:if space == "rs"
348 : ! do reduction using maximum accuracy
349 : IF (iop == 1) THEN
350 106540096 : total_fun = fun%pw_grid%dvol*accurate_sum(ABS(fun%array))
351 : ELSE
352 318501 : 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 101161 : 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 421678 : IF (fun%pw_grid%para%mode /= PW_MODE_LOCAL) THEN
367 416234 : CALL fun%pw_grid%para%group%sum(total_fun)
368 : END IF
369 :
370 421678 : IF (PRESENT(isign)) THEN
371 305732 : total_fun = total_fun*SIGN(1._dp, REAL(isign, dp))
372 : END IF
373 :
374 421678 : 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 1291688 : 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 1291688 : CALL timeset(routineN, handle)
418 :
419 1291688 : 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 1291688 : 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 1291688 : !$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 1291688 : CALL timestop(handle)
455 :
456 1291688 : END SUBROUTINE pw_gather_p_${kind}$
457 :
458 : ! **************************************************************************************************
459 : !> \brief ...
460 : !> \param pw ...
461 : !> \param c ...
462 : !> \param scale ...
463 : ! **************************************************************************************************
464 1324472 : 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 1324472 : CALL timeset(routineN, handle)
474 :
475 1324472 : 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 41162456642 : IF (.NOT. PRESENT(scale)) c = z_zero
483 :
484 2648944 : 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 1324472 : !$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 1324472 : 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 112444 : 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 112444 : !$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 1324472 : CALL timestop(handle)
546 :
547 1324472 : 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 2956486 : 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 2956486 : CALL timeset(routineN, handle)
579 :
580 2956486 : 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 610319 : IF (pw1%pw_grid%spherical) &
584 0 : CPABORT("Spherical grids only exist in reciprocal space!")
585 : #:endif
586 :
587 : #:if kind[1]=='1'
588 2346167 : IF (.NOT. ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) THEN
589 649830 : 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 649830 : ng1 = SIZE(pw1%array)
607 649830 : ng2 = SIZE(pw2%array)
608 649830 : ns = 2*MAX(ng1, ng2)
609 :
610 649830 : IF ((pw1%pw_grid%id_nr == pw2%pw_grid%reference)) THEN
611 649462 : IF (ng1 >= ng2) THEN
612 649462 : !$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 1696337 : !$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 2441276 : IF (ANY(SHAPE(pw2%array) /= SHAPE(pw1%array))) &
661 610319 : CPABORT("3D grids must be compatible!")
662 610319 : IF (pw1%pw_grid%spherical) &
663 0 : CPABORT("3D grids must not be spherical!")
664 610319 : !$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 2956486 : CALL timestop(handle)
670 :
671 2956486 : END SUBROUTINE pw_copy_${kind}$_${kind2}$_${space}$
672 :
673 : ! **************************************************************************************************
674 : !> \brief ...
675 : !> \param pw ...
676 : !> \param array ...
677 : ! **************************************************************************************************
678 1571627 : 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 1571627 : 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 1571627 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw, array)
694 : array(:, :, :) = ${type2type("pw%array(:, :, :)", kind, kind2)}$
695 : !$OMP END PARALLEL WORKSHARE
696 : #:endif
697 :
698 1571627 : CALL timestop(handle)
699 1571627 : END SUBROUTINE pw_copy_to_array_${kind}$_${kind2}$_${space}$
700 :
701 : ! **************************************************************************************************
702 : !> \brief ...
703 : !> \param pw ...
704 : !> \param array ...
705 : ! **************************************************************************************************
706 1620366 : 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 1620366 : CALL timeset(routineN, handle)
715 :
716 1620366 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw, array)
717 : pw%array = ${type2type("array", kind2, kind)}$
718 : !$OMP END PARALLEL WORKSHARE
719 :
720 1620366 : CALL timestop(handle)
721 1620366 : 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 1856728 : 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 1856728 : CALL timeset(routineN, handle)
756 :
757 1856728 : my_alpha = 1.0_dp
758 1856728 : IF (PRESENT(alpha)) my_alpha = alpha
759 :
760 1856728 : my_beta = 1.0_dp
761 1856728 : IF (PRESENT(beta)) my_beta = beta
762 :
763 1856728 : my_allow_noncompatible_grids = .FALSE.
764 1856728 : IF (PRESENT(allow_noncompatible_grids)) my_allow_noncompatible_grids = allow_noncompatible_grids
765 :
766 1856728 : IF (my_beta /= 1.0_dp) THEN
767 94749 : IF (my_beta == 0.0_dp) THEN
768 1588 : CALL pw_zero(pw2)
769 : ELSE
770 93161 : !$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 1355050 : IF (ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) THEN
778 :
779 638742 : IF (my_alpha == 1.0_dp) THEN
780 620646 : !$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 716308 : ELSE IF (pw_compatible(pw1%pw_grid, pw2%pw_grid) .OR. my_allow_noncompatible_grids) THEN
790 :
791 716308 : ng1 = SIZE(pw1%array)
792 716308 : ng2 = SIZE(pw2%array)
793 716308 : ng = MIN(ng1, ng2)
794 :
795 716308 : 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 716308 : 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 715210 : ELSE IF ((pw2%pw_grid%id_nr == pw1%pw_grid%reference)) THEN
844 715210 : 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 715210 : IF (my_alpha == 1.0_dp) THEN
862 715210 : !$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 501678 : IF (ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) THEN
882 501148 : IF (my_alpha == 1.0_dp) THEN
883 288586 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw1, pw2)
884 : pw2%array = pw2%array + ${type2type("pw1%array", kind, kind2)}$
885 : !$OMP END PARALLEL WORKSHARE
886 212562 : ELSE IF (my_alpha /= 0.0_dp) THEN
887 211128 : !$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 530 : 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 1856728 : CALL timestop(handle)
915 :
916 1856728 : 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 3945 : 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 3945 : CALL timeset(routineN, handle)
940 :
941 3945 : my_alpha = 1.0_dp
942 3945 : IF (PRESENT(alpha)) my_alpha = alpha
943 :
944 7890 : 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 3945 : IF (my_alpha == 1.0_dp) THEN
948 3881 : !$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 3945 : CALL timestop(handle)
958 :
959 3945 : END SUBROUTINE pw_multiply_${kind}$_${kind2}$_${space}$
960 :
961 : ! **************************************************************************************************
962 : !> \brief ...
963 : !> \param pw1 ...
964 : !> \param pw2 ...
965 : ! **************************************************************************************************
966 385014 : 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 385014 : CALL timeset(routineN, handle)
975 :
976 385014 : IF (.NOT. ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) &
977 0 : CPABORT("Incompatible grids!")
978 :
979 385014 : !$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 385014 : CALL timestop(handle)
984 :
985 385014 : 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 698811 : 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 698811 : CALL timeset(routineN, handle)
1014 :
1015 698811 : loc_sumtype = do_accurate_sum
1016 698811 : IF (PRESENT(sumtype)) loc_sumtype = sumtype
1017 :
1018 698811 : my_local_only = .FALSE.
1019 698811 : IF (PRESENT(local_only)) my_local_only = local_only
1020 :
1021 698811 : IF (.NOT. ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) THEN
1022 0 : CPABORT("Grids incompatible")
1023 : END IF
1024 :
1025 698811 : my_just_sum = .FALSE.
1026 698811 : IF (PRESENT(just_sum)) my_just_sum = just_sum
1027 :
1028 : ! do standard sum
1029 698811 : 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 8729104727 : *pw2%array, KIND=dp)) !? complex bit
1058 : #:endif
1059 :
1060 : END IF
1061 :
1062 698811 : 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 392511 : 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 392511 : IF (pw1%pw_grid%grid_span == HALFSPACE) THEN
1074 236628 : integral_value = 2.0_dp*integral_value
1075 236628 : IF (pw1%pw_grid%have_g0) integral_value = integral_value - &
1076 : #:if kind[0]=="c"
1077 128723 : 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 698811 : IF (.NOT. my_local_only .AND. pw1%pw_grid%para%mode == PW_MODE_DISTRIBUTED) &
1087 381910 : CALL pw1%pw_grid%para%group%sum(integral_value)
1088 :
1089 698811 : CALL timestop(handle)
1090 :
1091 698811 : 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 3193065 : 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 3193065 : COMPLEX(KIND=dp), DIMENSION(:, :), CONTIGUOUS, POINTER :: grays
1169 3193065 : 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 3193065 : INTEGER, DIMENSION(:), POINTER :: n
1179 : LOGICAL :: test
1180 : REAL(KIND=dp) :: norm
1181 :
1182 3193065 : CALL timeset(routineN, handle2)
1183 3193065 : out_unit = cp_logger_get_default_io_unit()
1184 : CALL timeset(routineN//"_"//TRIM(ADJUSTL(cp_to_string( &
1185 3193065 : CEILING(pw1%pw_grid%cutoff/10)*10))), handle)
1186 :
1187 3193065 : NULLIFY (c_in)
1188 3193065 : NULLIFY (c_out)
1189 :
1190 3193065 : IF (PRESENT(debug)) THEN
1191 1072 : test = debug
1192 : ELSE
1193 3191993 : test = .FALSE.
1194 : END IF
1195 :
1196 : !..check if grids are compatible
1197 3193065 : IF (.NOT. ASSOCIATED(pw1%pw_grid, pw2%pw_grid)) THEN
1198 0 : IF (pw1%pw_grid%dvol /= pw2%pw_grid%dvol) THEN
1199 0 : CPABORT("PW grids not compatible")
1200 : END IF
1201 0 : IF (pw1%pw_grid%para%group /= pw2%pw_grid%para%group) THEN
1202 0 : CPABORT("PW grids have not compatible MPI groups")
1203 : END IF
1204 : END IF
1205 :
1206 : !..prepare input
1207 : #:if space == "rs"
1208 1572163 : norm = 1.0_dp/pw1%pw_grid%ngpts
1209 : #:elif space == "gs"
1210 1620902 : norm = 1.0_dp
1211 : #:else
1212 : #:stop "Unknown space tag: "+space
1213 : #:endif
1214 :
1215 3193065 : n => pw1%pw_grid%npts
1216 :
1217 3193065 : IF (pw1%pw_grid%para%mode == PW_MODE_LOCAL) THEN
1218 :
1219 : !
1220 : !..replicated data, use local FFT
1221 : !
1222 :
1223 576905 : IF (test .AND. out_unit > 0) THEN
1224 0 : WRITE (out_unit, '(A)') " FFT Protocol "
1225 : #:if space=="rs"
1226 0 : WRITE (out_unit, '(A,T76,A)') " Transform direction ", "FWFFT"
1227 0 : WRITE (out_unit, '(A,T72,A)') " in space ", "REALSPACE"
1228 0 : WRITE (out_unit, '(A,T72,A)') " out space ", "REALSPACE"
1229 : #:else
1230 0 : WRITE (out_unit, '(A,T76,A)') " Transform direction ", "BWFFT"
1231 0 : WRITE (out_unit, '(A,T66,A)') " in space ", "RECIPROCALSPACE"
1232 0 : WRITE (out_unit, '(A,T66,A)') " out space ", "RECIPROCALSPACE"
1233 : #:endif
1234 0 : WRITE (out_unit, '(A,T66,E15.6)') " scale factor", norm
1235 : END IF
1236 :
1237 : #:if space=="rs"
1238 : #:if kind==kind2=="c3d"
1239 0 : c_in => pw1%array
1240 0 : c_out => pw2%array
1241 0 : CALL fft3d(FWFFT, n, c_in, c_out, scale=norm, debug=test)
1242 : #:elif kind=="r3d" and kind2=="c3d"
1243 0 : pw2%array = CMPLX(pw1%array, 0.0_dp, KIND=dp)
1244 0 : c_out => pw2%array
1245 0 : CALL fft3d(FWFFT, n, c_out, scale=norm, debug=test)
1246 : #:elif kind=="c3d" and kind2=="c1d"
1247 0 : c_in => pw1%array
1248 0 : ALLOCATE (c_out(n(1), n(2), n(3)))
1249 : ! transform
1250 0 : CALL fft3d(FWFFT, n, c_in, c_out, scale=norm, debug=test)
1251 : ! gather results
1252 0 : IF (test .AND. out_unit > 0) WRITE (out_unit, '(A)') " PW_GATHER : 3d -> 1d "
1253 0 : CALL pw_gather_s_${kind2}$_c3d(pw2, c_out)
1254 0 : DEALLOCATE (c_out)
1255 : #:elif kind=="r3d" and kind2=="c1d"
1256 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
1257 : CALL pw_gpu_r3dc1d_3d(pw1, pw2, scale=norm)
1258 : #elif defined (__PW_FPGA)
1259 : ALLOCATE (c_out(n(1), n(2), n(3)))
1260 : ! check if bitstream for the fft size is present
1261 : ! if not, perform fft3d in CPU
1262 : IF (pw_fpga_init_bitstream(n) == 1) THEN
1263 : CALL pw_copy_to_array(pw1, c_out)
1264 : #if (__PW_FPGA_SP && __PW_FPGA)
1265 : CALL pw_fpga_r3dc1d_3d_sp(n, c_out)
1266 : #else
1267 : CALL pw_fpga_r3dc1d_3d_dp(n, c_out)
1268 : #endif
1269 : CALL zdscal(n(1)*n(2)*n(3), norm, c_out, 1)
1270 : CALL pw_gather_s_${kind2}$_c3d(pw2, c_out)
1271 : ELSE
1272 : CALL pw_copy_to_array(pw1, c_out)
1273 : CALL fft3d(FWFFT, n, c_out, scale=norm, debug=test)
1274 : CALL pw_gather_s_${kind2}$_c3d(pw2, c_out)
1275 : END IF
1276 : DEALLOCATE (c_out)
1277 : #else
1278 1402375 : ALLOCATE (c_out(n(1), n(2), n(3)))
1279 2605564036 : c_out = 0.0_dp
1280 280475 : CALL pw_copy_to_array(pw1, c_out)
1281 280475 : CALL fft3d(FWFFT, n, c_out, scale=norm, debug=test)
1282 280475 : CALL pw_gather_s_${kind2}$_c3d(pw2, c_out)
1283 280475 : DEALLOCATE (c_out)
1284 : #endif
1285 : #:endif
1286 : #:else
1287 : #:if kind=="c3d" and kind2=="c3d"
1288 0 : c_in => pw1%array
1289 0 : c_out => pw2%array
1290 0 : CALL fft3d(BWFFT, n, c_in, c_out, scale=norm, debug=test)
1291 : #:elif kind=="c3d" and kind2=="r3d"
1292 0 : c_in => pw1%array
1293 0 : ALLOCATE (c_out(n(1), n(2), n(3)))
1294 0 : CALL fft3d(BWFFT, n, c_in, c_out, scale=norm, debug=test)
1295 : ! use real part only
1296 0 : IF (test .AND. out_unit > 0) WRITE (out_unit, '(A)') " REAL part "
1297 0 : pw2%array = REAL(c_out, KIND=dp)
1298 0 : DEALLOCATE (c_out)
1299 : #:elif kind=="c1d" and kind2=="c3d"
1300 0 : c_out => pw2%array
1301 0 : IF (test .AND. out_unit > 0) WRITE (out_unit, '(A)') " PW_SCATTER : 3d -> 1d "
1302 0 : CALL pw_scatter_s_${kind}$_c3d(pw1, c_out)
1303 0 : CALL fft3d(BWFFT, n, c_out, scale=norm, debug=test)
1304 : #:elif kind=="c1d" and kind2=="r3d"
1305 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
1306 : CALL pw_gpu_c1dr3d_3d(pw1, pw2, scale=norm)
1307 : #elif defined (__PW_FPGA)
1308 : ALLOCATE (c_out(n(1), n(2), n(3)))
1309 : ! check if bitstream for the fft size is present
1310 : ! if not, perform fft3d in CPU
1311 : IF (pw_fpga_init_bitstream(n) == 1) THEN
1312 : CALL pw_scatter_s_${kind}$_c3d(pw1, c_out)
1313 : ! transform using FPGA
1314 : #if (__PW_FPGA_SP && __PW_FPGA)
1315 : CALL pw_fpga_c1dr3d_3d_sp(n, c_out)
1316 : #else
1317 : CALL pw_fpga_c1dr3d_3d_dp(n, c_out)
1318 : #endif
1319 : CALL zdscal(n(1)*n(2)*n(3), norm, c_out, 1)
1320 : ! use real part only
1321 : CALL pw_copy_from_array(pw2, c_out)
1322 : ELSE
1323 : IF (test .AND. out_unit > 0) WRITE (out_unit, '(A)') " PW_SCATTER : 3d -> 1d "
1324 : CALL pw_scatter_s_${kind}$_c3d(pw1, c_out)
1325 : ! transform
1326 : CALL fft3d(BWFFT, n, c_out, scale=norm, debug=test)
1327 : ! use real part only
1328 : IF (test .AND. out_unit > 0) WRITE (out_unit, '(A)') " REAL part "
1329 : CALL pw_copy_from_array(pw2, c_out)
1330 : END IF
1331 : DEALLOCATE (c_out)
1332 : #else
1333 1482150 : ALLOCATE (c_out(n(1), n(2), n(3)))
1334 296430 : IF (test .AND. out_unit > 0) WRITE (out_unit, '(A)') " PW_SCATTER : 3d -> 1d "
1335 296430 : CALL pw_scatter_s_${kind}$_c3d(pw1, c_out)
1336 : ! transform
1337 296430 : CALL fft3d(BWFFT, n, c_out, scale=norm, debug=test)
1338 : ! use real part only
1339 296430 : IF (test .AND. out_unit > 0) WRITE (out_unit, '(A)') " REAL part "
1340 296430 : CALL pw_copy_from_array(pw2, c_out)
1341 296430 : DEALLOCATE (c_out)
1342 : #endif
1343 : #:endif
1344 : #:endif
1345 :
1346 576905 : IF (test .AND. out_unit > 0) WRITE (out_unit, '(A)') " End of FFT Protocol "
1347 :
1348 : ELSE
1349 :
1350 : !
1351 : !..parallel FFT
1352 : !
1353 :
1354 2616160 : IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0) THEN
1355 8 : WRITE (out_unit, '(A)') " FFT Protocol "
1356 : #:if space=="rs"
1357 4 : WRITE (out_unit, '(A,T76,A)') " Transform direction ", "FWFFT"
1358 4 : WRITE (out_unit, '(A,T72,A)') " in space ", "REALSPACE"
1359 4 : WRITE (out_unit, '(A,T72,A)') " out space ", "REALSPACE"
1360 : #:else
1361 4 : WRITE (out_unit, '(A,T76,A)') " Transform direction ", "BWFFT"
1362 4 : WRITE (out_unit, '(A,T66,A)') " in space ", "RECIPROCALSPACE"
1363 4 : WRITE (out_unit, '(A,T66,A)') " out space ", "RECIPROCALSPACE"
1364 : #:endif
1365 8 : WRITE (out_unit, '(A,T66,E15.6)') " scale factor", norm
1366 : END IF
1367 :
1368 2616160 : my_pos = pw1%pw_grid%para%my_pos
1369 2616160 : nrays = pw1%pw_grid%para%nyzray(my_pos)
1370 2616160 : grays => pw1%pw_grid%grays
1371 :
1372 : #:if space=="rs"
1373 : #:if kind=="c3d" and kind2=="c1d"
1374 : !..prepare input
1375 536 : c_in => pw1%array
1376 1918288 : grays = z_zero
1377 : !..transform
1378 536 : IF (pw1%pw_grid%para%ray_distribution) THEN
1379 : CALL fft3d(FWFFT, n, c_in, grays, pw1%pw_grid%para%group, &
1380 : pw1%pw_grid%para%rs_group, &
1381 : pw1%pw_grid%para%yzp, pw1%pw_grid%para%nyzray, &
1382 432 : pw1%pw_grid%para%bo, scale=norm, debug=test)
1383 : ELSE
1384 : CALL fft3d(FWFFT, n, c_in, grays, pw1%pw_grid%para%rs_group, &
1385 104 : pw1%pw_grid%para%bo, scale=norm, debug=test)
1386 : END IF
1387 : !..prepare output
1388 536 : IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0) &
1389 4 : WRITE (out_unit, '(A)') " PW_GATHER : 2d -> 1d "
1390 536 : CALL pw_gather_p_${kind2}$ (pw2, grays)
1391 : #:elif kind=="r3d" and kind2=="c1d"
1392 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
1393 : ! (no ray dist. is not efficient in CUDA)
1394 : use_pw_gpu = pw1%pw_grid%para%ray_distribution
1395 : IF (use_pw_gpu) THEN
1396 : CALL pw_gpu_r3dc1d_3d_ps(pw1, pw2, scale=norm)
1397 : ELSE
1398 : #endif
1399 : !.. prepare input
1400 5164608 : nloc = pw1%pw_grid%npts_local
1401 6455760 : ALLOCATE (c_in(nloc(1), nloc(2), nloc(3)))
1402 1291152 : CALL pw_copy_to_array(pw1, c_in)
1403 37106806021 : grays = z_zero
1404 : !..transform
1405 1291152 : IF (pw1%pw_grid%para%ray_distribution) THEN
1406 : CALL fft3d(FWFFT, n, c_in, grays, pw1%pw_grid%para%group, &
1407 : pw1%pw_grid%para%rs_group, &
1408 : pw1%pw_grid%para%yzp, pw1%pw_grid%para%nyzray, &
1409 1291152 : pw1%pw_grid%para%bo, scale=norm, debug=test)
1410 : ELSE
1411 : CALL fft3d(FWFFT, n, c_in, grays, pw1%pw_grid%para%rs_group, &
1412 0 : pw1%pw_grid%para%bo, scale=norm, debug=test)
1413 : END IF
1414 : !..prepare output
1415 1291152 : IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0) &
1416 0 : WRITE (out_unit, '(A)') " PW_GATHER : 2d -> 1d "
1417 1291152 : CALL pw_gather_p_${kind2}$ (pw2, grays)
1418 1291152 : DEALLOCATE (c_in)
1419 :
1420 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
1421 : END IF
1422 : #endif
1423 : #:endif
1424 : #:else
1425 : #:if kind=="c1d" and kind2=="c3d"
1426 : !..prepare input
1427 536 : IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0) &
1428 4 : WRITE (out_unit, '(A)') " PW_SCATTER : 2d -> 1d "
1429 1918288 : grays = z_zero
1430 536 : CALL pw_scatter_p_${kind}$ (pw1, grays)
1431 536 : c_in => pw2%array
1432 : !..transform
1433 536 : IF (pw1%pw_grid%para%ray_distribution) THEN
1434 : CALL fft3d(BWFFT, n, c_in, grays, pw1%pw_grid%para%group, &
1435 : pw1%pw_grid%para%rs_group, &
1436 : pw1%pw_grid%para%yzp, pw1%pw_grid%para%nyzray, &
1437 432 : pw1%pw_grid%para%bo, scale=norm, debug=test)
1438 : ELSE
1439 : CALL fft3d(BWFFT, n, c_in, grays, pw1%pw_grid%para%rs_group, &
1440 104 : pw1%pw_grid%para%bo, scale=norm, debug=test)
1441 : END IF
1442 : !..prepare output (nothing to do)
1443 : #:elif kind=="c1d" and kind2=="r3d"
1444 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
1445 : ! (no ray dist. is not efficient in CUDA)
1446 : use_pw_gpu = pw1%pw_grid%para%ray_distribution
1447 : IF (use_pw_gpu) THEN
1448 : CALL pw_gpu_c1dr3d_3d_ps(pw1, pw2, scale=norm)
1449 : ELSE
1450 : #endif
1451 : !.. prepare input
1452 1323936 : IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0) &
1453 0 : WRITE (out_unit, '(A)') " PW_SCATTER : 2d -> 1d "
1454 41160538354 : grays = z_zero
1455 1323936 : CALL pw_scatter_p_${kind}$ (pw1, grays)
1456 5295744 : nloc = pw2%pw_grid%npts_local
1457 6619680 : ALLOCATE (c_in(nloc(1), nloc(2), nloc(3)))
1458 : !..transform
1459 1323936 : IF (pw1%pw_grid%para%ray_distribution) THEN
1460 : CALL fft3d(BWFFT, n, c_in, grays, pw1%pw_grid%para%group, &
1461 : pw1%pw_grid%para%rs_group, &
1462 : pw1%pw_grid%para%yzp, pw1%pw_grid%para%nyzray, &
1463 1323936 : pw1%pw_grid%para%bo, scale=norm, debug=test)
1464 : ELSE
1465 : CALL fft3d(BWFFT, n, c_in, grays, pw1%pw_grid%para%rs_group, &
1466 0 : pw1%pw_grid%para%bo, scale=norm, debug=test)
1467 : END IF
1468 : !..prepare output
1469 1323936 : IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0) &
1470 0 : WRITE (out_unit, '(A)') " Real part "
1471 1323936 : CALL pw_copy_from_array(pw2, c_in)
1472 1323936 : 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 3193065 : IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0) THEN
1481 8 : WRITE (out_unit, '(A)') " End of FFT Protocol "
1482 : END IF
1483 :
1484 3193065 : CALL timestop(handle)
1485 3193065 : CALL timestop(handle2)
1486 :
1487 3193065 : 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 280475 : 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 280475 : 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 280475 : 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 280475 : !$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 280475 : CALL timestop(handle)
1560 :
1561 280475 : 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 296430 : 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 296430 : 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 2516498065 : IF (.NOT. PRESENT(scale)) c = 0.0_dp
1608 :
1609 592860 : 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 296430 : !$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 296430 : 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 296430 : CALL timestop(handle)
1661 :
1662 296430 : 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
1687 : REAL(KIND=dp) :: omega_2
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 WORKSHARE DEFAULT(NONE) SHARED(pw,omega_2)
1696 : pw%array = pw%array*EXP(-pw%pw_grid%gsq*omega_2)
1697 : !$OMP END PARALLEL WORKSHARE
1698 :
1699 3841 : CALL timestop(handle)
1700 :
1701 3841 : END SUBROUTINE pw_gauss_damp
1702 :
1703 : ! **************************************************************************************************
1704 : !> \brief Multiply all data points with the logarithmic derivative of a Gaussian
1705 : !> \param pw ...
1706 : !> \param omega ...
1707 : !> \note
1708 : !> Performs a Gaussian damping
1709 : ! **************************************************************************************************
1710 1329 : SUBROUTINE pw_log_deriv_gauss(pw, omega)
1711 :
1712 : TYPE(pw_c1d_gs_type), INTENT(IN) :: pw
1713 : REAL(KIND=dp), INTENT(IN) :: omega
1714 :
1715 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_log_deriv_gauss'
1716 :
1717 : INTEGER :: handle
1718 : REAL(KIND=dp) :: omega_2
1719 :
1720 1329 : CALL timeset(routineN, handle)
1721 1329 : CPASSERT(omega >= 0)
1722 :
1723 1329 : omega_2 = omega*omega
1724 1329 : omega_2 = 0.25_dp/omega_2
1725 :
1726 1329 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw,omega_2)
1727 : pw%array = pw%array*(1.0_dp + omega_2*pw%pw_grid%gsq)
1728 : !$OMP END PARALLEL WORKSHARE
1729 :
1730 1329 : CALL timestop(handle)
1731 1329 : END SUBROUTINE pw_log_deriv_gauss
1732 :
1733 : ! **************************************************************************************************
1734 : !> \brief Multiply all data points with a Gaussian damping factor
1735 : !> Needed for longrange Coulomb potential
1736 : !> V(\vec r)=erf(omega*r)/r
1737 : !> V(\vec g)=\frac{4*\pi}{g**2}*exp(-g**2/omega**2)
1738 : !> \param pw ...
1739 : !> \param omega ...
1740 : !> \par History
1741 : !> Frederick Stein (12-04-2019) created
1742 : !> \author Frederick Stein (12-Apr-2019)
1743 : !> \note
1744 : !> Performs a Gaussian damping
1745 : ! **************************************************************************************************
1746 0 : SUBROUTINE pw_compl_gauss_damp(pw, omega)
1747 :
1748 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: pw
1749 : REAL(KIND=dp), INTENT(IN) :: omega
1750 :
1751 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_compl_gauss_damp'
1752 :
1753 : INTEGER :: cnt, handle, i
1754 : REAL(KIND=dp) :: omega_2, tmp
1755 :
1756 0 : CALL timeset(routineN, handle)
1757 :
1758 0 : omega_2 = omega*omega
1759 0 : omega_2 = 0.25_dp/omega_2
1760 :
1761 0 : cnt = SIZE(pw%array)
1762 :
1763 0 : !$OMP PARALLEL DO PRIVATE(i, tmp) DEFAULT(NONE) SHARED(cnt, pw,omega_2)
1764 : DO i = 1, cnt
1765 : tmp = -omega_2*pw%pw_grid%gsq(i)
1766 : IF (ABS(tmp) > 1.0E-5_dp) THEN
1767 : pw%array(i) = pw%array(i)*(1.0_dp - EXP(tmp))
1768 : ELSE
1769 : pw%array(i) = pw%array(i)*(tmp + 0.5_dp*tmp*(tmp + (1.0_dp/3.0_dp)*tmp**2))
1770 : END IF
1771 : END DO
1772 : !$OMP END PARALLEL DO
1773 :
1774 0 : CALL timestop(handle)
1775 :
1776 0 : END SUBROUTINE pw_compl_gauss_damp
1777 :
1778 : ! **************************************************************************************************
1779 : !> \brief Multiply all data points with the logarithmic derivative of the complementary Gaussian damping factor
1780 : !> \param pw ...
1781 : !> \param omega ...
1782 : !> \note
1783 : ! **************************************************************************************************
1784 0 : SUBROUTINE pw_log_deriv_compl_gauss(pw, omega)
1785 :
1786 : TYPE(pw_c1d_gs_type), INTENT(IN) :: pw
1787 : REAL(KIND=dp), INTENT(IN) :: omega
1788 :
1789 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_log_deriv_compl_gauss'
1790 :
1791 : INTEGER :: handle, i
1792 : REAL(KIND=dp) :: omega_2, tmp
1793 :
1794 0 : CALL timeset(routineN, handle)
1795 :
1796 0 : omega_2 = omega*omega
1797 0 : omega_2 = 0.25_dp/omega_2
1798 :
1799 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,tmp) &
1800 0 : !$OMP SHARED(pw,omega_2)
1801 : DO i = 1, SIZE(pw%array)
1802 : tmp = omega_2*pw%pw_grid%gsq(i)
1803 : ! For too small arguments, use the Taylor polynomial to prevent division by zero
1804 : IF (ABS(tmp) >= 0.003_dp) THEN
1805 : pw%array(i) = pw%array(i)*(1.0_dp - tmp*EXP(-tmp)/(1.0_dp - EXP(-tmp)))
1806 : ELSE
1807 : pw%array(i) = pw%array(i)*(0.5_dp*tmp - tmp**2/12.0_dp)
1808 : END IF
1809 : END DO
1810 : !$OMP END PARALLEL DO
1811 :
1812 0 : CALL timestop(handle)
1813 :
1814 0 : END SUBROUTINE pw_log_deriv_compl_gauss
1815 :
1816 : ! **************************************************************************************************
1817 : !> \brief Multiply all data points with a Gaussian damping factor and mixes it with the original function
1818 : !> Needed for mixed longrange/Coulomb potential
1819 : !> V(\vec r)=(a+b*erf(omega*r))/r
1820 : !> V(\vec g)=\frac{4*\pi}{g**2}*(a+b*exp(-g**2/omega**2))
1821 : !> \param pw ...
1822 : !> \param omega ...
1823 : !> \param scale_coul ...
1824 : !> \param scale_long ...
1825 : !> \par History
1826 : !> Frederick Stein (16-Dec-2021) created
1827 : !> \author Frederick Stein (16-Dec-2021)
1828 : !> \note
1829 : !> Performs a Gaussian damping
1830 : ! **************************************************************************************************
1831 332 : SUBROUTINE pw_gauss_damp_mix(pw, omega, scale_coul, scale_long)
1832 :
1833 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: pw
1834 : REAL(KIND=dp), INTENT(IN) :: omega, scale_coul, scale_long
1835 :
1836 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_gauss_damp_mix'
1837 :
1838 : INTEGER :: handle
1839 : REAL(KIND=dp) :: omega_2
1840 :
1841 332 : CALL timeset(routineN, handle)
1842 :
1843 332 : omega_2 = omega*omega
1844 332 : omega_2 = 0.25_dp/omega_2
1845 :
1846 332 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw, omega_2, scale_coul, scale_long)
1847 : pw%array = pw%array*(scale_coul + scale_long*EXP(-pw%pw_grid%gsq*omega_2))
1848 : !$OMP END PARALLEL WORKSHARE
1849 :
1850 332 : CALL timestop(handle)
1851 :
1852 332 : END SUBROUTINE pw_gauss_damp_mix
1853 :
1854 : ! **************************************************************************************************
1855 : !> \brief Multiply all data points with the logarithmic derivative of the mixed longrange/Coulomb potential
1856 : !> Needed for mixed longrange/Coulomb potential
1857 : !> \param pw ...
1858 : !> \param omega ...
1859 : !> \param scale_coul ...
1860 : !> \param scale_long ...
1861 : !> \note
1862 : ! **************************************************************************************************
1863 249 : SUBROUTINE pw_log_deriv_mix_cl(pw, omega, scale_coul, scale_long)
1864 :
1865 : TYPE(pw_c1d_gs_type), INTENT(IN) :: pw
1866 : REAL(KIND=dp), INTENT(IN) :: omega, scale_coul, scale_long
1867 :
1868 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_log_deriv_mix_cl'
1869 :
1870 : INTEGER :: handle, i
1871 : REAL(KIND=dp) :: omega_2, tmp
1872 :
1873 249 : CALL timeset(routineN, handle)
1874 :
1875 249 : omega_2 = omega*omega
1876 249 : omega_2 = 0.25_dp/omega_2
1877 :
1878 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,tmp) &
1879 249 : !$OMP SHARED(pw,omega_2,scale_long,scale_coul)
1880 : DO i = 1, SIZE(pw%array)
1881 : tmp = omega_2*pw%pw_grid%gsq(i)
1882 : pw%array(i) = pw%array(i)*(1.0_dp + scale_long*tmp*EXP(-tmp)/(scale_coul + scale_long*EXP(-tmp)))
1883 : END DO
1884 : !$OMP END PARALLEL DO
1885 :
1886 249 : CALL timestop(handle)
1887 :
1888 249 : END SUBROUTINE pw_log_deriv_mix_cl
1889 :
1890 : ! **************************************************************************************************
1891 : !> \brief Multiply all data points with a complementary cosine
1892 : !> Needed for truncated Coulomb potential
1893 : !> V(\vec r)=1/r if r<rc, 0 else
1894 : !> V(\vec g)=\frac{4*\pi}{g**2}*(1-cos(g*rc))
1895 : !> \param pw ...
1896 : !> \param rcutoff ...
1897 : !> \par History
1898 : !> Frederick Stein (07-06-2021) created
1899 : !> \author Frederick Stein (07-Jun-2021)
1900 : !> \note
1901 : !> Multiplies by complementary cosine
1902 : ! **************************************************************************************************
1903 0 : SUBROUTINE pw_truncated(pw, rcutoff)
1904 :
1905 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: pw
1906 : REAL(KIND=dp), INTENT(IN) :: rcutoff
1907 :
1908 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_truncated'
1909 :
1910 : INTEGER :: handle, i
1911 : REAL(KIND=dp) :: tmp
1912 :
1913 0 : CALL timeset(routineN, handle)
1914 0 : CPASSERT(rcutoff >= 0)
1915 :
1916 0 : !$OMP PARALLEL DO PRIVATE(i,tmp) DEFAULT(NONE) SHARED(pw, rcutoff)
1917 : DO i = 1, SIZE(pw%array)
1918 : tmp = SQRT(pw%pw_grid%gsq(i))*rcutoff
1919 : IF (tmp >= 0.005_dp) THEN
1920 : pw%array(i) = pw%array(i)*(1.0_dp - COS(tmp))
1921 : ELSE
1922 : pw%array(i) = pw%array(i)*tmp**2/2.0_dp*(1.0 - tmp**2/12.0_dp)
1923 : END IF
1924 : END DO
1925 : !$OMP END PARALLEL DO
1926 :
1927 0 : CALL timestop(handle)
1928 :
1929 0 : END SUBROUTINE pw_truncated
1930 :
1931 : ! **************************************************************************************************
1932 : !> \brief Multiply all data points with the logarithmic derivative of the complementary cosine
1933 : !> This function is needed for virials using truncated Coulomb potentials
1934 : !> \param pw ...
1935 : !> \param rcutoff ...
1936 : !> \note
1937 : ! **************************************************************************************************
1938 0 : SUBROUTINE pw_log_deriv_trunc(pw, rcutoff)
1939 :
1940 : TYPE(pw_c1d_gs_type), INTENT(IN) :: pw
1941 : REAL(KIND=dp), INTENT(IN) :: rcutoff
1942 :
1943 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_log_deriv_trunc'
1944 :
1945 : INTEGER :: handle, i
1946 : REAL(KIND=dp) :: rchalf, tmp
1947 :
1948 0 : CALL timeset(routineN, handle)
1949 0 : CPASSERT(rcutoff >= 0)
1950 :
1951 0 : rchalf = 0.5_dp*rcutoff
1952 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,tmp) &
1953 0 : !$OMP SHARED(pw,rchalf)
1954 : DO i = 1, SIZE(pw%array)
1955 : tmp = rchalf*SQRT(pw%pw_grid%gsq(i))
1956 : ! For too small arguments, use the Taylor polynomial to prevent division by zero
1957 : IF (ABS(tmp) >= 0.0001_dp) THEN
1958 : pw%array(i) = pw%array(i)*(1.0_dp - tmp/TAN(tmp))
1959 : ELSE
1960 : pw%array(i) = pw%array(i)*tmp**2*(1.0_dp + tmp**2/15.0_dp)/3.0_dp
1961 : END IF
1962 : END DO
1963 : !$OMP END PARALLEL DO
1964 :
1965 0 : CALL timestop(handle)
1966 :
1967 0 : END SUBROUTINE pw_log_deriv_trunc
1968 :
1969 : ! **************************************************************************************************
1970 : !> \brief Calculate the derivative of a plane wave vector
1971 : !> \param pw ...
1972 : !> \param n ...
1973 : !> \par History
1974 : !> JGH (06-10-2002) allow only for inplace derivatives
1975 : !> \author JGH (25-Feb-2001)
1976 : !> \note
1977 : !> Calculate the derivative dx^n(1) dy^n(2) dz^n(3) PW
1978 : ! **************************************************************************************************
1979 1224399 : SUBROUTINE pw_derive(pw, n)
1980 :
1981 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: pw
1982 : INTEGER, DIMENSION(3), INTENT(IN) :: n
1983 :
1984 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_derive'
1985 :
1986 : COMPLEX(KIND=dp) :: im
1987 : INTEGER :: handle, m
1988 :
1989 1224399 : CALL timeset(routineN, handle)
1990 :
1991 4897596 : IF (ANY(n < 0)) &
1992 0 : CPABORT("Nonnegative exponents are not supported!")
1993 :
1994 4897596 : m = SUM(n)
1995 1224399 : im = CMPLX(0.0_dp, 1.0_dp, KIND=dp)**m
1996 :
1997 1224399 : IF (n(1) == 1) THEN
1998 408589 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw)
1999 : pw%array(:) = pw%array(:)*pw%pw_grid%g(1, :)
2000 : !$OMP END PARALLEL WORKSHARE
2001 815810 : ELSE IF (n(1) > 1) THEN
2002 132 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(n, pw)
2003 : pw%array(:) = pw%array(:)*(pw%pw_grid%g(1, :)**n(1))
2004 : !$OMP END PARALLEL WORKSHARE
2005 : END IF
2006 :
2007 1224399 : IF (n(2) == 1) THEN
2008 408133 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw)
2009 : pw%array(:) = pw%array(:)*pw%pw_grid%g(2, :)
2010 : !$OMP END PARALLEL WORKSHARE
2011 816266 : ELSE IF (n(2) > 1) THEN
2012 132 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(n, pw)
2013 : pw%array(:) = pw%array(:)*(pw%pw_grid%g(2, :)**n(2))
2014 : !$OMP END PARALLEL WORKSHARE
2015 : END IF
2016 :
2017 1224399 : IF (n(3) == 1) THEN
2018 407677 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw)
2019 : pw%array(:) = pw%array(:)*pw%pw_grid%g(3, :)
2020 : !$OMP END PARALLEL WORKSHARE
2021 816722 : ELSE IF (n(3) > 1) THEN
2022 132 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(n, pw)
2023 : pw%array(:) = pw%array(:)*(pw%pw_grid%g(3, :)**n(3))
2024 : !$OMP END PARALLEL WORKSHARE
2025 : END IF
2026 :
2027 : ! im can take the values 1, -1, i, -i
2028 : ! skip this if im == 1
2029 1224399 : IF (ABS(REAL(im, KIND=dp) - 1.0_dp) > 1.0E-10_dp) THEN
2030 1224399 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(im, pw)
2031 : pw%array(:) = im*pw%array(:)
2032 : !$OMP END PARALLEL WORKSHARE
2033 : END IF
2034 :
2035 1224399 : CALL timestop(handle)
2036 :
2037 1224399 : 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 2406 : 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 2406 : CALL timeset(routineN, handle)
2057 :
2058 2406 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw)
2059 : pw%array(:) = -pw%array(:)*pw%pw_grid%gsq(:)
2060 : !$OMP END PARALLEL WORKSHARE
2061 :
2062 2406 : CALL timestop(handle)
2063 :
2064 2406 : 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%array(ig)/pw%pw_grid%gsq(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 716838 : ELEMENTAL FUNCTION pw_compatible(grida, gridb) RESULT(compat)
2206 :
2207 : TYPE(pw_grid_type), INTENT(IN) :: grida, gridb
2208 : LOGICAL :: compat
2209 :
2210 716838 : compat = .FALSE.
2211 :
2212 716838 : IF (grida%id_nr == gridb%id_nr) THEN
2213 : compat = .TRUE.
2214 716838 : ELSE IF (grida%reference == gridb%id_nr) THEN
2215 : compat = .TRUE.
2216 1628 : ELSE IF (gridb%reference == grida%id_nr) THEN
2217 1098 : compat = .TRUE.
2218 : END IF
2219 :
2220 716838 : 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
|
