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 57242 : integer function fftcore_set_mixprec(wp) result(old_wp)
149 :
150 : !Arguments ------------------------------------
151 : !scalars
152 : integer,intent(in) :: wp
153 : ! *************************************************************************
154 :
155 57242 : old_wp = fftcore_mixprec
156 57242 : fftcore_mixprec = abs(wp)
157 :
158 57194 : select case (fftcore_mixprec)
159 : case (0)
160 57194 : 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 57242 : ABI_ERROR(sjoin("Wrong value for input wp:", itoa(fftcore_mixprec)))
165 : end select
166 :
167 57242 : 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 256 : 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 256 : ans = .False.
237 256 : fftalga = fftalg/100; fftalgb = mod(fftalg,100)/10; fftalgc = mod(fftalg,10)
238 :
239 256 : if (fftalga == FFT_FFTW3) ans = .True.
240 : !if (fftalga == FFT_DFTI) ans = .True.
241 256 : if (fftalga == FFT_SG2002) ans = .True.
242 :
243 256 : 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 22303 : 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 22303 : 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 22303 : fftalg = 512
287 22303 : 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 22303 : 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 8501 : pure function get_cache_kb()
357 :
358 : !Local variables-------------------------------
359 : !scalars
360 : integer :: get_cache_kb
361 : ! *************************************************************************
362 :
363 : ! Default value
364 154715 : get_cache_kb = 16
365 : !get_cache_kb = 32
366 : !get_cache_kb = 256
367 :
368 8501 : 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 146214 : 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 146214 : fftalg = 512
405 : #endif
406 :
407 584856 : ngfft(1:3) = n123
408 146214 : ngfft(4) = 2*(ngfft(1)/2)+1
409 146214 : ngfft(5) = 2*(ngfft(2)/2)+1
410 146214 : ngfft(6) = ngfft(3)
411 146214 : ngfft(7)= fftalg ! fftalg
412 146214 : ngfft(8)= get_cache_kb() ! cache_kb
413 146214 : ngfft(9)= 0 ! paral_fft_
414 146214 : ngfft(10)=1 ! nproc_fft
415 146214 : ngfft(11)=0 ! me_fft
416 146214 : ngfft(12)=0 ! n2proc
417 146214 : ngfft(13)=0 ! n3proc
418 877284 : ngfft(14:18)=0 ! not used
419 :
420 146214 : 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 20396 : 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 10198 : my_prtvol=0; if (PRESENT(prtvol)) my_prtvol=prtvol
459 :
460 10198 : msg = ch10//' ==== FFT mesh description (ngfft) ==== '
461 10198 : if (PRESENT(header)) msg=ch10//' ==== '//TRIM(ADJUSTL(header))//' ==== '
462 10198 : call wrtout(units, msg)
463 :
464 : write(msg,'(2(a,3i5,a),a,i5,2a,i5)')&
465 10198 : '- FFT mesh divisions ........................ ',ngfft(1),ngfft(2),ngfft(3),ch10,&
466 10198 : '- Augmented FFT divisions ................... ',ngfft(4),ngfft(5),ngfft(6),ch10,&
467 10198 : '- FFT algorithm ............................. ',ngfft(7),ch10,&
468 20396 : '- FFT cache size ............................ ',ngfft(8)
469 10198 : call wrtout(units, msg)
470 :
471 10198 : 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 10198 : 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 260098 : 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 260098 : plane=0
533 :
534 : !look at +/- g1 planes:
535 260098 : dsqmax=zero
536 260098 : dsqmin=dsq(ngfft(1)/2,-ngfft(2)/2,-ngfft(3)/2,gmet,kpt)+0.01_dp
537 3602984 : do i2=-ngfft(2)/2,ngfft(2)/2
538 81043142 : do i3=-ngfft(3)/2,ngfft(3)/2
539 77440158 : dsp = dsq(ngfft(1)/2, i2, i3,gmet,kpt)
540 77440158 : dsm = dsq( - ngfft(1)/2, i2, i3,gmet,kpt)
541 77440158 : if (dsp>dsqmax) dsqmax = dsp
542 77440158 : if (dsm>dsqmax) dsqmax = dsm
543 77440158 : if (dsp<dsqmin) then
544 4290521 : dsqmin = dsp
545 4290521 : i1min = ngfft(1)/2
546 4290521 : i2min = i2
547 4290521 : i3min = i3
548 4290521 : plane=1
549 : end if
550 80783044 : if (dsm<dsqmin) then
551 1177269 : dsqmin = dsm
552 1177269 : i1min = - ngfft(1)/2
553 1177269 : i2min = i2
554 1177269 : i3min = i3
555 1177269 : plane=1
556 : end if
557 : end do
558 : end do
559 : !
560 : !+/- g2 planes:
561 3851272 : do i1=-ngfft(1)/2,ngfft(1)/2
562 85565238 : do i3=-ngfft(3)/2,ngfft(3)/2
563 81713966 : dsp = dsq(i1,ngfft(2)/2,i3,gmet,kpt)
564 81713966 : dsm = dsq(i1,-ngfft(2)/2,i3,gmet,kpt)
565 81713966 : if (dsp>dsqmax) dsqmax = dsp
566 81713966 : if (dsm>dsqmax) dsqmax = dsm
567 81713966 : if (dsp<dsqmin) then
568 375291 : dsqmin = dsp
569 375291 : i1min = i1
570 375291 : i2min = ngfft(2)/2
571 375291 : i3min = i3
572 375291 : plane=2
573 : end if
574 85305140 : if (dsm<dsqmin) then
575 208842 : dsqmin = dsm
576 208842 : i1min = i1
577 208842 : i2min = - ngfft(2)/2
578 208842 : i3min = i3
579 208842 : plane=2
580 : end if
581 : end do
582 : end do
583 : !
584 : !+/- g3 planes:
585 3851272 : do i1=-ngfft(1)/2,ngfft(1)/2
586 83399162 : do i2=-ngfft(2)/2,ngfft(2)/2
587 79547890 : dsp = dsq(i1,i2,ngfft(3)/2,gmet,kpt)
588 79547890 : dsm = dsq(i1,i2,-ngfft(3)/2,gmet,kpt)
589 79547890 : if (dsp>dsqmax) dsqmax = dsp
590 79547890 : if (dsm>dsqmax) dsqmax = dsm
591 79547890 : if (dsp<dsqmin) then
592 294650 : dsqmin = dsp
593 294650 : i1min = i1
594 294650 : i2min = i2
595 294650 : i3min = ngfft(3)/2
596 294650 : plane=3
597 : end if
598 83139064 : if (dsm<dsqmin) then
599 141909 : dsqmin = dsm
600 141909 : i1min = i1
601 141909 : i2min = i2
602 141909 : i3min = - ngfft(3)/2
603 141909 : plane=3
604 : end if
605 : end do
606 : end do
607 :
608 260098 : 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 260098 : gbound(1)=i1min
620 260098 : gbound(2)=i2min
621 260098 : gbound(3)=i3min
622 :
623 : contains
624 :
625 477664126 : 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 477664126 : & +gmet(3,1)*(kpt(3)+dble(i3))*(kpt(1)+dble(i1)))
637 477664126 : 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 8392 : subroutine getng(boxcutmin, chksymtnons, ecut, gmet, kpt, me_fft, mgfft, nfft, ngfft, &
696 8392 : 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 8392 : ount = std_out; if (present(unit)) ount = unit
749 :
750 : !write(std_out,*)' m_fftcore/getng : enter'
751 :
752 8392 : 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 8392 : paral_fft_=paral_fft;if (nproc_fft==0) paral_fft_=0
760 :
761 8392 : if(first==1.or.paral_fft_ /= previous_paral_mode) then
762 2114 : first=0; previous_paral_mode=paral_fft_
763 2114 : srch(:,:)=0
764 :
765 : ! Factors of 2
766 2114 : srch(1,1)=1
767 35938 : do ii=1,maxpow2
768 35938 : srch(ii+1,1)=srch(ii,1)*2
769 : end do
770 :
771 : ! Factors of 3
772 14798 : index=maxpow2+1
773 : if(maxpow3>0)then
774 14798 : do ii=1,max(1,maxpow3)
775 446054 : 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 14798 : index=(maxpow3+1)*index
781 : if(maxpow5>0)then
782 14798 : do ii=1,max(1,maxpow5)
783 12684 : li_srch = 0
784 1522080 : li_srch(1+ii*index:(ii+1)*index)=5_i8b*srch(1+(ii-1)*index:ii*index,1)
785 10578456 : where (li_srch > huge(maxpow3)) li_srch = huge(maxpow3)
786 1524194 : 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 2114 : 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 2114 : call sort_int(mmsrch,srch(:,1),iperm)
807 :
808 566552 : do ii=1,mmsrch
809 566552 : if(srch(ii,1)>largest_ngfft)exit
810 : end do
811 2114 : 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 2114 : if(paral_fft_==1)then
815 1447 : msrch(2)=0
816 387796 : do ii=1,msrch(1)
817 387796 : if(modulo(srch(ii,1),nproc_fft)==0) then
818 384129 : msrch(2)=msrch(2)+1
819 384129 : 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 667 : msrch(2)=msrch(1)
828 556278 : srch(:,2)=srch(:,1)
829 : end if
830 :
831 : ! The second and third search list have the same constraint
832 2114 : msrch(3)=msrch(2)
833 1769354 : 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 33568 : ngsav(1:3) = ngfft(1:3)
845 :
846 : !As an initial guess for ngfft, use the provided coarse mesh grid
847 8392 : if (PRESENT(ngfftc)) then
848 4720 : ngfft(1:3)=ngfftc(1:3)
849 1180 : call wrtout(ount,' Using supplied coarse mesh as initial guess.')
850 : else
851 28848 : ngfft(1:3)=2
852 : end if
853 :
854 : !Enlarge the initial guess until the set of ngfft entirely comprises the sphere
855 240728 : do
856 :
857 240728 : 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 240728 : if (dsqmin>=(half*boxcutmin**2*ecut/pi**2)) exit
861 :
862 : ! Fix nearest boundary
863 1902558 : do ii=1,msrch(plane)-1
864 1894166 : if (srch(ii,plane)>=ngfft(plane)) then
865 : ! redefine ngfft(plane) to next higher choice
866 232336 : ngfft(plane)=srch(ii+1,plane)
867 232336 : exit ! Exit the loop over ii
868 : end if
869 :
870 1661830 : 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 8392 : 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 8392 : ' For input ecut=',ecut,' best grid ngfft=',ngfft(1:3),ch10,&
888 16784 : ' max ecut=',ecutmx
889 8392 : 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 8392 : divisor(:,:)=0; tobechecked=0
905 :
906 33568 : do ii=1,3
907 109096 : do jj=1,3
908 75528 : if(jj==ii)cycle
909 1249710 : do isym=1,nsym
910 1174182 : if(symrel(jj,ii,isym)==0 .or. divisor(jj,ii)==1 )cycle
911 30701 : tobechecked=1
912 30701 : element=abs(symrel(jj,ii,isym))
913 30701 : testdiv= ( divisor(jj,ii)==0 .or. divisor(jj,ii)==element .or. element==1)
914 75528 : if(testdiv)then
915 30701 : 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 1143481 : 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 8392 : testok=1
931 8392 : if(tobechecked==1)then
932 24512 : do ii=1,3
933 73536 : do jj=1,3
934 55152 : xx=divisor(jj,ii)*ngfft(jj)
935 55152 : yy=xx/ngfft(ii)
936 73536 : if(abs(yy-nint(yy))>tol8)testok=0
937 : end do
938 24512 : 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 6128 : 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 8392 : 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 586 : testok=0
1130 2344 : do mu=1,3
1131 23339 : do ii=1,msrch(mu)
1132 22753 : if(srch(ii,mu)==ngsav(mu))then
1133 1758 : testok=testok+1
1134 1758 : exit
1135 : end if
1136 : end do
1137 : end do
1138 586 : if(testok==3)then
1139 586 : write(msg,'(a,3(a,i1,a,i3),a)') ' input values of',&
1140 2930 : (' ngfft(',mu,') =',ngsav(mu),mu=1,3),' are alright and will be used'
1141 586 : call wrtout(ount,msg)
1142 2344 : do mu = 1,3
1143 2344 : 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 33568 : mgfft = maxval(ngfft(1:3))
1151 :
1152 8392 : if (paral_fft_==1) then
1153 : ! For the time being, one need ngfft(2) and ngfft(3) to be multiple of nproc_fft
1154 6677 : 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 6677 : 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 1715 : 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 33568 : 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 8392 : ngfft(4)=2*(ngfft(1)/2)+1
1179 8392 : ngfft(5)=2*(ngfft(2)/2)+1
1180 8392 : ngfft(6)=ngfft(3)
1181 8392 : 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 8392 : if (present(gpu_option)) then
1193 8384 : 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 50352 : ngfft(14:18)=0 ! ngfft(14) to be filled outside of getng
1201 :
1202 8392 : if (paral_fft_==0) then
1203 1715 : ngfft(9)=0 ! paral_fft_
1204 1715 : ngfft(10)=1 ! nproc_fft
1205 1715 : ngfft(11)=0 ! me_fft
1206 1715 : ngfft(12)=0 ! n2proc
1207 1715 : ngfft(13)=0 ! n3proc
1208 : else
1209 6677 : ngfft(9)=1 ! paral_fft_
1210 6677 : ngfft(10)=nproc_fft
1211 6677 : ngfft(11)=me_fft
1212 6677 : ngfft(12)=ngfft(2)/nproc_fft
1213 6677 : ngfft(13)=ngfft(3)/nproc_fft
1214 : end if
1215 :
1216 16784 : call print_ngfft([ount], ngfft, header="FFT mesh")
1217 :
1218 8392 : 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 3264335 : 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 609460001 : gbound(1,1)=minval(kg_k(1,:))
1287 609460001 : gbound(2,1)=maxval(kg_k(1,:))
1288 9793005 : gbound(1:2,2)=gbound(1:2,1)
1289 :
1290 : !Treat differently the fftalgc cases
1291 9793005 : do ii=1,2
1292 :
1293 6528670 : fftalgc=3-ii
1294 :
1295 6528670 : if(fftalgc/=2)then
1296 : dim_a=3
1297 : dim_b=2
1298 : else
1299 3264335 : dim_a=2
1300 3264335 : dim_b=1
1301 : end if
1302 :
1303 : ! Relevant boundaries
1304 1218920002 : gbound(3,ii)=minval(kg_k(dim_a,:))
1305 1218920002 : gbound(4,ii)=maxval(kg_k(dim_a,:))
1306 6528670 : gmin_a=gbound(3,ii)
1307 6528670 : 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 6528670 : if(istwf_k>=2 .and. fftalgc==1)then
1314 30755 : if( istwf_k==2 .or. istwf_k==3 .or. istwf_k==6 .or. istwf_k==7 )then
1315 24884 : gbound(4,2)=max(gmax_a,-gmin_a)
1316 24884 : gbound(3,2)=-gbound(4,2)
1317 5871 : else if( istwf_k==4 .or. istwf_k==5 .or. istwf_k==8 .or. istwf_k==9 )then
1318 5871 : gbound(4,2)=max(gmax_a,-gmin_a-1)
1319 5871 : gbound(3,2)=-gbound(4,2)-1
1320 : end if
1321 30755 : gmax_a=gbound(4,2) ; gmin_a=gbound(3,2)
1322 : end if
1323 :
1324 6528670 : igb=5
1325 :
1326 : ! Consider first every g_a in range 0 ... gmax_a, then gmin_a ... -1
1327 6528670 : gmin(1)=0 ; gmax(1)=gmax_a
1328 6528670 : gmin(2)=gmin_a ; gmax(2)=-1
1329 :
1330 22850345 : do iloop=1,2
1331 :
1332 19586010 : if( gmin(iloop) <= gmax(iloop) )then
1333 :
1334 58322884 : do g_a=gmin(iloop),gmax(iloop)
1335 :
1336 45306510 : if(istwf_k==1 .or. fftalgc/=1)then
1337 : ! Select the minimal and maximal values, in the selected plane
1338 44950196 : gmin_b=mgfft+1 ! Initialized with a value larger than all possible ones
1339 44950196 : gmax_b=-mgfft-1 ! Initialized with a value smaller than all possible ones
1340 13419018168 : do ipw=1,npw
1341 13419018168 : if(kg_k(dim_a,ipw)==g_a)then
1342 1189933933 : kgk=kg_k(dim_b,ipw)
1343 1189933933 : if(kgk<=gmin_b)gmin_b=kgk
1344 1189933933 : if(kgk>=gmax_b)gmax_b=kgk
1345 : end if
1346 : end do
1347 :
1348 356314 : 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 356314 : testp=0
1354 356314 : if(g_a<=gmax_a)then
1355 : ! Select the minimal and maximal values, in the selected plane
1356 356314 : gmin_b1=mgfft+1 ! Initialized with a value larger than all possible ones
1357 356314 : gmax_b1=-mgfft-1 ! Initialized with a value smaller than all possible ones
1358 642690710 : do ipw=1,npw
1359 642690710 : if(kg_k(dim_a,ipw)==g_a)then
1360 22457399 : kgk=kg_k(dim_b,ipw)
1361 22457399 : if(kgk<=gmin_b1)gmin_b1=kgk
1362 22457399 : 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 356314 : testm=0
1372 356314 : if( istwf_k==2 .or. istwf_k==3 .or. istwf_k==6 .or. istwf_k==7 )then
1373 :
1374 288530 : 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 288530 : gmax_b2=mgfft+1 ! Initialized with a value larger than all possible ones
1378 288530 : gmin_b2=-mgfft-1 ! Initialized with a value smaller than all possible ones
1379 505047560 : do ipw=1,npw
1380 505047560 : if(kg_k(dim_a,ipw)==-g_a)then
1381 18377624 : kgk=kg_k(dim_b,ipw)
1382 18377624 : if(kgk<=gmax_b2)gmax_b2=kgk
1383 18377624 : if(kgk>=gmin_b2)gmin_b2=kgk
1384 : end if
1385 : end do
1386 : testm=1
1387 : end if
1388 :
1389 67784 : else if( istwf_k==4 .or. istwf_k==5 .or. istwf_k==8 .or. istwf_k==9 )then
1390 :
1391 67784 : 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 67784 : gmax_b2=mgfft+1 ! Initialized with a value larger than all possible ones
1395 67784 : gmin_b2=-mgfft-1 ! Initialized with a value smaller than all possible ones
1396 137643150 : do ipw=1,npw
1397 137643150 : if(kg_k(dim_a,ipw)==-g_a-1)then
1398 4079775 : kgk=kg_k(dim_b,ipw)
1399 4079775 : if(kgk<=gmax_b2)gmax_b2=kgk
1400 4079775 : 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 356314 : gmin_b2=max(gmin_b2,-mgfft) ; gmax_b2=min(gmax_b2,mgfft)
1413 356314 : if(istwf_k<=5)then
1414 279544 : gmax_b2=-gmax_b2 ; gmin_b2=-gmin_b2
1415 : else
1416 76770 : gmax_b2=-gmax_b2-1 ; gmin_b2=-gmin_b2-1
1417 : end if
1418 : end if
1419 :
1420 356314 : if( testp==1 .and. testm==1)then
1421 356314 : 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 45306510 : 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 45306510 : gbound(igb,ii)=gmin_b
1441 45306510 : gbound(igb+1,ii)=gmax_b
1442 :
1443 45306510 : 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 23717 : if(istwf_k==2 .or. istwf_k==4)then
1446 18779 : gbound(igb+1,ii)=max(gmax_b,-gmin_b)
1447 18779 : gbound(igb,ii)=-gbound(igb+1,ii)
1448 4938 : else if(istwf_k==3 .or. istwf_k==5)then
1449 4938 : gbound(igb+1,ii)=max(gmax_b,-gmin_b-1)
1450 4938 : gbound(igb,ii)=-gbound(igb+1,ii)-1
1451 : end if
1452 :
1453 : end if
1454 :
1455 58322884 : 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 3264335 : 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 2549427 : 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 5098854 : 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 2549427 : if (istwf_k>=2) then
1566 54202 : if(istwf_k==2 .or. istwf_k==4 .or. istwf_k==6 .or. istwf_k==8)then
1567 21567 : i1inver(1)=1
1568 207639 : do i1=2,n1
1569 207639 : i1inver(i1)=n1+2-i1
1570 : end do
1571 : else
1572 344318 : do i1=1,n1
1573 344318 : i1inver(i1)=n1+1-i1
1574 : end do
1575 : end if
1576 54202 : if(istwf_k>=2 .and. istwf_k<=5)then
1577 34537 : i2inver(1)=1
1578 309049 : do i2=2,n2
1579 309049 : i2inver(i2)=n2+2-i2
1580 : end do
1581 : else
1582 216898 : do i2=1,n2
1583 216898 : i2inver(i2)=n2+1-i2
1584 : end do
1585 : end if
1586 54202 : if(istwf_k==2 .or. istwf_k==3 .or. istwf_k==6 .or. istwf_k==7)then
1587 35558 : i3inver(1)=1
1588 380351 : do i3=2,n3
1589 380351 : i3inver(i3)=n3+2-i3
1590 : end do
1591 : else
1592 236819 : do i3=1,n3
1593 236819 : i3inver(i3)=n3+1-i3
1594 : end do
1595 : end if
1596 : end if
1597 :
1598 2549427 : if (iflag==1 .or. iflag==2) then
1599 : ! Insert cg into cfft with extra 0 s around outside:
1600 19199044578 : cfft = zero
1601 :
1602 : ! Take care of each plane wave, and complete cfft if needed
1603 1286749 : if (istwf_k==1) then
1604 :
1605 1233290 : if (iflag==1) then
1606 : !$OMP PARALLEL DO PRIVATE(i1,i2,i3) IF (ndat>1)
1607 2259584 : do idat=1,ndat
1608 492645561 : do ipw=1,npw
1609 490385977 : i1=kg_k(1,ipw); if(i1<0) i1=i1+n1; i1=i1+1
1610 490385977 : i2=kg_k(2,ipw); if(i2<0) i2=i2+n2; i2=i2+1
1611 490385977 : i3=kg_k(3,ipw); if(i3<0) i3=i3+n3; i3=i3+1
1612 :
1613 490385977 : cfft(1,i1,i2,i3+n6*(idat-1))=cg(1,ipw+npw*(idat-1))
1614 491515793 : 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 1233290 : 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 53459 : else if (istwf_k>=2) then
1634 :
1635 53459 : npwmin=1
1636 53459 : if (istwf_k==2 .and. me_g0==1) then
1637 : ! If gamma point, then cfft must be completed
1638 19960 : do idat=1,ndat
1639 9982 : cfft(1,1,1,1+n6*(idat-1))=cg(1,1+npw*(idat-1))
1640 19960 : cfft(2,1,1,1+n6*(idat-1))=zero
1641 : end do
1642 : npwmin=2
1643 : end if
1644 :
1645 53459 : if (iflag==1) then
1646 : !$OMP PARALLEL DO PRIVATE(i1,i1inv,i2,i2inv,i3,i3inv) IF (ndat>1)
1647 103546 : do idat=1,ndat
1648 15752734 : do ipw=npwmin,npw
1649 15649188 : i1=kg_k(1,ipw); if(i1<0) i1=i1+n1; i1=i1+1
1650 15649188 : i2=kg_k(2,ipw); if(i2<0) i2=i2+n2; i2=i2+1
1651 15649188 : i3=kg_k(3,ipw); if(i3<0) i3=i3+n3; i3=i3+1
1652 : ! Construct the coordinates of -k-G
1653 15649188 : i1inv=i1inver(i1) ; i2inv=i2inver(i2) ; i3inv=i3inver(i3)
1654 :
1655 15649188 : cfft(1,i1,i2,i3+n6*(idat-1))=cg(1,ipw+npw*(idat-1))
1656 15649188 : cfft(2,i1,i2,i3+n6*(idat-1))=cg(2,ipw+npw*(idat-1))
1657 15649188 : cfft(1,i1inv,i2inv,i3inv+n6*(idat-1))= cg(1,ipw+npw*(idat-1))
1658 15700963 : 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 53459 : 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 1262678 : else if (iflag==-1 .or. iflag==-2) then
1685 : ! extract cg(output) from cfft(input)
1686 :
1687 1262678 : use_symmetry=0
1688 1262678 : identity(:,:)=0; identity(1,1)=1 ; identity(2,2)=1 ; identity(3,3)=1
1689 16414814 : if(sum((symrec(:,:)-identity(:,:))**2)/=0)use_symmetry=1
1690 5050712 : if(sum(shiftg(:)**2)/=0)use_symmetry=1
1691 :
1692 : ! Extract cg from cfft, ignoring components outside range of cg:
1693 1262678 : if (istwf_k==1) then
1694 :
1695 1261935 : if (use_symmetry==0) then
1696 158524 : if (iflag==-1) then
1697 : !$OMP PARALLEL DO PRIVATE(i1,i2,i3,ipwdat,i3dat) IF (ndat>1)
1698 143122 : do idat=1,ndat
1699 38434135 : do ipw=1,npw
1700 38291013 : i1=kg_k(1,ipw); if(i1<0) i1=i1+n1; i1=i1+1
1701 38291013 : i2=kg_k(2,ipw); if(i2<0) i2=i2+n2; i2=i2+1
1702 38291013 : i3=kg_k(3,ipw); if(i3<0) i3=i3+n3; i3=i3+1
1703 38291013 : ipwdat = ipw + (idat-1) * npw
1704 38291013 : i3dat = i3 + (idat-1) * n6
1705 :
1706 38291013 : cg(1,ipwdat)=cfft(1,i1,i2,i3dat)*xnorm
1707 38362576 : 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 2206822 : do idat=1,ndat
1731 464252121 : do ipw=1,npw
1732 462045299 : l1=kg_k(1,ipw)+shiftg(1)
1733 462045299 : l2=kg_k(2,ipw)+shiftg(2)
1734 462045299 : l3=kg_k(3,ipw)+shiftg(3)
1735 462045299 : j1=symrec(1,1)*l1+symrec(1,2)*l2+symrec(1,3)*l3
1736 462045299 : j2=symrec(2,1)*l1+symrec(2,2)*l2+symrec(2,3)*l3
1737 462045299 : j3=symrec(3,1)*l1+symrec(3,2)*l2+symrec(3,3)*l3
1738 462045299 : if(j1<0) j1=j1+n1; i1=j1+1
1739 462045299 : if(j2<0) j2=j2+n2; i2=j2+1
1740 462045299 : 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 462045299 : ipwdat = ipw + (idat-1) * npw
1745 462045299 : i3dat = i3 + (idat-1)*n6
1746 :
1747 : ! c(g) = cfft(S(g + shiftg))
1748 462045299 : cg(1,ipwdat)=cfft(1,i1,i2,i3dat)*xnorm
1749 463148710 : cg(2,ipwdat)=cfft(2,i1,i2,i3dat)*xnorm
1750 : end do
1751 : end do
1752 : end if
1753 :
1754 743 : else if (istwf_k>=2) then
1755 :
1756 743 : npwmin=1
1757 743 : 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 743 : 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 376 : id1=n1/2+2
1817 376 : id2=n2/2+2
1818 376 : 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 752 : do idat=1,ndat
1822 45834 : do ipw=npwmin,npw
1823 :
1824 45082 : i1=kg_k(1,ipw); if(i1<0) i1=i1+n1; i1=i1+1
1825 45082 : i2=kg_k(2,ipw); if(i2<0) i2=i2+n2; i2=i2+1
1826 45082 : i3=kg_k(3,ipw); if(i3<0) i3=i3+n3; i3=i3+1
1827 :
1828 45082 : i1inv=i1inver(i1) ; i2inv=i2inver(i2) ; i3inv=i3inver(i3)
1829 :
1830 45082 : l1=kg_k(1,ipw)+shiftg(1)
1831 45082 : l2=kg_k(2,ipw)+shiftg(2)
1832 45082 : l3=kg_k(3,ipw)+shiftg(3)
1833 45082 : j1=symrec(1,1)*l1+symrec(1,2)*l2+symrec(1,3)*l3
1834 45082 : j2=symrec(2,1)*l1+symrec(2,2)*l2+symrec(2,3)*l3
1835 45082 : j3=symrec(3,1)*l1+symrec(3,2)*l2+symrec(3,3)*l3
1836 45082 : if(j1<0)j1=j1+n1 ; i1=j1+1
1837 45082 : if(j2<0)j2=j2+n2 ; i2=j2+1
1838 45082 : if(j3<0)j3=j3+n3 ; i3=j3+1
1839 :
1840 : ! Construct the coordinates of -k-G
1841 45082 : l1=i1inv-(i1inv/id1)*n1-1+shiftg(1)
1842 45082 : l2=i2inv-(i2inv/id2)*n2-1+shiftg(2)
1843 45082 : l3=i3inv-(i3inv/id3)*n3-1+shiftg(3)
1844 45082 : j1=symrec(1,1)*l1+symrec(1,2)*l2+symrec(1,3)*l3
1845 45082 : j2=symrec(2,1)*l1+symrec(2,2)*l2+symrec(2,3)*l3
1846 45082 : j3=symrec(3,1)*l1+symrec(3,2)*l2+symrec(3,3)*l3
1847 45082 : if(j1<0)j1=j1+n1 ; i1inv=j1+1
1848 45082 : if(j2<0)j2=j2+n2 ; i2inv=j2+1
1849 45082 : if(j3<0)j3=j3+n3 ; i3inv=j3+1
1850 :
1851 45082 : ipwdat = ipw + (idat-1) * npw
1852 45082 : i3dat = i3 + (idat-1) * n6
1853 45082 : i3invdat = i3inv + (idat-1) * n6
1854 :
1855 : ! Here the time-reversal symmetry is used to project from cfft
1856 45082 : cg(1,ipwdat)=(cfft(1,i1,i2,i3dat) + cfft(1,i1inv,i2inv,i3invdat))*0.5d0*xnorm
1857 45458 : 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 2549427 : 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 696830 : 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 1247558986 : cfft = zero
2042 : !$OMP PARALLEL DO PRIVATE(i1,i2,i2_local,i3)
2043 1413984 : do idat=1,ndat
2044 139650584 : do ipw=1,npw
2045 138236600 : i1=kg_k(1,ipw); if(i1<0)i1=i1+n1; i1=i1+1
2046 138236600 : i2=kg_k(2,ipw); if(i2<0)i2=i2+n2; i2=i2+1
2047 138236600 : i3=kg_k(3,ipw); if(i3<0)i3=i3+n3; i3=i3+1
2048 138236600 : i2_local = tab_fftwf2_local(i2) + n6*(idat-1)
2049 138236600 : cfft(1,i1,i3,i2_local)=cg(1,ipw+npw*(idat-1))
2050 138953754 : cfft(2,i1,i3,i2_local)=cg(2,ipw+npw*(idat-1))
2051 : end do
2052 : end do
2053 :
2054 696830 : 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 1158611 : 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 17881710 : do j=1,n1dfft
2146 471276667 : do i=1,n2
2147 453394957 : zw(1,j,i)=zt(1,i,j)
2148 470118056 : zw(2,j,i)=zt(2,i,j)
2149 : end do
2150 : end do
2151 :
2152 1158611 : 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 10425179 : 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 216259320 : do j=1,n1dfft
2197 1514094604 : do i=1,max2+1
2198 1297835284 : zw(1,j,i)=zt(1,i,j)
2199 1503669425 : zw(2,j,i)=zt(2,i,j)
2200 : end do
2201 : end do
2202 :
2203 : ! Fill the center region with zeros
2204 146001875 : do i=max2+2,n2-m2+max2+1
2205 2752584651 : do j=1,n1dfft
2206 2606582776 : zw(1,j,i)=zero
2207 2742159472 : zw(2,j,i)=zero
2208 : end do
2209 : end do
2210 :
2211 : ! Here, negative frequencies
2212 10425179 : if (m2>=max2+2) then
2213 216259320 : do j=1,n1dfft
2214 1239573669 : do i=max2+2,m2
2215 1023314349 : zw(1,j,i+n2-m2)=zt(1,i,j)
2216 1229148490 : zw(2,j,i+n2-m2)=zt(2,i,j)
2217 : end do
2218 : end do
2219 : end if
2220 :
2221 10425179 : 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 5991720 : 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 149972498 : do i3=1,n3
2448 1685672926 : do i=0,n1dfft-1
2449 1535700428 : zmpi2(1,i1+i,j2,i3)=zw(1,i+1,i3)
2450 1679681206 : zmpi2(2,i1+i,j2,i3)=zw(2,i+1,i3)
2451 : end do
2452 : end do
2453 :
2454 5991720 : 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 1163090 : 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 33608291 : do i3=1,n3
2491 491344524 : do i1=1,n1dfft
2492 457736233 : zw(1,i1,i3)=zf(1,i1,i3)
2493 490181434 : zw(2,i1,i3)=zf(2,i1,i3)
2494 : end do
2495 : end do
2496 :
2497 1163090 : 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 4828630 : 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 32935356 : do i3=1,max3+1
2540 303796692 : do i1=1,n1dfft
2541 270861336 : zw(1,i1,i3)=zf(1,i1,i3)
2542 298968062 : zw(2,i1,i3)=zf(2,i1,i3)
2543 : end do
2544 : end do
2545 :
2546 : ! Fill the center region with zeros
2547 64163636 : do i3=max3+2,n3-m3+max3+1
2548 638200577 : do i1=1,n1dfft
2549 574036941 : zw(1,i1,i3)=zero
2550 633371947 : zw(2,i1,i3)=zero
2551 : end do
2552 : end do
2553 :
2554 : ! Here, negative frequencies
2555 28922475 : do i3=max3+2,m3
2556 261988393 : do i1=1,n1dfft
2557 233065918 : zw(1,i1,i3+n3-m3)=zf(1,i1,i3)
2558 257159763 : zw(2,i1,i3+n3-m3)=zf(2,i1,i3)
2559 : end do
2560 : end do
2561 :
2562 4828630 : 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 1205401 : 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 34634785 : do i3=1,n3
2598 524509437 : do i1=1,n1dfft
2599 489874652 : zf(1,i1,i3)=zw(1,i1,i3)
2600 523304036 : zf(2,i1,i3)=zw(2,i1,i3)
2601 : end do
2602 : end do
2603 :
2604 1205401 : 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 4423910 : 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 29728922 : do i3=1,max3+1
2644 273951362 : do i1=1,n1dfft
2645 244222440 : zf(1,i1,i3)=zw(1,i1,i3)
2646 269527452 : zf(2,i1,i3)=zw(2,i1,i3)
2647 : end do
2648 : end do
2649 :
2650 : ! Here, negative frequencies
2651 26044554 : do i3=max3+2,m3
2652 235512128 : do i1=1,n1dfft
2653 209467574 : zf(1,i1,i3)=zw(1,i1,i3+n3-m3)
2654 231088218 : zf(2,i1,i3)=zw(2,i1,i3+n3-m3)
2655 : end do
2656 : end do
2657 :
2658 4423910 : 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 1207862 : 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 1207862 : mfft=0
2688 1207862 : if (ioption == 2) then
2689 2347290 : do Jp2=Jp2st,nproc
2690 19093694 : do J2=J2st,nd2proc
2691 17954266 : mfft=mfft+1
2692 17954266 : if (mfft.gt.n1dfft) then
2693 275482 : Jp2st=Jp2
2694 275482 : J2st=J2
2695 275482 : return
2696 : end if
2697 508692864 : do I1=1,n1
2698 489874652 : zmpi1(1,I1,J2,j3,Jp2)=zw(1,mfft,I1)
2699 507553436 : zmpi1(2,I1,J2,j3,Jp2)=zw(2,mfft,I1)
2700 : end do
2701 : end do
2702 2071808 : 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 9179512 : 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 9179512 : mfft=0
2759 :
2760 9179512 : 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 25525700 : do Jp2=Jp2stf,nproc
2789 107067504 : do J2=J2stf,md2proc
2790 90721316 : mfft=mfft+1
2791 90721316 : if (mfft.gt.n1dfft) then
2792 320636 : Jp2stf=Jp2
2793 320636 : J2stf=J2
2794 320636 : return
2795 : end if
2796 90400680 : ind=(Jp2-1) * md2proc + J2
2797 90400680 : jj2=(ind-1)/nproc +1
2798 :
2799 : !jjp2=modulo(ind,nproc) +1
2800 90400680 : jjp2=modulo(ind-1,nproc)+1
2801 :
2802 : ! Here, zero and positive frequencies
2803 626085608 : do I1=1,max1+1
2804 535684928 : zmpi1(1,I1,Jj2,j3,Jjp2)=zw(1,mfft,I1)
2805 626085608 : zmpi1(2,I1,Jj2,j3,Jjp2)=zw(2,mfft,I1)
2806 : end do
2807 :
2808 : ! Here, negative frequencies
2809 542520640 : do I1=max1+2,m1
2810 435773772 : zmpi1(1,I1,Jj2,j3,Jjp2)=zw(1,mfft,I1+n1-m1)
2811 526174452 : zmpi1(2,I1,Jj2,j3,Jjp2)=zw(2,mfft,I1+n1-m1)
2812 : end do
2813 : end do
2814 25205064 : 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 5629311 : 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 139400815 : do i3=1,n3
2855 1600734167 : do i=0,n1dfft-1
2856 1461333352 : zw(1,i+1,i3)=zmpi2(1,i1+i,j2,i3)
2857 1595104856 : zw(2,i+1,i3)=zmpi2(2,i1+i,j2,i3)
2858 : end do
2859 : end do
2860 :
2861 5629311 : 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 1201940 : 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 18833000 : do j=1,n1dfft
2897 504686164 : do i=1,n2
2898 485853164 : zt(1,i,j)=zw(1,j,i)
2899 503484224 : zt(2,i,j)=zw(2,j,i)
2900 : end do
2901 : end do
2902 :
2903 1201940 : 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 9503572 : 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 198521796 : do j=1,n1dfft
2942 1376154088 : do i=1,max2+1
2943 1177632292 : zt(1,i,j)=zw(1,j,i)
2944 1366650516 : zt(2,i,j)=zw(2,j,i)
2945 : end do
2946 : end do
2947 :
2948 : ! Here, negative frequencies
2949 9503572 : if(m2>=max2+2)then
2950 198521796 : do j=1,n1dfft
2951 1124123960 : do i=max2+2,m2
2952 925602164 : zt(1,i,j)=zw(1,j,i+n2-m2)
2953 1114620388 : zt(2,i,j)=zw(2,j,i+n2-m2)
2954 : end do
2955 : end do
2956 : end if
2957 :
2958 9503572 : 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 1165091 : 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 1165091 : mfft=0
3108 :
3109 1165091 : if (ioption /= 1) then
3110 2277840 : do Jp2=Jp2st,nproc
3111 18124028 : do J2=J2st,nd2proc
3112 17011279 : mfft=mfft+1
3113 17011279 : if (mfft.gt.n1dfft) then
3114 244558 : Jp2st=Jp2
3115 244558 : J2st=J2
3116 244558 : return
3117 : end if
3118 475615703 : do I1=1,n1
3119 457736233 : zw(1,mfft,I1)=zmpi1(1,I1,J2,j3,Jp2)
3120 474502954 : zw(2,mfft,I1)=zmpi1(2,I1,J2,j3,Jp2)
3121 : end do
3122 : end do
3123 2033282 : 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 10021215 : pure subroutine mpiswitch_cent(j3,n1dfft,Jp2stb,J2stb,lot,max1,md1,m1,n1,md2proc,&
3180 10021215 : 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 10021215 : mfft=0
3196 :
3197 10021215 : 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 27644968 : do Jp2=Jp2stb,nproc
3234 116430980 : do J2=J2stb,md2proc
3235 :
3236 98807227 : mfft=mfft+1
3237 98807227 : if (mfft.gt.n1dfft) then
3238 397780 : Jp2stb=Jp2
3239 397780 : J2stb=J2
3240 : !ABI_WARNING("Returning from mpiswithc_cent")
3241 397780 : return
3242 : end if
3243 :
3244 98409447 : ind=(Jp2-1) * md2proc + J2
3245 98409447 : jj2=(ind-1)/nproc +1
3246 :
3247 : !jjp2=modulo(ind,nproc) +1
3248 98409447 : 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 689750779 : do I1=1,max1+1
3254 591341332 : zw(1,mfft,I1)=zmpi1(1,I1,Jj2,j3,Jjp2)
3255 689750779 : zw(2,mfft,I1)=zmpi1(2,I1,Jj2,j3,Jjp2)
3256 : end do
3257 :
3258 : ! Fill the center region with zeros
3259 1355559091 : do I1=max1+2,n1-m1+max1+1
3260 1257149644 : zw(1,mfft,I1)=zero
3261 1355559091 : zw(2,mfft,I1)=zero
3262 : end do
3263 :
3264 : ! Here, negative frequencies
3265 : ! In zmpi1, they are stored from 1 to m1half
3266 599339553 : do I1=max1+2,m1
3267 483306353 : zw(1,mfft,I1+n1-m1)=zmpi1(1,I1,Jj2,j3,Jjp2)
3268 581715800 : zw(2,mfft,I1+n1-m1)=zmpi1(2,I1,Jj2,j3,Jjp2)
3269 : end do
3270 :
3271 : end do
3272 27247188 : 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 67683 : 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 135374 : do idat=1,ndat
3308 1202910 : do i3=1,n3
3309 22358765 : do i2=1,n2
3310 22291074 : if (fftn2_distrib(i2) == me_fft) then
3311 15318706 : i2_local = ffti2_local(i2)
3312 15318706 : i2_ldat = i2_local + (idat-1) * nd2proc
3313 15318706 : fgbase= n1*(i2_local-1 + nd2proc*(i3-1)) + (idat-1) * nfft
3314 377619812 : do i1=1,n1
3315 362301106 : workf(1,i1,i3,i2_ldat)=fofg(1,i1+fgbase)
3316 377619812 : 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 67683 : 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 67437 : 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 67437 : xnorm=one/dble(n1*n2*n3)
3402 :
3403 : ! Transfer fft output to the original fft box
3404 134882 : do idat=1,ndat
3405 1236953 : do i2=1,n2
3406 1169516 : if( fftn2_distrib(i2) == me_fft) then
3407 895075 : i2_local = ffti2_local(i2)
3408 895075 : i2_ldat = i2_local + (idat-1) * nd2proc
3409 17511538 : do i3=1,n3
3410 16616463 : fgbase = n1*(i2_local - 1 + nd2proc*(i3-1)) + (idat - 1) * nfft
3411 440298291 : do i1=1,n1
3412 422786753 : fofg(1,i1+fgbase)=workf(1,i1,i3,i2_ldat)*xnorm
3413 439403216 : 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 67437 : 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 67683 : 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 135352 : do idat=1,ndat
3506 1202318 : do i3=1,n3
3507 1134644 : if( fftn3_distrib(i3) == me_fft) then
3508 874750 : i3_local = ffti3_local(i3)
3509 874750 : i3_ldat = i3_local + (idat - 1) * nd3proc
3510 16150796 : do i2=1,n2
3511 15276046 : frbase=n1*(i2-1+n2*(i3_local-1)) + (idat - 1) * nfft
3512 374780982 : do i1=1,n1
3513 373906232 : 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 67683 : fofr = huge(one)
3545 : end select
3546 :
3547 67683 : 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 67437 : 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 134860 : do idat=1,ndat
3654 1236637 : do i3=1,n3
3655 1169209 : if( me_fft == fftn3_distrib(i3) ) then
3656 894729 : i3_local = ffti3_local(i3)
3657 894729 : i3_ldat = i3_local + (idat-1) * nd3proc
3658 17468532 : do i2=1,n2
3659 16573803 : frbase=n1*(i2-1+n2*(i3_local-1)) + (idat-1) * nfft
3660 436584365 : do i1=1,n1
3661 419115833 : workr(1,i1,i2,i3_ldat)=fofr(i1+frbase)
3662 435689636 : 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 67437 : workr = huge(one)
3690 : end select
3691 :
3692 67437 : 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 117247689 : 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 117247689 : integer,allocatable :: index2d(:,:)
3806 : ! *************************************************************************
3807 :
3808 117247689 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
3809 :
3810 : !First, generate a 2d index for each column of data
3811 468990756 : ABI_MALLOC(index2d,(n1,n2))
3812 36318395895 : index2d(:,:)=0
3813 117247689 : index=1
3814 117247689 : igb=3
3815 613473582 : do g2=0,gbound(2) ! g2max
3816 2567075530 : do g1=0,gbound(igb+1) ! g1max
3817 2070849637 : index2d(g1+1,g2+1)=index
3818 2567075530 : index=index+1
3819 : end do
3820 496225893 : if(gbound(igb)<=-1)then ! g1min
3821 1460826527 : do g1=gbound(igb)+n1,n1-1
3822 1063715045 : index2d(g1+1,g2+1)=index
3823 1460826527 : index=index+1
3824 : end do
3825 : end if
3826 613473582 : igb=igb+2
3827 : end do
3828 :
3829 117247689 : if(gbound(1)<=-1)then ! g2min
3830 504282407 : do g2=gbound(1)+n2,n2-1
3831 1397552686 : do g1=0,gbound(igb+1)
3832 1007899391 : index2d(g1+1,g2+1)=index
3833 1397552686 : index=index+1
3834 : end do
3835 389653295 : if(gbound(igb)<=-1)then
3836 1702060799 : do g1=gbound(igb)+n1,n1-1
3837 1324378791 : index2d(g1+1,g2+1)=index
3838 1702060799 : index=index+1
3839 : end do
3840 : end if
3841 504282407 : igb=igb+2
3842 : end do
3843 : end if
3844 :
3845 117247689 : 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 31864812876 : do ipw=1,npw_k
3856 31747565187 : i1=kg_k(1,ipw); if(i1<0)i1=i1+n1 ; i1=i1+1
3857 31747565187 : i2=kg_k(2,ipw); if(i2<0)i2=i2+n2 ; i2=i2+1
3858 31747565187 : i3=kg_k(3,ipw); if(i3<0)i3=i3+n3 ; i3=i3+1
3859 31747565187 : indpw_k(4,ipw)=index2d(i1,i2)
3860 31864812876 : indpw_k(3,ipw)=i3
3861 : end do
3862 :
3863 2027861285 : do i1=1,n1
3864 36323241967 : do i2=1,n2
3865 34295380682 : index=index2d(i1,i2)
3866 36205994278 : if(index/=0)then
3867 5466842864 : indpw_k(1,index)=i1
3868 5466842864 : indpw_k(2,index)=i2
3869 : end if
3870 : end do
3871 : end do
3872 :
3873 117247689 : ABI_FREE(index2d)
3874 :
3875 117247689 : 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 1787939 : 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 1787939 : integer,allocatable :: array_ipw(:),ig1arr(:),ig2arr(:)
3938 1787939 : integer,allocatable :: ig3arr(:),kg_ind(:),kg_small(:,:)
3939 1787939 : integer, allocatable :: npw_gather(:),npw_disp(:) ,kg_ind_gather(:),kg_small_gather(:,:)
3940 : real(dp) :: kmax(3),minor(3),numer_(3),tsec(2)
3941 1787939 : real(dp),allocatable :: kg1arr(:),kg2arr(:),kg3arr(:)
3942 : ! *************************************************************************
3943 :
3944 : DBG_ENTER("COLL")
3945 :
3946 1787939 : call timab(23,1,tsec)
3947 1787939 : 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 1787939 : 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 1787939 : np_band=0
3965 1787939 : if (mpw>0) then
3966 918203 : np_band=1; if(mpi_enreg%paral_kgb==1) np_band=max(1,mpi_enreg%nproc_band)
3967 918203 : alloc_size=max(alloc_size,(mpw+1)*np_band)
3968 2754609 : ABI_MALLOC(kg_small,(3, alloc_size))
3969 2754609 : ABI_MALLOC(kg_ind,(alloc_size))
3970 326396457 : 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 1787939 : 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 1787939 : minor(1)=gmet(2,2)*gmet(3,3)-gmet(2,3)**2
3981 1787939 : 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 1787939 : minor(2)=gmet(1,1)*gmet(3,3)-gmet(1,3)**2
3983 1787939 : 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 1787939 : minor(3)=gmet(2,2)*gmet(1,1)-gmet(1,2)**2
3985 1787939 : 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 1787939 : gmet_trace=gmet(1,1)+gmet(2,2)+gmet(3,3)
3989 :
3990 7151756 : do ii=1,3
3991 5363817 : xx=gmet(ii,ii)*minor(ii)-numer_(ii)
3992 5363817 : 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 5363817 : kmax(ii)=sqrt(gscut*minor(ii)/xx)
3996 5363817 : nmax(ii)=floor(kmax(ii)-kpt(ii)+tol10)
3997 5363817 : nmin(ii)=ceiling(-kmax(ii)-kpt(ii)-tol10)
3998 7151756 : 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 1787939 : ig=0;ind=0
4002 1787939 : in=0
4003 1787939 : 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 5363817 : ABI_MALLOC(ig1arr,(ngrid(1)))
4010 5363817 : ABI_MALLOC(ig2arr,(ngrid(2)))
4011 5363817 : ABI_MALLOC(ig3arr,(ngrid(3)))
4012 5363817 : ABI_MALLOC(kg1arr,(ngrid(1)))
4013 5363817 : ABI_MALLOC(kg2arr,(ngrid(2)))
4014 5363817 : ABI_MALLOC(kg3arr,(ngrid(3)))
4015 :
4016 16686838 : do ig1p=1,ngrid(1)
4017 14898899 : ig1arr(ig1p)=ig1p-1
4018 14898899 : if (ig1p-1>nmax(1)) ig1arr(ig1p)=ig1p-ngrid(1)-1
4019 16686838 : 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 1787939 : ig2pmax=ngrid(2)
4025 1787939 : if(istwf_k>=2 .and. exchn2n3d==0) ig2pmax=nmax(2)+1
4026 5363817 : ABI_MALLOC(array_ipw,(-ig2pmax:ig2pmax))
4027 33266572 : array_ipw(:)=0
4028 16633286 : do ig2p=1,ig2pmax
4029 14845347 : ig2arr(ig2p)=ig2p-1
4030 14845347 : if (ig2p-1>nmax(2)) ig2arr(ig2p)=ig2p-ngrid(2)-1
4031 16633286 : 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 1787939 : ig3pmax=ngrid(3)
4037 1787939 : if (istwf_k>=2 .and. exchn2n3d==1) ig3pmax=nmax(3)+1
4038 :
4039 16933611 : do ig3p=1,ig3pmax
4040 15145672 : ig3arr(ig3p)=ig3p-1
4041 15145672 : if(ig3p-1>nmax(3)) ig3arr(ig3p)=ig3p-ngrid(3)-1
4042 16933611 : kg3arr(ig3p)=kpt(3)+dble(ig3arr(ig3p))
4043 : end do
4044 :
4045 : !Performs loop on all grid points.
4046 : !---------------------------------------------------------------------
4047 1787939 : igtot = 0
4048 1787939 : if(exchn2n3d==0)then
4049 1787937 : mpi_enreg%me_g0=0
4050 16933575 : do ig3p=1,ngrid(3)
4051 15145638 : ig3=ig3arr(ig3p)
4052 15145638 : v3=kg3arr(ig3p)
4053 15145638 : ig2pmin=1
4054 15145638 : if( istwf_k>=2 .and. istwf_k<=5 .and. ig3<0)then
4055 30246 : ig2pmin=2
4056 : end if
4057 : ! ig2pmax was initialized previously
4058 161625722 : do ig2p=ig2pmin,ig2pmax
4059 144692147 : ig2=ig2arr(ig2p)
4060 : ! PAY ATTENTION : the proc 0 must have me_g0=1
4061 144692147 : ipw_ok = .true.
4062 144692147 : if(mpi_enreg%paral_kgb==1 ) then
4063 504583 : n2 =mpi_enreg%distribfft%n2_coarse
4064 504583 : ipw_ok = ipw_ok.and.(mpi_enreg%me_fft == mpi_enreg%distribfft%tab_fftwf2_distrib( modulo(ig2,n2) + 1))
4065 : end if
4066 144692147 : if (ig2==0 .and. ipw_ok) mpi_enreg%me_g0=1
4067 144692147 : v2=kg2arr(ig2p)
4068 144692147 : gs_part=gmet(2,2)*v2*v2+gmet(3,3)*v3*v3+2.0_dp*gmet(2,3)*v2*v3
4069 144692147 : gs_fact=2.0_dp*(gmet(1,2)*v2+gmet(3,1)*v3)
4070 144692147 : ig1pmax=ngrid(1)
4071 144692147 : if( (istwf_k==2.or.istwf_k==3) .and. ig3p==1 .and. ig2p==1)ig1pmax=nmax(1)+1
4072 1763033375 : do ig1p=1,ig1pmax
4073 1603195590 : v1=kg1arr(ig1p)
4074 1603195590 : gmin=gs_part+v1*(gs_fact+v1*gmet11)
4075 : ! If inside sphere:
4076 1747887737 : if (gmin<=gscut) then
4077 623286029 : if (ipw_ok) then
4078 622831221 : ig=ig+1 ! inside sphere
4079 622831221 : igtot=igtot+1
4080 622831221 : if (mpw>0.and.ig<=alloc_size) then
4081 : ! Keep coords of pw:
4082 312218333 : kg_small(1,ig)=ig1arr(ig1p)
4083 312218333 : kg_small(2,ig)=ig2
4084 312218333 : kg_small(3,ig)=ig3
4085 312218333 : kg_ind(ig)=igtot
4086 : end if
4087 622831221 : array_ipw(ig2)=array_ipw(ig2)+1
4088 : else
4089 454808 : 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 1787939 : if (mpw*np_band > 0 .and. ig > alloc_size) then
4153 132 : npw = mpw*np_band
4154 : else
4155 1787807 : npw=ig
4156 : end if
4157 1787939 : alloc_size = max(alloc_size,npw)
4158 :
4159 1787939 : ABI_FREE(ig1arr)
4160 1787939 : ABI_FREE(ig2arr)
4161 1787939 : ABI_FREE(ig3arr)
4162 1787939 : ABI_FREE(kg1arr)
4163 1787939 : ABI_FREE(kg2arr)
4164 1787939 : ABI_FREE(kg3arr)
4165 :
4166 : !BandFFT: plane-wave load balancing
4167 1787939 : 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 2018 : np_fft=max(1,mpi_enreg%nproc_fft)
4170 6054 : ABI_MALLOC(npw_gather,(np_fft)) ! Count pw before balancing
4171 2018 : call xmpi_allgather(npw,npw_gather,mpi_enreg%comm_fft,ierr)
4172 16046 : gap_pw = 100._dp*(maxval(npw_gather(:))-minval(npw_gather))/(1.*sum(npw_gather(:))/np_fft)
4173 2018 : write(msg,'(a,f5.2)' ) ' Relative gap for number of plane waves between process (%): ',gap_pw
4174 2018 : call wrtout(std_out,msg)
4175 2018 : 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 4036 : ABI_FREE(npw_gather)
4213 : end if
4214 :
4215 : !BandFFT: band load balancing
4216 1787939 : if(mpi_enreg%paral_kgb==1.and.mpi_enreg%nproc_band>0) then
4217 6991 : np_band=max(1,mpi_enreg%nproc_band)
4218 6991 : npw_split=ig;npw=npw_split/np_band
4219 6991 : npw_remain=modulo(npw_split,np_band)
4220 6991 : if(mpi_enreg%me_band < npw_remain) npw=npw+1
4221 6991 : if(mpw > 0) then ! This is for the case when we only compute npw and put mpw=0
4222 3003 : npw_before=mpi_enreg%me_band*(npw_split/np_band)+min(npw_remain,mpi_enreg%me_band)
4223 2501483 : kg_small(:,1:npw)=kg_small(:,npw_before+1:npw_before+npw)
4224 627623 : kg_ind ( 1:npw)=kg_ind ( npw_before+1:npw_before+npw)
4225 : end if
4226 : end if
4227 1787939 : if(mpw > 0) then
4228 312762186 : do i1=1,npw
4229 1247375932 : kg(:,i1+ikg)=kg_small(:,i1)
4230 312762186 : if (allocated(mpi_enreg%my_kgtab)) then
4231 270755 : ikpt_this_proc=mpi_enreg%my_kpttab(ikpt)
4232 270755 : mpi_enreg%my_kgtab(i1,ikpt_this_proc) = kg_ind(i1)
4233 : end if
4234 : end do
4235 918203 : ABI_FREE(kg_small)
4236 918203 : ABI_FREE(kg_ind)
4237 : end if
4238 :
4239 1787939 : ABI_FREE(array_ipw)
4240 :
4241 : !Take care of the me_g0 flag
4242 1787939 : mpi_enreg%me_g0_fft=mpi_enreg%me_g0
4243 1787939 : if(mpi_enreg%paral_kgb==1.and.mpi_enreg%nproc_band>0) then
4244 6991 : 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 2334 : mpi_enreg%me_g0=1
4247 : else
4248 : ! All other cases
4249 4657 : mpi_enreg%me_g0=0
4250 : end if
4251 : end if
4252 :
4253 : !Check that npw is not zero
4254 1787939 : 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 1787939 : call timab(23,2,tsec)
4263 :
4264 : DBG_EXIT("COLL")
4265 :
4266 1787939 : 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 75 : 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 75 : gscut=0.5_dp*ecut*piinv**2
4320 75 : gmet_trace=gmet(1,1)+gmet(2,2)+gmet(3,3)
4321 75 : minor(1)=gmet(2,2)*gmet(3,3)-gmet(2,3)**2
4322 75 : minor(2)=gmet(1,1)*gmet(3,3)-gmet(1,3)**2
4323 75 : minor(3)=gmet(2,2)*gmet(1,1)-gmet(1,2)**2
4324 75 : 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 75 : 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 75 : 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 300 : ngmin(:)=1000000;ngmax(:)=0
4329 784 : do ikpt=1,nkpt
4330 709 : istwf_k=istwfk(ikpt)
4331 :
4332 2836 : do ii=1,3
4333 2127 : xx=gmet(ii,ii)*minor(ii)-numer_(ii)
4334 2127 : 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 2127 : kmax=sqrt(gscut*minor(ii)/xx)
4338 2127 : nmax(ii)=floor(kmax-kpt(ii,ikpt)+tol10)
4339 2127 : nmin=ceiling(-kmax-kpt(ii,ikpt)-tol10)
4340 2836 : ngrid(ii)=nmax(ii)-nmin+1
4341 : end do
4342 :
4343 709 : ng1=ngrid(1);if(istwf_k==2.or.istwf_k==3) ng1=nmax(1)+1
4344 709 : if(exchn2n3d==0)then
4345 709 : ng3=ngrid(3)
4346 709 : 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 709 : if(ng1<ngmin(1)) ngmin(1)=ng1
4355 709 : if(ng2<ngmin(2)) ngmin(2)=ng2
4356 709 : if(ng3<ngmin(3)) ngmin(3)=ng3
4357 709 : if(ng1>ngmax(1)) ngmax(1)=ng1
4358 709 : if(ng2>ngmax(2)) ngmax(2)=ng2
4359 784 : if(ng3>ngmax(3)) ngmax(3)=ng3
4360 :
4361 : end do
4362 :
4363 : DBG_EXIT("COLL")
4364 :
4365 75 : 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 395398 : 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 395398 : type(MPI_type) :: MPI_enreg_seq
4415 : !arrays
4416 : integer :: kg_dum(3, 0)
4417 395398 : integer,allocatable :: iwork(:,:)
4418 : ! *********************************************************************
4419 :
4420 395398 : call initmpi_seq(MPI_enreg_seq)
4421 :
4422 : ! Calculate the number of G-vectors for this k-point.
4423 395398 : 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 1186194 : ABI_MALLOC(kg_k,(3,npw_k))
4427 395398 : call kpgsph(ecut, exchn2n3d0, gmet, ikg0, 0, istwf_k, kg_k, kpoint, mkmem_, MPI_enreg_seq, npw_k, npw_k_test)
4428 :
4429 395398 : call destroy_mpi_enreg(MPI_enreg_seq)
4430 :
4431 395398 : 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 395398 : if (present(mpw)) mpw = max(mpw, npw_k)
4440 :
4441 395398 : 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 395398 : 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 135620 : 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 135620 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
4491 : !integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
4492 : ! *************************************************************************
4493 :
4494 135620 : 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 135620 : me_fft=mpi_enreg%me_fft
4501 135620 : nd2=(n2-1)/mpi_enreg%nproc_fft+1
4502 :
4503 135620 : if (n2== mpi_enreg%distribfft%n2_coarse) then
4504 135620 : fftn2_distrib => mpi_enreg%distribfft%tab_fftdp2_distrib
4505 135620 : 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 10269551 : do ig=1,npw_k
4516 10133931 : ig1=modulo(kg_k(1,ig),n1)
4517 10133931 : ig2=modulo(kg_k(2,ig),n2)
4518 10133931 : ig3=modulo(kg_k(3,ig),n3)
4519 10133931 : if(me_fft==fftn2_distrib(ig2+1)) then
4520 10133931 : ig2=ffti2_local(ig2+1) - 1
4521 10133931 : indpw_k(ig)=ig1+1+n1*(ig2+nd2*ig3)
4522 10133931 : mask(ig)=.true.
4523 : else
4524 0 : indpw_k(ig)=0
4525 0 : mask(ig)=.false.
4526 : end if
4527 81207068 : 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 135620 : 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 : !!***
|
