TABLE OF CONTENTS
ABINIT/getngrec [ 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