TABLE OF CONTENTS


ABINIT/getngrec [ Functions ]

[ Top ] [ Functions ]

NAME

 getngrec

FUNCTION

 This routine computes the fft box for the recursion method, accordingly to the troncation radius.
 It is quite similar to getng, but : 
     - there is no xboxcut and ecut consistency
     - ngfft (the initial fft box) is the maximum fft box

COPYRIGHT

 Copyright (C) 2008-2018 ABINIT group ( ).
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

  ngfft(18)=non truncated fft box
  mgfft=maximum of ngfft(1:3)
  inf_rmet=define the infinitesimal metric : rprimd*(transpose(rprimd)), divided by the number of discretisation point
  recrcut=truncating
  delta=to obtain radius of truncation

OUTPUT

  ngfftrec=truncated fft box
  nfftrec= truncated nfft 
  tronc=True if truncation is made 

SIDE EFFECTS

PARENTS

      m_rec

CHILDREN

      sort_int,timab

SOURCE

 41 #if defined HAVE_CONFIG_H
 42 #include "config.h"
 43 #endif
 44 
 45 #include "abi_common.h"
 46 
 47 subroutine getngrec(ngfft,rmet,ngfftrec,nfftrec,recrcut,delta,tronc)
 48 
 49  use defs_basis
 50  use m_profiling_abi
 51  use m_sort
 52 
 53 !This section has been created automatically by the script Abilint (TD).
 54 !Do not modify the following lines by hand.
 55 #undef ABI_FUNC
 56 #define ABI_FUNC 'getngrec'
 57  use interfaces_18_timing
 58 !End of the abilint section
 59 
 60 implicit none
 61 
 62 !Arguments -------------------------------
 63 !scalars
 64 real(dp),intent(in) :: recrcut,delta
 65 integer,intent(out) :: nfftrec
 66 logical,intent(out) :: tronc
 67 !arrays
 68 integer,intent(in)  :: ngfft(18)
 69 real(dp),intent(in) :: rmet(3,3)
 70 integer,intent(out) :: ngfftrec(18)
 71 
 72 !Local variables-------------------------------
 73 !scalars
 74 integer :: ii,iimin,index,jj,jjmin,kk,kkmin,largest_ngfftrec,maxpow11,maxpow2
 75 integer :: maxpow3,maxpow5,maxpow7,mmsrch,plane
 76 real(dp) :: dsm,dsp,dsqmin,rtroncat
 77 !arrays
 78 integer :: get_ngfftrec(3),maxsrch(3),minsrch(3)
 79 integer,allocatable :: iperm(:),srch(:)
 80 real(dp) :: tsec(2)
 81 real(dp) :: inf_rmet(3,3)
 82 ! *************************************************************************
 83 
 84  call timab(602,1,tsec)
 85 
 86  ngfftrec(:) = ngfft(:) 
 87 
 88  if (recrcut>tol14) then  !default value dtset%recrcut = zero means no troncation 
 89    rtroncat = recrcut+delta
 90    get_ngfftrec(:)=1
 91    plane = 1 
 92 
 93    do ii=1,3
 94      inf_rmet(ii,:) = rmet(ii,:)/(real(ngfft(1:3)*ngfft(ii),dp))
 95    end do
 96 
 97 
 98 !  minimum value of ngfftrec
 99    do ii = 1,3
