Branch data Line data Source code
1 : : !!****m* ABINIT/m_fftcore
2 : : !! NAME
3 : : !! m_fftcore
4 : : !!
5 : : !! FUNCTION
6 : : !! Low-level tools for FFT (sequential and MPI parallel version)
7 : : !! It also provides helper functions to set up the list of G vectors
8 : : !! inside a sphere or to count them.
9 : : !!
10 : : !! COPYRIGHT
11 : : !! Copyright (C) 2014-2026 ABINIT group (SG, XG, AR, MG, MT)
12 : : !! This file is distributed under the terms of the
13 : : !! GNU General Public License, see ~abinit/COPYING
14 : : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : : !!
16 : : !! TODO
17 : : !! 1) Pass distribfft instead of MPI_enreg to simplify the API and facilitate code-reuse.
18 : : !!
19 : : !! 2) Merge this module with m_distribfft
20 : : !!
21 : : !! 3) Get rid of paral_kgb and MPI_type! This is a low-level module that may be called by other
22 : : !! code in which paral_kgb is meaningless! FFT tables and a MPI communicator are sufficient.
23 : : !!
24 : : !! SOURCE
25 : :
26 : : #if defined HAVE_CONFIG_H
27 : : #include "config.h"
28 : : #endif
29 : :
30 : : #include "abi_common.h"
31 : :
32 : : module m_fftcore
33 : :
34 : : use defs_basis
35 : : use m_abicore
36 : : use m_errors
37 : : use m_xmpi
38 : : use m_sort
39 : :
40 : : use m_time, only : timab
41 : : use m_fstrings, only : itoa, sjoin
42 : : use m_geometry, only : normv
43 : : use defs_abitypes, only : MPI_type
44 : : use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq
45 : :
46 : : implicit none
47 : :
48 : : private
49 : :
50 : : public :: fftalg_isavailable ! True if the FFT library specified by fftalg is available.
51 : : public :: fftalg_has_mpi ! True if fftalg provides MPI-FFTs.
52 : : public :: fftalg_for_npfft ! Returns the default value for fftalg given the number of processors for the FFT.
53 : : public :: fftalg_info ! Returns strings with info on the FFT library specified by fftalg.
54 : : public :: get_cache_kb ! Returns the cache size in Kbs (based on CPP variables).
55 : : public :: ngfft_seq ! initialize ngfft(18) from the FFT divisions (assume sequential FFT)
56 : : public :: print_ngfft ! Print the content of ngfft(18) in explicative format.
57 : : public :: bound ! Find distance**2 to boundary point of fft box nearest to kpt
58 : : public :: getng ! From ecut and metric tensor in reciprocal space, computes recommended ngfft(1:3)
59 : : public :: sphereboundary ! Finds the boundary of the basis sphere of G vectors
60 : : public :: sphere
61 : : public :: sphere_fft ! Insert cg inside box.
62 : : public :: sphere_fft1 ! TODO: This one should be replaced by sphere_fft.
63 : : public :: change_istwfk ! Change the istwfk mode of a set of wavefunctions (sequential version, same k-point)
64 : : public :: kpgsph ! Set up the G vector list
65 : : public :: kpgcount ! Give the number of G vector in each direction
66 : : public :: get_kg ! Helper function to calculate the set of G-vectors at a given kpoint (no MPI FFT)
67 : : public :: kgindex ! Compute the index of each plane wave on a FFT grid.
68 : :
69 : : ! Low-level tools for MPI FFT
70 : : public :: switch
71 : : public :: switch_cent
72 : : public :: switchreal
73 : : public :: switchreal_cent
74 : : public :: scramble
75 : : public :: fill
76 : : public :: fill_cent
77 : : public :: unfill
78 : : public :: unfill_cent
79 : : public :: unmpiswitch
80 : : public :: unswitch
81 : : public :: unswitchreal_cent
82 : : public :: unmpiswitch_cent
83 : : public :: unscramble
84 : : public :: unswitch_cent
85 : : public :: unswitchreal
86 : : public :: mpiswitch
87 : : public :: mpiswitch_cent
88 : :
89 : : public :: mpifft_fg2dbox
90 : : public :: mpifft_fg2dbox_dpc
91 : : public :: mpifft_dbox2fg
92 : : public :: mpifft_dbox2fg_dpc
93 : : public :: mpifft_dbox2fr
94 : : public :: mpifft_dbox2fr_dpc
95 : : public :: mpifft_fr2dbox
96 : : public :: mpifft_fr2dbox_dpc
97 : : public :: mpifft_collect_datar ! Collect a real-space MPI-FFT distributed array on each proc.
98 : :
99 : : public :: indfftrisc
100 : : public :: addrho
101 : : public :: multpot
102 : :
103 : : ! 0 for double precision version (default), 1 for mixed precision FFTs
104 : : integer, public, save, protected :: fftcore_mixprec = 0
105 : : public :: fftcore_set_mixprec
106 : : ! *************************************************************************
107 : :
108 : : !----------------------------------------------------------------------
109 : : ! Private variables
110 : :
111 : : #define FFTALGA_SIZE 5
112 : : character(len=*),private,parameter :: fftalga2name(1:FFTALGA_SIZE)= &
113 : : (/"Goedecker ", &
114 : : "Vendor FFT ", &
115 : : "FFTW3 ", &
116 : : "Goedecker2002 ", &
117 : : "DFTI " /)
118 : :
119 : : #define FFTALGB_SIZE 1
120 : : character(len=*),private,parameter :: fftalgb2name(0:FFTALGB_SIZE)= &
121 : : (/"C2C",&
122 : : "R2C"/)
123 : :
124 : : #define FFTALGC_SIZE 2
125 : : character(len=*),private,parameter :: fftalgc2name(0:FFTALGC_SIZE)= &
126 : : (/"No pad ",&
127 : : "zero-pad ",&
128 : : "zero-pad+cache "/)
129 : :
130 : : contains
131 : : !!***
132 : :
133 : : !----------------------------------------------------------------------
134 : :
135 : : !!****f* m_fftcore/fftcore_set_mixprec
136 : : !! NAME
137 : : !! fftalg_set_precision
138 : : !!
139 : : !! FUNCTION
140 : : !! Set the precision to be used in the FFT routines: 0 for standard double precision,
141 : : !! 1 for mixed precision (dp input, sp for intermediate arrays passed to FFT libs)
142 : : !! Return old value.
143 : : !!
144 : : !! INPUTS
145 : : !!
146 : : !! SOURCE
147 : :
148 : 58101 : integer function fftcore_set_mixprec(wp) result(old_wp)
149 : :
150 : : !Arguments ------------------------------------
151 : : !scalars
152 : : integer,intent(in) :: wp
153 : : ! *************************************************************************
154 : :
155 : 58101 : old_wp = fftcore_mixprec
156 : 58101 : fftcore_mixprec = abs(wp)
157 : :
158 : 58053 : select case (fftcore_mixprec)
159 : : case (0)
160 [ - + ]: 58053 : if (old_wp /= fftcore_mixprec) call wrtout(std_out, " fftcore_mixprec 0 --> Using double-precision FFT", newlines=1)
161 : : case (1)
162 [ + + ]: 48 : if (old_wp /= fftcore_mixprec) call wrtout(std_out, " fftcore_mixprec 1 --> Using mixed precision FFT", newlines=1)
163 : : case default
164 [ + + - ]: 58101 : ABI_ERROR(sjoin("Wrong value for input wp:", itoa(fftcore_mixprec)))
165 : : end select
166 : :
167 : 58101 : end function fftcore_set_mixprec
168 : : !!***
169 : :
170 : : !----------------------------------------------------------------------
171 : :
172 : : !!****f* m_fftcore/fftalg_isavailable
173 : : !! NAME
174 : : !! fftalg_isavailable
175 : : !!
176 : : !! FUNCTION
177 : : !! Returns TRUE if the FFT library specified by fftalg (ngfft(7)) is available
178 : : !!
179 : : !! INPUTS
180 : : !! fftalg=Input variable.
181 : : !!
182 : : !! SOURCE
183 : :
184 : 98 : logical pure function fftalg_isavailable(fftalg) result(ans)
185 : :
186 : : !Arguments ------------------------------------
187 : : integer,intent(in) :: fftalg
188 : :
189 : : !Local variables-------------------------------
190 : : integer :: fftalga,fftalgb,fftalgc
191 : : ! *************************************************************************
192 : :
193 : 98 : ans = .TRUE.
194 : 98 : fftalga = fftalg/100
195 : 98 : fftalgb = mod(fftalg,100)/10
196 : 98 : fftalgc = mod(fftalg,10)
197 : :
198 : : ! Optional FFT libraries.
199 : : #ifndef HAVE_FFTW3
200 [ + + ]: 98 : if (fftalga == FFT_FFTW3) ans = .FALSE.
201 : : #endif
202 : :
203 : : #ifndef HAVE_DFTI
204 : : if (fftalga == FFT_DFTI) ans = .FALSE.
205 : : #endif
206 : :
207 : 98 : end function fftalg_isavailable
208 : : !!***
209 : :
210 : : !----------------------------------------------------------------------
211 : :
212 : : !!****f* m_fftcore/fftalg_has_mpi
213 : : !! NAME
214 : : !! fftalg_has_mpi
215 : : !!
216 : : !! FUNCTION
217 : : !! True if the FFT library specified by fftalg is available.
218 : : !!
219 : : !! INPUTS
220 : : !! fftalg=Input variable.
221 : : !!
222 : : !! SOURCE
223 : :
224 : 268 : pure function fftalg_has_mpi(fftalg) result(ans)
225 : :
226 : : !Arguments ------------------------------------
227 : : !scalars
228 : : integer,intent(in) :: fftalg
229 : : logical :: ans
230 : :
231 : : !Local variables-------------------------------
232 : : !scalars
233 : : integer :: fftalga,fftalgb,fftalgc
234 : : ! *************************************************************************
235 : :
236 : 268 : ans = .False.
237 : 268 : fftalga = fftalg/100; fftalgb = mod(fftalg,100)/10; fftalgc = mod(fftalg,10)
238 : :
239 [ + - ]: 268 : if (fftalga == FFT_FFTW3) ans = .True.
240 : : !if (fftalga == FFT_DFTI) ans = .True.
241 [ + - ]: 268 : if (fftalga == FFT_SG2002) ans = .True.
242 : :
243 : 268 : end function fftalg_has_mpi
244 : : !!***
245 : :
246 : : !----------------------------------------------------------------------
247 : :
248 : : !!****f* m_fftcore/fftalg_for_npfft
249 : : !! NAME
250 : : !! fftalg_for_npfft
251 : : !!
252 : : !! FUNCTION
253 : : !! Returns the default value of fftalg given the number of MPI nodes
254 : : !! to be used in the FFTs.
255 : : !!
256 : : !! INPUTS
257 : : !! nproc_fft=Number of processors used for MPI FFT
258 : : !! nthreads =Number of openMP threads
259 : : !!
260 : : !! OUTPUT
261 : : !! fftalg=Integer used to select the FFT library.
262 : : !!
263 : : !! SOURCE
264 : :
265 : 22507 : pure function fftalg_for_npfft(nproc_fft, forbid_threads) result(fftalg)
266 : :
267 : : !Arguments ------------------------------------
268 : : !scalars
269 : : integer,intent(in) :: nproc_fft
270 : : logical,intent(in),optional :: forbid_threads
271 : : integer :: fftalg
272 : : ! *************************************************************************
273 : :
274 : : ! Default for the sequential case.
275 : : fftalg = 112
276 : :
277 : : ! Use Goedecker2002 if fftalg does not support MPI or threads (e.g 112)
278 : 22507 : if (nproc_fft > 1) fftalg = 401
279 : : if (present(forbid_threads)) then
280 : : if (forbid_threads) fftalg = 401
281 : : endif
282 : :
283 : : #ifdef HAVE_FFTW3
284 : : fftalg = 312
285 : : #elif defined HAVE_DFTI
286 : 22507 : fftalg = 512
287 [ + + ]: 22507 : if (nproc_fft > 1) fftalg = 401 ! MPI-FFT with DFTI is not implemented yet
288 : : #endif
289 : :
290 : : !if (nproc_fft > 1) fftalg = 401 ! This is to revert to the old behavior.
291 : :
292 : 22507 : end function fftalg_for_npfft
293 : : !!***
294 : :
295 : : !----------------------------------------------------------------------
296 : :
297 : : !!****f* m_fftcore/fftalg_info
298 : : !! NAME
299 : : !! fftalg_info
300 : : !!
301 : : !! FUNCTION
302 : : !! Returns info on the FFT library specified by fftalg (ngfft(7))
303 : : !!
304 : : !! INPUTS
305 : : !! fftalg=Input variable.
306 : : !!
307 : : !! OUTPUT
308 : : !! library=String with the name of FFT library
309 : : !! cplex_mode= String defining whether the FFT library supports real<-->complex transforms.
310 : : !! padding_mode=Padding mode.
311 : : !!
312 : : !! SOURCE
313 : :
314 : 423 : subroutine fftalg_info(fftalg,library,cplex_mode,padding_mode)
315 : :
316 : : !Arguments ------------------------------------
317 : : !scalars
318 : : integer,intent(in) :: fftalg
319 : : character(len=*),intent(out) :: library,cplex_mode,padding_mode
320 : :
321 : : !Local variables-------------------------------
322 : : !scalars
323 : : integer :: fftalga,fftalgb,fftalgc
324 : : ! *************************************************************************
325 : :
326 [ + - + - : 423 : library = "Unknown"; cplex_mode = "Unknown"; padding_mode = "Unknown"
+ - + - +
- ]
327 : :
328 : 423 : fftalga=fftalg/100
329 [ + - + - : 423 : if (fftalga>0 .and. fftalga<=FFTALGA_SIZE) library = fftalga2name(fftalga)
+ - ]
330 : :
331 : 423 : fftalgb=mod(fftalg,100)/10
332 [ + - + - : 423 : if (fftalgb>=0 .and. fftalgb<=FFTALGB_SIZE) cplex_mode = fftalgb2name(fftalgb)
+ - ]
333 : :
334 : 423 : fftalgc=mod(fftalg,10)
335 [ + - + - : 423 : if (fftalgc>=0 .and. fftalgc<=FFTALGC_SIZE) padding_mode = fftalgc2name(fftalgc)
+ - ]
336 : :
337 [ + - ]: 423 : end subroutine fftalg_info
338 : : !!***
339 : :
340 : : !----------------------------------------------------------------------
341 : :
342 : : !!****f* m_fftcore/get_cache_kb
343 : : !! NAME
344 : : !! get_cache_kb
345 : : !!
346 : : !! FUNCTION
347 : : !! Returns the cache size in KB to be used for cache blocking algorithms in the FFT routines.
348 : : !! The value is derived from the values of the CPP options defined in config.h
349 : : !!
350 : : !! TODO
351 : : !! Use C to get the real cache size.
352 : : !! See http://stackoverflow.com/questions/12594208/c-program-to-determine-levels-size-of-cache
353 : : !!
354 : : !! SOURCE
355 : :
356 : 8585 : pure function get_cache_kb()
357 : :
358 : : !Local variables-------------------------------
359 : : !scalars
360 : : integer :: get_cache_kb
361 : : ! *************************************************************************
362 : :
363 : : ! Default value
364 : 154954 : get_cache_kb = 16
365 : : !get_cache_kb = 32
366 : : !get_cache_kb = 256
367 : :
368 : 8585 : end function get_cache_kb
369 : : !!***
370 : :
371 : : !----------------------------------------------------------------------
372 : :
373 : : !!****f* m_fftcore/ngfft_seq
374 : : !! NAME
375 : : !! ngfft_seq
376 : : !!
377 : : !! FUNCTION
378 : : !! Helper function used to initialize ngfft(18) from the FFT divisions
379 : : !! in the case of sequential execution.
380 : : !!
381 : : !! INPUTS
382 : : !! n123(3)=FFT divisions.
383 : : !!
384 : : !! OUTPUT
385 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft.
386 : : !!
387 : : !! SOURCE
388 : :
389 : 146369 : pure subroutine ngfft_seq(ngfft, n123)
390 : :
391 : : !Arguments ------------------------------------
392 : : integer,intent(in) :: n123(3)
393 : : integer,intent(out) :: ngfft(18)
394 : :
395 : : !Local variables-------------------------------
396 : : integer :: fftalg
397 : : ! *************************************************************************
398 : :
399 : : ! Default for sequential case.
400 : : fftalg = 112
401 : : #ifdef HAVE_FFTW3
402 : : fftalg = 312
403 : : #elif defined HAVE_DFTI
404 : 146369 : fftalg = 512
405 : : #endif
406 : :
407 [ + + ]: 585476 : ngfft(1:3) = n123
408 : 146369 : ngfft(4) = 2*(ngfft(1)/2)+1
409 : 146369 : ngfft(5) = 2*(ngfft(2)/2)+1
410 : 146369 : ngfft(6) = ngfft(3)
411 : 146369 : ngfft(7)= fftalg ! fftalg
412 : 146369 : ngfft(8)= get_cache_kb() ! cache_kb
413 : 146369 : ngfft(9)= 0 ! paral_fft_
414 : 146369 : ngfft(10)=1 ! nproc_fft
415 : 146369 : ngfft(11)=0 ! me_fft
416 : 146369 : ngfft(12)=0 ! n2proc
417 : 146369 : ngfft(13)=0 ! n3proc
418 [ + + ]: 878214 : ngfft(14:18)=0 ! not used
419 : :
420 : 146369 : end subroutine ngfft_seq
421 : : !!***
422 : :
423 : : !----------------------------------------------------------------------
424 : :
425 : : !!****f* m_fftcore/print_ngfft
426 : : !! NAME
427 : : !! print_ngfft
428 : : !!
429 : : !! FUNCTION
430 : : !! Print the content of ngfft(18) in explicative format.
431 : : !!
432 : : !! INPUTS
433 : : !! units: Unit numbers
434 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft.
435 : : !! [unit]=unit number for output (defaults to std_out).
436 : : !! [prtvol]=verbosity level (defaults to 0).
437 : : !!
438 : : !! OUTPUT
439 : : !! Only writing
440 : : !!
441 : : !! SOURCE
442 : :
443 : 20672 : subroutine print_ngfft(units, ngfft, header, prtvol)
444 : :
445 : : !Arguments ------------------------------------
446 : : !scalars
447 : : integer,intent(in) :: units(:)
448 : : integer,intent(in),optional :: prtvol
449 : : character(len=*),intent(in),optional :: header
450 : : !arrays
451 : : integer,intent(in) :: ngfft(18)
452 : :
453 : : !Local variables-------------------------------
454 : : integer :: my_prtvol
455 : : character(len=500) :: msg
456 : : ! *************************************************************************
457 : :
458 [ + + ]: 10336 : my_prtvol=0; if (PRESENT(prtvol)) my_prtvol=prtvol
459 : :
460 : 10336 : msg = ch10//' ==== FFT mesh description (ngfft) ==== '
461 [ + - + - : 10336 : if (PRESENT(header)) msg=ch10//' ==== '//TRIM(ADJUSTL(header))//' ==== '
+ - ]
462 : 10336 : call wrtout(units, msg)
463 : :
464 : : write(msg,'(2(a,3i5,a),a,i5,2a,i5)')&
465 : 10336 : '- FFT mesh divisions ........................ ',ngfft(1),ngfft(2),ngfft(3),ch10,&
466 : 10336 : '- Augmented FFT divisions ................... ',ngfft(4),ngfft(5),ngfft(6),ch10,&
467 : 10336 : '- FFT algorithm ............................. ',ngfft(7),ch10,&
468 : 20672 : '- FFT cache size ............................ ',ngfft(8)
469 : 10336 : call wrtout(units, msg)
470 : :
471 [ - + ]: 10336 : if (my_prtvol > 0) then
472 : : write(msg,'(6(a,i5,a),a,4i5)')&
473 : 0 : '- FFT parallelization level ................. ',ngfft(9),ch10,&
474 : 0 : '- Number of processors in my FFT group ...... ',ngfft(10),ch10,&
475 : 0 : '- Index of me in my FFT group ............... ',ngfft(11),ch10,&
476 : 0 : '- No of xy planes in R space treated by me .. ',ngfft(12),ch10,&
477 : 0 : '- No of xy planes in G space treated by me .. ',ngfft(13),ch10,&
478 : 0 : '- MPI communicator for FFT .................. ',ngfft(14),ch10,&
479 : 0 : '- Value of ngfft(15:18) ..................... ',ngfft(15:18)
480 : 0 : call wrtout(units, msg)
481 : : end if
482 : :
483 [ - + ]: 10336 : end subroutine print_ngfft
484 : : !!***
485 : :
486 : : !!****f* m_fftcore/bound
487 : : !! NAME
488 : : !! bound
489 : : !!
490 : : !! FUNCTION
491 : : !! For given kpt, ngfft, and gmet,
492 : : !! Find distance**2 to boundary point of fft box nearest to kpt
493 : : !! Find distance**2 to boundary point of fft box farthest to kpt
494 : : !!
495 : : !! INPUTS
496 : : !! kpt(3)=real input k vector (reduced coordinates)
497 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
498 : : !! gmet(3,3)=reciprocal space metric (currently in Bohr**-2)
499 : : !!
500 : : !! OUTPUT
501 : : !! dsqmax=maximum distance**2 from k to boundary in Bohr**-2.
502 : : !! dsqmin=minimum distance**2 from k to boundary in Bohr**-2.
503 : : !! gbound(3)=coords of G on boundary (corresponding to gsqmin)
504 : : !! plane=which plane min occurs in (1,2, or 3 for G1,etc).
505 : : !!
506 : : !! NOTES
507 : : !! Potential trouble: this routine was written assuming kpt lies inside
508 : : !! first Brillouin zone. No measure is taken to fold input kpt back
509 : : !! into first zone. Given arbitrary kpt, this will cause trouble.
510 : : !!
511 : : !! SOURCE
512 : :
513 : 263234 : subroutine bound(dsqmax,dsqmin,gbound,gmet,kpt,ngfft,plane)
514 : :
515 : : !Arguments ------------------------------------
516 : : !scalars
517 : : integer,intent(out) :: plane
518 : : real(dp),intent(out) :: dsqmax,dsqmin
519 : : !arrays
520 : : integer,intent(in) :: ngfft(18)
521 : : integer,intent(out) :: gbound(3)
522 : : real(dp),intent(in) :: gmet(3,3),kpt(3)
523 : :
524 : : !Local variables-------------------------------
525 : : !scalars
526 : : integer :: i1,i1min,i2,i2min,i3,i3min
527 : : real(dp) :: dsm,dsp
528 : : character(len=500) :: msg
529 : : ! *************************************************************************
530 : :
531 : : !Set plane to impossible value
532 : 263234 : plane=0
533 : :
534 : : !look at +/- g1 planes:
535 : 263234 : dsqmax=zero
536 : 263234 : dsqmin=dsq(ngfft(1)/2,-ngfft(2)/2,-ngfft(3)/2,gmet,kpt)+0.01_dp
537 [ + + ]: 3673920 : do i2=-ngfft(2)/2,ngfft(2)/2
538 [ + + ]: 83941602 : do i3=-ngfft(3)/2,ngfft(3)/2
539 : 80267682 : dsp = dsq(ngfft(1)/2, i2, i3,gmet,kpt)
540 : 80267682 : dsm = dsq( - ngfft(1)/2, i2, i3,gmet,kpt)
541 [ + + ]: 80267682 : if (dsp>dsqmax) dsqmax = dsp
542 [ + + ]: 80267682 : if (dsm>dsqmax) dsqmax = dsm
543 [ + + ]: 80267682 : if (dsp<dsqmin) then
544 : 4423863 : dsqmin = dsp
545 : 4423863 : i1min = ngfft(1)/2
546 : 4423863 : i2min = i2
547 : 4423863 : i3min = i3
548 : 4423863 : plane=1
549 : : end if
550 [ + + ]: 83678368 : if (dsm<dsqmin) then
551 : 1191830 : dsqmin = dsm
552 : 1191830 : i1min = - ngfft(1)/2
553 : 1191830 : i2min = i2
554 : 1191830 : i3min = i3
555 : 1191830 : plane=1
556 : : end if
557 : : end do
558 : : end do
559 : : !
560 : : !+/- g2 planes:
561 [ + + ]: 3924208 : do i1=-ngfft(1)/2,ngfft(1)/2
562 [ + + ]: 88515534 : do i3=-ngfft(3)/2,ngfft(3)/2
563 : 84591326 : dsp = dsq(i1,ngfft(2)/2,i3,gmet,kpt)
564 : 84591326 : dsm = dsq(i1,-ngfft(2)/2,i3,gmet,kpt)
565 [ - + ]: 84591326 : if (dsp>dsqmax) dsqmax = dsp
566 [ - + ]: 84591326 : if (dsm>dsqmax) dsqmax = dsm
567 [ + + ]: 84591326 : if (dsp<dsqmin) then
568 : 383366 : dsqmin = dsp
569 : 383366 : i1min = i1
570 : 383366 : i2min = ngfft(2)/2
571 : 383366 : i3min = i3
572 : 383366 : plane=2
573 : : end if
574 [ + + ]: 88252300 : if (dsm<dsqmin) then
575 : 212241 : dsqmin = dsm
576 : 212241 : i1min = i1
577 : 212241 : i2min = - ngfft(2)/2
578 : 212241 : i3min = i3
579 : 212241 : plane=2
580 : : end if
581 : : end do
582 : : end do
583 : : !
584 : : !+/- g3 planes:
585 [ + + ]: 3924208 : do i1=-ngfft(1)/2,ngfft(1)/2
586 [ + + ]: 86414166 : do i2=-ngfft(2)/2,ngfft(2)/2
587 : 82489958 : dsp = dsq(i1,i2,ngfft(3)/2,gmet,kpt)
588 : 82489958 : dsm = dsq(i1,i2,-ngfft(3)/2,gmet,kpt)
589 [ - + ]: 82489958 : if (dsp>dsqmax) dsqmax = dsp
590 [ - + ]: 82489958 : if (dsm>dsqmax) dsqmax = dsm
591 [ + + ]: 82489958 : if (dsp<dsqmin) then
592 : 302716 : dsqmin = dsp
593 : 302716 : i1min = i1
594 : 302716 : i2min = i2
595 : 302716 : i3min = ngfft(3)/2
596 : 302716 : plane=3
597 : : end if
598 [ + + ]: 86150932 : if (dsm<dsqmin) then
599 : 144233 : dsqmin = dsm
600 : 144233 : i1min = i1
601 : 144233 : i2min = i2
602 : 144233 : i3min = - ngfft(3)/2
603 : 144233 : plane=3
604 : : end if
605 : : end do
606 : : end do
607 : :
608 [ - + ]: 263234 : if (plane==0) then
609 : : ! Trouble: missed boundary somehow
610 : : write(msg, '(a,a,a,3f9.4,a,3(i0,1x),a,a,a,a,a)' )&
611 : 0 : 'Trouble finding boundary of G sphere for',ch10,&
612 : 0 : 'kpt=',kpt(:),' and ng=',ngfft(1:3),ch10,&
613 : 0 : 'Action : check that kpt lies',&
614 : 0 : 'reasonably within first Brillouin zone; ',ch10,&
615 : 0 : 'else code bug, contact ABINIT group.'
616 : 0 : ABI_BUG(msg)
617 : : end if
618 : :
619 : 263234 : gbound(1)=i1min
620 : 263234 : gbound(2)=i2min
621 : 263234 : gbound(3)=i3min
622 : :
623 : : contains
624 : :
625 : 494961166 : function dsq(i1,i2,i3,gmet,kpt)
626 : :
627 : : integer :: i1,i2,i3
628 : : real(dp) :: dsq
629 : : real(dp) :: kpt(3),gmet(3,3)
630 : :
631 : : dsq=gmet(1,1)*(kpt(1)+dble(i1))**2&
632 : : & +gmet(2,2)*(kpt(2)+dble(i2))**2&
633 : : & +gmet(3,3)*(kpt(3)+dble(i3))**2&
634 : : & +2._dp*(gmet(1,2)*(kpt(1)+dble(i1))*(kpt(2)+dble(i2))&
635 : : & +gmet(2,3)*(kpt(2)+dble(i2))*(kpt(3)+dble(i3))&
636 : 494961166 : & +gmet(3,1)*(kpt(3)+dble(i3))*(kpt(1)+dble(i1)))
637 : 494961166 : end function dsq
638 : :
639 : : end subroutine bound
640 : : !!***
641 : :
642 : : !!****f* m_fftcore/getng
643 : : !! NAME
644 : : !! getng
645 : : !!
646 : : !! FUNCTION
647 : : !! From ecut and metric tensor in reciprocal space, computes recommended ngfft(1:3)
648 : : !! Also computes the recommended value of nfft and mgfft
649 : : !! Pay attention that the FFT grid must be compatible with the symmetry operations (see irrzg.f).
650 : : !!
651 : : !! INPUTS
652 : : !! boxcutmin=minimum value of boxcut admitted (boxcut is the ratio
653 : : !! between the radius of the sphere contained in the FFT box, and the
654 : : !! radius of the planewave sphere): usually 2.0.
655 : : !! chksymtnons= if==3, will impose the FFT grid to be invariant under the spatial symmetries.
656 : : !! ecut=energy cutoff in Hartrees
657 : : !! gmet(3,3)=reciprocal space metric (bohr**-2).
658 : : !! kpt(3)=input k vector in terms of reciprocal lattice primitive translations
659 : : !! me_fft=index of the processor in the FFT set (use 0 if sequential)
660 : : !! nproc_fft=number of processors in the FFT set (use 1 if sequential)
661 : : !! nsym=number of symmetry elements in group
662 : : !! paral_fft=0 if no FFT parallelisation; 1 if FFT parallelisation
663 : : !! symrel(3,3,nsym)=symmetry matrices in real space (integers)
664 : : !! tnons(3,nsym)=nonsymmorphic translations associated to symrel
665 : : !!
666 : : !! OUTPUT
667 : : !! mgfft= max(ngfft(1),ngfft(2),ngfft(3))
668 : : !! nfft=number of points in the FFT box=ngfft(1)*ngfft(2)*ngfft(3)/nproc_fft
669 : : !!
670 : : !! SIDE EFFECTS
671 : : !! Input/Output
672 : : !! ngfft(1:18)=integer array with FFT box dimensions and other information on FFTs.
673 : : !! On input ngfft(1:3) contains optional trial values. If ngfft(1:3)/minbox is greater than value
674 : : !! calculated to avoid wrap-around error and ngfft obeys constraint placed by the FFT routine that is used
675 : : !! then ngfft(1:3) is left unchanged. Otherwise set to value computed in now.
676 : : !!
677 : : !! Note that there is the possibility of an undetected error if we
678 : : !! are dealing with a cubic unit cell and ngfft(1), ngfft(2) and ngfft(3)
679 : : !! are different. In the future we should handle this case.
680 : : !!
681 : : !! ngfft(4),ngfft(5),ngfft(6)= modified values to avoid cache trashing,
682 : : !! presently: ngfft(4)=ngfft(1)+1 if ngfft(1) is even;
683 : : !! ngfft(5)=ngfft(2)+1 if ngfft(2) is even.
684 : : !! in the other cases, ngfft(4:6)=ngfft(1:3).
685 : : !! Other choices may better, but this is left for the future.
686 : : !! ngfft(7)=choice for FFT algorithm, see the input variable fftalg
687 : : !! ngfft(8)=size of the cache, in bytes (not used here presently).!!
688 : : !! other ngfft slots are used for parallelism see ~abinit/doc/variables/vargs.htm#ngfft
689 : : !! [ngfftc(1:18)]= -optional- value of ngfft for the "coarse" grid
690 : : !! [unit] = -optional- output unit number (DEFAULT std_out)
691 : : !! [gpu_option] = GPU implementation to use, i.e. cuda, openMP, ... (0=not using GPU)
692 : : !!
693 : : !! SOURCE
694 : :
695 : 8494 : subroutine getng(boxcutmin, chksymtnons, ecut, gmet, kpt, me_fft, mgfft, nfft, ngfft, &
696 : 8494 : nproc_fft, nsym, paral_fft, symrel, tnons, &
697 : : ngfftc, unit, gpu_option) ! optional
698 : :
699 : : use defs_fftdata, only : mg
700 : :
701 : : !Arguments ------------------------------------
702 : : !scalars
703 : : integer,intent(in) :: chksymtnons,me_fft,nproc_fft,nsym,paral_fft
704 : : integer,intent(out) :: mgfft,nfft
705 : : integer,optional,intent(in) :: unit,gpu_option
706 : : real(dp),intent(in) :: boxcutmin,ecut
707 : : !arrays
708 : : integer,intent(in) :: symrel(3,3,nsym)
709 : : integer,intent(in),optional :: ngfftc(3)
710 : : integer,intent(inout) :: ngfft(18)
711 : : real(dp),intent(in) :: gmet(3,3),kpt(3)
712 : : real(dp),intent(in) :: tnons(3,nsym)
713 : :
714 : : !Local variables-------------------------------
715 : : !scalars
716 : : integer,save :: first=1,msrch(3),previous_paral_mode=0
717 : : integer :: element,ifactor,ii,index,ipower,isrch,isrch1,isrch2,isrch3,isym,jj,mu,paral_fft_
718 : : integer :: plane,testok,tobechecked,ount,fftalga,nn,ngdiv,valpow
719 : : real(dp),parameter :: minbox=0.75_dp
720 : : real(dp) :: dsqmax,dsqmin,ecutmx,prodcurrent,prodtrial,tnscaled,xx,yy
721 : : logical :: testdiv
722 : : character(len=500) :: msg
723 : : integer,parameter :: largest_ngfft=mg ! Goedecker FFT: any powers of 2, 3, and 5 - must be coherent with defs_fftdata.F90
724 : : integer,parameter :: maxpow2 =16 ! int(log(largest_ngfft+half)/log(two))
725 : : integer,parameter :: maxpow3 =6 ! int(log(largest_ngfft+half)/log(three))
726 : : integer,parameter :: maxpow5 =6 ! int(log(largest_ngfft+half)/log(five))
727 : : !#if defined HAVE_FFTW3 || defined HAVE_DFTI
728 : : !#ifdef _GMATTEO_WHISH_LIST
729 : : #if 0
730 : : integer,parameter :: maxpow7 =5 ! FFTW3 and DFTI support powers of 7 and 11
731 : : integer,parameter :: maxpow11=4
732 : : #else
733 : : integer,parameter :: maxpow7 =0
734 : : integer,parameter :: maxpow11=0
735 : : #endif
736 : : integer,parameter :: mmsrch=(maxpow2+1)*(maxpow3+1)*(maxpow5+1)*(maxpow7+1)*(maxpow11+1)
737 : : integer,parameter :: nfactor=10, mpower=5
738 : : !Arrays
739 : : integer,save :: iperm(mmsrch),srch(mmsrch,3)
740 : : integer(i8b) :: li_srch(mmsrch)
741 : : integer :: divisor(3,3),gbound(3),imax(3),imin(3),ngcurrent(3)
742 : : integer :: ngmax(3),ngsav(3),ngtrial(3)
743 : : integer :: npower(3,mpower)
744 : : integer,parameter :: factor(10) = (/1,2,3,4,5,6,8,9,10,12/)
745 : : integer,parameter :: power(5) = (/2,3,5,7,11/)
746 : : ! *************************************************************************
747 : :
748 [ + + ]: 8494 : ount = std_out; if (present(unit)) ount = unit
749 : :
750 : : !write(std_out,*)' m_fftcore/getng : enter'
751 : :
752 : 8494 : fftalga = ngfft(7)/100
753 : :
754 : : !If not yet done, compute recommended (boxcut>=2) fft grid dimensions
755 : : !In case we switch for paral to sequential mode, recompute srch.
756 : : !This is the case e.g. when computing ngfftdiel in sequential mode
757 : : !after an initial computation of ngfft in parallel
758 : :
759 [ - + ]: 8494 : paral_fft_=paral_fft;if (nproc_fft==0) paral_fft_=0
760 : :
761 [ + + + + ]: 8494 : if(first==1.or.paral_fft_ /= previous_paral_mode) then
762 : 2148 : first=0; previous_paral_mode=paral_fft_
763 : 2148 : srch(:,:)=0
764 : :
765 : : ! Factors of 2
766 : 2148 : srch(1,1)=1
767 [ + + ]: 36516 : do ii=1,maxpow2
768 : 36516 : srch(ii+1,1)=srch(ii,1)*2
769 : : end do
770 : :
771 : : ! Factors of 3
772 : 15036 : index=maxpow2+1
773 : : if(maxpow3>0)then
774 [ + + ]: 15036 : do ii=1,max(1,maxpow3)
775 [ + - + + : 453228 : srch(1+ii*index:(ii+1)*index,1)=3*srch(1+(ii-1)*index:ii*index,1)
+ + ]
776 : : end do
777 : : end if
778 : :
779 : : ! Factors of 5
780 : 15036 : index=(maxpow3+1)*index
781 : : if(maxpow5>0)then
782 [ + + ]: 15036 : do ii=1,max(1,maxpow5)
783 : 12888 : li_srch = 0
784 [ + + ]: 1546560 : li_srch(1+ii*index:(ii+1)*index)=5_i8b*srch(1+(ii-1)*index:ii*index,1)
785 [ + + + + ]: 10748592 : where (li_srch > huge(maxpow3)) li_srch = huge(maxpow3)
786 [ + + ]: 1548708 : srch(1+ii*index:(ii+1)*index,1)=li_srch(1+ii*index:(ii+1)*index)
787 : : end do
788 : : end if
789 : :
790 : : ! Factors of 7
791 : : index=(maxpow5+1)*index
792 : : if(maxpow7>0)then
793 : : do ii=1,max(1,maxpow7)
794 : : srch(1+ii*index:(ii+1)*index,1)=7*srch(1+(ii-1)*index:ii*index,1)
795 : : end do
796 : : end if
797 : :
798 : : ! Factors of 11
799 : 2148 : index=(maxpow7+1)*index
800 : : if(maxpow11>0)then
801 : : do ii=1,max(1,maxpow11)
802 : : srch(1+ii*index:(ii+1)*index,1)=11*srch(1+(ii-1)*index:ii*index,1)
803 : : end do
804 : : end if
805 : :
806 : 2148 : call sort_int(mmsrch,srch(:,1),iperm)
807 : :
808 [ + - ]: 575664 : do ii=1,mmsrch
809 [ + + ]: 575664 : if(srch(ii,1)>largest_ngfft)exit
810 : : end do
811 : 2148 : msrch(1)=ii-1
812 : :
813 : : ! In case of FFT parallelism, one need ngfft(2) and ngfft(3) to be multiple of nproc_fft
814 [ + + ]: 2148 : if(paral_fft_==1)then
815 : 1472 : msrch(2)=0
816 [ + + ]: 394496 : do ii=1,msrch(1)
817 [ + + ]: 394496 : if(modulo(srch(ii,1),nproc_fft)==0) then
818 : 390360 : msrch(2)=msrch(2)+1
819 : 390360 : srch(msrch(2),2)=srch(ii,1)
820 : : end if
821 : : end do
822 : : !write(msg,'(a,i0,a,i0,2a,i0)')&
823 : : ! 'The second and third dimension of the FFT grid: ',ngfft(2),", ",ngfft(3),ch10,&
824 : : ! 'were imposed to be multiple of the number of processors for the FFT: ', nproc_fft
825 : : !if (ount /= dev_null) ABI_COMMENT(msg)
826 : : else
827 : 676 : msrch(2)=msrch(1)
828 [ + + ]: 563784 : srch(:,2)=srch(:,1)
829 : : end if
830 : :
831 : : ! The second and third search list have the same constraint
832 : 2148 : msrch(3)=msrch(2)
833 [ + + ]: 1797778 : srch(:,3)=srch(:,2)
834 : :
835 : : ! The set of allowed ngfft values has been found
836 : : end if ! first==1
837 : :
838 : : !=============================================================================================
839 : : !
840 : : ! Determination of sufficient values of ngfft with ngfft(2) and ngfft(3) taken inside
841 : : ! sets of values that take into account the constraint on nproc_fft
842 : :
843 : : !Save input values of ngfft
844 [ + + ]: 33976 : ngsav(1:3) = ngfft(1:3)
845 : :
846 : : !As an initial guess for ngfft, use the provided coarse mesh grid
847 [ + + ]: 8494 : if (PRESENT(ngfftc)) then
848 [ + + ]: 4888 : ngfft(1:3)=ngfftc(1:3)
849 : 1222 : call wrtout(ount,' Using supplied coarse mesh as initial guess.')
850 : : else
851 [ + + ]: 29088 : ngfft(1:3)=2
852 : : end if
853 : :
854 : : !Enlarge the initial guess until the set of ngfft entirely comprises the sphere
855 : 243694 : do
856 : :
857 : 243694 : call bound(dsqmax,dsqmin,gbound,gmet,kpt,ngfft,plane)
858 : :
859 : : ! Exit the infinite do-loop if the sphere is inside the FFT box
860 [ + + ]: 243694 : if (dsqmin>=(half*boxcutmin**2*ecut/pi**2)) exit
861 : :
862 : : ! Fix nearest boundary
863 [ - + ]: 1932620 : do ii=1,msrch(plane)-1
864 [ + + ]: 1924126 : if (srch(ii,plane)>=ngfft(plane)) then
865 : : ! redefine ngfft(plane) to next higher choice
866 : 235200 : ngfft(plane)=srch(ii+1,plane)
867 : 235200 : exit ! Exit the loop over ii
868 : : end if
869 : :
870 [ - + ]: 1688926 : if (ii==msrch(plane)-1)then
871 : : ! Here, we are in trouble
872 : : write(msg, '(a,i12,5a)' ) &
873 : 0 : 'ngfft is bigger than allowed value =',ngfft(plane),'.',ch10,&
874 : 0 : 'This indicates that desired ngfft is larger than getng',ch10,&
875 : 0 : 'can handle. The code has to be changed and compiled.'
876 : 0 : ABI_BUG(msg)
877 : : end if
878 : : end do
879 : :
880 : : end do ! End of the infinite do-loop : will either "exit", or abort
881 : :
882 : : !ecutmx=maximum ecut consistent with chosen ngfft
883 : 8494 : ecutmx=0.5_dp*pi**2*dsqmin
884 : :
885 : : !Print results
886 : : write(msg, '(a,1p,e14.6,a,3i8,a,a,e14.6)' ) &
887 : 8494 : ' For input ecut=',ecut,' best grid ngfft=',ngfft(1:3),ch10,&
888 : 16988 : ' max ecut=',ecutmx
889 : 8494 : call wrtout(ount,msg)
890 : :
891 : : ! The FFT grid is compatible with the symmetries if for each
892 : : ! symmetry isym, each ii and each jj, the quantity
893 : : ! (ngfft(jj)*symrel(jj,ii,isym))/ngfft(ii) is an integer.
894 : : ! This relation is immediately verified for diagonal elements, since
895 : : ! symrel is an integer. It is also verified if symrel(ii,jj,isym) is zero.
896 : : ! Moreover, to cope with the non-symmorphic translation vectors, at least the
897 : : ! origin must be sent to a point of the FFT grid. Hence, the non-symmorphic
898 : : ! translations tnons(i) multiplied by ngfft(i) must be an integer.
899 : : ! The latter condition is however imposed only when chksymtnons=3.
900 : : ! Indeed, ABINIT will be able to reimpose the symmetry at the level of the density and potential.
901 : : ! This might be a problem for GW calculations, though ...
902 : :
903 : : !Compute the biggest (positive) common divisor of each off-diagonal element of the symmetry matrices
904 : 8494 : divisor(:,:)=0; tobechecked=0
905 : :
906 [ + + ]: 33976 : do ii=1,3
907 [ + + ]: 110422 : do jj=1,3
908 [ + + ]: 76446 : if(jj==ii)cycle
909 [ + + ]: 1255416 : do isym=1,nsym
910 [ + + + + ]: 1178970 : if(symrel(jj,ii,isym)==0 .or. divisor(jj,ii)==1 )cycle
911 : 30845 : tobechecked=1
912 : 30845 : element=abs(symrel(jj,ii,isym))
913 [ + + - + : 30845 : testdiv= ( divisor(jj,ii)==0 .or. divisor(jj,ii)==element .or. element==1)
- - ]
914 : 76446 : if(testdiv)then
915 : 30845 : divisor(jj,ii)=element
916 : : else
917 : : ! Must evaluate common divisor between non-trivial numbers
918 : : do
919 [ # # ]: 0 : if(divisor(jj,ii)<element)element=element-divisor(jj,ii)
920 [ # # ]: 0 : if(divisor(jj,ii)>element)divisor(jj,ii)=divisor(jj,ii)-element
921 [ - - ]: 1148125 : if(divisor(jj,ii)==element)exit
922 : : end do
923 : : end if
924 : : end do
925 : : end do
926 : : end do
927 : :
928 : : !Check whether there is a problem: the grid must be invariant
929 : : !with respect to point symmetry operations, and spatial symmetry operations if chksymtnons==3
930 : 8494 : testok=1
931 [ + + ]: 8494 : if(tobechecked==1)then
932 [ + + ]: 24608 : do ii=1,3
933 [ + + ]: 73824 : do jj=1,3
934 : 55368 : xx=divisor(jj,ii)*ngfft(jj)
935 : 55368 : yy=xx/ngfft(ii)
936 [ + + ]: 73824 : if(abs(yy-nint(yy))>tol8)testok=0
937 : : end do
938 [ + + ]: 24608 : if(chksymtnons==3)then
939 [ + + ]: 5877 : do isym=1,nsym
940 : 5682 : tnscaled=tnons(ii,isym)*ngfft(ii)
941 [ + + ]: 5877 : if(abs(tnscaled-nint(tnscaled))>tol8)testok=0
942 : : enddo
943 : : endif
944 : : end do
945 : : end if
946 : :
947 : : !DEBUG
948 : : !write(std_out,*)' m_fftcore/getng : chksymtnons,testok,nproc_fft=',chksymtnons,testok,nproc_fft
949 : : !ENDDEBUG
950 : :
951 : : !If there is a problem
952 [ + + ]: 6152 : if(testok==0)then
953 : : ! Find the powers of 2, 3, 5 (possibly 7 and 11) that are needed or will be enough in ngfft,
954 : : ! taking both the constraint on nproc_nfft and the constraint on tnons.
955 : : ! First decompose nproc_fft, providing divisors of ngfft(2:3)
956 : 88 : nn=nproc_fft
957 : 88 : npower(:,:)=0
958 [ + + ]: 528 : do ipower=1,mpower
959 : 440 : valpow=power(ipower)
960 [ - + ]: 528 : do while (mod(nn,valpow)==0)
961 : 0 : nn=nn/valpow
962 : 0 : npower(3,ipower)=npower(3,ipower)+1
963 : : end do
964 : : enddo
965 [ - + ]: 88 : if(nn/=1)then
966 : 0 : ABI_ERROR(sjoin("nproc_fft: ", itoa(nproc_fft), "is not a multiple of 2, 3, 5, 7 or 11"))
967 : : endif
968 [ + + ]: 528 : npower(2,:)=npower(3,:)
969 : :
970 : : ! Then examine tnons
971 [ + + ]: 88 : if(chksymtnons==3)then
972 [ + + ]: 56 : do ii=1,3
973 [ + - ]: 90 : do ifactor=1,nfactor
974 : 90 : testok=1
975 [ + + ]: 2322 : do isym=1,nsym
976 : 2232 : tnscaled=factor(ifactor)*tnons(ii,isym)
977 [ + + ]: 2322 : if(abs(tnscaled-nint(tnscaled))>tol8)testok=0
978 : : enddo
979 [ + + ]: 90 : if(testok==1)exit
980 : : enddo
981 [ + - ]: 56 : if(testok==1)then
982 [ + + ]: 42 : if(ifactor/=1)then
983 [ + + ]: 36 : if(ifactor==2)npower(ii,1)=max(npower(ii,1),1) ! At least one power of 2
984 [ + + ]: 36 : if(ifactor==3)npower(ii,2)=max(npower(ii,2),1) ! At least one power of 3
985 [ - + ]: 36 : if(ifactor==4)npower(ii,1)=max(npower(ii,1),2) ! At least two powers of 2
986 [ - + ]: 36 : if(ifactor==5)npower(ii,3)=max(npower(ii,3),1) ! At least one power of 5
987 [ + + ]: 36 : if(ifactor==6)then
988 : 2 : npower(ii,1)=max(npower(ii,1),1) ! At least one power of 2
989 : 2 : npower(ii,2)=max(npower(ii,2),1) ! At least one power of 3
990 : : endif
991 [ - + ]: 36 : if(ifactor==7)npower(ii,1)=max(npower(ii,1),3) ! At least three powers of 2
992 [ - + ]: 36 : if(ifactor==8)npower(ii,2)=max(npower(ii,2),2) ! At least two powers of 3
993 [ - + ]: 36 : if(ifactor==9)then
994 : 0 : npower(ii,1)=max(npower(ii,1),1) ! At least one power of 2
995 : 0 : npower(ii,3)=max(npower(ii,3),1) ! At least one power of 5
996 : : end if
997 [ - + ]: 36 : if(ifactor==10)then
998 : 0 : npower(ii,1)=max(npower(ii,1),2) ! At least two powers of 2
999 : 0 : npower(ii,2)=max(npower(ii,2),1) ! At least one power of 3
1000 : : end if
1001 : : endif
1002 : : else
1003 : : write(msg, '(5a,i12,2a,9i12,2a,3f10.7,2a)' ) &
1004 : 0 : 'Chksymtnons=1 . Found potentially symmetry-breaking value of tnons, ', ch10,&
1005 : 0 : ' which is neither a rational fraction in 1/8th nor in 1/12th (1/9th and 1/10th are tolerated also) :', ch10,&
1006 : 0 : ' for the symmetry number ',isym,ch10,&
1007 : 0 : ' symrel is ',symrel(1:3,1:3,isym),ch10,&
1008 : 0 : ' tnons is ',tnons(1:3,isym),ch10,&
1009 : 0 : 'This problem should have been caught earlier.'
1010 : 0 : ABI_BUG(msg)
1011 : : endif
1012 : : enddo ! ii
1013 : : endif ! chksymtnons
1014 : :
1015 : : ! Get minimal search indices, simply those of the current ngfft
1016 [ + + ]: 352 : do ii=1,3
1017 [ + + ]: 70840 : do isrch=1,msrch(ii)
1018 : 70488 : index=srch(isrch,ii)
1019 [ + + ]: 70752 : if(index==ngfft(ii))imin(ii)=isrch
1020 : : end do
1021 : : end do
1022 : :
1023 : : ! Get maximal search indices : the ngtrial values must be identical (to fulfill the constraint induced by the
1024 : : ! off diagonal elements of symrel), but also must contain sufficient powers of basic primes (2, 3, 5, 7, 11),
1025 : : ! and be bigger than all current ngfft components. This should guarantee that such a triplet fulfills all constraints.
1026 : : ! Determine the divisor of allowed ngmax
1027 : : ngdiv=1
1028 [ + + ]: 528 : do ipower=1,mpower
1029 [ + + ]: 1848 : ngdiv=ngdiv*power(ipower)**(maxval(npower(:,ipower)))
1030 : : enddo
1031 [ + + ]: 352 : ngmax(1)=ngdiv*(maxval(ngfft(1:3)-1)/ngdiv+1)
1032 [ + + ]: 352 : ngmax(1:3)=ngmax(1)
1033 [ + + ]: 352 : do ii=1,3
1034 [ + - ]: 3847 : do isrch=1,msrch(ii)
1035 : 3759 : index=srch(isrch,ii)
1036 [ + + + + ]: 3759 : if(mod(index,ngdiv)==0 .and. index>=ngmax(ii))then
1037 : 264 : imax(ii)=isrch
1038 : 264 : ngmax(ii)=index
1039 : 264 : exit
1040 : : endif
1041 : : end do
1042 : : end do
1043 : : ! This gives a tentative symmetric triplet
1044 : : ngcurrent(1:3)=ngmax(1:3)
1045 : 88 : prodcurrent=ngmax(1)*ngmax(2)*ngmax(3)+1.0d-3
1046 : : ! However, it is perhaps possible to do better, by asymmetric triplets, still giving lower prodcurrent !
1047 : 88 : ngmax(1)=min(int(prodcurrent/(ngfft(2)*ngfft(3))),srch(msrch(1),1))
1048 : 88 : ngmax(2)=min(int(prodcurrent/(ngfft(1)*ngfft(3))),srch(msrch(2),2))
1049 : 88 : ngmax(3)=min(int(prodcurrent/(ngfft(1)*ngfft(2))),srch(msrch(3),3))
1050 [ + + ]: 352 : do ii=1,3
1051 [ + + ]: 70840 : do isrch=1,msrch(ii)
1052 : 70488 : index=srch(isrch,ii)
1053 : : ! One cannot suppose that ngmax belongs to the allowed list,
1054 : : ! so must use <= instead of == , to determine largest index
1055 [ + + ]: 70752 : if(index<=ngmax(ii))imax(ii)=isrch
1056 : : end do
1057 : : end do
1058 : :
1059 : : !write(std_out,*)' ngmin(1:3)=',srch(imin(1),1),srch(imin(2),2),srch(imin(3),3)
1060 : : !write(std_out,*)' ngmax(1:3)=',ngmax(1:3)
1061 : :
1062 : 88 : ngcurrent(1:3)=ngmax(1:3)
1063 : 88 : prodcurrent=ngmax(1)*ngmax(2)*ngmax(3)+1.0d-3
1064 : :
1065 : : ! Now, start brute force search
1066 [ + + ]: 968 : do isrch1=imin(1),imax(1)
1067 : 880 : ngtrial(1)=srch(isrch1,1)
1068 [ + + ]: 16009 : do isrch2=imin(2),imax(2)
1069 : 15041 : ngtrial(2)=srch(isrch2,2)
1070 [ + + ]: 130881 : do isrch3=imin(3),imax(3)
1071 : 128354 : ngtrial(3)=srch(isrch3,3)
1072 : 128354 : prodtrial=real(ngtrial(1))*real(ngtrial(2))*real(ngtrial(3))+1.0d-3
1073 [ + + ]: 128354 : if(prodtrial>prodcurrent-1.0d-4)exit
1074 : : ! The trial product is lower or equal to the current product,
1075 : : ! so now, checks whether the symmetry constraints are OK
1076 : : testok=1
1077 [ + + ]: 459840 : do ii=1,3
1078 [ + + ]: 1379520 : do jj=1,3
1079 : 1034640 : xx=divisor(jj,ii)*ngtrial(jj)
1080 : 1034640 : yy=xx/ngtrial(ii)
1081 [ + + ]: 1379520 : if(abs(yy-nint(yy))>tol8)testok=0
1082 : : end do
1083 [ + + ]: 459840 : if(chksymtnons==3)then
1084 [ + + ]: 12828 : do isym=1,nsym
1085 : 12240 : tnscaled=tnons(ii,isym)*ngtrial(ii)
1086 [ + + ]: 12828 : if(abs(tnscaled-nint(tnscaled))>tol8)testok=0
1087 : : enddo
1088 : : endif
1089 : : end do
1090 : : ! DEBUG
1091 : : ! write(ount,'(a,3i6,a,i3,a,es16.6)' )' getng : current trial triplet',ngtrial(1:3),&
1092 : : ! & ' testok=',testok,' prodtrial=',prodtrial
1093 : : ! ENDDEBUG
1094 [ + + ]: 114960 : if(testok==0)cycle
1095 : : ! When one arrives here, the symmetry constraints are fulfilled, so update current values
1096 : : ! Then continues the search, in hope of a better value.
1097 : 92 : ngcurrent(1:3)=ngtrial(1:3)
1098 : 130001 : prodcurrent=prodtrial
1099 : : end do
1100 : : end do
1101 : : end do
1102 : :
1103 [ + + ]: 352 : ngfft(1:3)=ngcurrent(1:3)
1104 : 88 : call bound(dsqmax,dsqmin,gbound,gmet,kpt,ngfft,plane)
1105 : : ! ecutmx=maximum ecut consistent with chosen ngfft
1106 : 88 : ecutmx=0.5_dp*pi**2*dsqmin
1107 : : ! Give results
1108 : : write(msg, '(a,3i8,a,a,e14.6)' ) &
1109 : 88 : ' However, must be changed due to symmetry =>',ngfft(1:3),ch10,&
1110 : 176 : ' with max ecut=',ecutmx
1111 : 88 : call wrtout(ount,msg)
1112 : :
1113 [ - + ]: 88 : if (prodcurrent>huge(ii)) then
1114 : : write(msg, '(5a)' )&
1115 : 0 : 'The best FFT grid will lead to indices larger',ch10,&
1116 : 0 : 'than the largest representable integer on this machine.',ch10,&
1117 : 0 : 'Action: try to deal with smaller problems. Also contact ABINIT group.'
1118 : 0 : ABI_ERROR(msg)
1119 : : end if
1120 : :
1121 : : end if ! testok==0
1122 : :
1123 : : !Possibly use the input values of ngfft
1124 : : if (int( dble(ngsav(1)) / minbox ) >= ngfft(1) .and.&
1125 [ + + + - : 8494 : int( dble(ngsav(2)) / minbox ) >= ngfft(2) .and.&
+ - ]
1126 : : int( dble(ngsav(3)) / minbox ) >= ngfft(3) ) then
1127 : :
1128 : : ! Must check whether the values are in the allowed list
1129 : 588 : testok=0
1130 [ + + ]: 2352 : do mu=1,3
1131 [ + - ]: 23479 : do ii=1,msrch(mu)
1132 [ + + ]: 22891 : if(srch(ii,mu)==ngsav(mu))then
1133 : 1764 : testok=testok+1
1134 : 1764 : exit
1135 : : end if
1136 : : end do
1137 : : end do
1138 [ + - ]: 588 : if(testok==3)then
1139 : 588 : write(msg,'(a,3(a,i1,a,i3),a)') ' input values of',&
1140 [ + - + + ]: 2940 : (' ngfft(',mu,') =',ngsav(mu),mu=1,3),' are alright and will be used'
1141 : 588 : call wrtout(ount,msg)
1142 [ + + ]: 2352 : do mu = 1,3
1143 : 2352 : ngfft(mu) = ngsav(mu)
1144 : : end do
1145 : : end if
1146 : :
1147 : : end if
1148 : :
1149 : : !mgfft needs to be set to the maximum of ngfft(1),ngfft(2),ngfft(3)
1150 [ + + ]: 33976 : mgfft = maxval(ngfft(1:3))
1151 : :
1152 [ + + ]: 8494 : if (paral_fft_==1) then
1153 : : ! For the time being, one need ngfft(2) and ngfft(3) to be multiple of nproc_fft
1154 [ - + ]: 6755 : if(modulo(ngfft(2),nproc_fft)/=0)then
1155 : : write(msg,'(4a,i5,a,i5)')&
1156 : 0 : 'The second dimension of the FFT grid, ngfft(2), should be ',&
1157 : 0 : 'a multiple of the number of processors for the FFT, nproc_fft.',ch10,&
1158 : 0 : 'However, ngfft(2)=',ngfft(2),' and nproc_fft=',nproc_fft
1159 : 0 : ABI_BUG(msg)
1160 : : end if
1161 [ - + ]: 6755 : if(modulo(ngfft(3),nproc_fft)/=0)then
1162 : : write(msg,'(4a,i5,a,i5)')&
1163 : 0 : 'The third dimension of the FFT grid, ngfft(3), should be ',&
1164 : 0 : 'a multiple of the number of processors for the FFT, nproc_fft.',ch10,&
1165 : 0 : 'However, ngfft(3)=',ngfft(3),' and nproc_fft=',nproc_fft
1166 : 0 : ABI_BUG(msg)
1167 : : end if
1168 : :
1169 [ - + ]: 1739 : else if (paral_fft_/=0) then
1170 : 0 : write(msg,'(a,i0)')'paral_fft_ should be 0 or 1, but its value is ',paral_fft_
1171 : 0 : ABI_BUG(msg)
1172 : : end if
1173 : :
1174 : : ! Compute effective number of FFT points (for this MPI node if parall FFT)
1175 [ + + ]: 33976 : nfft=product(ngfft(1:3))/max(1,nproc_fft)
1176 : :
1177 : : !Set up fft array dimensions ngfft(4,5,6) to avoid cache conflicts
1178 : 8494 : ngfft(4)=2*(ngfft(1)/2)+1
1179 : 8494 : ngfft(5)=2*(ngfft(2)/2)+1
1180 : 8494 : ngfft(6)=ngfft(3)
1181 : 8494 : if (any(fftalga == [FFT_FFTW3, FFT_DFTI])) then
1182 : : ! FFTW3 supports leading dimensions but at the price of a larger number of FFTs
1183 : : ! to be executed along z when the zero-padded version is used.
1184 : : ! One should check whether the augmentation is beneficial for FFTW3.
1185 : : ngfft(4)=2*(ngfft(1)/2)+1
1186 : : ngfft(5)=2*(ngfft(2)/2)+1
1187 : : !ngfft(4)=ngfft(1)
1188 : : !ngfft(5)=ngfft(2)
1189 : : ngfft(6)=ngfft(3)
1190 : : end if
1191 : :
1192 [ + + ]: 8494 : if (present(gpu_option)) then
1193 [ - + ]: 8486 : if (gpu_option/=ABI_GPU_DISABLED) then
1194 : 0 : ngfft(4)=ngfft(1)
1195 : 0 : ngfft(5)=ngfft(2)
1196 : 0 : ngfft(6)=ngfft(3)
1197 : : end if
1198 : : end if
1199 : :
1200 [ + + ]: 50964 : ngfft(14:18)=0 ! ngfft(14) to be filled outside of getng
1201 : :
1202 [ + + ]: 8494 : if (paral_fft_==0) then
1203 : 1739 : ngfft(9)=0 ! paral_fft_
1204 : 1739 : ngfft(10)=1 ! nproc_fft
1205 : 1739 : ngfft(11)=0 ! me_fft
1206 : 1739 : ngfft(12)=0 ! n2proc
1207 : 1739 : ngfft(13)=0 ! n3proc
1208 : : else
1209 : 6755 : ngfft(9)=1 ! paral_fft_
1210 : 6755 : ngfft(10)=nproc_fft
1211 : 6755 : ngfft(11)=me_fft
1212 : 6755 : ngfft(12)=ngfft(2)/nproc_fft
1213 : 6755 : ngfft(13)=ngfft(3)/nproc_fft
1214 : : end if
1215 : :
1216 [ + + ]: 16988 : call print_ngfft([ount], ngfft, header="FFT mesh")
1217 : :
1218 : 8494 : end subroutine getng
1219 : : !!***
1220 : :
1221 : : !!****f* m_fftcore/sphereboundary
1222 : : !! NAME
1223 : : !! sphereboundary
1224 : : !!
1225 : : !! FUNCTION
1226 : : !! Finds the boundary of the basis sphere of G vectors (at a given
1227 : : !! k point) for use in improved zero padding of ffts in 3 dimensions.
1228 : : !! Provides data to be used by subroutine fourwf, in the form of
1229 : : !! an array gbound(2*mgfft+8,2).
1230 : : !!
1231 : : !! The first component (for use when mod(fftalg,10)==2))
1232 : : !! provides integer values g1min,g1max,g2min,g2max
1233 : : !! and then for g2 in the
1234 : : !! sequence g2=0,1,2,...,g2max,g2min,g2min+1,...,-1, provides g1min, g1max.
1235 : : !!
1236 : : !! The second component (for use when mod(fftalg,10)==1))
1237 : : !! provides integer values g1min,g1max,g3min,g3max,
1238 : : !! where g3min and g3max have been corrected in case of time-reversal
1239 : : !! and then for g3 in the sequence
1240 : : !! g3=0,1,2,...,g3max,g3min,g3min+1,...,-1, provides g2min, g2max.
1241 : : !! (also corrected in case of time-reversal)
1242 : : !!
1243 : : !! These values are stored in the above order in array gbound.
1244 : : !! Debug mode, if fftalg is between 000 and 099
1245 : : !!
1246 : : !! INPUTS
1247 : : !! istwf_k=option parameter that describes the storage of wfs
1248 : : !! kg_k(3,npw)=integer coordinates of G vectors in basis sphere
1249 : : !! mgfft=maximum size of 1D FFTs (only for dimensioning purposes)
1250 : : !! npw=number of G vectors in basis at this k point
1251 : : !!
1252 : : !! OUTPUT
1253 : : !! gbound(2*mgfft+8,2)=defined above
1254 : : !!
1255 : : !! SOURCE
1256 : :
1257 : 3270680 : subroutine sphereboundary(gbound, istwf_k, kg_k, mgfft, npw)
1258 : :
1259 : : !Arguments ------------------------------------
1260 : : !scalars
1261 : : integer,intent(in) :: istwf_k,mgfft,npw
1262 : : !arrays
1263 : : integer,intent(in) :: kg_k(3,npw)
1264 : : integer,intent(out) :: gbound(2*mgfft+8,2)
1265 : :
1266 : : !Local variables-------------------------------
1267 : : !scalars
1268 : : integer :: dim_a,dim_b,fftalgc,g_a,gmax_a,gmax_b,gmax_b1,gmax_b2,gmin_a,gmin_b
1269 : : integer :: gmin_b1,gmin_b2,igb,ii,iloop,ipw,testm,testp,kgk
1270 : : character(len=500) :: msg
1271 : : !arrays
1272 : : integer :: gmax(2),gmin(2)
1273 : : ! *************************************************************************
1274 : : !
1275 : : !DEBUG
1276 : : !write(std_out,*)' sphereboundary : enter'
1277 : : !write(std_out, '(a)' )' sphereboundary : list of plane waves coordinates for k point '
1278 : : !write(std_out, '(a)' )' ipw kg_k(1:3,ipw) '
1279 : : !do ipw=1,npw
1280 : : !write(std_out, '(i10,a,3i6)' )ipw,' ',kg_k(1:3,ipw)
1281 : : !end do
1282 : : !gbound=-999
1283 : : !ENDDEBUG
1284 : :
1285 : : !Determine cube boundaries
1286 [ + + ]: 615611345 : gbound(1,1)=minval(kg_k(1,:))
1287 [ + + ]: 615611345 : gbound(2,1)=maxval(kg_k(1,:))
1288 [ + + ]: 9812040 : gbound(1:2,2)=gbound(1:2,1)
1289 : :
1290 : : !Treat differently the fftalgc cases
1291 [ + + ]: 9812040 : do ii=1,2
1292 : :
1293 : 6541360 : fftalgc=3-ii
1294 : :
1295 [ + + ]: 6541360 : if(fftalgc/=2)then
1296 : : dim_a=3
1297 : : dim_b=2
1298 : : else
1299 : 3270680 : dim_a=2
1300 : 3270680 : dim_b=1
1301 : : end if
1302 : :
1303 : : ! Relevant boundaries
1304 [ + + ]: 1231222690 : gbound(3,ii)=minval(kg_k(dim_a,:))
1305 [ + + ]: 1231222690 : gbound(4,ii)=maxval(kg_k(dim_a,:))
1306 : 6541360 : gmin_a=gbound(3,ii)
1307 : 6541360 : gmax_a=gbound(4,ii)
1308 : :
1309 : : ! Must complete the sphere for fftalgc==1 and special storage modes.
1310 : : ! Explanation : sg_fftpad is not able to take into account
1311 : : ! the time-reversal symmetry, so that the boundaries will not be delimited
1312 : : ! by the kg_k set, but by their symmetric also.
1313 [ + + + + ]: 6541360 : if(istwf_k>=2 .and. fftalgc==1)then
1314 [ + + + + ]: 31405 : if( istwf_k==2 .or. istwf_k==3 .or. istwf_k==6 .or. istwf_k==7 )then
1315 : 25450 : gbound(4,2)=max(gmax_a,-gmin_a)
1316 : 25450 : gbound(3,2)=-gbound(4,2)
1317 [ + + + - ]: 5955 : else if( istwf_k==4 .or. istwf_k==5 .or. istwf_k==8 .or. istwf_k==9 )then
1318 : 5955 : gbound(4,2)=max(gmax_a,-gmin_a-1)
1319 : 5955 : gbound(3,2)=-gbound(4,2)-1
1320 : : end if
1321 : 31405 : gmax_a=gbound(4,2) ; gmin_a=gbound(3,2)
1322 : : end if
1323 : :
1324 : 6541360 : igb=5
1325 : :
1326 : : ! Consider first every g_a in range 0 ... gmax_a, then gmin_a ... -1
1327 : 6541360 : gmin(1)=0 ; gmax(1)=gmax_a
1328 : 6541360 : gmin(2)=gmin_a ; gmax(2)=-1
1329 : :
1330 [ + + ]: 22894760 : do iloop=1,2
1331 : :
1332 [ + + ]: 19624080 : if( gmin(iloop) <= gmax(iloop) )then
1333 : :
1334 [ + + ]: 58487573 : do g_a=gmin(iloop),gmax(iloop)
1335 : :
1336 [ + + + + ]: 45446495 : if(istwf_k==1 .or. fftalgc/=1)then
1337 : : ! Select the minimal and maximal values, in the selected plane
1338 : 45082537 : gmin_b=mgfft+1 ! Initialized with a value larger than all possible ones
1339 : 45082537 : gmax_b=-mgfft-1 ! Initialized with a value smaller than all possible ones
1340 [ + + ]: 13747146737 : do ipw=1,npw
1341 [ + + ]: 13747146737 : if(kg_k(dim_a,ipw)==g_a)then
1342 : 1202103678 : kgk=kg_k(dim_b,ipw)
1343 : 1202103678 : if(kgk<=gmin_b)gmin_b=kgk
1344 : 1202103678 : if(kgk>=gmax_b)gmax_b=kgk
1345 : : end if
1346 : : end do
1347 : :
1348 [ + - ]: 363958 : else if(istwf_k>=2 .and. fftalgc==1)then
1349 : :
1350 : : ! Here, must take into account time-reversal symmetry explicitly
1351 : :
1352 : : ! Determine the boundaries for the plane g_a
1353 : 363958 : testp=0
1354 [ + - ]: 363958 : if(g_a<=gmax_a)then
1355 : : ! Select the minimal and maximal values, in the selected plane
1356 : 363958 : gmin_b1=mgfft+1 ! Initialized with a value larger than all possible ones
1357 : 363958 : gmax_b1=-mgfft-1 ! Initialized with a value smaller than all possible ones
1358 [ + + ]: 644259775 : do ipw=1,npw
1359 [ + + ]: 644259775 : if(kg_k(dim_a,ipw)==g_a)then
1360 : 22577652 : kgk=kg_k(dim_b,ipw)
1361 : 22577652 : if(kgk<=gmin_b1)gmin_b1=kgk
1362 : 22577652 : if(kgk>=gmax_b1)gmax_b1=kgk
1363 : : end if
1364 : : end do
1365 : :
1366 : :
1367 : : testp=1
1368 : : end if
1369 : :
1370 : : ! Determine the boundaries for the plane -g_a or -g_a-1
1371 : 363958 : testm=0
1372 [ + + + + ]: 363958 : if( istwf_k==2 .or. istwf_k==3 .or. istwf_k==6 .or. istwf_k==7 )then
1373 : :
1374 [ + - ]: 294854 : if(-g_a>=gmin_a)then
1375 : : ! Select the minimal and maximal values, in the selected plane
1376 : : ! Warning : there is an inversion of search (might be confusing)
1377 : 294854 : gmax_b2=mgfft+1 ! Initialized with a value larger than all possible ones
1378 : 294854 : gmin_b2=-mgfft-1 ! Initialized with a value smaller than all possible ones
1379 [ + + ]: 506293405 : do ipw=1,npw
1380 [ + + ]: 506293405 : if(kg_k(dim_a,ipw)==-g_a)then
1381 : 18477543 : kgk=kg_k(dim_b,ipw)
1382 : 18477543 : if(kgk<=gmax_b2)gmax_b2=kgk
1383 : 18477543 : if(kgk>=gmin_b2)gmin_b2=kgk
1384 : : end if
1385 : : end do
1386 : : testm=1
1387 : : end if
1388 : :
1389 [ + + + - ]: 69104 : else if( istwf_k==4 .or. istwf_k==5 .or. istwf_k==8 .or. istwf_k==9 )then
1390 : :
1391 [ + - ]: 69104 : if(-g_a-1>=gmin_a)then
1392 : : ! Select the minimal and maximal values, in the selected plane
1393 : : ! Warning : there is an inversion of search (might be confusing)
1394 : 69104 : gmax_b2=mgfft+1 ! Initialized with a value larger than all possible ones
1395 : 69104 : gmin_b2=-mgfft-1 ! Initialized with a value smaller than all possible ones
1396 [ + + ]: 137966370 : do ipw=1,npw
1397 [ + + ]: 137966370 : if(kg_k(dim_a,ipw)==-g_a-1)then
1398 : 4100109 : kgk=kg_k(dim_b,ipw)
1399 : 4100109 : if(kgk<=gmax_b2)gmax_b2=kgk
1400 : 4100109 : if(kgk>=gmin_b2)gmin_b2=kgk
1401 : : end if
1402 : : end do
1403 : : testm=1
1404 : : end if
1405 : :
1406 : : end if
1407 : :
1408 : : ! Must invert the boundaries, to use them for plane g_a
1409 : : if(testm==1)then
1410 : : ! This is needed to avoid border effect
1411 : : ! if the search did not lead to any element
1412 : 363958 : gmin_b2=max(gmin_b2,-mgfft) ; gmax_b2=min(gmax_b2,mgfft)
1413 [ + + ]: 363958 : if(istwf_k<=5)then
1414 : 284139 : gmax_b2=-gmax_b2 ; gmin_b2=-gmin_b2
1415 : : else
1416 : 79819 : gmax_b2=-gmax_b2-1 ; gmin_b2=-gmin_b2-1
1417 : : end if
1418 : : end if
1419 : :
1420 [ + - ]: 363958 : if( testp==1 .and. testm==1)then
1421 : 363958 : gmin_b=min(gmin_b1,gmin_b2) ; gmax_b=max(gmax_b1,gmax_b2)
1422 [ # # ]: 0 : else if( testp==1 )then
1423 : : gmin_b=gmin_b1 ; gmax_b=gmax_b1
1424 [ # # ]: 0 : else if( testm==1 )then
1425 : 0 : gmin_b=gmin_b2 ; gmax_b=gmax_b2
1426 : : end if
1427 : :
1428 : : end if ! Endif take into account time-reversal symmetry
1429 : :
1430 [ - + ]: 45446495 : if (igb+1>2*mgfft+4) then
1431 : : write(msg, '(2a, 4(a,3(i0,1x),a))' )&
1432 : 0 : "About to overwrite gbound array (FFT mesh too small) ",ch10, &
1433 : 0 : " iloop, igb, mgb = ",iloop,igb,2*mgfft+4, ch10, &
1434 : 0 : " istwfk, mgfft, npw = ",istwf_k, mgfft, npw, ch10, &
1435 [ # # ]: 0 : " minval(kg_k) = ",minval(kg_k, dim=2), ch10, &
1436 [ # # ]: 0 : " maxval(kg_k) = ",maxval(kg_k, dim=2), ch10
1437 : 0 : ABI_BUG(msg)
1438 : : end if
1439 : :
1440 : 45446495 : gbound(igb,ii)=gmin_b
1441 : 45446495 : gbound(igb+1,ii)=gmax_b
1442 : :
1443 [ + + + + : 45446495 : if( iloop==1 .and. istwf_k>=2 .and. istwf_k<=5 .and. fftalgc==2 .and. g_a==0)then
+ + + + ]
1444 : : ! If k_y=0 , for fftalgc==2, the g_a==0 plane must be completed
1445 [ + + ]: 24144 : if(istwf_k==2 .or. istwf_k==4)then
1446 : 19065 : gbound(igb+1,ii)=max(gmax_b,-gmin_b)
1447 : 19065 : gbound(igb,ii)=-gbound(igb+1,ii)
1448 [ + - ]: 5079 : else if(istwf_k==3 .or. istwf_k==5)then
1449 : 5079 : gbound(igb+1,ii)=max(gmax_b,-gmin_b-1)
1450 : 5079 : gbound(igb,ii)=-gbound(igb+1,ii)-1
1451 : : end if
1452 : :
1453 : : end if
1454 : :
1455 : 58487573 : igb=igb+2
1456 : :
1457 : : end do ! g_a
1458 : : end if
1459 : : end do ! iloop
1460 : : end do ! ii (fftalgc)
1461 : :
1462 : : !DEBUG
1463 : : !write(std_out,'(a)')' sphereoundary : list of plane waves coordinates for 1st k point '
1464 : : !write(std_out,'(a)')' ipw kg_k(1:3,ipw) '
1465 : : !do ipw=1,npw
1466 : : !write(std_out, '(i10,a,3i6)' )ipw,' ',kg_k(1:3,ipw)
1467 : : !end do
1468 : : !write(std_out, '(a)' )' sphereboundary : list of boundaries '
1469 : : !do igb=1,2*mgfft+8
1470 : : !write(std_out, '(i10,a,2i6)' )igb,' ',gbound(igb,1),gbound(igb,2)
1471 : : !end do
1472 : : !write(std_out,*)' sphereboundary : exit '
1473 : : !ENDDEBUG
1474 : :
1475 : 3270680 : end subroutine sphereboundary
1476 : : !!***
1477 : :
1478 : : !----------------------------------------------------------------------
1479 : :
1480 : : !!****f* m_fftcore/sphere
1481 : : !! NAME
1482 : : !! sphere
1483 : : !!
1484 : : !! FUNCTION
1485 : : !! Array cg is defined in sphere with npw g-vectors.
1486 : : !! Insert cg inside FFT box of n1*n2*n3 points to define array cfft for
1487 : : !! rest of cfft is filled with 0 s.
1488 : : !!
1489 : : !! iflag=1 ==> insert cg into cfft.
1490 : : !! iflag=2 ==> insert cg into cfft, where the second and third dimension
1491 : : !! have been switched (needed for new 2002 SGoedecker FFT)
1492 : : !! iflag=-1==> extract cg from cfft.
1493 : : !! iflag=-2==> extract cg from cfft, where the second and third dimension
1494 : : !! have been switched (needed for new 2002 SGoedecker FFT)
1495 : : !! WARNING: iflag=-2 cannot use symmetry operations.
1496 : : !!
1497 : : !! There is also the possibility to apply a symmetry operation,
1498 : : !! as well as to make a shift in reciprocal space, or to multiply
1499 : : !! by a constant factor, in the case iflag=-1.
1500 : : !! Multiplication by a constant factor is also possible in the case iflag=-2.
1501 : : !!
1502 : : !! INPUTS
1503 : : !! cg(2,npw*ndat)= contains values for npw G vectors in basis sphere
1504 : : !! ndat=number of wavefunctions
1505 : : !! npw=number of G vectors in basis at this k point
1506 : : !! cfft(2,n4,n5,n6*ndat) = array in FFT box
1507 : : !! n1,n2,n3=physical dimension of the box (cfft)
1508 : : !! n4,n5,n6=memory dimension of cfft
1509 : : !! kg_k(3,npw)=integer coordinates of G vectors in basis sphere
1510 : : !! istwf_k=option parameter that describes the storage of wfs
1511 : : !! iflag=option parameter. Possible values: -1, -2, 1, 2
1512 : : !! me_g0=1 if this node has G=0.
1513 : : !! shiftg(3)=The shift in reciprocal space.
1514 : : !! symrec(3,3)=symmetry operation in reciprocal space to be applied (symrec)
1515 : : !! xnorm=Normalization factor.
1516 : : !!
1517 : : !! SIDE EFFECTS
1518 : : !! Input/Output
1519 : : !! iflag=1 and 2, insert cg(input) into cfft(output)
1520 : : !! iflag=-1 and -2, extract cg(output) from cfft(input)
1521 : : !!
1522 : : !! NOTES
1523 : : !! cg and cfft are assumed to be of type COMPLEX, although this routine treats
1524 : : !! them as real of twice the length to avoid nonstandard complex*16.
1525 : : !! If istwf_k differs from 1, then special storage modes must be taken
1526 : : !! into account, for symmetric wavefunctions coming from k=(0 0 0) or other special k points.
1527 : : !!
1528 : : !! TODO
1529 : : !! 1) Order arguments
1530 : : !! 2) Split the two cases to avoid breaking intent: from and to sphere (merge with cg_box2gpsh and cg_gsph2box?)
1531 : : !! 3) If symmetries are used with or without shiftg, it might happen that the FFT mesh
1532 : : !! is not large enough to accommodate the rotated G, in this case one should return ierr /= 0
1533 : : !!
1534 : : !! SOURCE
1535 : :
1536 : 2553219 : subroutine sphere(cg, ndat, npw, cfft, n1, n2, n3, n4, n5, n6, kg_k, istwf_k, iflag, me_g0, shiftg, symrec, xnorm)
1537 : :
1538 : : !Arguments ------------------------------------
1539 : : !scalars
1540 : : integer,intent(in) :: iflag,istwf_k,n1,n2,n3,n4,n5,n6,ndat,npw,me_g0
1541 : : real(dp),intent(in) :: xnorm
1542 : : !arrays
1543 : : integer,intent(in) :: kg_k(3,npw),shiftg(3),symrec(3,3)
1544 : : real(dp),intent(inout) :: cfft(2,n4,n5,n6*ndat),cg(2,npw*ndat)
1545 : :
1546 : : !Local variables-------------------------------
1547 : : !scalars
1548 : : integer :: i1,i1inv,i2,i2inv,i3,i3inv,id1,id2,id3,idat,ipw
1549 : : integer :: j1,j2,j3,l1,l2,l3,npwmin,use_symmetry,i3dat,i3invdat,i2invdat,ipwdat,i2dat
1550 : : !character(len=500) :: msg
1551 : : !arrays
1552 : : integer :: identity(3,3)
1553 : 5106438 : integer :: i1inver(n1),i2inver(n2),i3inver(n3)
1554 : : ! *************************************************************************
1555 : :
1556 : : DBG_ENTER("COLL")
1557 : :
1558 : : ! In the case of special k-points, invariant under time-reversal,
1559 : : ! but not Gamma, initialize the inverse coordinates.
1560 : : ! Remember that:
1561 : : !
1562 : : ! u_k(G) = u_{k+G0}(G-G0); u_{-k}(G) = u_k(G)^* and therefore:
1563 : : ! u_{G0/2}(G) = u_{G0/2}(-G-G0)^*.
1564 : :
1565 [ + + ]: 2553219 : if (istwf_k>=2) then
1566 [ + + ]: 55506 : if(istwf_k==2 .or. istwf_k==4 .or. istwf_k==6 .or. istwf_k==8)then
1567 : 21767 : i1inver(1)=1
1568 [ + + ]: 210519 : do i1=2,n1
1569 : 210519 : i1inver(i1)=n1+2-i1
1570 : : end do
1571 : : else
1572 [ + + ]: 362262 : do i1=1,n1
1573 : 362262 : i1inver(i1)=n1+1-i1
1574 : : end do
1575 : : end if
1576 [ + + ]: 55506 : if(istwf_k>=2 .and. istwf_k<=5)then
1577 : 35197 : i2inver(1)=1
1578 [ + + ]: 318981 : do i2=2,n2
1579 : 318981 : i2inver(i2)=n2+2-i2
1580 : : end do
1581 : : else
1582 [ + + ]: 227330 : do i2=1,n2
1583 : 227330 : i2inver(i2)=n2+1-i2
1584 : : end do
1585 : : end if
1586 [ + + + + ]: 55506 : if(istwf_k==2 .or. istwf_k==3 .or. istwf_k==6 .or. istwf_k==7)then
1587 : 36666 : i3inver(1)=1
1588 [ + + ]: 396627 : do i3=2,n3
1589 : 396627 : i3inver(i3)=n3+2-i3
1590 : : end do
1591 : : else
1592 [ + + ]: 240219 : do i3=1,n3
1593 : 240219 : i3inver(i3)=n3+1-i3
1594 : : end do
1595 : : end if
1596 : : end if
1597 : :
1598 [ + + ]: 2553219 : if (iflag==1 .or. iflag==2) then
1599 : : ! Insert cg into cfft with extra 0 s around outside:
1600 [ + + + + : 19214889626 : cfft = zero
+ + + + ]
1601 : :
1602 : : ! Take care of each plane wave, and complete cfft if needed
1603 [ + + ]: 1288645 : if (istwf_k==1) then
1604 : :
1605 [ + + ]: 1233902 : if (iflag==1) then
1606 : : !$OMP PARALLEL DO PRIVATE(i1,i2,i3) IF (ndat>1)
1607 [ + + ]: 2260808 : do idat=1,ndat
1608 [ + + ]: 492825609 : do ipw=1,npw
1609 [ + + ]: 490564801 : i1=kg_k(1,ipw); if(i1<0) i1=i1+n1; i1=i1+1
1610 [ + + ]: 490564801 : i2=kg_k(2,ipw); if(i2<0) i2=i2+n2; i2=i2+1
1611 [ + + ]: 490564801 : i3=kg_k(3,ipw); if(i3<0) i3=i3+n3; i3=i3+1
1612 : :
1613 : 490564801 : cfft(1,i1,i2,i3+n6*(idat-1))=cg(1,ipw+npw*(idat-1))
1614 : 491695229 : cfft(2,i1,i2,i3+n6*(idat-1))=cg(2,ipw+npw*(idat-1))
1615 : : end do
1616 : : end do
1617 : : end if
1618 : :
1619 [ + + ]: 1233902 : if (iflag==2) then
1620 : : !$OMP PARALLEL DO PRIVATE(i1,i2,i3) IF (ndat>1)
1621 [ + + ]: 239643 : do idat=1,ndat
1622 [ + + ]: 48993833 : do ipw=1,npw
1623 [ + + ]: 48754190 : i1=kg_k(1,ipw); if(i1<0) i1=i1+n1; i1=i1+1
1624 [ + + ]: 48754190 : i2=kg_k(2,ipw); if(i2<0) i2=i2+n2; i2=i2+1
1625 [ + + ]: 48754190 : i3=kg_k(3,ipw); if(i3<0) i3=i3+n3; i3=i3+1
1626 : :
1627 : 48754190 : cfft(1,i1,i3,i2+n6*(idat-1))=cg(1,ipw+npw*(idat-1))
1628 : 48890311 : cfft(2,i1,i3,i2+n6*(idat-1))=cg(2,ipw+npw*(idat-1))
1629 : : end do
1630 : : end do
1631 : : end if
1632 : :
1633 [ + - ]: 54743 : else if (istwf_k>=2) then
1634 : :
1635 : 54743 : npwmin=1
1636 [ + + + - ]: 54743 : if (istwf_k==2 .and. me_g0==1) then
1637 : : ! If gamma point, then cfft must be completed
1638 [ + + ]: 20336 : do idat=1,ndat
1639 : 10170 : cfft(1,1,1,1+n6*(idat-1))=cg(1,1+npw*(idat-1))
1640 : 20336 : cfft(2,1,1,1+n6*(idat-1))=zero
1641 : : end do
1642 : : npwmin=2
1643 : : end if
1644 : :
1645 [ + + ]: 54743 : if (iflag==1) then
1646 : : !$OMP PARALLEL DO PRIVATE(i1,i1inv,i2,i2inv,i3,i3inv) IF (ndat>1)
1647 [ + + ]: 106114 : do idat=1,ndat
1648 [ + + ]: 16057466 : do ipw=npwmin,npw
1649 [ + + ]: 15951352 : i1=kg_k(1,ipw); if(i1<0) i1=i1+n1; i1=i1+1
1650 [ - + ]: 15951352 : i2=kg_k(2,ipw); if(i2<0) i2=i2+n2; i2=i2+1
1651 [ + + ]: 15951352 : i3=kg_k(3,ipw); if(i3<0) i3=i3+n3; i3=i3+1
1652 : : ! Construct the coordinates of -k-G
1653 : 15951352 : i1inv=i1inver(i1) ; i2inv=i2inver(i2) ; i3inv=i3inver(i3)
1654 : :
1655 : 15951352 : cfft(1,i1,i2,i3+n6*(idat-1))=cg(1,ipw+npw*(idat-1))
1656 : 15951352 : cfft(2,i1,i2,i3+n6*(idat-1))=cg(2,ipw+npw*(idat-1))
1657 : 15951352 : cfft(1,i1inv,i2inv,i3inv+n6*(idat-1))= cg(1,ipw+npw*(idat-1))
1658 : 16004411 : cfft(2,i1inv,i2inv,i3inv+n6*(idat-1))=-cg(2,ipw+npw*(idat-1))
1659 : : end do
1660 : : end do
1661 : : end if
1662 : :
1663 [ + + ]: 54743 : if (iflag==2) then
1664 : : !$OMP PARALLEL DO PRIVATE(i1,i1inv,i2,i2inv,i3,i3inv) IF (ndat>1)
1665 [ + + ]: 3376 : do idat=1,ndat
1666 [ + + ]: 5483232 : do ipw=npwmin,npw
1667 [ + + ]: 5479856 : i1=kg_k(1,ipw); if(i1<0) i1=i1+n1; i1=i1+1
1668 [ - + ]: 5479856 : i2=kg_k(2,ipw); if(i2<0) i2=i2+n2; i2=i2+1
1669 [ + + ]: 5479856 : i3=kg_k(3,ipw); if(i3<0) i3=i3+n3; i3=i3+1
1670 : :
1671 : : ! Construct the coordinates of -k-G
1672 : 5479856 : i1inv=i1inver(i1) ; i2inv=i2inver(i2) ; i3inv=i3inver(i3)
1673 : :
1674 : 5479856 : cfft(1,i1,i3,i2+n6*(idat-1))=cg(1,ipw+npw*(idat-1))
1675 : 5479856 : cfft(2,i1,i3,i2+n6*(idat-1))=cg(2,ipw+npw*(idat-1))
1676 : 5479856 : cfft(1,i1inv,i3inv,i2inv+n6*(idat-1))= cg(1,ipw+npw*(idat-1))
1677 : 5481544 : cfft(2,i1inv,i3inv,i2inv+n6*(idat-1))=-cg(2,ipw+npw*(idat-1))
1678 : : end do
1679 : : end do
1680 : : end if
1681 : :
1682 : : end if
1683 : :
1684 [ + - ]: 1264574 : else if (iflag==-1 .or. iflag==-2) then
1685 : : ! extract cg(output) from cfft(input)
1686 : :
1687 : 1264574 : use_symmetry=0
1688 : 1264574 : identity(:,:)=0; identity(1,1)=1 ; identity(2,2)=1 ; identity(3,3)=1
1689 [ + + + + : 16439462 : if(sum((symrec(:,:)-identity(:,:))**2)/=0)use_symmetry=1
+ + ]
1690 [ + + + + ]: 5058296 : if(sum(shiftg(:)**2)/=0)use_symmetry=1
1691 : :
1692 : : ! Extract cg from cfft, ignoring components outside range of cg:
1693 [ + + ]: 1264574 : if (istwf_k==1) then
1694 : :
1695 [ + + ]: 1263811 : if (use_symmetry==0) then
1696 [ + + ]: 159268 : if (iflag==-1) then
1697 : : !$OMP PARALLEL DO PRIVATE(i1,i2,i3,ipwdat,i3dat) IF (ndat>1)
1698 [ + + ]: 144610 : do idat=1,ndat
1699 [ + + ]: 38793659 : do ipw=1,npw
1700 [ + + ]: 38649049 : i1=kg_k(1,ipw); if(i1<0) i1=i1+n1; i1=i1+1
1701 [ + + ]: 38649049 : i2=kg_k(2,ipw); if(i2<0) i2=i2+n2; i2=i2+1
1702 [ + + ]: 38649049 : i3=kg_k(3,ipw); if(i3<0) i3=i3+n3; i3=i3+1
1703 : 38649049 : ipwdat = ipw + (idat-1) * npw
1704 : 38649049 : i3dat = i3 + (idat-1) * n6
1705 : :
1706 : 38649049 : cg(1,ipwdat)=cfft(1,i1,i2,i3dat)*xnorm
1707 : 38721356 : cg(2,ipwdat)=cfft(2,i1,i2,i3dat)*xnorm
1708 : : end do
1709 : : end do
1710 : : else
1711 : : ! iflag==-2
1712 : : !$OMP PARALLEL DO PRIVATE(i1,i2,i3,ipwdat,i2dat) IF (ndat>1)
1713 [ + + ]: 206509 : do idat=1,ndat
1714 [ + + ]: 40573243 : do ipw=1,npw
1715 [ + + ]: 40366734 : i1=kg_k(1,ipw); if(i1<0) i1=i1+n1; i1=i1+1
1716 [ + + ]: 40366734 : i2=kg_k(2,ipw); if(i2<0) i2=i2+n2; i2=i2+1
1717 [ + + ]: 40366734 : i3=kg_k(3,ipw); if(i3<0) i3=i3+n3; i3=i3+1
1718 : :
1719 : 40366734 : ipwdat = ipw + (idat-1) * npw
1720 : 40366734 : i2dat = i2 + (idat-1) * n6
1721 : :
1722 : 40366734 : cg(1,ipwdat)=cfft(1,i1,i3,i2dat)*xnorm
1723 : 40486278 : cg(2,ipwdat)=cfft(2,i1,i3,i2dat)*xnorm
1724 : : end do
1725 : : end do
1726 : : end if
1727 : : else
1728 : : ! use_symmetry == 1
1729 : : !$OMP PARALLEL DO PRIVATE(i1,i2,i3,j1,j2,j3,l1,l2,l3,ipwdat,i3dat) IF (ndat>1)
1730 [ + + ]: 2209086 : do idat=1,ndat
1731 [ + + ]: 464674009 : do ipw=1,npw
1732 : 462464923 : l1=kg_k(1,ipw)+shiftg(1)
1733 : 462464923 : l2=kg_k(2,ipw)+shiftg(2)
1734 : 462464923 : l3=kg_k(3,ipw)+shiftg(3)
1735 : 462464923 : j1=symrec(1,1)*l1+symrec(1,2)*l2+symrec(1,3)*l3
1736 : 462464923 : j2=symrec(2,1)*l1+symrec(2,2)*l2+symrec(2,3)*l3
1737 : 462464923 : j3=symrec(3,1)*l1+symrec(3,2)*l2+symrec(3,3)*l3
1738 [ + + ]: 462464923 : if(j1<0) j1=j1+n1; i1=j1+1
1739 [ + + ]: 462464923 : if(j2<0) j2=j2+n2; i2=j2+1
1740 [ + + ]: 462464923 : if(j3<0) j3=j3+n3; i3=j3+1
1741 : : ! [i1, i2, i3] are the indices of S(g + g0) in the FFT box.
1742 : : ! while ipw is the index of g in kg_k
1743 : :
1744 : 462464923 : ipwdat = ipw + (idat-1) * npw
1745 : 462464923 : i3dat = i3 + (idat-1)*n6
1746 : :
1747 : : ! c(g) = cfft(S(g + shiftg))
1748 : 462464923 : cg(1,ipwdat)=cfft(1,i1,i2,i3dat)*xnorm
1749 : 463569466 : cg(2,ipwdat)=cfft(2,i1,i2,i3dat)*xnorm
1750 : : end do
1751 : : end do
1752 : : end if
1753 : :
1754 [ + - ]: 763 : else if (istwf_k>=2) then
1755 : :
1756 : 763 : npwmin=1
1757 [ + + + - ]: 763 : if (istwf_k==2 .and. me_g0==1) then
1758 : : ! Extract cg from cfft, in a way that projects on a
1759 : : ! wavefunction with time-reversal symmetry
1760 [ + + ]: 274 : do idat=1,ndat
1761 : 137 : ipwdat = 1 + (idat-1) * npw
1762 : 137 : i3dat = 1 + (idat-1)*n6
1763 : 137 : cg(1,ipwdat)=cfft(1,1,1,i3dat)*xnorm
1764 : 274 : cg(2,ipwdat)=zero
1765 : : end do
1766 : : npwmin=2
1767 : : end if
1768 : :
1769 [ + + ]: 763 : if (use_symmetry==0) then
1770 : :
1771 [ + + ]: 367 : if (iflag==-1) then
1772 : : !$OMP PARALLEL DO PRIVATE(i1,i1inv,i2,i2inv,i3,i3inv,ipwdat,i3dat,i3invdat) IF (ndat>1)
1773 [ + + ]: 434 : do idat=1,ndat
1774 [ + + ]: 42206 : do ipw=npwmin,npw
1775 [ + + ]: 41772 : i1=kg_k(1,ipw); if(i1<0) i1=i1+n1; i1=i1+1
1776 [ - + ]: 41772 : i2=kg_k(2,ipw); if(i2<0) i2=i2+n2; i2=i2+1
1777 [ + + ]: 41772 : i3=kg_k(3,ipw); if(i3<0) i3=i3+n3; i3=i3+1
1778 : :
1779 : : ! Construct the coordinates of -k-G
1780 : 41772 : i1inv=i1inver(i1); i2inv=i2inver(i2); i3inv=i3inver(i3)
1781 : :
1782 : 41772 : ipwdat = ipw + (idat-1) * npw
1783 : 41772 : i3dat = i3 + (idat-1) * n6
1784 : 41772 : i3invdat = i3inv + (idat-1) * n6
1785 : :
1786 : : ! Here the time-reversal symmetry is used to project from cfft
1787 : 41772 : cg(1,ipwdat)=(cfft(1,i1,i2,i3dat) + cfft(1,i1inv,i2inv,i3invdat))*0.5d0*xnorm
1788 : 41989 : cg(2,ipwdat)=(cfft(2,i1,i2,i3dat) - cfft(2,i1inv,i2inv,i3invdat))*0.5d0*xnorm
1789 : : end do
1790 : : end do
1791 : :
1792 : : else
1793 : : !$OMP PARALLEL DO PRIVATE(i1,i1inv,i2,i2inv,i3,i3inv,ipwdat,i2dat,i2invdat) IF (ndat>1)
1794 [ + + ]: 300 : do idat=1,ndat
1795 [ + + ]: 1745174 : do ipw=npwmin,npw
1796 [ + + ]: 1744874 : i1=kg_k(1,ipw); if(i1<0) i1=i1+n1; i1=i1+1
1797 [ - + ]: 1744874 : i2=kg_k(2,ipw); if(i2<0) i2=i2+n2; i2=i2+1
1798 [ + + ]: 1744874 : i3=kg_k(3,ipw); if(i3<0) i3=i3+n3; i3=i3+1
1799 : :
1800 : : ! Construct the coordinates of -k-G
1801 : 1744874 : i1inv=i1inver(i1) ; i2inv=i2inver(i2) ; i3inv=i3inver(i3)
1802 : :
1803 : 1744874 : ipwdat = ipw + (idat-1) * npw
1804 : 1744874 : i2dat = i2 + (idat-1) * n6
1805 : 1744874 : i2invdat = i2inv + (idat-1) * n6
1806 : :
1807 : : ! Here the time-reversal symmetry is used to project from cfft
1808 : 1744874 : cg(1,ipwdat)=(cfft(1,i1,i3,i2dat) + cfft(1,i1inv,i3inv,i2invdat))*0.5d0*xnorm
1809 : 1745024 : cg(2,ipwdat)=(cfft(2,i1,i3,i2dat) - cfft(2,i1inv,i3inv,i2invdat))*0.5d0*xnorm
1810 : : end do
1811 : : end do
1812 : : end if
1813 : :
1814 : : else
1815 : : ! Use symmetry
1816 : 396 : id1=n1/2+2
1817 : 396 : id2=n2/2+2
1818 : 396 : id3=n3/2+2
1819 : :
1820 : : !$OMP PARALLEL DO PRIVATE(i1,i1inv,i2,i2inv,i3,i3inv,j1,j2,j3,l1,l2,l3,ipwdat,i3dat,i3invdat) IF (ndat>1)
1821 [ + + ]: 792 : do idat=1,ndat
1822 [ + + ]: 48714 : do ipw=npwmin,npw
1823 : :
1824 [ + + ]: 47922 : i1=kg_k(1,ipw); if(i1<0) i1=i1+n1; i1=i1+1
1825 [ - + ]: 47922 : i2=kg_k(2,ipw); if(i2<0) i2=i2+n2; i2=i2+1
1826 [ + + ]: 47922 : i3=kg_k(3,ipw); if(i3<0) i3=i3+n3; i3=i3+1
1827 : :
1828 : 47922 : i1inv=i1inver(i1) ; i2inv=i2inver(i2) ; i3inv=i3inver(i3)
1829 : :
1830 : 47922 : l1=kg_k(1,ipw)+shiftg(1)
1831 : 47922 : l2=kg_k(2,ipw)+shiftg(2)
1832 : 47922 : l3=kg_k(3,ipw)+shiftg(3)
1833 : 47922 : j1=symrec(1,1)*l1+symrec(1,2)*l2+symrec(1,3)*l3
1834 : 47922 : j2=symrec(2,1)*l1+symrec(2,2)*l2+symrec(2,3)*l3
1835 : 47922 : j3=symrec(3,1)*l1+symrec(3,2)*l2+symrec(3,3)*l3
1836 [ + + ]: 47922 : if(j1<0)j1=j1+n1 ; i1=j1+1
1837 [ + + ]: 47922 : if(j2<0)j2=j2+n2 ; i2=j2+1
1838 [ + + ]: 47922 : if(j3<0)j3=j3+n3 ; i3=j3+1
1839 : :
1840 : : ! Construct the coordinates of -k-G
1841 : 47922 : l1=i1inv-(i1inv/id1)*n1-1+shiftg(1)
1842 : 47922 : l2=i2inv-(i2inv/id2)*n2-1+shiftg(2)
1843 : 47922 : l3=i3inv-(i3inv/id3)*n3-1+shiftg(3)
1844 : 47922 : j1=symrec(1,1)*l1+symrec(1,2)*l2+symrec(1,3)*l3
1845 : 47922 : j2=symrec(2,1)*l1+symrec(2,2)*l2+symrec(2,3)*l3
1846 : 47922 : j3=symrec(3,1)*l1+symrec(3,2)*l2+symrec(3,3)*l3
1847 [ + + ]: 47922 : if(j1<0)j1=j1+n1 ; i1inv=j1+1
1848 [ + + ]: 47922 : if(j2<0)j2=j2+n2 ; i2inv=j2+1
1849 [ + + ]: 47922 : if(j3<0)j3=j3+n3 ; i3inv=j3+1
1850 : :
1851 : 47922 : ipwdat = ipw + (idat-1) * npw
1852 : 47922 : i3dat = i3 + (idat-1) * n6
1853 : 47922 : i3invdat = i3inv + (idat-1) * n6
1854 : :
1855 : : ! Here the time-reversal symmetry is used to project from cfft
1856 : 47922 : cg(1,ipwdat)=(cfft(1,i1,i2,i3dat) + cfft(1,i1inv,i2inv,i3invdat))*0.5d0*xnorm
1857 : 48318 : cg(2,ipwdat)=(cfft(2,i1,i2,i3dat) - cfft(2,i1inv,i2inv,i3invdat))*0.5d0*xnorm
1858 : : end do
1859 : : end do
1860 : : end if
1861 : :
1862 : : end if
1863 : :
1864 : : else
1865 : 0 : ABI_BUG(sjoin('iflag: ', itoa(iflag),' not acceptable.'))
1866 : : end if
1867 : :
1868 : : DBG_EXIT("COLL")
1869 : :
1870 : 2553219 : end subroutine sphere
1871 : : !!***
1872 : :
1873 : : !----------------------------------------------------------------------
1874 : :
1875 : : !!****f* m_fftcore/sphere_fft
1876 : : !! NAME
1877 : : !! sphere_fft
1878 : : !!
1879 : : !! FUNCTION
1880 : : !! Array cg is defined in sphere with npw points. Insert cg inside box
1881 : : !! of n1*n2*n3 points to define array cfft for fft box.
1882 : : !! corresponds to given element in cg. rest of cfft is filled with 0 s.
1883 : : !!
1884 : : !! iflag=1==>insert cg into cfft.
1885 : : !! iflag=2==>insert cg into cfft, where the second and third dimension
1886 : : !! have been switched (needed for new 2002 SGoedecker FFT)
1887 : : !! iflag=-1==> extract cg from cfft.
1888 : : !! iflag=-2==> extract cg from cfft, where the second and third dimension
1889 : : !! have been switched (needed for new 2002 SGoedecker FFT)
1890 : : !! (WARNING : iflag=-2 cannot use symmetry operations)
1891 : : !!
1892 : : !! There is also the possibility to apply a symmetry operation,
1893 : : !! as well as to make a shift in reciprocal space, or to multiply
1894 : : !! by a constant factor, in the case iflag=-1.
1895 : : !! Multiplication by a constant factor is also possible in the case iflag=-2.
1896 : : !!
1897 : : !! INPUTS
1898 : : !! cg(2,npw)= contains values for npw G vectors in basis sphere
1899 : : !! ndat=number of FFT to do in //
1900 : : !! npw=number of G vectors in basis at this k point
1901 : : !! cfft(2,n4,n5,n6) = fft box
1902 : : !! n1,n2,n3=physical dimension of the box (cfft)
1903 : : !! n4,n5,n6=memory dimension of cfft
1904 : : !! kg_k(3,npw)=integer coordinates of G vectors in basis sphere
1905 : : !! mpi_enreg=information about MPI parallelization
1906 : : !! tab_fftwf2_local(n2)=local i2 indices in fourwf
1907 : : !! nd2proc TO BE DESCRIBED SB 090831
1908 : : !! iflag=option parameter. Possible values: -1, -2, 1, 2 ; this is used only in debug option
1909 : : !!
1910 : : !! OUTPUT
1911 : : !! (see side effects)
1912 : : !!
1913 : : !! SIDE EFFECTS
1914 : : !! Input/Output
1915 : : !! iflag=1 and 2, insert cg(input) into cfft(output)
1916 : : !! iflag=-1 and -2, extract cg(output) from cfft(input)
1917 : : !!
1918 : : !! NOTES
1919 : : !! cg and cfft are assumed to be of type COMPLEX, although this routine treats
1920 : : !! them as real of twice the length to avoid nonstandard complex*16.
1921 : : !!
1922 : : !! WARNING
1923 : : !! NO CHECK is DONE over iflag.
1924 : : !!
1925 : : !! TODO
1926 : : !! Order arguments
1927 : : !!
1928 : : !! SOURCE
1929 : :
1930 : 0 : subroutine sphere_fft(cg,ndat,npw,cfft,n1,n2,n3,n4,n5,kg_k,tab_fftwf2_local,nd2proc)
1931 : :
1932 : : !Arguments ------------------------------------
1933 : : !scalars
1934 : : integer,intent(in) :: n1,n2,n3,n4,n5,nd2proc,ndat,npw
1935 : : integer,intent(in) :: tab_fftwf2_local(n2)
1936 : : !arrays
1937 : : integer,intent(in) :: kg_k(3,npw)
1938 : : real(dp),intent(in) :: cg(2,npw*ndat)
1939 : : real(dp),intent(out) :: cfft(2,n4,n5,nd2proc*ndat)
1940 : :
1941 : : !Local variables-------------------------------
1942 : : !scalars
1943 : : integer :: i1,i2,i2_local,i3,idat,ipw
1944 : : ! *************************************************************************
1945 : :
1946 : : !Insert cg into cfft with extra 0 s around outside:
1947 [ # # # # : 0 : cfft = zero
# # # # ]
1948 : :
1949 : : !$OMP PARALLEL DO PRIVATE(i1,i2,i2_local,i3)
1950 [ # # ]: 0 : do ipw=1,npw
1951 [ # # ]: 0 : i1=kg_k(1,ipw); if(i1<0)i1=i1+n1; i1=i1+1
1952 [ # # ]: 0 : i2=kg_k(2,ipw); if(i2<0)i2=i2+n2; i2=i2+1
1953 [ # # ]: 0 : i3=kg_k(3,ipw); if(i3<0)i3=i3+n3; i3=i3+1
1954 : 0 : i2_local = tab_fftwf2_local(i2)
1955 [ # # ]: 0 : do idat=1,ndat
1956 : 0 : cfft(1,i1,i3,i2_local + nd2proc*(idat-1))=cg(1,ipw+npw*(idat-1))
1957 : 0 : cfft(2,i1,i3,i2_local + nd2proc*(idat-1))=cg(2,ipw+npw*(idat-1))
1958 : : end do
1959 : : end do
1960 : :
1961 : 0 : end subroutine sphere_fft
1962 : : !!***
1963 : :
1964 : : !----------------------------------------------------------------------
1965 : :
1966 : : !!****f* m_fftcore/sphere_fft1
1967 : : !! NAME
1968 : : !! sphere_fft1
1969 : : !!
1970 : : !! FUNCTION
1971 : : !! Array cg is defined in sphere with npw points. Insert cg inside box
1972 : : !! of n1*n2*n3 points to define array cfft for fft box.
1973 : : !! corresponds to given element in cg. rest of cfft is filled with 0 s.
1974 : : !!
1975 : : !! iflag=1==>insert cg into cfft.
1976 : : !! iflag=2==>insert cg into cfft, where the second and third dimension
1977 : : !! have been switched (needed for new 2002 SGoedecker FFT)
1978 : : !! iflag=-1==> extract cg from cfft.
1979 : : !! iflag=-2==> extract cg from cfft, where the second and third dimension
1980 : : !! have been switched (needed for new 2002 SGoedecker FFT)
1981 : : !! (WARNING : iflag=-2 cannot use symmetry operations)
1982 : : !!
1983 : : !! There is also the possibility to apply a symmetry operation,
1984 : : !! as well as to make a shift in reciprocal space, or to multiply
1985 : : !! by a constant factor, in the case iflag=-1.
1986 : : !! Multiplication by a constant factor is also possible in the case iflag=-2.
1987 : : !!
1988 : : !! INPUTS
1989 : : !! cg(2,npw)= contains values for npw G vectors in basis sphere
1990 : : !! ndat=number of FFT to do in //
1991 : : !! npw=number of G vectors in basis at this k point
1992 : : !! cfft(2,n4,n5,n6) = fft box
1993 : : !! n1,n2,n3=physical dimension of the box (cfft)
1994 : : !! n4,n5,n6=memory dimension of cfft
1995 : : !! kg_k(3,npw)=integer coordinates of G vectors in basis sphere
1996 : : !! nd2proc TO BE DESCRIBED SB 090831
1997 : : !! iflag=option parameter. Possible values: -1, -2, 1, 2
1998 : : !! tab_fftwf2_local(n2)=local i2 indices in fourwf
1999 : : !!
2000 : : !! OUTPUT
2001 : : !! (see side effects)
2002 : : !!
2003 : : !! SIDE EFFECTS
2004 : : !! Input/Output
2005 : : !! iflag=1 and 2, insert cg(input) into cfft(output)
2006 : : !! iflag=-1 and -2, extract cg(output) from cfft(input)
2007 : : !!
2008 : : !! NOTES
2009 : : !! cg and cfft are assumed to be of type COMPLEX, although this routine treats
2010 : : !! them as real of twice the length to avoid nonstandard complex*16.
2011 : : !!
2012 : : !! WARNING
2013 : : !! NO CHECK is DONE over iflag.
2014 : : !!
2015 : : !! TODO
2016 : : !! Order arguments
2017 : : !! sphere_fft1 is similar to sphere_fft, the only difference being that ndat > 1 is not supported.
2018 : : !! Why? Should merge the two APIs.
2019 : : !!
2020 : : !! SOURCE
2021 : :
2022 : 700212 : subroutine sphere_fft1(cg,ndat,npw,cfft,n1,n2,n3,n4,n5,n6,kg_k,tab_fftwf2_local)
2023 : :
2024 : :
2025 : : !Arguments ------------------------------------
2026 : : !scalars
2027 : : integer,intent(in) :: n1,n2,n3,n4,n5,n6,ndat,npw
2028 : : !arrays
2029 : : integer,intent(in) :: kg_k(3,npw)
2030 : : integer,intent(in) :: tab_fftwf2_local(n2)
2031 : : real(dp),intent(in) :: cg(2,npw*ndat)
2032 : : real(dp),intent(inout) :: cfft(2,n4,n5,n6*ndat)
2033 : :
2034 : : !Local variables-------------------------------
2035 : : !scalars
2036 : : integer :: i1,i2,i2_local,i3,idat,ipw
2037 : : ! *************************************************************************
2038 : :
2039 : : !Insert cg into cfft with extra 0 s around outside:
2040 : :
2041 [ + + + + : 1631896820 : cfft = zero
+ + + + ]
2042 : : !$OMP PARALLEL DO PRIVATE(i1,i2,i2_local,i3)
2043 [ + + ]: 1426840 : do idat=1,ndat
2044 [ + + ]: 196417943 : do ipw=1,npw
2045 [ + + ]: 194991103 : i1=kg_k(1,ipw); if(i1<0)i1=i1+n1; i1=i1+1
2046 [ + + ]: 194991103 : i2=kg_k(2,ipw); if(i2<0)i2=i2+n2; i2=i2+1
2047 [ + + ]: 194991103 : i3=kg_k(3,ipw); if(i3<0)i3=i3+n3; i3=i3+1
2048 : 194991103 : i2_local = tab_fftwf2_local(i2) + n6*(idat-1)
2049 : 194991103 : cfft(1,i1,i3,i2_local)=cg(1,ipw+npw*(idat-1))
2050 : 195717731 : cfft(2,i1,i3,i2_local)=cg(2,ipw+npw*(idat-1))
2051 : : end do
2052 : : end do
2053 : :
2054 : 700212 : end subroutine sphere_fft1
2055 : : !!***
2056 : :
2057 : : !----------------------------------------------------------------------
2058 : :
2059 : : !!****f* m_fftcore/change_istwfk
2060 : : !! NAME
2061 : : !! change_istwfk
2062 : : !!
2063 : : !! FUNCTION
2064 : : !! This function allows one to change the time-reversal storage mode (istwfk)
2065 : : !! of a *full* set of u(G). It does not support MPI-FFT!
2066 : : !!
2067 : : !! INPUTS
2068 : : !! from_npw=number of G vectors in input from_cg
2069 : : !! from_kg_k(3,npw)=integer coordinates of the G vectors of from_cg
2070 : : !! from_istwfk=option parameter that describes the storage in from_cg
2071 : : !! to_npw=number of G vectors in output to_cg
2072 : : !! to_kg_k(3,npw)=integer coordinates of the G vectors in to_cg
2073 : : !! to_istwfk=option parameter that describes the storage in to_cg
2074 : : !! n1,n2,n3=physical dimension of the box (must be large enough to contain the sphere, no check is done)
2075 : : !! ndat=number of wavefunctions
2076 : : !! from_cg(2,from_npw*ndat)= Input u(g) values
2077 : : !!
2078 : : !! OUTPUTS
2079 : : !! to_cg(2,to_npw*ndat)= Output u(g) defined on the list of vectors to_kg_k with time-reversal mode to_istwfk
2080 : : !!
2081 : : !! SOURCE
2082 : :
2083 : 4 : subroutine change_istwfk(from_npw,from_kg,from_istwfk,to_npw,to_kg,to_istwfk,n1,n2,n3,ndat,from_cg,to_cg)
2084 : :
2085 : : !Arguments ------------------------------------
2086 : : !scalars
2087 : : integer,intent(in) :: from_npw,from_istwfk,to_npw,to_istwfk,n1,n2,n3,ndat
2088 : : !arrays
2089 : : integer,intent(in) :: from_kg(3,from_npw),to_kg(3,to_npw)
2090 : : real(dp),intent(inout) :: from_cg(2,from_npw*ndat) ! out due to sphere!
2091 : : real(dp),intent(inout) :: to_cg(2,to_npw*ndat)
2092 : :
2093 : : !Local variables-------------------------------
2094 : : !scalars
2095 : : integer :: n4,n5,n6
2096 : : real(dp),parameter :: xnorm1=one
2097 : : !arrays
2098 : : integer,parameter :: shiftg0(3)=0,me_g0=1
2099 : : integer,parameter :: symmE(3,3)=reshape([1,0,0,0,1,0,0,0,1],[3,3])
2100 : 4 : real(dp),allocatable :: cfft(:,:,:,:)
2101 : : ! *************************************************************************
2102 : :
2103 : 4 : n4=2*(n1/2)+1
2104 : 4 : n5=2*(n2/2)+1
2105 : 4 : n6=2*(n3/2)+1
2106 : :
2107 [ - + + - : 20 : ABI_MALLOC(cfft, (2,n4,n5,n6*ndat))
- + + - -
+ + - + -
+ - - + -
+ ]
2108 : :
2109 : : ! iflag=1 ==> insert from_cg into cfft.
2110 : 4 : call sphere(from_cg,ndat,from_npw,cfft,n1,n2,n3,n4,n5,n6,from_kg,from_istwfk,+1,me_g0,shiftg0,symmE,xnorm1)
2111 : :
2112 : : ! iflag=-1 ==> extract to_cg from cfft.
2113 : 4 : call sphere(to_cg,ndat,to_npw,cfft,n1,n2,n3,n4,n5,n6,to_kg,to_istwfk,-1,me_g0,shiftg0,symmE,xnorm1)
2114 : :
2115 [ - + ]: 4 : ABI_FREE(cfft)
2116 : :
2117 : 4 : end subroutine change_istwfk
2118 : : !!***
2119 : :
2120 : : !----------------------------------------------------------------------
2121 : :
2122 : : !!****f* m_fftcore/switch
2123 : : !! NAME
2124 : : !! switch
2125 : : !!
2126 : : !! FUNCTION
2127 : : !!
2128 : : !! INPUTS
2129 : : !!
2130 : : !! OUTPUT
2131 : : !!
2132 : : !! SOURCE
2133 : :
2134 : 1625379 : pure subroutine switch(n1dfft,n2,lot,n1,lzt,zt,zw)
2135 : :
2136 : : !Arguments ------------------------------------
2137 : : integer,intent(in) :: n1dfft,n2,lot,n1,lzt
2138 : : real(dp),intent(in) :: zt(2,lzt,n1)
2139 : : real(dp),intent(inout) :: zw(2,lot,n2)
2140 : :
2141 : : !Local variables-------------------------------
2142 : : integer :: i,j
2143 : : ! *************************************************************************
2144 : :
2145 [ + + ]: 23327006 : do j=1,n1dfft
2146 [ + + ]: 926458027 : do i=1,n2
2147 : 903131021 : zw(1,j,i)=zt(1,i,j)
2148 : 924832648 : zw(2,j,i)=zt(2,i,j)
2149 : : end do
2150 : : end do
2151 : :
2152 : 1625379 : end subroutine switch
2153 : : !!***
2154 : :
2155 : : !----------------------------------------------------------------------
2156 : :
2157 : : !!****f* m_fftcore/switch_cent
2158 : : !! NAME
2159 : : !! switch_cent
2160 : : !!
2161 : : !! FUNCTION
2162 : : !! Perform the rotation:
2163 : : !!
2164 : : !! input: I2,i1,j3,(jp3)
2165 : : !! output: i1,I2,j3,(jp3)
2166 : : !!
2167 : : !! and pad the signal with zeros.
2168 : : !!
2169 : : !! INPUTS
2170 : : !! n1dfft=Number of 1D FFTs to perform
2171 : : !! max2=Max G_y in the small box enclosing the G-sphere.
2172 : : !! m2=Size of the small box enclosing the G-sphere along y
2173 : : !! n2=Dimension of the transform along y
2174 : : !! lot=Cache blocking factor.
2175 : : !! n1=Dimension of the transform along x
2176 : : !! lzt=Second dimension of z
2177 : : !! zt(2,lzt,n1)
2178 : : !!
2179 : : !! OUTPUT
2180 : : !! zw(2,lot,n2)=Cache working array
2181 : : !!
2182 : : !! SOURCE
2183 : :
2184 : 11686163 : pure subroutine switch_cent(n1dfft,max2,m2,n2,lot,n1,lzt,zt,zw)
2185 : :
2186 : : !Arguments ------------------------------------
2187 : : integer,intent(in) :: n1dfft,max2,m2,n2,lot,n1,lzt
2188 : : real(dp),intent(in) :: zt(2,lzt,n1)
2189 : : real(dp),intent(inout) :: zw(2,lot,n2)
2190 : :
2191 : : !Local variables-------------------------------
2192 : : integer :: i,j
2193 : : ! *************************************************************************
2194 : :
2195 : : ! Here, zero and positive frequencies
2196 [ + + ]: 233206512 : do j=1,n1dfft
2197 [ + + ]: 1815774148 : do i=1,max2+1
2198 : 1582567636 : zw(1,j,i)=zt(1,i,j)
2199 : 1804087985 : zw(2,j,i)=zt(2,i,j)
2200 : : end do
2201 : : end do
2202 : :
2203 : : ! Fill the center region with zeros
2204 [ + + ]: 191549483 : do i=max2+2,n2-m2+max2+1
2205 [ + + ]: 3336224547 : do j=1,n1dfft
2206 : 3144675064 : zw(1,j,i)=zero
2207 : 3324538384 : zw(2,j,i)=zero
2208 : : end do
2209 : : end do
2210 : :
2211 : : ! Here, negative frequencies
2212 [ + - ]: 11686163 : if (m2>=max2+2) then
2213 [ + + ]: 233206512 : do j=1,n1dfft
2214 [ + + ]: 1509880797 : do i=max2+2,m2
2215 : 1276674285 : zw(1,j,i+n2-m2)=zt(1,i,j)
2216 : 1498194634 : zw(2,j,i+n2-m2)=zt(2,i,j)
2217 : : end do
2218 : : end do
2219 : : end if
2220 : :
2221 : 11686163 : end subroutine switch_cent
2222 : : !!***
2223 : :
2224 : : !----------------------------------------------------------------------
2225 : :
2226 : : !!****f* m_fftcore/switchreal
2227 : : !! NAME
2228 : : !! switchreal
2229 : : !!
2230 : : !! FUNCTION
2231 : : !! Perform the rotation:
2232 : : !!
2233 : : !! input: I2,i1,j3,(jp3)
2234 : : !! output: i1,I2,j3,(jp3)
2235 : : !!
2236 : : !! and pad the signal with zeros.
2237 : : !! Used for real wavefunctions.
2238 : : !!
2239 : : !! INPUTS
2240 : : !! includelast
2241 : : !! n1dfft=Number of 1D FFTs to perform
2242 : : !! n2=Dimension of the transform along y
2243 : : !! n2eff
2244 : : !! lot=Cache blocking factor.
2245 : : !! n1zt
2246 : : !! lzt
2247 : : !! zt(2,lzt,n1zt)
2248 : : !!
2249 : : !! OUTPUT
2250 : : !! zw(2,lot,n2)
2251 : : !!
2252 : : !! SOURCE
2253 : :
2254 : 6480 : pure subroutine switchreal(includelast,n1dfft,n2,n2eff,lot,n1zt,lzt,zt,zw)
2255 : :
2256 : : !Arguments ------------------------------------
2257 : : integer,intent(in) :: includelast,n1dfft,n2,n2eff,lot,n1zt,lzt
2258 : : real(dp),intent(in) :: zt(2,lzt,n1zt)
2259 : : real(dp),intent(inout) :: zw(2,lot,n2)
2260 : :
2261 : : !Local variables-------------------------------
2262 : : integer :: i,j
2263 : : ! *************************************************************************
2264 : :
2265 [ + - ]: 6480 : if (includelast==1) then
2266 : :
2267 : : ! Compute symmetric and antisymmetric combinations
2268 [ + + ]: 74196 : do j=1,n1dfft
2269 : 67716 : zw(1,j,1)=zt(1,1,2*j-1)
2270 : 74196 : zw(2,j,1)=zt(1,1,2*j )
2271 : : end do
2272 [ + + ]: 220698 : do i=2,n2eff
2273 [ + + ]: 2323620 : do j=1,n1dfft
2274 : 2102922 : zw(1,j,i)= zt(1,i,2*j-1)-zt(2,i,2*j)
2275 : 2102922 : zw(2,j,i)= zt(2,i,2*j-1)+zt(1,i,2*j)
2276 : 2102922 : zw(1,j,n2+2-i)= zt(1,i,2*j-1)+zt(2,i,2*j)
2277 : 2317140 : zw(2,j,n2+2-i)=-zt(2,i,2*j-1)+zt(1,i,2*j)
2278 : : end do
2279 : : end do
2280 : :
2281 : : else
2282 : :
2283 : : ! An odd number of FFTs
2284 : : ! Compute symmetric and antisymmetric combinations
2285 [ # # ]: 0 : do j=1,n1dfft-1
2286 : 0 : zw(1,j,1)=zt(1,1,2*j-1)
2287 : 0 : zw(2,j,1)=zt(1,1,2*j )
2288 : : end do
2289 : 0 : zw(1,n1dfft,1)=zt(1,1,2*n1dfft-1)
2290 : 0 : zw(2,n1dfft,1)=zero
2291 : :
2292 [ # # ]: 0 : do i=2,n2eff
2293 [ # # ]: 0 : do j=1,n1dfft-1
2294 : 0 : zw(1,j,i)= zt(1,i,2*j-1)-zt(2,i,2*j)
2295 : 0 : zw(2,j,i)= zt(2,i,2*j-1)+zt(1,i,2*j)
2296 : 0 : zw(1,j,n2+2-i)= zt(1,i,2*j-1)+zt(2,i,2*j)
2297 : 0 : zw(2,j,n2+2-i)=-zt(2,i,2*j-1)+zt(1,i,2*j)
2298 : : end do
2299 : 0 : zw(1,n1dfft,i)= zt(1,i,2*n1dfft-1)
2300 : 0 : zw(2,n1dfft,i)= zt(2,i,2*n1dfft-1)
2301 : 0 : zw(1,n1dfft,n2+2-i)= zt(1,i,2*n1dfft-1)
2302 : 0 : zw(2,n1dfft,n2+2-i)=-zt(2,i,2*n1dfft-1)
2303 : : end do
2304 : : end if
2305 : :
2306 : 6480 : end subroutine switchreal
2307 : : !!***
2308 : :
2309 : : !----------------------------------------------------------------------
2310 : :
2311 : : !!****f* m_fftcore/switchreal_cent
2312 : : !! NAME
2313 : : !! switchreal_cent
2314 : : !!
2315 : : !! FUNCTION
2316 : : !! Perform the rotation:
2317 : : !!
2318 : : !! input: I2,i1,j3,(jp3)
2319 : : !! output: i1,I2,j3,(jp3)
2320 : : !!
2321 : : !! and pad the signal with zeros.
2322 : : !! Used for the Fourier transform of real wavefunctions.
2323 : : !!
2324 : : !! INPUTS
2325 : : !! includelast
2326 : : !! n1dfft=Number of 1D FFTs to perform
2327 : : !! max2
2328 : : !! n2=Dimension of the transform along y
2329 : : !! lot=Cache blocking factor.
2330 : : !! n1zt, lzt=Dimensions of zt.
2331 : : !! zt(2,lzt,n1zt)
2332 : : !!
2333 : : !! OUTPUT
2334 : : !! zw(2,lot,n2)
2335 : : !!
2336 : : !! SOURCE
2337 : :
2338 : 22764 : pure subroutine switchreal_cent(includelast,n1dfft,max2,n2,lot,n1zt,lzt,zt,zw)
2339 : :
2340 : : !Arguments ------------------------------------
2341 : : integer,intent(in) :: includelast,n1dfft,max2,n2,lot,n1zt,lzt
2342 : : real(dp),intent(in) :: zt(2,lzt,n1zt)
2343 : : real(dp),intent(inout) :: zw(2,lot,n2)
2344 : :
2345 : : !Local variables-------------------------------
2346 : : integer :: i,j
2347 : : ! *************************************************************************
2348 : :
2349 [ + + ]: 22764 : if (includelast==1) then
2350 : :
2351 : : ! Compute symmetric and antisymmetric combinations
2352 [ + + ]: 244872 : do j=1,n1dfft
2353 : 224508 : zw(1,j,1)=zt(1,1,2*j-1)
2354 : 244872 : zw(2,j,1)=zt(1,1,2*j )
2355 : : end do
2356 : :
2357 [ + + ]: 459984 : do i=2,max2+1
2358 [ + + ]: 5171724 : do j=1,n1dfft
2359 : 4711740 : zw(1,j,i)= zt(1,i,2*j-1)-zt(2,i,2*j)
2360 : 4711740 : zw(2,j,i)= zt(2,i,2*j-1)+zt(1,i,2*j)
2361 : 4711740 : zw(1,j,n2+2-i)= zt(1,i,2*j-1)+zt(2,i,2*j)
2362 : 5151360 : zw(2,j,n2+2-i)=-zt(2,i,2*j-1)+zt(1,i,2*j)
2363 : : end do
2364 : : end do
2365 : :
2366 [ + - ]: 20364 : if(max2+1<n2-max2)then
2367 [ + + ]: 979896 : do i=max2+2,n2-max2
2368 [ + + ]: 11345700 : do j=1,n1dfft
2369 : 10365804 : zw(1,j,i)=zero
2370 : 11325336 : zw(2,j,i)=zero
2371 : : end do
2372 : : end do
2373 : : end if
2374 : :
2375 : : else
2376 : : ! Compute symmetric and antisymmetric combinations
2377 [ + + ]: 24000 : do j=1,n1dfft-1
2378 : 21600 : zw(1,j,1)=zt(1,1,2*j-1)
2379 : 24000 : zw(2,j,1)=zt(1,1,2*j )
2380 : : end do
2381 : :
2382 : 2400 : zw(1,n1dfft,1)=zt(1,1,2*n1dfft-1)
2383 : 2400 : zw(2,n1dfft,1)=zero
2384 [ + + ]: 40800 : do i=2,max2+1
2385 [ + + ]: 384000 : do j=1,n1dfft-1
2386 : 345600 : zw(1,j,i)= zt(1,i,2*j-1)-zt(2,i,2*j)
2387 : 345600 : zw(2,j,i)= zt(2,i,2*j-1)+zt(1,i,2*j)
2388 : 345600 : zw(1,j,n2+2-i)= zt(1,i,2*j-1)+zt(2,i,2*j)
2389 : 384000 : zw(2,j,n2+2-i)=-zt(2,i,2*j-1)+zt(1,i,2*j)
2390 : : end do
2391 : 38400 : zw(1,n1dfft,i)= zt(1,i,2*n1dfft-1)
2392 : 38400 : zw(2,n1dfft,i)= zt(2,i,2*n1dfft-1)
2393 : 38400 : zw(1,n1dfft,n2+2-i)= zt(1,i,2*n1dfft-1)
2394 : 40800 : zw(2,n1dfft,n2+2-i)=-zt(2,i,2*n1dfft-1)
2395 : : end do
2396 : :
2397 [ + - ]: 2400 : if(max2+1<n2-max2)then
2398 [ + + ]: 96000 : do i=max2+2,n2-max2
2399 [ + + ]: 1032000 : do j=1,n1dfft
2400 : 936000 : zw(1,j,i)=zero
2401 : 1029600 : zw(2,j,i)=zero
2402 : : end do
2403 : : end do
2404 : : end if
2405 : : end if
2406 : :
2407 : 22764 : end subroutine switchreal_cent
2408 : : !!***
2409 : :
2410 : : !----------------------------------------------------------------------
2411 : :
2412 : : !!****f* m_fftcore/scramble
2413 : : !! NAME
2414 : : !! scramble
2415 : : !!
2416 : : !! FUNCTION
2417 : : !! This routine performs the local rotation
2418 : : !!
2419 : : !! input: G1,R3,G2,(Gp2)
2420 : : !! output: G1,G2,R3,(Gp2)
2421 : : !!
2422 : : !! INPUTS
2423 : : !! i1=Index of x in the small box enclosing the G-sphere.
2424 : : !! j2
2425 : : !! lot=Cache blocking factor
2426 : : !! n1dfft=Number of 1D FFTs performed.
2427 : : !! md1,md2proc,nnd3=Used to dimension zmpi2
2428 : : !! n3=Dimension of the transform along z.
2429 : : !! zw(2,lot,n3): zw(:,1:n1dfft,n3) contains the lines transformed along z
2430 : : !!
2431 : : !! OUTPTU
2432 : : !! zmpi2(2,md1,md2proc,nnd3)
2433 : : !!
2434 : : !! SOURCE
2435 : :
2436 : 6785284 : pure subroutine scramble(i1,j2,lot,n1dfft,md1,n3,md2proc,nnd3,zw,zmpi2)
2437 : :
2438 : : !Arguments ------------------------------------
2439 : : integer,intent(in) :: i1,j2,lot,n1dfft,md1,n3,md2proc,nnd3
2440 : : real(dp),intent(in) :: zw(2,lot,n3)
2441 : : real(dp),intent(inout) :: zmpi2(2,md1,md2proc,nnd3)
2442 : :
2443 : : !Local variables-------------------------------
2444 : : integer :: i3,i
2445 : : ! *************************************************************************
2446 : :
2447 [ + + ]: 216269518 : do i3=1,n3
2448 [ + + ]: 2462909050 : do i=0,n1dfft-1
2449 : 2246639532 : zmpi2(1,i1+i,j2,i3)=zw(1,i+1,i3)
2450 : 2456123766 : zmpi2(2,i1+i,j2,i3)=zw(2,i+1,i3)
2451 : : end do
2452 : : end do
2453 : :
2454 : 6785284 : end subroutine scramble
2455 : : !!***
2456 : :
2457 : : !----------------------------------------------------------------------
2458 : :
2459 : : !!****f* m_fftcore/fill
2460 : : !! NAME
2461 : : !! fill
2462 : : !!
2463 : : !! FUNCTION
2464 : : !! Receives a set of z-lines in reciprocal space,
2465 : : !! insert the values in the cache work array zw (no padding)
2466 : : !!
2467 : : !! INPUTS
2468 : : !! nd1,nd3=Dimensions of the input array zf.
2469 : : !! lot=Cache blocking factor.
2470 : : !! n1dfft=Number of 1D FFTs to perform
2471 : : !! n3=Dimension of the transform along z
2472 : : !! zf(2,nd1,nd3)=Input array
2473 : : !!
2474 : : !! OUTPUT
2475 : : !! zw(2,lot,n3)=Cache work array with the z-lines.
2476 : : !!
2477 : : !! SOURCE
2478 : :
2479 : 1629858 : pure subroutine fill(nd1,nd3,lot,n1dfft,n3,zf,zw)
2480 : :
2481 : : !Arguments ------------------------------------
2482 : : integer,intent(in) :: nd1,nd3,lot,n1dfft,n3
2483 : : real(dp),intent(in) :: zf(2,nd1,nd3)
2484 : : real(dp),intent(inout) :: zw(2,lot,n3)
2485 : :
2486 : : ! local variables
2487 : : integer :: i1,i3
2488 : : ! *************************************************************************
2489 : :
2490 [ + + ]: 77158003 : do i3=1,n3
2491 [ + + ]: 984630300 : do i1=1,n1dfft
2492 : 907472297 : zw(1,i1,i3)=zf(1,i1,i3)
2493 : 983000442 : zw(2,i1,i3)=zf(2,i1,i3)
2494 : : end do
2495 : : end do
2496 : :
2497 : 1629858 : end subroutine fill
2498 : : !!***
2499 : :
2500 : : !----------------------------------------------------------------------
2501 : :
2502 : : !!****f* m_fftcore/fill_cent
2503 : : !! NAME
2504 : : !! fill_cent
2505 : : !!
2506 : : !! FUNCTION
2507 : : !! Receives a set of z-lines in reciprocal space,
2508 : : !! insert the values in cache work array defined on the FFT box
2509 : : !! and pads the central of the frequency region with zeros.
2510 : : !!
2511 : : !! INPUTS
2512 : : !! md1,md3=Leading dimension of zf along x and z
2513 : : !! lot=second dimension of zw (cache blocking factor)
2514 : : !! n1dfft=Number of 1d transforms to be performed along z.
2515 : : !! max3=Max G_z in the small box enclosing the G-sphere.
2516 : : !! m3=Number of points in the *small* box enclosing the G-sphere
2517 : : !! n3=Dimension of the FFT transform along z
2518 : : !! zf(2,md1,md3)=x-z planes in reciprocal space
2519 : : !!
2520 : : !! OUTPUT
2521 : : !! zw(2,lot,n3)= Filled cache work array.
2522 : : !! zw(:,1:n1dfft,n3) contains the lines to be transformed along.
2523 : : !!
2524 : : !! SOURCE
2525 : :
2526 : 5155426 : pure subroutine fill_cent(md1,md3,lot,n1dfft,max3,m3,n3,zf,zw)
2527 : :
2528 : : !Arguments ------------------------------------
2529 : : integer,intent(in) :: md1,md3,lot,n1dfft,max3,m3,n3
2530 : : real(dp),intent(in) :: zf(2,md1,md3)
2531 : : real(dp),intent(inout) :: zw(2,lot,n3)
2532 : :
2533 : : !Local variables-------------------------------
2534 : : !scalars
2535 : : integer :: i1,i3
2536 : : ! *************************************************************************
2537 : :
2538 : : ! Here, zero and positive frequencies
2539 [ + + ]: 38867280 : do i3=1,max3+1
2540 [ + + ]: 375029376 : do i1=1,n1dfft
2541 : 336162096 : zw(1,i1,i3)=zf(1,i1,i3)
2542 : 369873950 : zw(2,i1,i3)=zf(2,i1,i3)
2543 : : end do
2544 : : end do
2545 : :
2546 : : ! Fill the center region with zeros
2547 [ + + ]: 76027484 : do i3=max3+2,n3-m3+max3+1
2548 [ + + ]: 784463165 : do i1=1,n1dfft
2549 : 708435681 : zw(1,i1,i3)=zero
2550 : 779307739 : zw(2,i1,i3)=zero
2551 : : end do
2552 : : end do
2553 : :
2554 : : ! Here, negative frequencies
2555 [ + + ]: 34527603 : do i3=max3+2,m3
2556 [ + + ]: 329097061 : do i1=1,n1dfft
2557 : 294569458 : zw(1,i1,i3+n3-m3)=zf(1,i1,i3)
2558 : 323941635 : zw(2,i1,i3+n3-m3)=zf(2,i1,i3)
2559 : : end do
2560 : : end do
2561 : :
2562 : 5155426 : end subroutine fill_cent
2563 : : !!***
2564 : :
2565 : : !----------------------------------------------------------------------
2566 : :
2567 : : !!****f* m_fftcore/unfill
2568 : : !! NAME
2569 : : !! unfill
2570 : : !!
2571 : : !! FUNCTION
2572 : : !! Move data from the cache work array to zf
2573 : : !!
2574 : : !! INPUTS
2575 : : !! nd1,nd3=Dimensions of the input array zf.
2576 : : !! lot=Cache blocking factor.
2577 : : !! n1dfft=Number of 1D FFTs to perform
2578 : : !! n3=Dimension of the transform along z
2579 : : !! zw(2,lot,n3)=Cache work array with the z-lines.
2580 : : !!
2581 : : !! OUTPUT
2582 : : !! zf(2,nd1,nd3)= zf(:,1:n1dfft,:1:n3) is filled with the results stored in zw
2583 : : !!
2584 : : !! SOURCE
2585 : :
2586 : 1797673 : pure subroutine unfill(nd1,nd3,lot,n1dfft,n3,zw,zf)
2587 : :
2588 : : !Arguments ------------------------------------
2589 : : integer,intent(in) :: nd1,nd3,lot,n1dfft,n3
2590 : : real(dp),intent(in) :: zw(2,lot,n3)
2591 : : real(dp),intent(inout) :: zf(2,nd1,nd3)
2592 : :
2593 : : !Local variables-------------------------------
2594 : : integer :: i1,i3
2595 : : ! *************************************************************************
2596 : :
2597 [ + + ]: 92554993 : do i3=1,n3
2598 [ + + ]: 1165947981 : do i1=1,n1dfft
2599 : 1073392988 : zf(1,i1,i3)=zw(1,i1,i3)
2600 : 1164150308 : zf(2,i1,i3)=zw(2,i1,i3)
2601 : : end do
2602 : : end do
2603 : :
2604 : 1797673 : end subroutine unfill
2605 : : !!***
2606 : :
2607 : : !----------------------------------------------------------------------
2608 : :
2609 : : !!****f* m_fftcore/unfill_cent
2610 : : !! NAME
2611 : : !! unfill_cent
2612 : : !!
2613 : : !! FUNCTION
2614 : : !! Transfer data from the cache work array to zf.
2615 : : !! Takes into account zero padding (only the non-zero entries are moved)
2616 : : !!
2617 : : !! INPUTS
2618 : : !! md1,md3=Leading dimension of zf along x and z
2619 : : !! lot=Cache blocking factor.
2620 : : !! n1dfft=Number of 1d transforms performed along z.
2621 : : !! max3=Max index of G_z the small box enclosing the G-sphere.
2622 : : !! m3=Number of points in the *small* box enclosing the G-sphere
2623 : : !! n3=Dimension of the FFT transform along z
2624 : : !! zw(2,lot,n3)=Cache work array
2625 : : !!
2626 : : !! OUTPUT
2627 : : !! zf(2,md1,md3)= zf(:,1:n1dfft,1:m3) is filled with the non-zero components.
2628 : : !!
2629 : : !! SOURCE
2630 : :
2631 : 4737806 : pure subroutine unfill_cent(md1,md3,lot,n1dfft,max3,m3,n3,zw,zf)
2632 : :
2633 : : !Arguments ------------------------------------
2634 : : integer,intent(in) :: md1,md3,lot,n1dfft,max3,m3,n3
2635 : : real(dp),intent(in) :: zw(2,lot,n3)
2636 : : real(dp),intent(inout) :: zf(2,md1,md3)
2637 : :
2638 : : !Local variables-------------------------------
2639 : : integer :: i1,i3
2640 : : ! *************************************************************************
2641 : :
2642 : : ! Here, zero and positive frequencies
2643 [ + + ]: 35440946 : do i3=1,max3+1
2644 [ + + ]: 342557546 : do i1=1,n1dfft
2645 : 307116600 : zf(1,i1,i3)=zw(1,i1,i3)
2646 : 337819740 : zf(2,i1,i3)=zw(2,i1,i3)
2647 : : end do
2648 : : end do
2649 : :
2650 : : ! Here, negative frequencies
2651 [ + + ]: 31442682 : do i3=max3+2,m3
2652 [ + + ]: 300156296 : do i1=1,n1dfft
2653 : 268713614 : zf(1,i1,i3)=zw(1,i1,i3+n3-m3)
2654 : 295418490 : zf(2,i1,i3)=zw(2,i1,i3+n3-m3)
2655 : : end do
2656 : : end do
2657 : :
2658 : 4737806 : end subroutine unfill_cent
2659 : : !!***
2660 : :
2661 : : !----------------------------------------------------------------------
2662 : :
2663 : : !!****f* m_fftcore/unmpiswitch
2664 : : !! NAME
2665 : : !! unmpiswitch
2666 : : !!
2667 : : !! FUNCTION
2668 : : !!
2669 : : !! INPUTS
2670 : : !!
2671 : : !! OUTPUT
2672 : : !!
2673 : : !! SOURCE
2674 : :
2675 : 1800134 : pure subroutine unmpiswitch(j3,n1dfft,Jp2st,J2st,lot,n1,nd2proc,nd3proc,nproc,ioption,zw,zmpi1)
2676 : :
2677 : : !Arguments ------------------------------------
2678 : : integer,intent(in) :: j3,n1dfft,lot,n1,nd2proc,nd3proc,nproc,ioption
2679 : : integer,intent(inout) :: Jp2st,J2st
2680 : : real(dp),intent(in) :: zw(2,lot,n1)
2681 : : real(dp),intent(inout) :: zmpi1(2,n1,nd2proc,nd3proc,nproc)
2682 : :
2683 : : !Local variables-------------------------------
2684 : : integer :: i1,jp2,j2,ind,jjp2,mfft,jj2
2685 : : ! *************************************************************************
2686 : :
2687 : 1800134 : mfft=0
2688 [ + - ]: 1800134 : if (ioption == 2) then
2689 [ + + ]: 3073050 : do Jp2=Jp2st,nproc
2690 [ + + ]: 25870502 : do J2=J2st,nd2proc
2691 : 24597586 : mfft=mfft+1
2692 [ + + ]: 24597586 : if (mfft.gt.n1dfft) then
2693 : 801010 : Jp2st=Jp2
2694 : 801010 : J2st=J2
2695 : 801010 : return
2696 : : end if
2697 [ + + ]: 1098462480 : do I1=1,n1
2698 : 1073392988 : zmpi1(1,I1,J2,j3,Jp2)=zw(1,mfft,I1)
2699 : 1097189564 : zmpi1(2,I1,J2,j3,Jp2)=zw(2,mfft,I1)
2700 : : end do
2701 : : end do
2702 : 2272040 : J2st=1
2703 : : end do
2704 : :
2705 : : else
2706 [ # # ]: 0 : do Jp2=Jp2st,nproc
2707 [ # # ]: 0 : do J2=J2st,nd2proc
2708 : 0 : mfft=mfft+1
2709 [ # # ]: 0 : if (mfft.gt.n1dfft) then
2710 : 0 : Jp2st=Jp2
2711 : 0 : J2st=J2
2712 : 0 : return
2713 : : end if
2714 : 0 : ind=(Jp2-1) * nd2proc + J2
2715 : 0 : jj2=(ind-1)/nproc +1
2716 : :
2717 : : !jjp2=modulo(ind,nproc) +1
2718 : 0 : jjp2=modulo(ind-1,nproc)+1
2719 : :
2720 [ # # ]: 0 : do I1=1,n1
2721 : 0 : zmpi1(1,I1,jj2,j3,jjp2)=zw(1,mfft,I1)
2722 : 0 : zmpi1(2,I1,jj2,j3,jjp2)=zw(2,mfft,I1)
2723 : : end do
2724 : : end do
2725 : 0 : J2st=1
2726 : : end do
2727 : : end if
2728 : :
2729 : : end subroutine unmpiswitch
2730 : : !!***
2731 : :
2732 : : !----------------------------------------------------------------------
2733 : :
2734 : : !!****f* m_fftcore/unmpiswitch_cent
2735 : : !! NAME
2736 : : !! unmpiswitch_cent
2737 : : !!
2738 : : !! FUNCTION
2739 : : !!
2740 : : !! INPUTS
2741 : : !!
2742 : : !! OUTPUT
2743 : : !!
2744 : : !! SOURCE
2745 : :
2746 : 9807304 : pure subroutine unmpiswitch_cent(j3,n1dfft,Jp2stf,J2stf,lot,max1,md1,m1,n1,md2proc,nd3proc,nproc,ioption,zw,zmpi1)
2747 : :
2748 : : !Arguments ------------------------------------
2749 : : integer,intent(in) :: j3,n1dfft,lot,max1,md1,m1,n1,md2proc,nd3proc,nproc,ioption
2750 : : integer,intent(inout) :: Jp2stf,J2stf
2751 : : real(dp),intent(inout) :: zmpi1(2,md1,md2proc,nd3proc,nproc)
2752 : : real(dp),intent(in) :: zw(2,lot,n1)
2753 : :
2754 : : !Local variables-------------------------------
2755 : : integer :: mfft,Jp2,J2,I1,ind,jj2,jjp2
2756 : : ! *************************************************************************
2757 : :
2758 : 9807304 : mfft=0
2759 : :
2760 [ - + ]: 9807304 : if (ioption == 2) then
2761 [ # # ]: 0 : do Jp2=Jp2stf,nproc
2762 [ # # ]: 0 : do J2=J2stf,md2proc
2763 : 0 : mfft=mfft+1
2764 : :
2765 [ # # ]: 0 : if (mfft.gt.n1dfft) then
2766 : 0 : Jp2stf=Jp2
2767 : 0 : J2stf=J2
2768 : 0 : return
2769 : : end if
2770 : :
2771 : : ! Here, zero and positive frequencies
2772 [ # # ]: 0 : do I1=1,max1+1
2773 : 0 : zmpi1(1,I1,J2,j3,Jp2)=zw(1,mfft,I1)
2774 : 0 : zmpi1(2,I1,J2,j3,Jp2)=zw(2,mfft,I1)
2775 : : end do
2776 : :
2777 : : ! Here, negative frequencies
2778 [ # # ]: 0 : do I1=max1+2,m1
2779 : 0 : zmpi1(1,I1,J2,j3,Jp2)=zw(1,mfft,I1+n1-m1)
2780 : 0 : zmpi1(2,I1,J2,j3,Jp2)=zw(2,mfft,I1+n1-m1)
2781 : : end do
2782 : :
2783 : : end do
2784 : 0 : J2stf=1
2785 : : end do
2786 : :
2787 : : else
2788 [ + + ]: 26628020 : do Jp2=Jp2stf,nproc
2789 [ + + ]: 115466064 : do J2=J2stf,md2proc
2790 : 98645348 : mfft=mfft+1
2791 [ + + ]: 98645348 : if (mfft.gt.n1dfft) then
2792 : 711164 : Jp2stf=Jp2
2793 : 711164 : J2stf=J2
2794 : 711164 : return
2795 : : end if
2796 : 97934184 : ind=(Jp2-1) * md2proc + J2
2797 : 97934184 : jj2=(ind-1)/nproc +1
2798 : :
2799 : : !jjp2=modulo(ind,nproc) +1
2800 : 97934184 : jjp2=modulo(ind-1,nproc)+1
2801 : :
2802 : : ! Here, zero and positive frequencies
2803 [ + + ]: 763174184 : do I1=1,max1+1
2804 : 665240000 : zmpi1(1,I1,Jj2,j3,Jjp2)=zw(1,mfft,I1)
2805 : 763174184 : zmpi1(2,I1,Jj2,j3,Jjp2)=zw(2,mfft,I1)
2806 : : end do
2807 : :
2808 : : ! Here, negative frequencies
2809 [ + + ]: 672550240 : do I1=max1+2,m1
2810 : 557795340 : zmpi1(1,I1,Jj2,j3,Jjp2)=zw(1,mfft,I1+n1-m1)
2811 : 655729524 : zmpi1(2,I1,Jj2,j3,Jjp2)=zw(2,mfft,I1+n1-m1)
2812 : : end do
2813 : : end do
2814 : 25916856 : J2stf=1
2815 : : end do
2816 : : end if
2817 : :
2818 : : end subroutine unmpiswitch_cent
2819 : : !!***
2820 : :
2821 : : !----------------------------------------------------------------------
2822 : :
2823 : : !!****f* m_fftcore/unscramble
2824 : : !! NAME
2825 : : !! unscramble
2826 : : !!
2827 : : !! FUNCTION
2828 : : !!
2829 : : !! INPUTS
2830 : : !! i1
2831 : : !! j2
2832 : : !! lot=Cache blocking factor.
2833 : : !! n1dfft=Number of 1D FFTs to perform
2834 : : !! md1,n3,md2proc,nnd3
2835 : : !! zmpi2(2,md1,md2proc,nnd3)
2836 : : !!
2837 : : !! OUTPUT
2838 : : !! zw(2,lot,n3)= cache work array
2839 : : !!
2840 : : !! SOURCE
2841 : :
2842 : 6535479 : pure subroutine unscramble(i1,j2,lot,n1dfft,md1,n3,md2proc,nnd3,zmpi2,zw)
2843 : :
2844 : : !Arguments ------------------------------------
2845 : : integer,intent(in) :: i1,j2,lot,n1dfft,md1,n3,md2proc,nnd3
2846 : : real(dp),intent(in) :: zmpi2(2,md1,md2proc,nnd3)
2847 : : real(dp),intent(inout) :: zw(2,lot,n3)
2848 : :
2849 : : !Local variables-------------------------------
2850 : : !scalars
2851 : : integer :: i,i3
2852 : : ! *************************************************************************
2853 : :
2854 [ + + ]: 219227431 : do i3=1,n3
2855 [ + + ]: 2515655759 : do i=0,n1dfft-1
2856 : 2296428328 : zw(1,i+1,i3)=zmpi2(1,i1+i,j2,i3)
2857 : 2509120280 : zw(2,i+1,i3)=zmpi2(2,i1+i,j2,i3)
2858 : : end do
2859 : : end do
2860 : :
2861 : 6535479 : end subroutine unscramble
2862 : : !!***
2863 : :
2864 : : !----------------------------------------------------------------------
2865 : :
2866 : : !!****f* m_fftcore/unswitch
2867 : : !! NAME
2868 : : !! unswitch
2869 : : !!
2870 : : !! FUNCTION
2871 : : !!
2872 : : !! INPUTS
2873 : : !! n1dfft=Number of 1D FFTs
2874 : : !! n2=Dimension of the transform along y
2875 : : !! lot=Cache blocking factor.
2876 : : !! n1=Dimension of the transform along x.
2877 : : !! lzt
2878 : : !! zw(2,lot,n2)=Cache work array
2879 : : !!
2880 : : !! OUTPUT
2881 : : !! zt(2,lzt,n1)
2882 : : !!
2883 : : !! SOURCE
2884 : :
2885 : 1794212 : pure subroutine unswitch(n1dfft,n2,lot,n1,lzt,zw,zt)
2886 : :
2887 : : !Arguments ------------------------------------
2888 : : integer,intent(in) :: n1dfft,n2,lot,n1,lzt
2889 : : real(dp),intent(in) :: zw(2,lot,n2)
2890 : : real(dp),intent(inout) :: zt(2,lzt,n1)
2891 : :
2892 : : !Local variables-------------------------------
2893 : : integer :: i,j
2894 : : ! *************************************************************************
2895 : :
2896 [ + + ]: 25543064 : do j=1,n1dfft
2897 [ + + ]: 1094914564 : do i=1,n2
2898 : 1069371500 : zt(1,i,j)=zw(1,j,i)
2899 : 1093120352 : zt(2,i,j)=zw(2,j,i)
2900 : : end do
2901 : : end do
2902 : :
2903 : 1794212 : end subroutine unswitch
2904 : : !!***
2905 : :
2906 : : !----------------------------------------------------------------------
2907 : :
2908 : : !!****f* m_fftcore/unswitch_cent
2909 : : !! NAME
2910 : : !! unswitch_cent
2911 : : !!
2912 : : !! FUNCTION
2913 : : !!
2914 : : !! INPUTS
2915 : : !! n1dfft=Number of 1D FFTs to perform
2916 : : !! max2=Max G_y in the small box enclosing the G-sphere.
2917 : : !! m2=Size of the small box enclosing the G-sphere along y
2918 : : !! n2=Dimension of the transform along y
2919 : : !! lot=Cache blocking factor.
2920 : : !! n1=Dimension of the transform along x
2921 : : !! lzt
2922 : : !! zw(2,lot,n2)=Cache working array
2923 : : !!
2924 : : !! OUTPUT
2925 : : !! zt(2,lzt,n1)
2926 : : !!
2927 : : !! SOURCE
2928 : :
2929 : 10717156 : pure subroutine unswitch_cent(n1dfft,max2,m2,n2,lot,n1,lzt,zw,zt)
2930 : :
2931 : : !Arguments ------------------------------------
2932 : : integer,intent(in) :: n1dfft,max2,m2,n2,lot,n1,lzt
2933 : : real(dp),intent(in) :: zw(2,lot,n2)
2934 : : real(dp),intent(inout) :: zt(2,lzt,n1)
2935 : :
2936 : : !Local variables-------------------------------
2937 : : integer :: i,j
2938 : : ! *************************************************************************
2939 : :
2940 : : ! Here, zero and positive frequencies
2941 [ + + ]: 214802388 : do j=1,n1dfft
2942 [ + + ]: 1666611832 : do i=1,max2+1
2943 : 1451809444 : zt(1,i,j)=zw(1,j,i)
2944 : 1655894676 : zt(2,i,j)=zw(2,j,i)
2945 : : end do
2946 : : end do
2947 : :
2948 : : ! Here, negative frequencies
2949 [ + - ]: 10717156 : if(m2>=max2+2)then
2950 [ + + ]: 214802388 : do j=1,n1dfft
2951 [ + + ]: 1384447688 : do i=max2+2,m2
2952 : 1169645300 : zt(1,i,j)=zw(1,j,i+n2-m2)
2953 : 1373730532 : zt(2,i,j)=zw(2,j,i+n2-m2)
2954 : : end do
2955 : : end do
2956 : : end if
2957 : :
2958 : 10717156 : end subroutine unswitch_cent
2959 : : !!***
2960 : :
2961 : : !----------------------------------------------------------------------
2962 : :
2963 : : !!****f* m_fftcore/unswitchreal
2964 : : !! NAME
2965 : : !! unswitchreal
2966 : : !!
2967 : : !! FUNCTION
2968 : : !!
2969 : : !! INPUTS
2970 : : !! n1dfft=Number of 1D FFTs to perform
2971 : : !! n2=Dimension of the transform along y
2972 : : !! n2eff=
2973 : : !! lot=Cache blocking factor.
2974 : : !! n1zt
2975 : : !! lzt
2976 : : !! zw(2,lot,n2)=Cache working array
2977 : : !!
2978 : : !! OUTPUT
2979 : : !! zt(2,lzt,n1)
2980 : : !!
2981 : : !! SOURCE
2982 : :
2983 : 5922 : pure subroutine unswitchreal(n1dfft,n2,n2eff,lot,n1zt,lzt,zw,zt)
2984 : :
2985 : : !Arguments ------------------------------------
2986 : : integer,intent(in) :: n1dfft,n2,n2eff,lot,n1zt,lzt
2987 : : real(dp),intent(in) :: zw(2,lot,n2)
2988 : : real(dp),intent(inout) :: zt(2,lzt,n1zt)
2989 : :
2990 : : !Local variables-------------------------------
2991 : : integer :: i,j
2992 : : ! *************************************************************************
2993 : :
2994 : : ! Decompose symmetric and antisymmetric parts
2995 [ + + ]: 67320 : do j=1,n1dfft
2996 : 61398 : zt(1,1,2*j-1)=zw(1,j,1)
2997 : 61398 : zt(2,1,2*j-1)=zero
2998 : 61398 : zt(1,1,2*j) =zw(2,j,1)
2999 : 67320 : zt(2,1,2*j) =zero
3000 : : end do
3001 : :
3002 [ + + ]: 205668 : do i=2,n2eff
3003 [ + + ]: 2155014 : do j=1,n1dfft
3004 : 1949346 : zt(1,i,2*j-1)= (zw(1,j,i)+zw(1,j,n2+2-i))*half
3005 : 1949346 : zt(2,i,2*j-1)= (zw(2,j,i)-zw(2,j,n2+2-i))*half
3006 : 1949346 : zt(1,i,2*j) = (zw(2,j,i)+zw(2,j,n2+2-i))*half
3007 : 2149092 : zt(2,i,2*j) =-(zw(1,j,i)-zw(1,j,n2+2-i))*half
3008 : : end do
3009 : : end do
3010 : :
3011 : 5922 : end subroutine unswitchreal
3012 : : !!***
3013 : :
3014 : : !----------------------------------------------------------------------
3015 : :
3016 : : !!****f* m_fftcore/unswitchreal_cent
3017 : : !! NAME
3018 : : !! unswitchreal_cent
3019 : : !!
3020 : : !! FUNCTION
3021 : : !!
3022 : : !! INPUTS
3023 : : !! n1dfft=Number of 1D FFTs to perform
3024 : : !! max2=Max G_y in the small box enclosing the G-sphere.
3025 : : !! n2=Dimension of the transform along y
3026 : : !! lot=Cache blocking factor.
3027 : : !! n1zt
3028 : : !! lzt
3029 : : !! zw(2,lot,n2)=Cache working array
3030 : : !!
3031 : : !! OUTPUT
3032 : : !! zt(2,lzt,n1)
3033 : : !!
3034 : : !! SOURCE
3035 : :
3036 : 15052 : pure subroutine unswitchreal_cent(n1dfft,max2,n2,lot,n1zt,lzt,zw,zt)
3037 : :
3038 : : !Arguments ------------------------------------
3039 : : integer,intent(in) :: n1dfft,max2,n2,lot,n1zt,lzt
3040 : : real(dp),intent(in) :: zw(2,lot,n2)
3041 : : real(dp),intent(inout) :: zt(2,lzt,n1zt)
3042 : :
3043 : : !Local variables-------------------------------
3044 : : integer :: i,j
3045 : : ! *************************************************************************
3046 : :
3047 [ + + ]: 179056 : do j=1,n1dfft
3048 : 164004 : zt(1,1,2*j-1)=zw(1,j,1)
3049 : 164004 : zt(2,1,2*j-1)=zero
3050 : 164004 : zt(1,1,2*j) =zw(2,j,1)
3051 : 179056 : zt(2,1,2*j) =zero
3052 : : end do
3053 : :
3054 [ + + ]: 333112 : do i=2,max2+1
3055 [ + + ]: 3721932 : do j=1,n1dfft
3056 : 3388820 : zt(1,i,2*j-1)= (zw(1,j,i)+zw(1,j,n2+2-i))*half
3057 : 3388820 : zt(2,i,2*j-1)= (zw(2,j,i)-zw(2,j,n2+2-i))*half
3058 : 3388820 : zt(1,i,2*j) = (zw(2,j,i)+zw(2,j,n2+2-i))*half
3059 : 3706880 : zt(2,i,2*j) =-(zw(1,j,i)-zw(1,j,n2+2-i))*half
3060 : : end do
3061 : : end do
3062 : :
3063 : : ! Here, zero and positive frequencies
3064 : : ! do 90,j=1,n1dfft
3065 : : ! do 90,i=1,max2+1
3066 : : ! zt(1,i,j)=zw(1,j,i)
3067 : : ! zt(2,i,j)=zw(2,j,i)
3068 : : !90 continue
3069 : :
3070 : : ! Here, negative frequencies
3071 : : ! if(m2>=max2+2)then
3072 : : ! do 110,j=1,n1dfft
3073 : : ! do 110,i=max2+2,m2
3074 : : ! zt(1,i,j)=zw(1,j,i+n2-m2)
3075 : : ! zt(2,i,j)=zw(2,j,i+n2-m2)
3076 : : !110 continue
3077 : : ! end if
3078 : :
3079 : 15052 : end subroutine unswitchreal_cent
3080 : : !!***
3081 : :
3082 : : !----------------------------------------------------------------------
3083 : :
3084 : : !!****f* m_fftcore/mpiswitch
3085 : : !! NAME
3086 : : !! mpiswitch
3087 : : !!
3088 : : !! FUNCTION
3089 : : !!
3090 : : !! INPUTS
3091 : : !!
3092 : : !! OUTPUT
3093 : : !!
3094 : : !! SOURCE
3095 : :
3096 : 1631859 : pure subroutine mpiswitch(j3,n1dfft,Jp2st,J2st,lot,n1,nd2proc,nd3proc,nproc,ioption,zmpi1,zw)
3097 : :
3098 : : !Arguments ------------------------------------
3099 : : integer,intent(in) :: j3,n1dfft,lot,n1,nd2proc,nd3proc,nproc,ioption
3100 : : integer,intent(inout) :: Jp2st,J2st
3101 : : real(dp),intent(in) :: zmpi1(2,n1,nd2proc,nd3proc,nproc)
3102 : : real(dp),intent(inout) :: zw(2,lot,n1)
3103 : :
3104 : : !Local variables-------------------------------
3105 : : integer :: Jp2,J2,I1,ind,jj2,mfft,jjp2
3106 : : ! *************************************************************************
3107 : 1631859 : mfft=0
3108 : :
3109 [ + - ]: 1631859 : if (ioption /= 1) then
3110 [ + + ]: 2860560 : do Jp2=Jp2st,nproc
3111 [ + + ]: 23627300 : do J2=J2st,nd2proc
3112 : 22398599 : mfft=mfft+1
3113 [ + + ]: 22398599 : if (mfft.gt.n1dfft) then
3114 : 653350 : Jp2st=Jp2
3115 : 653350 : J2st=J2
3116 : 653350 : return
3117 : : end if
3118 [ + + ]: 930446247 : do I1=1,n1
3119 : 907472297 : zw(1,mfft,I1)=zmpi1(1,I1,J2,j3,Jp2)
3120 : 929217546 : zw(2,mfft,I1)=zmpi1(2,I1,J2,j3,Jp2)
3121 : : end do
3122 : : end do
3123 : 2207210 : J2st=1
3124 : : end do
3125 : :
3126 : : else
3127 [ # # ]: 0 : do Jp2=Jp2st,nproc
3128 [ # # ]: 0 : do J2=J2st,nd2proc
3129 : 0 : mfft=mfft+1
3130 [ # # ]: 0 : if (mfft.gt.n1dfft) then
3131 : 0 : Jp2st=Jp2
3132 : 0 : J2st=J2
3133 : 0 : return
3134 : : end if
3135 : 0 : ind=(Jp2-1) * nd2proc + J2
3136 : 0 : jj2=(ind-1)/nproc +1
3137 : :
3138 : : !jjp2=modulo(ind,nproc) +1
3139 : 0 : jjp2=modulo(ind-1,nproc)+1
3140 : :
3141 : : !in other words: mfft=(jj2-1)*nproc+jjp2 (modulo case)
3142 : : !instead of mfft=(Jjp2-1) * nd2proc + Jj2 (slice case)
3143 : : !with 1<=jjp2<=nproc, jj2=1,nd2proc
3144 [ # # ]: 0 : do I1=1,n1
3145 : : ! zw(1,mfft,I1)=zmpi1(1,I1,J2,j3,Jp2)
3146 : : ! zw(2,mfft,I1)=zmpi1(2,I1,J2,j3,Jp2)
3147 : 0 : zw(1,mfft,I1)=zmpi1(1,I1,jj2,j3,jjp2)
3148 : 0 : zw(2,mfft,I1)=zmpi1(2,I1,jj2,j3,jjp2)
3149 : : end do
3150 : : end do
3151 : 0 : J2st=1
3152 : : end do
3153 : : end if
3154 : :
3155 : : end subroutine mpiswitch
3156 : : !!***
3157 : :
3158 : : !----------------------------------------------------------------------
3159 : :
3160 : : !!****f* m_fftcore/mpiswitch_cent
3161 : : !! NAME
3162 : : !! mpiswitch_cent
3163 : : !!
3164 : : !! FUNCTION
3165 : : !! Perform the local rotation
3166 : : !!
3167 : : !! input: I1,J2,j3,Jp2,(jp3)
3168 : : !! output: J2,Jp2,I1,j3,(jp3)
3169 : : !!
3170 : : !! and fill the central region of the frequency spectrum with zeros
3171 : : !!
3172 : : !! INPUTS
3173 : : !!
3174 : : !! OUTPUT
3175 : : !!
3176 : : !! SOURCE
3177 : :
3178 : :
3179 : 10674807 : pure subroutine mpiswitch_cent(j3,n1dfft,Jp2stb,J2stb,lot,max1,md1,m1,n1,md2proc,&
3180 : 10674807 : nd3proc,nproc,ioption,zmpi1,zw,max2,m2,n2)
3181 : :
3182 : : !Arguments ------------------------------------
3183 : : integer,intent(in) :: j3,n1dfft,lot,max1,md1,m1,n1,md2proc,nd3proc,nproc,ioption
3184 : : integer,intent(in) :: m2,max2,n2
3185 : : integer,intent(inout) :: Jp2stb,J2stb
3186 : : real(dp),intent(in) :: zmpi1(2,md1,md2proc,nd3proc,nproc)
3187 : : real(dp),intent(inout) :: zw(2,lot,n1)
3188 : :
3189 : : !Local variables-------------------------------
3190 : : integer :: mfft,jp2,j2,jjp2,jj2,i1,ind
3191 : : ! *************************************************************************
3192 : :
3193 : : ABI_UNUSED((/m2,max2,n2/))
3194 : :
3195 : 10674807 : mfft=0
3196 : :
3197 [ - + ]: 10674807 : if (ioption /= 1) then
3198 [ # # ]: 0 : do Jp2=Jp2stb,nproc
3199 [ # # ]: 0 : do J2=J2stb,md2proc
3200 : :
3201 : 0 : mfft=mfft+1
3202 [ # # ]: 0 : if (mfft.gt.n1dfft) then
3203 : 0 : Jp2stb=Jp2
3204 : 0 : J2stb=J2
3205 : : !ABI_WARNING("Returning from mpiswithc_cent")
3206 : 0 : return
3207 : : end if
3208 : :
3209 : : ! Here, zero and positive frequencies
3210 : : ! In zmpi1, they are stored from 1 to max1+1
3211 [ # # ]: 0 : do I1=1,max1+1
3212 : 0 : zw(1,mfft,I1)=zmpi1(1,I1,J2,j3,Jp2)
3213 : 0 : zw(2,mfft,I1)=zmpi1(2,I1,J2,j3,Jp2)
3214 : : end do
3215 : :
3216 : : ! Fill the center region with zeros
3217 [ # # ]: 0 : do I1=max1+2,n1-m1+max1+1
3218 : 0 : zw(1,mfft,I1)=zero
3219 : 0 : zw(2,mfft,I1)=zero
3220 : : end do
3221 : :
3222 : : ! Here, negative frequencies
3223 : : ! In zmpi1, they are stored from 1 to m1half
3224 [ # # ]: 0 : do I1=max1+2,m1
3225 : 0 : zw(1,mfft,I1+n1-m1)=zmpi1(1,I1,J2,j3,Jp2)
3226 : 0 : zw(2,mfft,I1+n1-m1)=zmpi1(2,I1,J2,j3,Jp2)
3227 : : end do
3228 : : end do
3229 : 0 : J2stb=1
3230 : : end do
3231 : :
3232 : : else
3233 [ + + ]: 28795888 : do Jp2=Jp2stb,nproc
3234 [ + + ]: 125176340 : do J2=J2stb,md2proc
3235 : :
3236 : 107055259 : mfft=mfft+1
3237 [ + + ]: 107055259 : if (mfft.gt.n1dfft) then
3238 : 802708 : Jp2stb=Jp2
3239 : 802708 : J2stb=J2
3240 : : !ABI_WARNING("Returning from mpiswithc_cent")
3241 : 802708 : return
3242 : : end if
3243 : :
3244 : 106252551 : ind=(Jp2-1) * md2proc + J2
3245 : 106252551 : jj2=(ind-1)/nproc +1
3246 : :
3247 : : !jjp2=modulo(ind,nproc) +1
3248 : 106252551 : jjp2=modulo(ind-1,nproc)+1
3249 : :
3250 : : ! I gather consecutive I2 indexes in mfft in the modulo case
3251 : : ! Here, zero and positive frequencies
3252 : : ! In zmpi1, they are stored from 1 to max1+1
3253 [ + + ]: 832116955 : do I1=1,max1+1
3254 : 725864404 : zw(1,mfft,I1)=zmpi1(1,I1,Jj2,j3,Jjp2)
3255 : 832116955 : zw(2,mfft,I1)=zmpi1(2,I1,Jj2,j3,Jjp2)
3256 : : end do
3257 : :
3258 : : ! Fill the center region with zeros
3259 [ + + ]: 1640291443 : do I1=max1+2,n1-m1+max1+1
3260 : 1534038892 : zw(1,mfft,I1)=zero
3261 : 1640291443 : zw(2,mfft,I1)=zero
3262 : : end do
3263 : :
3264 : : ! Here, negative frequencies
3265 : : ! In zmpi1, they are stored from 1 to m1half
3266 [ + + ]: 734359953 : do I1=max1+2,m1
3267 : 609986321 : zw(1,mfft,I1+n1-m1)=zmpi1(1,I1,Jj2,j3,Jjp2)
3268 : 716238872 : zw(2,mfft,I1+n1-m1)=zmpi1(2,I1,Jj2,j3,Jjp2)
3269 : : end do
3270 : :
3271 : : end do
3272 : 27993180 : J2stb=1
3273 : : end do
3274 : : end if
3275 : :
3276 : : end subroutine mpiswitch_cent
3277 : : !!***
3278 : :
3279 : : !----------------------------------------------------------------------
3280 : :
3281 : : !!****f* m_fftcore/mpifft_fg2dbox
3282 : : !! NAME
3283 : : !! mpifft_fg2dbox
3284 : : !!
3285 : : !! FUNCTION
3286 : : !!
3287 : : !! INPUTS
3288 : : !!
3289 : : !! OUTPUT
3290 : : !!
3291 : : !! SOURCE
3292 : :
3293 : 69143 : pure subroutine mpifft_fg2dbox(nfft,ndat,fofg,n1,n2,n3,n4,nd2proc,n6,fftn2_distrib,ffti2_local,me_fft,workf)
3294 : :
3295 : : !Arguments ------------------------------------
3296 : : !scalars
3297 : : integer,intent(in) :: nfft,ndat,n1,n2,n3,n4,nd2proc,n6,me_fft
3298 : : !arrays
3299 : : integer,intent(in) :: fftn2_distrib(n2),ffti2_local(n2)
3300 : : real(dp),intent(in) :: fofg(2,nfft*ndat)
3301 : : real(dp),intent(inout) :: workf(2,n4,n6,nd2proc*ndat)
3302 : :
3303 : : !Local variables-------------------------------
3304 : : integer :: idat,i1,i2,i3,i2_local,i2_ldat,fgbase
3305 : : ! *************************************************************************
3306 : :
3307 [ + + ]: 138294 : do idat=1,ndat
3308 [ + + ]: 1321782 : do i3=1,n3
3309 [ + + ]: 32433233 : do i2=1,n2
3310 [ + + ]: 32364082 : if (fftn2_distrib(i2) == me_fft) then
3311 : 20297234 : i2_local = ffti2_local(i2)
3312 : 20297234 : i2_ldat = i2_local + (idat-1) * nd2proc
3313 : 20297234 : fgbase= n1*(i2_local-1 + nd2proc*(i3-1)) + (idat-1) * nfft
3314 [ + + ]: 832334404 : do i1=1,n1
3315 : 812037170 : workf(1,i1,i3,i2_ldat)=fofg(1,i1+fgbase)
3316 : 832334404 : workf(2,i1,i3,i2_ldat)=fofg(2,i1+fgbase)
3317 : : end do
3318 : : end if
3319 : : end do
3320 : : end do
3321 : : end do
3322 : :
3323 : 69143 : end subroutine mpifft_fg2dbox
3324 : : !!***
3325 : :
3326 : : !----------------------------------------------------------------------
3327 : :
3328 : : !!****f* m_fftcore/mpifft_fg2dbox_dpc
3329 : : !! NAME
3330 : : !! mpifft_fg2dbox_dpc
3331 : : !!
3332 : : !! FUNCTION
3333 : : !!
3334 : : !! INPUTS
3335 : : !!
3336 : : !! OUTPUT
3337 : : !!
3338 : : !! SOURCE
3339 : :
3340 : 0 : pure subroutine mpifft_fg2dbox_dpc(nfft,ndat,fofg,n1,n2,n3,n4,nd2proc,n6,fftn2_distrib,ffti2_local,me_fft,workf)
3341 : :
3342 : : !Arguments ------------------------------------
3343 : : !scalars
3344 : : integer,intent(in) :: nfft,ndat,n1,n2,n3,n4,nd2proc,n6,me_fft
3345 : : !arrays
3346 : : integer,intent(in) :: fftn2_distrib(n2),ffti2_local(n2)
3347 : : real(dp),intent(in) :: fofg(2,nfft*ndat)
3348 : : complex(dp),intent(inout) :: workf(n4,n6,nd2proc*ndat)
3349 : :
3350 : : !Local variables-------------------------------
3351 : : integer :: idat,i1,i2,i3,i2_local,i2_ldat,fgbase
3352 : : ! *************************************************************************
3353 : :
3354 [ # # ]: 0 : do idat=1,ndat
3355 [ # # ]: 0 : do i3=1,n3
3356 [ # # ]: 0 : do i2=1,n2
3357 [ # # ]: 0 : if (fftn2_distrib(i2) == me_fft) then
3358 : 0 : i2_local = ffti2_local(i2)
3359 : 0 : i2_ldat = i2_local + (idat-1) * nd2proc
3360 : 0 : fgbase= n1*(i2_local-1 + nd2proc*(i3-1)) + (idat-1) * nfft
3361 [ # # ]: 0 : do i1=1,n1
3362 : 0 : workf(i1,i3,i2_ldat)=CMPLX(fofg(1,i1+fgbase), fofg(2,i1+fgbase), kind=dp)
3363 : : end do
3364 : : end if
3365 : : end do
3366 : : end do
3367 : : end do
3368 : :
3369 : 0 : end subroutine mpifft_fg2dbox_dpc
3370 : : !!***
3371 : :
3372 : : !----------------------------------------------------------------------
3373 : :
3374 : : !!****f* m_fftcore/mpifft_dbox2fg
3375 : : !! NAME
3376 : : !! mpifft_dbox2fg
3377 : : !!
3378 : : !! FUNCTION
3379 : : !!
3380 : : !! INPUTS
3381 : : !!
3382 : : !! OUTPUT
3383 : : !!
3384 : : !! SOURCE
3385 : :
3386 : 69005 : pure subroutine mpifft_dbox2fg(n1,n2,n3,n4,nd2proc,n6,ndat,fftn2_distrib,ffti2_local,me_fft,workf,nfft,fofg)
3387 : :
3388 : : !Arguments ------------------------------------
3389 : : !scalars
3390 : : integer,intent(in) :: n1,n2,n3,n4,nd2proc,n6,ndat,me_fft,nfft
3391 : : !arrays
3392 : : integer,intent(in) :: fftn2_distrib(n2),ffti2_local(n2)
3393 : : real(dp),intent(in) :: workf(2,n4,n6,nd2proc*ndat)
3394 : : real(dp),intent(out) :: fofg(2,nfft*ndat)
3395 : :
3396 : : !Local variables-------------------------------
3397 : : integer :: idat,i1,i2,i3,i2_local,i2_ldat,fgbase
3398 : : real(dp) :: xnorm
3399 : : ! *************************************************************************
3400 : :
3401 : 69005 : xnorm=one/dble(n1*n2*n3)
3402 : :
3403 : : ! Transfer fft output to the original fft box
3404 [ + + ]: 138018 : do idat=1,ndat
3405 [ + + ]: 1373577 : do i2=1,n2
3406 [ + + ]: 1304572 : if( fftn2_distrib(i2) == me_fft) then
3407 : 961819 : i2_local = ffti2_local(i2)
3408 : 961819 : i2_ldat = i2_local + (idat-1) * nd2proc
3409 [ + + ]: 23696074 : do i3=1,n3
3410 : 22734255 : fgbase = n1*(i2_local - 1 + nd2proc*(i3-1)) + (idat - 1) * nfft
3411 [ + + ]: 1030001163 : do i1=1,n1
3412 : 1006305089 : fofg(1,i1+fgbase)=workf(1,i1,i3,i2_ldat)*xnorm
3413 : 1029039344 : fofg(2,i1+fgbase)=workf(2,i1,i3,i2_ldat)*xnorm
3414 : : end do
3415 : : end do
3416 : : end if
3417 : : end do
3418 : : end do
3419 : :
3420 : 69005 : end subroutine mpifft_dbox2fg
3421 : : !!***
3422 : :
3423 : : !----------------------------------------------------------------------
3424 : :
3425 : : !!****f* m_fftcore/mpifft_dbox2fg_dpc
3426 : : !! NAME
3427 : : !! mpifft_dbox2fg_dpc
3428 : : !!
3429 : : !! FUNCTION
3430 : : !!
3431 : : !! INPUTS
3432 : : !!
3433 : : !! OUTPUT
3434 : : !!
3435 : : !! SOURCE
3436 : :
3437 : 0 : pure subroutine mpifft_dbox2fg_dpc(n1,n2,n3,n4,nd2proc,n6,ndat,fftn2_distrib,ffti2_local,me_fft,workf,nfft,fofg)
3438 : :
3439 : : !Arguments ------------------------------------
3440 : : !scalars
3441 : : integer,intent(in) :: n1,n2,n3,n4,nd2proc,n6,ndat,me_fft,nfft
3442 : : !arrays
3443 : : integer,intent(in) :: fftn2_distrib(n2),ffti2_local(n2)
3444 : : complex(dp),intent(in) :: workf(n4,n6,nd2proc*ndat)
3445 : : real(dp),intent(out) :: fofg(2,nfft*ndat)
3446 : :
3447 : : !Local variables-------------------------------
3448 : : integer :: idat,i1,i2,i3,i2_local,i2_ldat,fgbase
3449 : : real(dp) :: xnorm
3450 : : ! *************************************************************************
3451 : :
3452 : 0 : xnorm=one/dble(n1*n2*n3)
3453 : :
3454 : : ! Transfer fft output to the original fft box
3455 [ # # ]: 0 : do idat=1,ndat
3456 [ # # ]: 0 : do i2=1,n2
3457 [ # # ]: 0 : if( fftn2_distrib(i2) == me_fft) then
3458 : 0 : i2_local = ffti2_local(i2)
3459 : 0 : i2_ldat = i2_local + (idat-1) * nd2proc
3460 [ # # ]: 0 : do i3=1,n3
3461 : 0 : fgbase = n1*(i2_local - 1 + nd2proc*(i3-1)) + (idat - 1) * nfft
3462 [ # # ]: 0 : do i1=1,n1
3463 : 0 : fofg(1,i1+fgbase)=REAL (workf(i1,i3,i2_ldat))*xnorm
3464 : 0 : fofg(2,i1+fgbase)=AIMAG(workf(i1,i3,i2_ldat))*xnorm
3465 : : end do
3466 : : end do
3467 : : end if
3468 : : end do
3469 : : end do
3470 : :
3471 : 0 : end subroutine mpifft_dbox2fg_dpc
3472 : : !!***
3473 : :
3474 : : !----------------------------------------------------------------------
3475 : :
3476 : : !!****f* m_fftcore/mpifft_dbox2fr
3477 : : !! NAME
3478 : : !! mpifft_dbox2fr
3479 : : !!
3480 : : !! FUNCTION
3481 : : !!
3482 : : !! INPUTS
3483 : : !!
3484 : : !! OUTPUT
3485 : : !!
3486 : : !! SOURCE
3487 : :
3488 : 69143 : pure subroutine mpifft_dbox2fr(n1,n2,n3,n4,n5,nd3proc,ndat,fftn3_distrib,ffti3_local,me_fft,workr,cplex,nfft,fofr)
3489 : :
3490 : : !Arguments ------------------------------------
3491 : : !scalars
3492 : : integer,intent(in) :: n1,n2,n3,n4,n5,nd3proc,ndat,me_fft,nfft,cplex
3493 : : !!arrays
3494 : : integer,intent(in) :: fftn3_distrib(n3),ffti3_local(n3)
3495 : : real(dp),intent(in) :: workr(2,n4,n5,nd3proc*ndat)
3496 : : real(dp),intent(out) :: fofr(cplex*nfft*ndat)
3497 : :
3498 : : !Local variables-------------------------------
3499 : : integer :: idat,i1,i2,i3,i3_local,i3_ldat,frbase
3500 : : ! *************************************************************************
3501 : :
3502 : : select case (cplex)
3503 : : case (1)
3504 : :
3505 [ + + ]: 138272 : do idat=1,ndat
3506 [ + + ]: 1321190 : do i3=1,n3
3507 [ + + ]: 1252056 : if( fftn3_distrib(i3) == me_fft) then
3508 : 932726 : i3_local = ffti3_local(i3)
3509 : 932726 : i3_ldat = i3_local + (idat - 1) * nd3proc
3510 [ + + ]: 21187300 : do i2=1,n2
3511 : 20254574 : frbase=n1*(i2-1+n2*(i3_local-1)) + (idat - 1) * nfft
3512 [ + + ]: 829553550 : do i1=1,n1
3513 : 828620824 : fofr(i1+frbase)=workr(1,i1,i2,i3_ldat)
3514 : : end do
3515 : : end do
3516 : : end if
3517 : : end do
3518 : : end do
3519 : :
3520 : : case (2)
3521 : :
3522 [ + + ]: 22 : do idat=1,ndat
3523 [ + + ]: 592 : do i3=1,n3
3524 [ + - ]: 583 : if (fftn3_distrib(i3) == me_fft) then
3525 : 570 : i3_local = ffti3_local(i3)
3526 : 570 : i3_ldat = i3_local + (idat - 1) * nd3proc
3527 [ + + ]: 43230 : do i2=1,n2
3528 : 42660 : frbase=2*n1*(i2-1+n2*(i3_local-1)) + (idat - 1) * cplex * nfft
3529 : : !if (frbase > cplex*nfft*ndat - 2*n1) then
3530 : : ! write(std_out,*)i2,i3_local,frbase,cplex*nfft*ndat
3531 : : ! ABI_ERROR("frbase")
3532 : : !end if
3533 [ + + ]: 3714150 : do i1=1,n1
3534 : 3670920 : fofr(2*i1-1+frbase)=workr(1,i1,i2,i3_ldat)
3535 : 3713580 : fofr(2*i1 +frbase)=workr(2,i1,i2,i3_ldat)
3536 : : end do
3537 : : end do
3538 : : end if
3539 : : end do
3540 : : end do
3541 : :
3542 : : case default
3543 : : !ABI_BUG("Wrong cplex")
3544 [ + + - - : 69143 : fofr = huge(one)
- ]
3545 : : end select
3546 : :
3547 : 69143 : end subroutine mpifft_dbox2fr
3548 : : !!***
3549 : :
3550 : : !----------------------------------------------------------------------
3551 : :
3552 : : !!****f* m_fftcore/mpifft_dbox2fr_dpc
3553 : : !! NAME
3554 : : !! mpifft_dbox2fr_dpc
3555 : : !!
3556 : : !! FUNCTION
3557 : : !!
3558 : : !! INPUTS
3559 : : !!
3560 : : !! OUTPUT
3561 : : !!
3562 : : !! SOURCE
3563 : :
3564 : 0 : pure subroutine mpifft_dbox2fr_dpc(n1,n2,n3,n4,n5,nd3proc,ndat,fftn3_distrib,ffti3_local,me_fft,workr,cplex,nfft,fofr)
3565 : :
3566 : : !Arguments ------------------------------------
3567 : : !scalars
3568 : : integer,intent(in) :: n1,n2,n3,n4,n5,nd3proc,ndat,me_fft,nfft,cplex
3569 : : !!arrays
3570 : : integer,intent(in) :: fftn3_distrib(n3),ffti3_local(n3)
3571 : : complex(dp),intent(in) :: workr(n4,n5,nd3proc*ndat)
3572 : : real(dp),intent(out) :: fofr(cplex*nfft*ndat)
3573 : :
3574 : : !Local variables-------------------------------
3575 : : integer :: idat,i1,i2,i3,i3_local,i3_ldat,frbase
3576 : : ! *************************************************************************
3577 : :
3578 : : select case (cplex)
3579 : : case (1)
3580 : :
3581 [ # # ]: 0 : do idat=1,ndat
3582 [ # # ]: 0 : do i3=1,n3
3583 [ # # ]: 0 : if( fftn3_distrib(i3) == me_fft) then
3584 : 0 : i3_local = ffti3_local(i3)
3585 : 0 : i3_ldat = i3_local + (idat - 1) * nd3proc
3586 [ # # ]: 0 : do i2=1,n2
3587 : 0 : frbase=n1*(i2-1+n2*(i3_local-1)) + (idat - 1) * nfft
3588 [ # # ]: 0 : do i1=1,n1
3589 : 0 : fofr(i1+frbase)=REAL(workr(i1,i2,i3_ldat))
3590 : : end do
3591 : : end do
3592 : : end if
3593 : : end do
3594 : : end do
3595 : :
3596 : : case (2)
3597 : :
3598 [ # # ]: 0 : do idat=1,ndat
3599 [ # # ]: 0 : do i3=1,n3
3600 [ # # ]: 0 : if (fftn3_distrib(i3) == me_fft) then
3601 : 0 : i3_local = ffti3_local(i3)
3602 : 0 : i3_ldat = i3_local + (idat - 1) * nd3proc
3603 [ # # ]: 0 : do i2=1,n2
3604 : 0 : frbase=2*n1*(i2-1+n2*(i3_local-1)) + (idat - 1) * cplex * nfft
3605 [ # # ]: 0 : do i1=1,n1
3606 : 0 : fofr(2*i1-1+frbase)=REAL (workr(i1,i2,i3_ldat))
3607 : 0 : fofr(2*i1 +frbase)=AIMAG(workr(i1,i2,i3_ldat))
3608 : : end do
3609 : : end do
3610 : : end if
3611 : : end do
3612 : : end do
3613 : :
3614 : : case default
3615 : : !ABI_BUG("Wrong cplex")
3616 [ # # # # : 0 : fofr = huge(one)
# ]
3617 : : end select
3618 : :
3619 : 0 : end subroutine mpifft_dbox2fr_dpc
3620 : : !!***
3621 : :
3622 : : !----------------------------------------------------------------------
3623 : :
3624 : : !!****f* m_fftcore/mpifft_fr2dbox
3625 : : !! NAME
3626 : : !! mpifft_fr2dbox
3627 : : !!
3628 : : !! FUNCTION
3629 : : !!
3630 : : !! INPUTS
3631 : : !!
3632 : : !! OUTPUT
3633 : : !!
3634 : : !! SOURCE
3635 : :
3636 : 69005 : pure subroutine mpifft_fr2dbox(cplex,nfft,ndat,fofr,n1,n2,n3,n4,n5,nd3proc,fftn3_distrib,ffti3_local,me_fft,workr)
3637 : :
3638 : : !Arguments ------------------------------------
3639 : : !scalars
3640 : : integer,intent(in) :: cplex,nfft,ndat,n1,n2,n3,n4,n5,nd3proc,me_fft
3641 : : !!arrays
3642 : : integer,intent(in) :: fftn3_distrib(n3),ffti3_local(n3)
3643 : : real(dp),intent(in) :: fofr(cplex*nfft*ndat)
3644 : : real(dp),intent(inout) :: workr(2,n4,n5,nd3proc*ndat)
3645 : :
3646 : : !Local variables-------------------------------
3647 : : integer :: idat,i1,i2,i3,i3_local,i3_ldat,frbase
3648 : : ! *************************************************************************
3649 : :
3650 : : select case (cplex)
3651 : : case (1)
3652 : :
3653 [ + + ]: 137996 : do idat=1,ndat
3654 [ + + ]: 1373261 : do i3=1,n3
3655 [ + + ]: 1304265 : if( me_fft == fftn3_distrib(i3) ) then
3656 : 961473 : i3_local = ffti3_local(i3)
3657 : 961473 : i3_ldat = i3_local + (idat-1) * nd3proc
3658 [ + + ]: 23653068 : do i2=1,n2
3659 : 22691595 : frbase=n1*(i2-1+n2*(i3_local-1)) + (idat-1) * nfft
3660 [ + + ]: 1026287237 : do i1=1,n1
3661 : 1002634169 : workr(1,i1,i2,i3_ldat)=fofr(i1+frbase)
3662 : 1025325764 : workr(2,i1,i2,i3_ldat)=zero
3663 : : end do
3664 : : end do
3665 : : end if
3666 : : end do
3667 : : end do
3668 : :
3669 : : case (2)
3670 : :
3671 [ + + ]: 22 : do idat=1,ndat
3672 [ + + ]: 592 : do i3=1,n3
3673 [ + - ]: 583 : if( me_fft == fftn3_distrib(i3) ) then
3674 : 570 : i3_local = ffti3_local(i3)
3675 : 570 : i3_ldat = i3_local + (idat-1) * nd3proc
3676 [ + + ]: 43230 : do i2=1,n2
3677 : 42660 : frbase=2*n1*(i2-1+n2*(i3_local-1)) + (idat-1) * cplex * nfft
3678 [ + + ]: 3714150 : do i1=1,n1
3679 : 3670920 : workr(1,i1,i2,i3_ldat)=fofr(2*i1-1+frbase)
3680 : 3713580 : workr(2,i1,i2,i3_ldat)=fofr(2*i1 +frbase)
3681 : : end do
3682 : : end do
3683 : : end if
3684 : : end do
3685 : : end do
3686 : :
3687 : : case default
3688 : : !ABI_BUG("Wrong cplex")
3689 [ + + - - : 69005 : workr = huge(one)
- - - - -
- - ]
3690 : : end select
3691 : :
3692 : 69005 : end subroutine mpifft_fr2dbox
3693 : : !!***
3694 : :
3695 : : !----------------------------------------------------------------------
3696 : :
3697 : : !!****f* m_fftcore/mpifft_fr2dbox_dpc
3698 : : !! NAME
3699 : : !! mpifft_fr2dbox_dpc
3700 : : !!
3701 : : !! FUNCTION
3702 : : !!
3703 : : !! INPUTS
3704 : : !!
3705 : : !! OUTPUT
3706 : : !!
3707 : : !! SOURCE
3708 : :
3709 : 0 : pure subroutine mpifft_fr2dbox_dpc(cplex,nfft,ndat,fofr,n1,n2,n3,n4,n5,nd3proc,fftn3_distrib,ffti3_local,me_fft,workr)
3710 : :
3711 : : !Arguments ------------------------------------
3712 : : !scalars
3713 : : integer,intent(in) :: cplex,nfft,ndat,n1,n2,n3,n4,n5,nd3proc,me_fft
3714 : : !!arrays
3715 : : integer,intent(in) :: fftn3_distrib(n3),ffti3_local(n3)
3716 : : real(dp),intent(in) :: fofr(cplex*nfft*ndat)
3717 : : complex(dp),intent(inout) :: workr(n4,n5,nd3proc*ndat)
3718 : :
3719 : : !Local variables-------------------------------
3720 : : integer :: idat,i1,i2,i3,i3_local,i3_ldat,frbase
3721 : : ! *************************************************************************
3722 : :
3723 : : select case (cplex)
3724 : : case (1)
3725 : :
3726 [ # # ]: 0 : do idat=1,ndat
3727 [ # # ]: 0 : do i3=1,n3
3728 [ # # ]: 0 : if( me_fft == fftn3_distrib(i3) ) then
3729 : 0 : i3_local = ffti3_local(i3)
3730 : 0 : i3_ldat = i3_local + (idat-1) * nd3proc
3731 [ # # ]: 0 : do i2=1,n2
3732 : 0 : frbase=n1*(i2-1+n2*(i3_local-1)) + (idat-1) * nfft
3733 [ # # ]: 0 : do i1=1,n1
3734 : 0 : workr(i1,i2,i3_ldat)=CMPLX(fofr(i1+frbase), zero, kind=dp)
3735 : : end do
3736 : : end do
3737 : : end if
3738 : : end do
3739 : : end do
3740 : :
3741 : : case (2)
3742 : :
3743 [ # # ]: 0 : do idat=1,ndat
3744 [ # # ]: 0 : do i3=1,n3
3745 [ # # ]: 0 : if( me_fft == fftn3_distrib(i3) ) then
3746 : 0 : i3_local = ffti3_local(i3)
3747 : 0 : i3_ldat = i3_local + (idat-1) * nd3proc
3748 [ # # ]: 0 : do i2=1,n2
3749 : 0 : frbase=2*n1*(i2-1+n2*(i3_local-1)) + (idat-1) * cplex * nfft
3750 [ # # ]: 0 : do i1=1,n1
3751 : 0 : workr(i1,i2,i3_ldat)=CMPLX(fofr(2*i1-1+frbase), fofr(2*i1 +frbase), kind=dp)
3752 : : end do
3753 : : end do
3754 : : end if
3755 : : end do
3756 : : end do
3757 : :
3758 : : case default
3759 : : !ABI_BUG("Wrong cplex")
3760 [ # # # # : 0 : workr = huge(one)
# # # #
# ]
3761 : : end select
3762 : :
3763 : 0 : end subroutine mpifft_fr2dbox_dpc
3764 : : !!***
3765 : :
3766 : : !----------------------------------------------------------------------
3767 : :
3768 : : !!****f* m_fftcore/indfftrisc
3769 : : !!
3770 : : !! NAME
3771 : : !! indfftrisc
3772 : : !!
3773 : : !! FUNCTION
3774 : : !! Take the data for sphere boundary and list of planewave in sphere (kg_k), manipulate them
3775 : : !! for convenient use in fourwf, and output them in indpw_k
3776 : : !!
3777 : : !! INPUTS
3778 : : !! gbound(2*mgfft+4)=sphere boundary data
3779 : : !! kg_k(3,npw_k)=reduced planewave coordinates
3780 : : !! mgfft=maximum size of 1D FFTs
3781 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
3782 : : !! npw_k=number of G vectors in basis at this k point
3783 : : !!
3784 : : !! OUTPUT
3785 : : !! indpw_k(4,npw_k)=array which gives fft box index for given basis sphere
3786 : : !! in a representation that is directly usable by sg_fftrisc.f
3787 : : !! ngb=number of FFTs along z
3788 : : !!
3789 : : !! SOURCE
3790 : :
3791 : 118554342 : subroutine indfftrisc(gbound,indpw_k,kg_k,mgfft,ngb,ngfft,npw_k)
3792 : :
3793 : : !Arguments ------------------------------------
3794 : : !scalars
3795 : : integer,intent(in) :: mgfft,npw_k
3796 : : integer,intent(out) :: ngb
3797 : : !arrays
3798 : : integer,intent(in) :: gbound(2*mgfft+4),kg_k(3,npw_k),ngfft(18)
3799 : : integer,intent(out) :: indpw_k(4,npw_k)
3800 : :
3801 : : !Local variables-------------------------------
3802 : : !scalars
3803 : : integer :: g1,g2,i1,i2,i3,igb,index,ipw,n1,n2,n3
3804 : : !arrays
3805 : 118554342 : integer,allocatable :: index2d(:,:)
3806 : : ! *************************************************************************
3807 : :
3808 : 118554342 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
3809 : :
3810 : : !First, generate a 2d index for each column of data
3811 [ - + + - : 474217368 : ABI_MALLOC(index2d,(n1,n2))
- + + - +
- + - - +
- + ]
3812 [ + + + + ]: 37198458290 : index2d(:,:)=0
3813 : 118554342 : index=1
3814 : 118554342 : igb=3
3815 [ + + ]: 622368045 : do g2=0,gbound(2) ! g2max
3816 [ + + ]: 2616828106 : do g1=0,gbound(igb+1) ! g1max
3817 : 2113014403 : index2d(g1+1,g2+1)=index
3818 : 2616828106 : index=index+1
3819 : : end do
3820 [ + + ]: 503813703 : if(gbound(igb)<=-1)then ! g1min
3821 [ + + ]: 1481350316 : do g1=gbound(igb)+n1,n1-1
3822 : 1078966992 : index2d(g1+1,g2+1)=index
3823 : 1481350316 : index=index+1
3824 : : end do
3825 : : end if
3826 : 622368045 : igb=igb+2
3827 : : end do
3828 : :
3829 [ + + ]: 118554342 : if(gbound(1)<=-1)then ! g2min
3830 [ + + ]: 511562979 : do g2=gbound(1)+n2,n2-1
3831 [ + + ]: 1417774160 : do g1=0,gbound(igb+1)
3832 : 1022012451 : index2d(g1+1,g2+1)=index
3833 : 1417774160 : index=index+1
3834 : : end do
3835 [ + + ]: 395761709 : if(gbound(igb)<=-1)then
3836 [ + + ]: 1738962267 : do g1=gbound(igb)+n1,n1-1
3837 : 1355171845 : index2d(g1+1,g2+1)=index
3838 : 1738962267 : index=index+1
3839 : : end do
3840 : : end if
3841 : 511562979 : igb=igb+2
3842 : : end do
3843 : : end if
3844 : :
3845 : 118554342 : ngb=index-1
3846 : :
3847 : :
3848 : : !The 2d index has been generated
3849 : : !Now, contract indpw_k(1,ipw) and indpw_k(2,ipw) into indpw_k(4,ipw)
3850 : : !indpw_k(1,ipw) and indpw_k(2,ipw) are used to hold inverse of index2d,
3851 : : !and for them, the second index does not fill 1:npw . It is only
3852 : : !the number of z-transform FFTs.
3853 : :
3854 : : !$OMP PARALLEL DO PRIVATE(i1,i2,i3)
3855 [ + + ]: 32485821382 : do ipw=1,npw_k
3856 [ + + ]: 32367267040 : i1=kg_k(1,ipw); if(i1<0)i1=i1+n1 ; i1=i1+1
3857 [ + + ]: 32367267040 : i2=kg_k(2,ipw); if(i2<0)i2=i2+n2 ; i2=i2+1
3858 [ + + ]: 32367267040 : i3=kg_k(3,ipw); if(i3<0)i3=i3+n3 ; i3=i3+1
3859 : 32367267040 : indpw_k(4,ipw)=index2d(i1,i2)
3860 : 32485821382 : indpw_k(3,ipw)=i3
3861 : : end do
3862 : :
3863 [ + + ]: 2061271339 : do i1=1,n1
3864 [ + + ]: 37203304362 : do i2=1,n2
3865 : 35142033023 : index=index2d(i1,i2)
3866 [ + + ]: 37084750020 : if(index/=0)then
3867 : 5569165691 : indpw_k(1,index)=i1
3868 : 5569165691 : indpw_k(2,index)=i2
3869 : : end if
3870 : : end do
3871 : : end do
3872 : :
3873 [ - + ]: 118554342 : ABI_FREE(index2d)
3874 : :
3875 : 118554342 : end subroutine indfftrisc
3876 : : !!***
3877 : :
3878 : : !----------------------------------------------------------------------
3879 : :
3880 : : !!****f* m_fftcore/kpgsph
3881 : : !! NAME
3882 : : !! kpgsph
3883 : : !!
3884 : : !! FUNCTION
3885 : : !! Use reciprocal space metric gmet(3,3) to set up the list
3886 : : !! of G vectors inside a sphere out to $ (1/2)*(2*\pi*(k+G))^2=ecut $.
3887 : : !! If mkmem=0 and mpw=0, then only count the number of planewaves
3888 : : !!
3889 : : !! INPUTS
3890 : : !! ecut=planewave kinetic energy cutoff (hartrees)
3891 : : !! exchn2n3d=if 1, n2 and n3 are exchanged
3892 : : !! gmet(3,3)=reciprocal space metric (bohr^-2)
3893 : : !! ikg=shift to be given to the location of the output data in the array kg
3894 : : !! ikpt=number of the k-point
3895 : : !! istwf_k=option parameter that describes the storage of wfs
3896 : : !! kpt(3)=reduced coords of k point (in terms of recip latt vecs)
3897 : : !! mkmem =maximum number of k points which can fit in core memory
3898 : : !! mpi_enreg=information about MPI parallelization
3899 : : !! mpw=maximum number of planewaves as dimensioned in calling routine
3900 : : !!
3901 : : !! OUTPUT
3902 : : !! kg(3,mpw*mkmem)=dimensionless coords of resulting G vecs (integer)
3903 : : !! npw=resulting number of planewaves inside ecut centered at kpt
3904 : : !!
3905 : : !! SIDE EFFECTS
3906 : : !! mpi_enreg
3907 : : !! %me_g0=if 1, the plane wave G(0 0 0) is in the set of plane waves (and is the first)
3908 : : !! TODO: other SIDE EFFECTS on mpi_enreg should be described !!!
3909 : : !!
3910 : : !! NOTES
3911 : : !! Must take into account the time-reversal symmetry when istwf_k is not 1.
3912 : : !!
3913 : : !! SOURCE
3914 : :
3915 : 1789409 : subroutine kpgsph(ecut, exchn2n3d, gmet, ikg, ikpt, istwf_k, kg, kpt, mkmem, mpi_enreg, mpw, npw)
3916 : :
3917 : : !Arguments ------------------------------------
3918 : : !scalars
3919 : : integer,intent(in) :: exchn2n3d,ikg,ikpt,istwf_k,mkmem,mpw
3920 : : integer,intent(out) :: npw
3921 : : real(dp),intent(in) :: ecut
3922 : : type(MPI_type),intent(inout) :: mpi_enreg
3923 : : !arrays
3924 : : integer,intent(inout) :: kg(3,mpw*mkmem)
3925 : : real(dp),intent(in) :: gmet(3,3),kpt(3)
3926 : :
3927 : : !Local variables-------------------------------
3928 : : !scalars
3929 : : integer :: i1,ig,ig1p,ig1pmax,ig2,ig2p,ig2pmax,ig2pmin,ig3,ig3p,ig3pmax
3930 : : integer :: ig3pmin,igtot,ii,ikpt_this_proc,in,ind,np_band,np_fft,npw_before,npw_remain,npw_split
3931 : : integer, save :: alloc_size=0
3932 : : real(dp) :: gap_pw,gmet11,gmet_trace,gmin,gs_fact,gs_part,gscut,v1,v2,v3,xx
3933 : : logical :: ipw_ok
3934 : : character(len=500) :: msg
3935 : : !arrays
3936 : : integer :: ngrid(3),nmax(3),nmin(3),n2,ierr
3937 : 1789409 : integer,allocatable :: array_ipw(:),ig1arr(:),ig2arr(:)
3938 : 1789409 : integer,allocatable :: ig3arr(:),kg_ind(:),kg_small(:,:)
3939 : 1789409 : integer, allocatable :: npw_gather(:),npw_disp(:) ,kg_ind_gather(:),kg_small_gather(:,:)
3940 : : real(dp) :: kmax(3),minor(3),numer_(3),tsec(2)
3941 : 1789409 : real(dp),allocatable :: kg1arr(:),kg2arr(:),kg3arr(:)
3942 : : ! *************************************************************************
3943 : :
3944 : : DBG_ENTER("COLL")
3945 : :
3946 : 1789409 : call timab(23,1,tsec)
3947 [ - + ]: 1789409 : if(istwf_k<1 .or. istwf_k>9)then
3948 : : write(msg,'(3a,i0,a)' )&
3949 : 0 : 'The variable istwf_k must be between 1 and 9, while',ch10,&
3950 : 0 : 'the argument of the routine istwf_k =',istwf_k,'.'
3951 : 0 : ABI_BUG(msg)
3952 : : end if
3953 : :
3954 [ - + ]: 1789409 : if(ikg+mpw>mkmem*mpw)then
3955 : : write(msg,'(5a,i0,a,i0,a,i0,4a)')&
3956 : 0 : 'The variables ikg, mkmem, and mpw must satisfy ikg<=(mkmem-1)*mpw,',ch10,&
3957 : 0 : 'while the arguments of the routine are',ch10,&
3958 : 0 : 'ikg =',ikg,', mkmem =',mkmem,', mpw =',mpw,ch10,&
3959 : 0 : 'Probable cause: Known error in invars1 for parallel spin-polarized case.',ch10,&
3960 : 0 : 'Temporary solution: Change the number of parallel processes.'
3961 : 0 : ABI_BUG(msg)
3962 : : end if
3963 : :
3964 : 1789409 : np_band=0
3965 [ + + ]: 1789409 : if (mpw>0) then
3966 [ + + ]: 918841 : np_band=1; if(mpi_enreg%paral_kgb==1) np_band=max(1,mpi_enreg%nproc_band)
3967 : 918841 : alloc_size=max(alloc_size,(mpw+1)*np_band)
3968 [ - + + - : 2756523 : ABI_MALLOC(kg_small,(3, alloc_size))
+ - - + -
+ ]
3969 [ - + + - : 2756523 : ABI_MALLOC(kg_ind,(alloc_size))
+ - - + -
+ ]
3970 [ + + ]: 327297840 : kg_ind(:)=0
3971 : : end if
3972 : :
3973 : : !A larger array, that will be split on the correct processor
3974 : : !G**2 cutoff, gscut=Ecut/2 /Pi^2
3975 : :
3976 : 1789409 : gscut=0.5_dp*ecut*piinv**2
3977 : :
3978 : : !In reduced coordinates, determine maximal value of k+G and G for each direction
3979 : :
3980 : 1789409 : minor(1)=gmet(2,2)*gmet(3,3)-gmet(2,3)**2
3981 : 1789409 : numer_(1)=gmet(1,2)**2*gmet(3,3)-2.0_dp*gmet(1,2)*gmet(1,3)*gmet(2,3) +gmet(1,3)**2*gmet(2,2)
3982 : 1789409 : minor(2)=gmet(1,1)*gmet(3,3)-gmet(1,3)**2
3983 : 1789409 : numer_(2)=gmet(2,3)**2*gmet(1,1)-2.0_dp*gmet(1,2)*gmet(1,3)*gmet(2,3) +gmet(2,1)**2*gmet(3,3)
3984 : 1789409 : minor(3)=gmet(2,2)*gmet(1,1)-gmet(1,2)**2
3985 : 1789409 : numer_(3)=gmet(3,2)**2*gmet(1,1)-2.0_dp*gmet(1,2)*gmet(1,3)*gmet(2,3) +gmet(1,3)**2*gmet(2,2)
3986 : :
3987 : : !Take the trace of the gmet tensor as dimensional reference
3988 : 1789409 : gmet_trace=gmet(1,1)+gmet(2,2)+gmet(3,3)
3989 : :
3990 [ + + ]: 7157636 : do ii=1,3
3991 : 5368227 : xx=gmet(ii,ii)*minor(ii)-numer_(ii)
3992 [ + - - + ]: 5368227 : if(xx<tol10*gmet_trace**3 .or. minor(ii)<tol10*gmet_trace**2)then
3993 : 0 : ABI_BUG('The metric tensor seem incorrect')
3994 : : end if
3995 : 5368227 : kmax(ii)=sqrt(gscut*minor(ii)/xx)
3996 [ - + ]: 5368227 : nmax(ii)=floor(kmax(ii)-kpt(ii)+tol10)
3997 [ - + ]: 5368227 : nmin(ii)=ceiling(-kmax(ii)-kpt(ii)-tol10)
3998 : 7157636 : ngrid(ii)=nmax(ii)-nmin(ii)+1
3999 : : end do
4000 : : !perform looping over fft box grid of size ngfft(1)*ngfft(2)*ngfft(3):
4001 : 1789409 : ig=0;ind=0
4002 : 1789409 : in=0
4003 : 1789409 : gmet11=gmet(1,1)
4004 : :
4005 : : !Set up standard search sequence for grid points, in standard storage mode :
4006 : : !0 1 2 3 ... nmax nmin ... -1
4007 : : !If the mode is not standard, then some part of the FFT grid must be selected
4008 : : !
4009 [ - + + - : 5368227 : ABI_MALLOC(ig1arr,(ngrid(1)))
+ - + - -
+ - + ]
4010 [ - + + - : 5368227 : ABI_MALLOC(ig2arr,(ngrid(2)))
+ - + - -
+ - + ]
4011 [ - + + - : 5368227 : ABI_MALLOC(ig3arr,(ngrid(3)))
+ - + - -
+ - + ]
4012 [ - + + - : 5368227 : ABI_MALLOC(kg1arr,(ngrid(1)))
+ - + - -
+ - + ]
4013 [ - + + - : 5368227 : ABI_MALLOC(kg2arr,(ngrid(2)))
+ - + - -
+ - + ]
4014 [ - + + - : 5368227 : ABI_MALLOC(kg3arr,(ngrid(3)))
+ - + - -
+ - + ]
4015 : :
4016 [ + + ]: 16704284 : do ig1p=1,ngrid(1)
4017 : 14914875 : ig1arr(ig1p)=ig1p-1
4018 [ + + ]: 14914875 : if (ig1p-1>nmax(1)) ig1arr(ig1p)=ig1p-ngrid(1)-1
4019 : 16704284 : kg1arr(ig1p)=kpt(1)+dble(ig1arr(ig1p))
4020 : : end do
4021 : :
4022 : : !For the second direction, the number of points might depend on istwf_k
4023 : : !---------------------------------------------------------------------
4024 : 1789409 : ig2pmax=ngrid(2)
4025 [ + + + - ]: 1789409 : if(istwf_k>=2 .and. exchn2n3d==0) ig2pmax=nmax(2)+1
4026 [ - + + - : 5368227 : ABI_MALLOC(array_ipw,(-ig2pmax:ig2pmax))
+ - + - -
+ - + ]
4027 [ + + ]: 33300436 : array_ipw(:)=0
4028 [ + + ]: 16650218 : do ig2p=1,ig2pmax
4029 : 14860809 : ig2arr(ig2p)=ig2p-1
4030 [ + + ]: 14860809 : if (ig2p-1>nmax(2)) ig2arr(ig2p)=ig2p-ngrid(2)-1
4031 : 16650218 : kg2arr(ig2p)=kpt(2)+dble(ig2arr(ig2p))
4032 : : end do
4033 : :
4034 : : !For the third direction, the number of points might depend on istwf_k
4035 : : !---------------------------------------------------------------------
4036 : 1789409 : ig3pmax=ngrid(3)
4037 [ + + - + ]: 1789409 : if (istwf_k>=2 .and. exchn2n3d==1) ig3pmax=nmax(3)+1
4038 : :
4039 [ + + ]: 16951031 : do ig3p=1,ig3pmax
4040 : 15161622 : ig3arr(ig3p)=ig3p-1
4041 [ + + ]: 15161622 : if(ig3p-1>nmax(3)) ig3arr(ig3p)=ig3p-ngrid(3)-1
4042 : 16951031 : kg3arr(ig3p)=kpt(3)+dble(ig3arr(ig3p))
4043 : : end do
4044 : :
4045 : : !Performs loop on all grid points.
4046 : : !---------------------------------------------------------------------
4047 : 1789409 : igtot = 0
4048 [ + + ]: 1789409 : if(exchn2n3d==0)then
4049 : 1789407 : mpi_enreg%me_g0=0
4050 [ + + ]: 16950995 : do ig3p=1,ngrid(3)
4051 : 15161588 : ig3=ig3arr(ig3p)
4052 : 15161588 : v3=kg3arr(ig3p)
4053 : 15161588 : ig2pmin=1
4054 [ + + + + ]: 15161588 : if( istwf_k>=2 .and. istwf_k<=5 .and. ig3<0)then
4055 : 30578 : ig2pmin=2
4056 : : end if
4057 : : ! ig2pmax was initialized previously
4058 [ + + ]: 161861250 : do ig2p=ig2pmin,ig2pmax
4059 : 144910255 : ig2=ig2arr(ig2p)
4060 : : ! PAY ATTENTION : the proc 0 must have me_g0=1
4061 : 144910255 : ipw_ok = .true.
4062 [ + + ]: 144910255 : if(mpi_enreg%paral_kgb==1 ) then
4063 : 516319 : n2 =mpi_enreg%distribfft%n2_coarse
4064 [ + + ]: 516319 : ipw_ok = ipw_ok.and.(mpi_enreg%me_fft == mpi_enreg%distribfft%tab_fftwf2_distrib( modulo(ig2,n2) + 1))
4065 : : end if
4066 [ + + ]: 144910255 : if (ig2==0 .and. ipw_ok) mpi_enreg%me_g0=1
4067 : 144910255 : v2=kg2arr(ig2p)
4068 : 144910255 : gs_part=gmet(2,2)*v2*v2+gmet(3,3)*v3*v3+2.0_dp*gmet(2,3)*v2*v3
4069 : 144910255 : gs_fact=2.0_dp*(gmet(1,2)*v2+gmet(3,1)*v3)
4070 : 144910255 : ig1pmax=ngrid(1)
4071 [ + + + + ]: 144910255 : if( (istwf_k==2.or.istwf_k==3) .and. ig3p==1 .and. ig2p==1)ig1pmax=nmax(1)+1
4072 [ + + ]: 1767679991 : do ig1p=1,ig1pmax
4073 : 1607608148 : v1=kg1arr(ig1p)
4074 : 1607608148 : gmin=gs_part+v1*(gs_fact+v1*gmet11)
4075 : : ! If inside sphere:
4076 [ + + ]: 1752518403 : if (gmin<=gscut) then
4077 [ + + ]: 625190544 : if (ipw_ok) then
4078 : 624642308 : ig=ig+1 ! inside sphere
4079 : 624642308 : igtot=igtot+1
4080 [ + + + + ]: 624642308 : if (mpw>0.and.ig<=alloc_size) then
4081 : : ! Keep coords of pw:
4082 : 313054579 : kg_small(1,ig)=ig1arr(ig1p)
4083 : 313054579 : kg_small(2,ig)=ig2
4084 : 313054579 : kg_small(3,ig)=ig3
4085 : 313054579 : kg_ind(ig)=igtot
4086 : : end if
4087 : 624642308 : array_ipw(ig2)=array_ipw(ig2)+1
4088 : : else
4089 : 548236 : igtot=igtot+1
4090 : : end if
4091 : : end if
4092 : : end do !ig1p
4093 : : end do !ig2p
4094 : : end do !ig3p
4095 : :
4096 : : else ! if (exchn2n3d/=0)
4097 : :
4098 : : ! ig2pmax was initialized previously
4099 : 2 : mpi_enreg%me_g0=0
4100 [ + + ]: 68 : do ig2p=1,ngrid(2)
4101 : 66 : ig2=ig2arr(ig2p)
4102 : : ! PAY ATTENTION : the proc 0 must have me_g0=1
4103 : 66 : ipw_ok = .true.
4104 [ + - ]: 66 : if(mpi_enreg%paral_kgb==1 ) then
4105 : 0 : n2 =mpi_enreg%distribfft%n2_coarse
4106 [ # # ]: 0 : ipw_ok = ipw_ok.and.(mpi_enreg%me_fft == mpi_enreg%distribfft%tab_fftwf2_distrib( modulo(ig2,n2) + 1))
4107 : : end if
4108 [ + + - + : 66 : if(ig2==0 .and. istwf_k>=2 .and. ipw_ok) mpi_enreg%me_g0=1
- - ]
4109 : 66 : v2 =kg2arr(ig2p)
4110 : 66 : ig3pmin=1
4111 [ + - - + : 66 : if( (istwf_k==2 .or. istwf_k==3 .or. istwf_k==6 .or. istwf_k==7) .and. ig2<0)then
- - ]
4112 : 0 : ig3pmin=2
4113 : : end if
4114 [ + + ]: 1190 : do ig3p=ig3pmin,ig3pmax
4115 : 1122 : ig3=ig3arr(ig3p)
4116 : 1122 : v3=kg3arr(ig3p)
4117 : 1122 : gs_part=gmet(2,2)*v2*v2+gmet(3,3)*v3*v3+2.0_dp*gmet(2,3)*v2*v3
4118 : 1122 : gs_fact=2.0_dp*(gmet(1,2)*v2+gmet(3,1)*v3)
4119 : 1122 : ig1pmax=ngrid(1)
4120 [ - + - - ]: 1122 : if( (istwf_k==2.or.istwf_k==3) .and. ig3p==1 .and. ig2p==1)ig1pmax=nmax(1)+1
4121 [ + + ]: 32604 : do ig1p=1,ig1pmax
4122 : 31416 : v1=kg1arr(ig1p)
4123 : 31416 : gmin=gs_part+v1*(gs_fact+v1*gmet11)
4124 : : ! If inside sphere:
4125 [ + + ]: 32538 : if (gmin<=gscut) then
4126 [ + - ]: 16078 : if (ipw_ok) then
4127 : 16078 : ig=ig+1 ! inside sphere
4128 : 16078 : igtot=igtot+1
4129 : : ! Make sure not to overrun array, or simply do not store if mpw=0
4130 [ + + - + ]: 16078 : if (mpw>0.and.ig<=alloc_size) then
4131 : : ! Keep coords of pw:
4132 : 8039 : kg_small(1,ig)=ig1arr(ig1p)
4133 : 8039 : kg_small(2,ig)=ig2
4134 : 8039 : kg_small(3,ig)=ig3
4135 : 8039 : kg_ind(ig)=igtot
4136 : : end if
4137 : : else
4138 : 0 : igtot=igtot+1
4139 : : end if
4140 : : end if
4141 : :
4142 : : end do ! ig1p
4143 : : end do ! ig3p
4144 : : ! end if ! if the ig2 plane is to be treated by this processor
4145 : : end do ! ig2p
4146 : :
4147 : : end if ! exchn2n3d==0 or ==1
4148 : :
4149 : : !Total number of G vectors at this k point is assigned: npw
4150 : : !when getcell = 1 it can be that ig exceeds mpw, the bound on kp_small
4151 : : !here is a workaround:
4152 [ + + + + ]: 1789409 : if (mpw*np_band > 0 .and. ig > alloc_size) then
4153 : 132 : npw = mpw*np_band
4154 : : else
4155 : 1789277 : npw=ig
4156 : : end if
4157 : 1789409 : alloc_size = max(alloc_size,npw)
4158 : :
4159 [ - + ]: 1789409 : ABI_FREE(ig1arr)
4160 [ - + ]: 1789409 : ABI_FREE(ig2arr)
4161 [ - + ]: 1789409 : ABI_FREE(ig3arr)
4162 [ - + ]: 1789409 : ABI_FREE(kg1arr)
4163 [ - + ]: 1789409 : ABI_FREE(kg2arr)
4164 [ - + ]: 1789409 : ABI_FREE(kg3arr)
4165 : :
4166 : : !BandFFT: plane-wave load balancing
4167 [ + + + + : 1789409 : if (mpi_enreg%paral_kgb==1.and.mpi_enreg%nproc_fft>1.and. mpi_enreg%pw_unbal_thresh>zero.and. istwf_k==1) then
+ - + + ]
4168 : : ! Check for reequilibration
4169 : 2042 : np_fft=max(1,mpi_enreg%nproc_fft)
4170 [ - + + - : 6126 : ABI_MALLOC(npw_gather,(np_fft)) ! Count pw before balancing
+ - - + -
+ ]
4171 : 2042 : call xmpi_allgather(npw,npw_gather,mpi_enreg%comm_fft,ierr)
4172 [ + + + + : 16214 : gap_pw = 100._dp*(maxval(npw_gather(:))-minval(npw_gather))/(1.*sum(npw_gather(:))/np_fft)
+ + ]
4173 : 2042 : write(msg,'(a,f5.2)' ) ' Relative gap for number of plane waves between process (%): ',gap_pw
4174 : 2042 : call wrtout(std_out,msg)
4175 [ - + ]: 2042 : if(gap_pw > mpi_enreg%pw_unbal_thresh) then ! Effective reequilibration
4176 : : write(msg,'(a,f5.2,a,i4,a,f5.2,a)') &
4177 : 0 : 'Plane-wave unbalancing (',gap_pw,'%) for kpt ',ikpt,' is higher than threshold (',&
4178 : 0 : mpi_enreg%pw_unbal_thresh,'%); a plane-wave balancing procedure is activated!'
4179 : 0 : call wrtout(std_out,msg)
4180 : : !Get optimal number
4181 [ # # ]: 0 : npw_split=sum(npw_gather(:))
4182 : 0 : npw=npw_split/np_fft
4183 : 0 : npw_remain=modulo(npw_split,np_fft)
4184 [ # # ]: 0 : if(mpi_enreg%me_fft < npw_remain) npw=npw+1
4185 : 0 : ig=npw
4186 : : !write(msg,*) 'New npw_fft = ', npw
4187 : : !call wrtout(std_out,msg)
4188 : 0 : alloc_size = max(alloc_size,npw)
4189 [ # # ]: 0 : if(mpw>0 ) then !Step for requilibration between fft process
4190 [ # # # # : 0 : ABI_MALLOC(npw_disp,(np_fft))
# # # # ]
4191 [ # # ]: 0 : npw_disp=0
4192 [ # # ]: 0 : do i1=2,np_fft
4193 : 0 : npw_disp(i1) = npw_disp(i1-1) + npw_gather(i1-1)
4194 : : end do
4195 [ # # # # : 0 : ABI_MALLOC(kg_ind_gather,(npw_split))
# # # # #
# # # ]
4196 [ # # # # : 0 : ABI_MALLOC(kg_small_gather,(3,npw_split))
# # # # #
# # # ]
4197 : : call xmpi_allgatherv(kg_ind, npw_gather(mpi_enreg%me_fft+1) , &
4198 : 0 : & kg_ind_gather,npw_gather,npw_disp,mpi_enreg%comm_fft,ierr)
4199 : : call xmpi_allgatherv(kg_small,3*npw_gather(mpi_enreg%me_fft+1), &
4200 [ # # # # : 0 : & kg_small_gather,3*npw_gather, 3*npw_disp,mpi_enreg%comm_fft,ierr)
# # ]
4201 : 0 : npw_before=mpi_enreg%me_fft*(npw_split/np_fft)+min(npw_remain,mpi_enreg%me_fft)
4202 [ # # # # ]: 0 : kg_small(:,1:npw)=kg_small_gather(:,npw_before+1:npw_before+npw)
4203 [ # # ]: 0 : kg_ind( 1:npw)=kg_ind_gather(npw_before+1:npw_before+npw)
4204 : : #ifdef DEBUG_MODE
4205 : : call wrtout(std_out,"keeping values done")
4206 : : #endif
4207 [ # # ]: 0 : ABI_FREE(npw_disp)
4208 [ # # ]: 0 : ABI_FREE(kg_ind_gather)
4209 [ # # ]: 0 : ABI_FREE(kg_small_gather)
4210 : : end if
4211 : : end if!End of reequilibration step for paral KGB
4212 [ - + ]: 4084 : ABI_FREE(npw_gather)
4213 : : end if
4214 : :
4215 : : !BandFFT: band load balancing
4216 [ + + + - ]: 1789409 : if(mpi_enreg%paral_kgb==1.and.mpi_enreg%nproc_band>0) then
4217 : 7015 : np_band=max(1,mpi_enreg%nproc_band)
4218 : 7015 : npw_split=ig;npw=npw_split/np_band
4219 : 7015 : npw_remain=modulo(npw_split,np_band)
4220 [ + + ]: 7015 : if(mpi_enreg%me_band < npw_remain) npw=npw+1
4221 [ + + ]: 7015 : if(mpw > 0) then ! This is for the case when we only compute npw and put mpw=0
4222 : 3015 : npw_before=mpi_enreg%me_band*(npw_split/np_band)+min(npw_remain,mpi_enreg%me_band)
4223 [ + - + + : 2688351 : kg_small(:,1:npw)=kg_small(:,npw_before+1:npw_before+npw)
+ + + + +
+ ]
4224 [ + - + + : 674349 : kg_ind ( 1:npw)=kg_ind ( npw_before+1:npw_before+npw)
+ + ]
4225 : : end if
4226 : : end if
4227 [ + + ]: 1789409 : if(mpw > 0) then
4228 [ + + ]: 313575713 : do i1=1,npw
4229 [ + + ]: 1250627488 : kg(:,i1+ikg)=kg_small(:,i1)
4230 [ + + ]: 313575713 : if (allocated(mpi_enreg%my_kgtab)) then
4231 : 294112 : ikpt_this_proc=mpi_enreg%my_kpttab(ikpt)
4232 : 294112 : mpi_enreg%my_kgtab(i1,ikpt_this_proc) = kg_ind(i1)
4233 : : end if
4234 : : end do
4235 [ - + ]: 918841 : ABI_FREE(kg_small)
4236 [ - + ]: 918841 : ABI_FREE(kg_ind)
4237 : : end if
4238 : :
4239 [ - + ]: 1789409 : ABI_FREE(array_ipw)
4240 : :
4241 : : !Take care of the me_g0 flag
4242 : 1789409 : mpi_enreg%me_g0_fft=mpi_enreg%me_g0
4243 [ + + + - ]: 1789409 : if(mpi_enreg%paral_kgb==1.and.mpi_enreg%nproc_band>0) then
4244 [ + + + + ]: 7015 : if(mpi_enreg%me_band==0.and.mpi_enreg%me_g0==1) then
4245 : : ! In this case, the processors had the 0 G vector before the new distribution, and still keeps it
4246 : 2340 : mpi_enreg%me_g0=1
4247 : : else
4248 : : ! All other cases
4249 : 4675 : mpi_enreg%me_g0=0
4250 : : end if
4251 : : end if
4252 : :
4253 : : !Check that npw is not zero
4254 [ + + - + ]: 1789409 : if(mpi_enreg%paral_kgb==1.and.npw==0) then
4255 : : write(msg,'(5a)' )&
4256 : 0 : 'Please decrease the number of npband*npfft MPI processes!',ch10,&
4257 : 0 : 'One of the MPI process has no plane-wave to handle.',ch10,&
4258 : 0 : 'Action: decrease npband and/or npfft.'
4259 : 0 : ABI_ERROR(msg)
4260 : : endif
4261 : :
4262 : 1789409 : call timab(23,2,tsec)
4263 : :
4264 : : DBG_EXIT("COLL")
4265 : :
4266 [ - + - + : 1789409 : end subroutine kpgsph
- + - + -
+ - + ]
4267 : : !!***
4268 : :
4269 : : !----------------------------------------------------------------------
4270 : :
4271 : : !!****f* m_fftcore/kpgcount
4272 : : !! NAME
4273 : : !! kpgcount
4274 : : !!
4275 : : !! FUNCTION
4276 : : !! Give the minimum and maximum number of G vectors in each direction:
4277 : : !! for each k-point, compute the number of G vectors in each direction,
4278 : : !! then store the min and max value over the set of k-points.
4279 : : !!
4280 : : !! INPUTS
4281 : : !! ecut=planewave kinetic energy cutoff (hartrees)
4282 : : !! exchn2n3d=if 1, n2 and n3 are exchanged
4283 : : !! gmet(3,3)=reciprocal space metric (bohr^-2)
4284 : : !! istwfk(nkpt)=option parameter that describes the storage of wfs
4285 : : !! kpt(3,nkpt)=reduced coords of k point (in terms of recip latt vecs)
4286 : : !! nkpt=number of k points
4287 : : !!
4288 : : !! OUTPUT
4289 : : !! ngmax(3)=maximum number of G vectors in each direction (x,y,z)
4290 : : !! ngmin(3)=minimum number of G vectors in each direction (x,y,z)
4291 : : !!
4292 : : !! NOTES
4293 : : !! This routine has been extracted from kpgsph...
4294 : : !!
4295 : : !! SOURCE
4296 : :
4297 : 79 : subroutine kpgcount(ecut,exchn2n3d,gmet,istwfk,kpt,ngmax,ngmin,nkpt)
4298 : :
4299 : : !Arguments ------------------------------------
4300 : : !scalars
4301 : : integer,intent(in) :: exchn2n3d,nkpt
4302 : : integer,intent(out) :: ngmax(3),ngmin(3)
4303 : : real(dp),intent(in) :: ecut
4304 : : !arrays
4305 : : integer,intent(in) :: istwfk(nkpt)
4306 : : real(dp),intent(in) :: gmet(3,3),kpt(3,nkpt)
4307 : :
4308 : : !Local variables-------------------------------
4309 : : !scalars
4310 : : integer :: ii,ikpt,istwf_k,kmax,ng1,ng2,ng3,nmin
4311 : : real(dp) :: gmet_trace,gscut,xx
4312 : : !arrays
4313 : : integer :: ngrid(3),nmax(3)
4314 : : real(dp) :: minor(3),numer_(3)
4315 : : ! *************************************************************************
4316 : :
4317 : : DBG_ENTER("COLL")
4318 : :
4319 : 79 : gscut=0.5_dp*ecut*piinv**2
4320 : 79 : gmet_trace=gmet(1,1)+gmet(2,2)+gmet(3,3)
4321 : 79 : minor(1)=gmet(2,2)*gmet(3,3)-gmet(2,3)**2
4322 : 79 : minor(2)=gmet(1,1)*gmet(3,3)-gmet(1,3)**2
4323 : 79 : minor(3)=gmet(2,2)*gmet(1,1)-gmet(1,2)**2
4324 : 79 : numer_(1)=gmet(1,2)**2*gmet(3,3)-2.0_dp*gmet(1,2)*gmet(1,3)*gmet(2,3)+gmet(1,3)**2*gmet(2,2)
4325 : 79 : numer_(2)=gmet(2,3)**2*gmet(1,1)-2.0_dp*gmet(1,2)*gmet(1,3)*gmet(2,3)+gmet(2,1)**2*gmet(3,3)
4326 : 79 : numer_(3)=gmet(3,2)**2*gmet(1,1)-2.0_dp*gmet(1,2)*gmet(1,3)*gmet(2,3)+gmet(1,3)**2*gmet(2,2)
4327 : :
4328 [ + + ]: 316 : ngmin(:)=1000000;ngmax(:)=0
4329 [ + + ]: 792 : do ikpt=1,nkpt
4330 : 713 : istwf_k=istwfk(ikpt)
4331 : :
4332 [ + + ]: 2852 : do ii=1,3
4333 : 2139 : xx=gmet(ii,ii)*minor(ii)-numer_(ii)
4334 [ + - - + ]: 2139 : if(xx<tol10*gmet_trace**3.or.minor(ii)<tol10*gmet_trace**2)then
4335 : 0 : ABI_BUG('The metric tensor seem incorrect')
4336 : : end if
4337 : 2139 : kmax=sqrt(gscut*minor(ii)/xx)
4338 [ - + ]: 2139 : nmax(ii)=floor(kmax-kpt(ii,ikpt)+tol10)
4339 [ - + ]: 2139 : nmin=ceiling(-kmax-kpt(ii,ikpt)-tol10)
4340 : 2852 : ngrid(ii)=nmax(ii)-nmin+1
4341 : : end do
4342 : :
4343 [ + + ]: 713 : ng1=ngrid(1);if(istwf_k==2.or.istwf_k==3) ng1=nmax(1)+1
4344 [ + - ]: 713 : if(exchn2n3d==0)then
4345 : 713 : ng3=ngrid(3)
4346 [ + + ]: 713 : ng2=ngrid(2);if(istwf_k>=2) ng2=nmax(2)+1
4347 [ + - ]: 4 : if(istwf_k>=2.and.istwf_k<=5) ng2=ng2-1
4348 : : else
4349 : 0 : ng2=ngrid(2)
4350 [ # # ]: 0 : ng3=ngrid(3);if(istwf_k>=2) ng3=nmax(3)+1
4351 [ # # # # ]: 0 : if(istwf_k==2.or.istwf_k==3.or.istwf_k==6.or.istwf_k==7) ng3=ng3-1
4352 : : end if
4353 : :
4354 [ + + ]: 713 : if(ng1<ngmin(1)) ngmin(1)=ng1
4355 [ + + ]: 713 : if(ng2<ngmin(2)) ngmin(2)=ng2
4356 [ + + ]: 713 : if(ng3<ngmin(3)) ngmin(3)=ng3
4357 [ + + ]: 713 : if(ng1>ngmax(1)) ngmax(1)=ng1
4358 [ + + ]: 713 : if(ng2>ngmax(2)) ngmax(2)=ng2
4359 [ + + ]: 792 : if(ng3>ngmax(3)) ngmax(3)=ng3
4360 : :
4361 : : end do
4362 : :
4363 : : DBG_EXIT("COLL")
4364 : :
4365 : 79 : end subroutine kpgcount
4366 : : !!***
4367 : :
4368 : : !----------------------------------------------------------------------
4369 : :
4370 : : !!****f* m_fftcore/get_kg
4371 : : !! NAME
4372 : : !! get_kg
4373 : : !!
4374 : : !! FUNCTION
4375 : : !! Helper function to calculate the set of G-vectors at a given kpoint.
4376 : : !! without taking advantage of FFT parallelism and G-vector distributions.
4377 : : !!
4378 : : !! INPUTS
4379 : : !! kpoint(3)=The k-point in reduced coordinates.
4380 : : !! ecut=Cutoff energy for planewave basis set.
4381 : : !! gmet(3,3)=reciprocal space metric ($\textrm{bohr}^{-2}$).
4382 : : !! istwfk=Options defining if time-reversal is used to decrease the number of G"s.
4383 : : !! [kin_sorted]=True if output g-vectors should be sorted by |k+g|^2/2. Default: False.
4384 : : !!
4385 : : !! OUTPUT
4386 : : !! npw_k=Total number of G-vectors in the full G-sphere.
4387 : : !! kg_k(3,npw_k) list of G-vectors allocated by the routine.
4388 : : !!
4389 : : !! SIDE EFFECTS
4390 : : !! [mpw]: Used to to compute the maximum number of PWs when looping over multiple k-points.
4391 : : !! [gmax(3)]: Max G-component when looping over multiple k-points.
4392 : : !!
4393 : : !! SOURCE
4394 : :
4395 : 395653 : subroutine get_kg(kpoint, istwf_k, ecut, gmet, npw_k, kg_k, &
4396 : : kin_sorted, mpw, gmax) ! optional
4397 : :
4398 : : !Arguments ------------------------------------
4399 : : !scalars
4400 : : integer,intent(in) :: istwf_k
4401 : : integer,intent(out) :: npw_k
4402 : : real(dp),intent(in) :: ecut
4403 : : integer,optional,intent(inout) :: mpw
4404 : : !arrays
4405 : : integer,allocatable,intent(out) :: kg_k(:,:)
4406 : : real(dp),intent(in) :: gmet(3,3),kpoint(3)
4407 : : logical,optional,intent(in) :: kin_sorted
4408 : : integer,optional,intent(inout) :: gmax(3)
4409 : :
4410 : : !Local variables-------------------------------
4411 : : !scalars
4412 : : integer,parameter :: mkmem_ = 1, exchn2n3d0 = 0, ikg0 = 0
4413 : : integer :: npw_k_test, ipw, ii
4414 : 395653 : type(MPI_type) :: MPI_enreg_seq
4415 : : !arrays
4416 : : integer :: kg_dum(3, 0)
4417 : 395653 : integer,allocatable :: iwork(:,:)
4418 : : ! *********************************************************************
4419 : :
4420 : 395653 : call initmpi_seq(MPI_enreg_seq)
4421 : :
4422 : : ! Calculate the number of G-vectors for this k-point.
4423 : 395653 : call kpgsph(ecut, exchn2n3d0, gmet, ikg0, 0, istwf_k, kg_dum, kpoint, 0, MPI_enreg_seq, 0, npw_k)
4424 : :
4425 : : ! Allocate and calculate the set of G-vectors.
4426 [ - + + - : 1186959 : ABI_MALLOC(kg_k,(3,npw_k))
+ - + - -
+ - + -
+ ]
4427 : 395653 : call kpgsph(ecut, exchn2n3d0, gmet, ikg0, 0, istwf_k, kg_k, kpoint, mkmem_, MPI_enreg_seq, npw_k, npw_k_test)
4428 : :
4429 : 395653 : call destroy_mpi_enreg(MPI_enreg_seq)
4430 : :
4431 [ - + ]: 395653 : if (present(kin_sorted)) then
4432 [ # # ]: 0 : if (kin_sorted) then
4433 : 0 : call sort_gvecs(npw_k, kpoint, gmet, kg_k, out_gvec=iwork)
4434 [ # # # # : 0 : kg_k = iwork
# # # # #
# # # # #
# # ]
4435 [ # # ]: 0 : ABI_FREE(iwork)
4436 : : end if
4437 : : end if
4438 : :
4439 [ + + ]: 395653 : if (present(mpw)) mpw = max(mpw, npw_k)
4440 : :
4441 [ + + ]: 395653 : if (present(gmax)) then
4442 [ + + ]: 71743650 : do ipw=1,npw_k
4443 [ + + ]: 286738737 : do ii=1,3
4444 : 286660116 : gmax(ii) = max(gmax(ii), abs(kg_k(ii,ipw)))
4445 : : end do
4446 : : end do
4447 : : end if
4448 : :
4449 [ - + - + : 395653 : end subroutine get_kg
- + - + -
+ - + - +
- + - + -
+ - + -
+ ]
4450 : : !!***
4451 : :
4452 : : !!****f* m_fftcore/kgindex
4453 : : !! NAME
4454 : : !! kgindex
4455 : : !!
4456 : : !! FUNCTION
4457 : : !! Compute the index of each plane wave on a FFT grid.
4458 : : !!
4459 : : !! INPUTS
4460 : : !! kg_k(3,npw_k)=dimensionless coords of G vecs (integer)
4461 : : !! mpi_enreg=information about MPI parallelization
4462 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
4463 : : !! npw_k=number of planewaves
4464 : : !!
4465 : : !! OUTPUT
4466 : : !! indpw_k(npw_k)=linear list number (in fft box) of given G vector for the current processor (local address)
4467 : : !! =0 if kg_k(ipw) is not treated by this processor
4468 : : !! mask(npw_k)=True if kg_k(ipw) belongs to this processor, false otherwise.
4469 : : !!
4470 : : !! NOTES
4471 : : !! mpi_enreg is not necessary in this case (the info is also in ngfft), but much more easy to read...
4472 : : !!
4473 : : !! SOURCE
4474 : :
4475 : 135738 : subroutine kgindex(indpw_k, kg_k, mask, mpi_enreg, ngfft, npw_k)
4476 : :
4477 : : !Arguments ------------------------------------
4478 : : !scalars
4479 : : integer,intent(in) :: npw_k
4480 : : type(MPI_type),intent(in) :: mpi_enreg
4481 : : !arrays
4482 : : integer,intent(in) :: kg_k(3,npw_k),ngfft(18)
4483 : : integer,intent(out) :: indpw_k(npw_k)
4484 : : logical,intent(out) :: mask(npw_k)
4485 : : !Local variables-------------------------------
4486 : : !scalars
4487 : : integer :: ig,ig1,ig2,ig3,me_fft,n1,n2,n3,nd2
4488 : : character(len=500) :: msg
4489 : : !arrays
4490 : 135738 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
4491 : : !integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
4492 : : ! *************************************************************************
4493 : :
4494 : 135738 : n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
4495 : :
4496 : : ! Use the following indexing (N means ngfft of the adequate direction)
4497 : : ! 0 1 2 3 ... N/2 -(N-1)/2 ... -1 <= kg
4498 : : ! 1 2 3 4 ....N/2+1 N/2+2 ... N <= index
4499 : :
4500 : 135738 : me_fft=mpi_enreg%me_fft
4501 : 135738 : nd2=(n2-1)/mpi_enreg%nproc_fft+1
4502 : :
4503 : 135738 : if (n2== mpi_enreg%distribfft%n2_coarse) then
4504 : 135738 : fftn2_distrib => mpi_enreg%distribfft%tab_fftdp2_distrib
4505 : 135738 : ffti2_local => mpi_enreg%distribfft%tab_fftdp2_local
4506 [ # # ]: 0 : else if (n2 == mpi_enreg%distribfft%n2_fine) then
4507 : 0 : fftn2_distrib => mpi_enreg%distribfft%tab_fftdp2dg_distrib
4508 : 0 : ffti2_local => mpi_enreg%distribfft%tab_fftdp2dg_local
4509 : : else
4510 : 0 : ABI_BUG("Unable to find an allocated distrib for this fft grid")
4511 : : end if
4512 : :
4513 : : !call ptabs_fourwf(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
4514 : :
4515 [ + + ]: 10274371 : do ig=1,npw_k
4516 : 10138633 : ig1=modulo(kg_k(1,ig),n1)
4517 : 10138633 : ig2=modulo(kg_k(2,ig),n2)
4518 : 10138633 : ig3=modulo(kg_k(3,ig),n3)
4519 [ + - ]: 10138633 : if(me_fft==fftn2_distrib(ig2+1)) then
4520 : 10138633 : ig2=ffti2_local(ig2+1) - 1
4521 : 10138633 : indpw_k(ig)=ig1+1+n1*(ig2+nd2*ig3)
4522 : 10138633 : mask(ig)=.true.
4523 : : else
4524 : 0 : indpw_k(ig)=0
4525 : 0 : mask(ig)=.false.
4526 : : end if
4527 [ + + + - : 81244802 : if (any(kg_k(:,ig) > ngfft(1:3)/2) .or. any(kg_k(:,ig) < -(ngfft(1:3)-1)/2) ) then
+ + + - -
+ ]
4528 : 0 : write(msg,'(a,3(i0,1x),a)')" The G-vector: ",kg_k(:, ig)," falls outside the FFT box. Increase boxcutmin (?)"
4529 : 0 : ABI_ERROR(msg)
4530 : : end if
4531 : : end do
4532 : :
4533 [ + - ]: 135738 : end subroutine kgindex
4534 : : !!***
4535 : :
4536 : : !----------------------------------------------------------------------
4537 : :
4538 : : !!****f* m_fftcore/addrho
4539 : : !! NAME
4540 : : !! addrho
4541 : : !!
4542 : : !! FUNCTION
4543 : : !! Add the contribution to the density generated by n1dfft x-y planes
4544 : : !!
4545 : : !! INPUTS
4546 : : !! icplexwf=1 if u(r) is real, 2 otherwise.
4547 : : !! includelast
4548 : : !! nd1,nd2=Leading dimensions of rhopart.
4549 : : !! n2=FFT dimension along y.
4550 : : !! lot=2nd Leading dimension of zw (cache blocking factor).
4551 : : !! n1dfft=Number of 1D FFTs along y performed.
4552 : : !! zw(2,lot,n2)=Array with the x-y planes (wavefunction in real space).
4553 : : !! weight=Weight factor for the density.
4554 : : !!
4555 : : !! SIDE EFFECTS
4556 : : !! rhopart(nd1,nd2)=density in the x-y plane, accumulated in output.
4557 : : !!
4558 : : !! SOURCE
4559 : :
4560 : 53744 : pure subroutine addrho(icplexwf,includelast,nd1,nd2,n2,lot,n1dfft,zw,rhopart,weight_r,weight_i)
4561 : :
4562 : :
4563 : : !Arguments ------------------------------------
4564 : : integer,intent(in) :: icplexwf,includelast,nd1,nd2,n2,lot,n1dfft
4565 : : real(dp),intent(in) :: zw(2,lot,n2)
4566 : : real(dp),intent(inout) :: rhopart(nd1,nd2)
4567 : : real(dp),intent(in) :: weight_i,weight_r
4568 : :
4569 : : !Local variables-------------------------------
4570 : : integer :: i2,j
4571 : : ! *************************************************************************
4572 : :
4573 [ + + ]: 53744 : if (icplexwf==2) then
4574 : : ! Complex wavefunction in real space.
4575 [ + - + + ]: 2326428 : do i2=1,n2-1,2
4576 [ + + ]: 25289120 : do j=1,n1dfft
4577 : 22962692 : rhopart(j,i2) = rhopart(j,i2) + weight_r*zw(1,j,i2)**2+weight_i*zw(2,j,i2)**2
4578 : 25239232 : rhopart(j,i2+1) = rhopart(j,i2+1) + weight_r*zw(1,j,i2+1)**2+weight_i*zw(2,j,i2+1)**2
4579 : : end do
4580 : : end do
4581 : :
4582 [ + + ]: 49888 : if (2*(n2/2)/=n2) then
4583 [ + + ]: 6720 : do j=1,n1dfft
4584 : 6720 : rhopart(j,n2 )=rhopart(j,n2 )+weight_r*zw(1,j,n2 )**2+weight_i*zw(2,j,n2 )**2
4585 : : end do
4586 : : end if
4587 : : else
4588 : : ! The wavefunction is real, in real space
4589 [ + + ]: 3856 : if (includelast==1) then
4590 [ + + ]: 314800 : do i2=1,n2
4591 [ + + ]: 3670448 : do j=1,n1dfft
4592 : 3355648 : rhopart(2*j-1,i2)=rhopart(2*j-1,i2)+zw(1,j,i2)**2*weight_r
4593 : 3666992 : rhopart(2*j ,i2)=rhopart(2*j ,i2)+zw(2,j,i2)**2*weight_i
4594 : : end do
4595 : : end do
4596 : : else
4597 [ + + ]: 29200 : do i2=1,n2
4598 [ + + ]: 288000 : do j=1,n1dfft-1
4599 : 259200 : rhopart(2*j-1,i2)=rhopart(2*j-1,i2)+zw(1,j,i2)**2*weight_r
4600 : 288000 : rhopart(2*j ,i2)=rhopart(2*j ,i2)+zw(2,j,i2)**2*weight_i
4601 : : end do
4602 : 29200 : rhopart(2*n1dfft-1,i2)=rhopart(2*n1dfft-1,i2)+zw(1,n1dfft,i2)**2*weight_r
4603 : : end do
4604 : : end if
4605 : :
4606 : : end if
4607 : :
4608 : 53744 : end subroutine addrho
4609 : : !!***
4610 : :
4611 : : !----------------------------------------------------------------------
4612 : :
4613 : : !!****f* m_fftcore/multpot
4614 : : !! NAME
4615 : : !! multpot
4616 : : !!
4617 : : !! FUNCTION
4618 : : !!
4619 : : !! INPUTS
4620 : : !! icplexwf=1 if u(r) is real, 2 otherwise.
4621 : : !! icplex=1 if v(r) is real, 2 otherwise.
4622 : : !! includelast
4623 : : !! nd1,nd2=Leading dimensions of pot(icplex*nd1,nd2)
4624 : : !! n2
4625 : : !! lot
4626 : : !! n1dfft
4627 : : !!
4628 : : !! OUTPUT
4629 : : !!
4630 : : !! SOURCE
4631 : :
4632 : 134610 : subroutine multpot(icplexwf,icplex,includelast,nd1,nd2,n2,lot,n1dfft,pot,zw)
4633 : :
4634 : : !Arguments ------------------------------------
4635 : : integer,intent(in) :: icplexwf,icplex,includelast,nd1,nd2,n2,lot,n1dfft
4636 : : real(dp),intent(in) :: pot(icplex*nd1,nd2)
4637 : : real(dp),intent(inout) :: zw(2,lot,n2)
4638 : :
4639 : : !Local variables-------------------------------
4640 : : integer :: i2,j
4641 : : real(dp) :: rew,imw
4642 : : ! *************************************************************************
4643 : :
4644 [ + + ]: 134610 : if (icplexwf==1) then
4645 : : ! Real u(r)
4646 : :
4647 [ - + ]: 3826 : if (icplex==2) then
4648 : 0 : ABI_BUG('multpot: icplexwf=1 and icplex=2')
4649 : : else
4650 : : ! TO BE SPEEDED UP : should use the same trick as Stefan
4651 [ + + ]: 3826 : if(includelast==1)then
4652 [ + + ]: 314050 : do i2=1,n2
4653 [ + + ]: 3658898 : do j=1,n1dfft
4654 : 3344848 : zw(1,j,i2)=zw(1,j,i2)*pot(2*j-1,i2)
4655 : 3655472 : zw(2,j,i2)=zw(2,j,i2)*pot(2*j ,i2)
4656 : : end do
4657 : : end do
4658 : : else
4659 [ + + ]: 29200 : do i2=1,n2
4660 [ + + ]: 288000 : do j=1,n1dfft-1
4661 : 259200 : zw(1,j,i2)=zw(1,j,i2)*pot(2*j-1,i2)
4662 : 288000 : zw(2,j,i2)=zw(2,j,i2)*pot(2*j ,i2)
4663 : : end do
4664 : 29200 : zw(1,n1dfft,i2)=zw(1,n1dfft,i2)*pot(2*n1dfft-1,i2)
4665 : : end do
4666 : : end if
4667 : : end if
4668 : :
4669 [ + - ]: 130784 : else if (icplexwf==2) then
4670 : : ! Complex u(r)
4671 : :
4672 [ + + ]: 130784 : if (icplex==1) then
4673 : :
4674 [ + - + + ]: 4235358 : do i2=1,n2-1,2
4675 [ + + ]: 47725952 : do j=1,n1dfft
4676 : 43490594 : zw(1,j,i2)=zw(1,j,i2)*pot(j,i2)
4677 : 43490594 : zw(2,j,i2)=zw(2,j,i2)*pot(j,i2)
4678 : 43490594 : zw(1,j,i2+1)=zw(1,j,i2+1)*pot(j,i2+1)
4679 : 47608168 : zw(2,j,i2+1)=zw(2,j,i2+1)*pot(j,i2+1)
4680 : : end do
4681 : : end do
4682 : :
4683 [ + + ]: 117784 : if (2*(n2/2)/=n2) then
4684 [ + + ]: 7200 : do j=1,n1dfft
4685 : 6750 : zw(1,j,n2)=zw(1,j,n2)*pot(j,n2)
4686 : 7200 : zw(2,j,n2)=zw(2,j,n2)*pot(j,n2)
4687 : : end do
4688 : : end if
4689 : :
4690 : : else
4691 : :
4692 [ + - + + ]: 663000 : do i2=1,n2-1,2
4693 [ + + ]: 7163000 : do j=1,n1dfft
4694 : 6500000 : rew = zw(1,j,i2); imw = zw(2,j,i2)
4695 : 6500000 : zw(1,j,i2) = rew*pot(2*j-1,i2) - imw*pot(2*j,i2)
4696 : 6500000 : zw(2,j,i2) = imw*pot(2*j-1,i2) + rew*pot(2*j,i2)
4697 : :
4698 : 6500000 : rew = zw(1,j,i2+1); imw = zw(2,j,i2+1)
4699 : 6500000 : zw(1,j,i2+1) = rew*pot(2*j-1,i2+1) - imw*pot(2*j,i2+1)
4700 : 7150000 : zw(2,j,i2+1) = imw*pot(2*j-1,i2+1) + rew*pot(2*j,i2+1)
4701 : : end do
4702 : : end do
4703 : :
4704 [ - + ]: 13000 : if (2*(n2/2)/=n2) then
4705 [ # # ]: 0 : do j=1,n1dfft
4706 : 0 : rew = zw(1,j,n2); imw = zw(2,j,n2)
4707 : 0 : zw(1,j,n2) = rew*pot(2*j-1,n2) - imw*pot(2*j,n2)
4708 : 0 : zw(2,j,n2) = imw*pot(2*j-1,n2) + rew*pot(2*j,n2)
4709 : : end do
4710 : : end if
4711 : :
4712 : : end if
4713 : : end if
4714 : :
4715 : 134610 : end subroutine multpot
4716 : : !!***
4717 : :
4718 : : !----------------------------------------------------------------------
4719 : :
4720 : : !!****f* m_fftcore/mpifft_collect_datar
4721 : : !! NAME
4722 : : !! mpifft_collect_datar
4723 : : !!
4724 : : !! FUNCTION
4725 : : !! Collect a real-space MPI-FFT distributed array on each proc.
4726 : : !!
4727 : : !! INPUTS
4728 : : !! ngfft(18)=contain all needed information about 3D FFT (see NOTES at beginning of scfcv)
4729 : : !! cplex=1 if real array, 2 for complex
4730 : : !! nfft=Number of FFT points treated by this MPI proc
4731 : : !! nspden=Second dimension of rhor
4732 : : !! rhor(cplex*nfft,nspden)=Array in real space (MPI-FFT distributed)
4733 : : !! fftn3_distrib(n3)=rank of the processors which own fft planes in 3rd dimension.
4734 : : !! fftn3_local(n3)=local i3 indices
4735 : : !! comm_fft=MPI-FFT communicator
4736 : : !! [master]=MPI rank, Optional. If present, the global array is available only on master node.
4737 : : !!
4738 : : !! OUTPUT
4739 : : !! rhor_glob(cplex*nfft_tot,nspden)=Global array
4740 : : !!
4741 : : !! SOURCE
4742 : :
4743 [ # # ]: 0 : subroutine mpifft_collect_datar(ngfft,cplex,nfft,nspden,rhor,comm_fft,fftn3_distrib,ffti3_local,rhor_glob,master)
4744 : :
4745 : : !Arguments ------------------------------------
4746 : : !scalars
4747 : : integer,intent(in) :: cplex,nfft,nspden,comm_fft
4748 : : integer,optional,intent(in) :: master
4749 : : !arrays
4750 : : integer,intent(in) :: ngfft(18),fftn3_distrib(ngfft(3)),ffti3_local(ngfft(3))
4751 : : real(dp),intent(in) :: rhor(cplex*nfft,nspden)
4752 : : real(dp),intent(out) :: rhor_glob(cplex*product(ngfft(1:3)),nspden)
4753 : :
4754 : : !Local variables-------------------------------
4755 : : integer :: ispden,i1,i2,i3,me_fft,i3_local,my_fftbase,glob_fftbase
4756 : : integer :: n1,n2,n3,ierr,nfft_tot
4757 : : ! *************************************************************************
4758 : :
4759 [ # # ]: 0 : nfft_tot = product(ngfft(1:3)); me_fft = xmpi_comm_rank(comm_fft)
4760 : 0 : n1 = ngfft(1); n2 = ngfft(2); n3 = ngfft(3)
4761 : :
4762 [ # # ]: 0 : if (nfft_tot == nfft) then
4763 : : ! full rhor on each node, just do a copy
4764 [ # # # # ]: 0 : rhor_glob = rhor
4765 : : else
4766 : : ! if MPI-FFT we have to gather the full rhor on each node.
4767 [ # # # # ]: 0 : rhor_glob = zero
4768 [ # # ]: 0 : do ispden=1,nspden
4769 [ # # ]: 0 : do i3=1,n3
4770 [ # # ]: 0 : if (me_fft == fftn3_distrib(i3)) then
4771 : 0 : i3_local = ffti3_local(i3)
4772 [ # # ]: 0 : do i2=1,n2
4773 : 0 : my_fftbase = cplex * ( (i2-1)*n1 + (i3_local-1)*n1*n2 )
4774 : 0 : glob_fftbase = cplex * ( (i2-1)*n1 + (i3-1)*n1*n2 )
4775 [ # # ]: 0 : do i1=1,cplex * n1
4776 : 0 : rhor_glob(i1+glob_fftbase,ispden) = rhor(i1+my_fftbase,ispden)
4777 : : end do
4778 : : end do
4779 : : end if
4780 : : end do
4781 : : end do
4782 [ # # ]: 0 : if (present(master)) then
4783 : 0 : call xmpi_sum_master(rhor_glob,master,comm_fft,ierr)
4784 : : else
4785 : 0 : call xmpi_sum(rhor_glob,comm_fft,ierr)
4786 : : end if
4787 : : end if
4788 : :
4789 : 0 : end subroutine mpifft_collect_datar
4790 : : !!***
4791 : :
4792 : : END MODULE m_fftcore
4793 : : !!***
|
