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