100      ngfftrec(ii)=floor(2*rtroncat/sqrt(inf_rmet(ii,ii)))+1  !minimum value
101      if(ngfftrec(ii)>=ngfft(ii))then
102        ngfftrec(ii)=ngfft(ii)
103        get_ngfftrec(ii)=0
104      end if
105    end do
106 
107 
108    if(sum(get_ngfftrec)/=0)then
109      largest_ngfftrec=maxval(ngfft(1:3))
110      maxpow2=int(log(largest_ngfftrec+0.5d0)/log(two))
111      maxpow3=int(log(largest_ngfftrec+0.5d0)/log(three))
112      maxpow5=int(log(largest_ngfftrec+0.5d0)/log(five))
113      maxpow7=0
114      maxpow11=0
115      mmsrch=(maxpow2+1)*(maxpow3+1)*(maxpow5+1)*(maxpow7+1)*(maxpow11+1)
116      ABI_ALLOCATE(srch,(mmsrch))
117      ABI_ALLOCATE(iperm,(mmsrch))
118 !    Factors of 2
119      srch(1)=1
120      do ii=1,maxpow2
121        srch(ii+1)=srch(ii)*2
122      end do
123 !    Factors of 3
124      index=maxpow2+1
125      if(maxpow3>0)then
126        do ii=1,maxpow3
127          srch(1+ii*index:(ii+1)*index)=3*srch(1+(ii-1)*index:ii*index)
128        end do
129      end if
130 !    Factors of 5
131      index=(maxpow3+1)*index
132      if(maxpow5>0)then
133        do ii=1,maxpow5
134          srch(1+ii*index:(ii+1)*index)=5*srch(1+(ii-1)*index:ii*index)
135        end do
136      end if
137 !    Factors of 7
138      index=(maxpow5+1)*index
139      if(maxpow7>0)then
140        do ii=1,maxpow7
141          srch(1+ii*index:(ii+1)*index)=7*srch(1+(ii-1)*index:ii*index)
142        end do
143      end if
144 !    Factors of 11
145      index=(maxpow7+1)*index
146      if(maxpow11>0)then
147        do ii=1,maxpow11
148          srch(1+ii*index:(ii+1)*index)=11*srch(1+(ii-1)*index:ii*index)
149        end do
150      end if
151 !    srch is the set of allowed ngfftrec values
152 
153      call sort_int(mmsrch,srch,iperm)
154      ABI_DEALLOCATE(iperm)
155 
156      do ii=1,3
157        if(get_ngfftrec(ii)==1)then
158          do jj=1,mmsrch
159            if(srch(jj)>=ngfftrec(ii))then
160              minsrch(ii)=jj
161              ngfftrec(ii)=srch(jj)
162              exit
163            end if
164          end do
165          do jj=minsrch(ii),mmsrch
166            if(srch(jj)>ngfft(ii))then 
167 !            since ngfftrec(ii)<ngfft(ii) for get_ngfftrec(ii)==1, 
168 !            and srch(mmsrch)maxval(ngfft(1:3)), 
169 !            that will appens in the range minsrch(ii),mmsrch
170              maxsrch(ii)=jj-1
171              exit
172            end if
173          end do
174        end if
175 !      since ngfft(ii) is in srch, we have here srch(maxsrch(ii))=ngfft(ii)
176 !      minsrch(ii), maxsrch(ii) is the range of index of srch in which we can 
177 !      search ngfftrec(ii)
178 
179        if(ngfftrec(ii)>=ngfft(ii))then
180          ngfftrec(ii)=ngfft(ii)
181          get_ngfftrec(ii)=0
182        end if
183      end do
184    end if
185 
186 !  verify that the entiere truncation sphere is in the fft box ; 
187 !  but only in the dimension in which we do not consider the entiere fft box
188    do while(sum(get_ngfftrec)/=0)  !again...
189 
190 !    determining the minimum distance between 0 and the boundary 
191 !    of the fft box
192 !    quite similar to the subroutine "bound", but only over the plane which 
193 !    are not the whole fft box
194      dsqmin=dsq_rec(ngfftrec(1)/2,-ngfftrec(2)/2,-ngfftrec(3)/2,inf_rmet)+0.01d0
195 
196      if(get_ngfftrec(1)/=0)then
197 !      look at +/- g1 planes:
198        do jj=-ngfftrec(2)/2,ngfftrec(2)/2
199          do kk=-ngfftrec(3)/2,ngfftrec(3)/2
200            dsp = dsq_rec(ngfftrec(1)/2, jj, kk,inf_rmet)
201            dsm = dsq_rec( - ngfftrec(1)/2, jj, kk,inf_rmet)
202            if (dsp<dsqmin) then
203              dsqmin = dsp
204              iimin = ngfftrec(1)/2
205              jjmin = jj
206              kkmin = kk
207              plane=1
208            end if
209            if (dsm<dsqmin) then
210              dsqmin = dsm
211              iimin =  - ngfftrec(1)/2
212              jjmin = jj
213              kkmin = kk
214              plane=1
215            end if
216          end do
217        end do
218      end if
219 
220      if(get_ngfftrec(2)/=0)then
221 !      +/- g2 planes:
222        do ii=-ngfftrec(1)/2,ngfftrec(1)/2
223          do kk=-ngfftrec(3)/2,ngfftrec(3)/2
224            dsp = dsq_rec(ii,ngfftrec(2)/2,kk,inf_rmet)
225            dsm = dsq_rec(ii,-ngfftrec(2)/2,kk,inf_rmet)
226            if (dsp<dsqmin) then
227              dsqmin = dsp
228              iimin = ii
229              jjmin = ngfftrec(2)/2
230              kkmin = kk
231              plane=2
232            end if
233            if (dsm<dsqmin) then
234              dsqmin = dsm
235              iimin = ii
236              jjmin =  - ngfftrec(2)/2
237              kkmin = kk
238              plane=2
239            end if
240          end do
241        end do
242      end if
243 
244      if(get_ngfftrec(3)/=0)then
245 !      +/- g3 planes:
246        do ii=-ngfftrec(1)/2,ngfftrec(1)/2
247          do jj=-ngfftrec(2)/2,ngfftrec(2)/2
248            dsp = dsq_rec(ii,jj,ngfftrec(3)/2,inf_rmet)
249            dsm = dsq_rec(ii,jj,-ngfftrec(3)/2,inf_rmet)
250            if (dsp<dsqmin) then
251              dsqmin = dsp
252              iimin = ii
253              jjmin = jj
254              kkmin = ngfftrec(3)/2
255              plane=3
256            end if
257            if (dsm<dsqmin) then
258              dsqmin = dsm
259              iimin = ii
260              jjmin = jj
261              kkmin =  - ngfftrec(3)/2
262              plane=3
263            end if
264          end do
265        end do
266      end if
267 
268      if(dsqmin>=rtroncat)then
269        get_ngfftrec=0
270        exit
271      end if
272 
273 !    Fix nearest boundary
274      do ii=minsrch(plane),maxsrch(plane)
275        if (srch(ii)>=ngfftrec(plane)) then
276 !        redefine ngfft(plane) to next higher choice
277          ngfftrec(plane)=srch(ii+1)
278 !        verify if we cover the whole box
279          if(ngfftrec(plane)>=ngfft(plane))then
280            ngfftrec(plane)=ngfft(plane)
281            get_ngfftrec(plane)=0
282          end if
283 !        Exit the loop over ii
284          exit
285        end if
286      end do
287 
288    end do
289 
290    if (allocated(srch)) then
291      ABI_DEALLOCATE(srch)
292    end if
293 
294 !  if(mod(ngfftrec(1),16)/=0) then
295 !  ngfftrec(1) = ngfftrec(1)+(16-mod(ngfftrec(1),16))
296 !  ngfftrec(2:3) = ngfftrec(1)
297 !  endif
298 
299    ngfftrec(4)=2*(ngfftrec(1)/2)+1
300    ngfftrec(5)=2*(ngfftrec(2)/2)+1
301    ngfftrec(6)=ngfftrec(3)
302 
303 !  --algorithm
304    ngfftrec(7)=ngfft(7)   ! to be improved for a better non-parallel algorithm - here it is automatically 401
305    ngfftrec(8)=ngfft(8)
306 
307  end if
308 
309 !--For now, recursion method doesn't use paralelism on FFT - which would require a great number of processors 
310  nfftrec = product(ngfftrec(1:3))
311  ngfftrec(9:11) = (/0,1,0/)   !--(/ paral, nproc, %me \)
312  ngfftrec(12:13) = ngfftrec(2:3)   ! n2proc ! n3proc
313 
314  tronc  = all(ngfftrec(:3)/=ngfft(:3))
315  call timab(602,2,tsec)
316 
317  contains
318 
319    function dsq_rec(ii,jj,kk,inf_rmet)
320 
321 
322 !This section has been created automatically by the script Abilint (TD).
323 !Do not modify the following lines by hand.
324 #undef ABI_FUNC
325 #define ABI_FUNC 'dsq_rec'
326 !End of the abilint section
327 
328    real(dp) :: dsq_rec
329    integer,intent(in) :: ii,jj,kk
330    real(dp),intent(in) :: inf_rmet(3,3)
331    dsq_rec=sqrt(&
332 &   inf_rmet(1,1)*dble(ii**2)&
333 &   +inf_rmet(2,2)*dble(jj**2)&
334 &   +inf_rmet(3,3)*dble(kk**2)&
335 &   +two*(inf_rmet(1,2)*dble(ii*jj)&
336 &   +inf_rmet(2,3)*dble(jj*kk)&
337 &   +inf_rmet(3,1)*dble(kk*ii)))
338  end function dsq_rec
339 
340 
341 end subroutine getngrec