TABLE OF CONTENTS
ABINIT/sphereboundary [ Functions ]
NAME
sphereboundary
FUNCTION
Finds the boundary of the basis sphere of G vectors (at a given k point) for use in improved zero padding of ffts in 3 dimensions. Provides data to be used by subroutine fourwf, in the form of an array gbound(2*mgfft+8,2). The first component (for use when mod(fftalg,10)==2)) provides integer values g1min,g1max,g2min,g2max and then for g2 in the sequence g2=0,1,2,...,g2max,g2min,g2min+1,...,-1, provides g1min, g1max. The second component (for use when mod(fftalg,10)==1)) provides integer values g1min,g1max,g3min,g3max, where g3min and g3max have been corrected in case of time-reversal and then for g3 in the sequence g3=0,1,2,...,g3max,g3min,g3min+1,...,-1, provides g2min, g2max. (also corrected in case of time-reversal) These values are stored in the above order in array gbound. Debug mode, if fftalg is between 000 and 099
COPYRIGHT
Copyright (C) 2002-2018 ABINIT group (DCA, XG) 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
istwf_k=option parameter that describes the storage of wfs kg_k(3,npw)=integer coordinates of G vectors in basis sphere mgfft=maximum size of 1D FFTs (only for dimensioning purposes) npw=number of G vectors in basis at this k point
OUTPUT
gbound(2*mgfft+8,2)=defined above
PARENTS
dfpt_eltfrkin,dfpt_mkrho,fock_getghc,m_bandfft_kpt,m_cut3d,m_epjdos m_fft,m_fft_prof,m_fock,m_gsphere,m_hamiltonian,m_wfd,mkrho,mlwfovlp pawmkaewf,posdoppler,scfcv,spin_current,suscep_stat,susk,tddft,wfconv
CHILDREN
SOURCE
52 #if defined HAVE_CONFIG_H 53 #include "config.h" 54 #endif 55 56 #include "abi_common.h" 57 58 59 subroutine sphereboundary(gbound,istwf_k,kg_k,mgfft,npw) 60 61 use defs_basis 62 use m_errors 63 use m_profiling_abi 64 65 !This section has been created automatically by the script Abilint (TD). 66 !Do not modify the following lines by hand. 67 #undef ABI_FUNC 68 #define ABI_FUNC 'sphereboundary' 69 !End of the abilint section 70 71 implicit none 72 73 !Arguments ------------------------------------ 74 !scalars 75 integer,intent(in) :: istwf_k,mgfft,npw 76 !arrays 77 integer,intent(in) :: kg_k(3,npw) 78 integer,intent(out) :: gbound(2*mgfft+8,2) 79 80 !Local variables------------------------------- 81 !scalars 82 integer :: dim_a,dim_b,fftalgc,g_a,gmax_a,gmax_b,gmax_b1,gmax_b2,gmin_a,gmin_b 83 integer :: gmin_b1,gmin_b2,igb,ii,iloop,ipw,testm,testp,kgk 84 character(len=500) :: message 85 !arrays 86 integer :: gmax(2),gmin(2) 87 88 ! ************************************************************************* 89 ! 90 !DEBUG 91 !write(std_out,*)' sphereboundary : enter' 92 !write(std_out, '(a)' )' sphereboundary : list of plane waves coordinates for k point ' 93 !write(std_out, '(a)' )' ipw kg_k(1:3,ipw) ' 94 !do ipw=1,npw 95 !write(std_out, '(i10,a,3i6)' )ipw,' ',kg_k(1:3,ipw) 96 !end do 97 !gbound=-999 98 !ENDDEBUG 99 100 !Determine cube boundaries 101 gbound(1,1)=minval(kg_k(1,:)) 102 gbound(2,1)=maxval(kg_k(1,:)) 103 gbound(1:2,2)=gbound(1:2,1) 104 105 !Treat differently the fftalgc cases 106 do ii=1,2 107 108 fftalgc=3-ii 109 110 if(fftalgc/=2)then 111 dim_a=3 112 dim_b=2 113 else 114 dim_a=2 115 dim_b=1 116 end if 117 118 ! Relevant boundaries 119 gbound(3,ii)=minval(kg_k(dim_a,:)) 120 gbound(4,ii)=maxval(kg_k(dim_a,:)) 121 gmin_a=gbound(3,ii) 122 gmax_a=gbound(4,ii) 123 124 ! Must complete the sphere for fftalgc==1 and special storage modes. 125 ! Explanation : sg_fftpad is not able to take into account 126 ! the time-reversal symmetry, so that the boundaries will not be delimited 127 ! by the kg_k set, but by their symmetric also. 128 if(istwf_k>=2 .and. fftalgc==1)then 129 if( istwf_k==2 .or. istwf_k==3 .or. istwf_k==6 .or. istwf_k==7 )then 130 gbound(4,2)=max(gmax_a,-gmin_a) 131 gbound(3,2)=-gbound(4,2) 132 else if( istwf_k==4 .or. istwf_k==5 .or. istwf_k==8 .or. istwf_k==9 )then 133 gbound(4,2)=max(gmax_a,-gmin_a-1) 134 gbound(3,2)=-gbound(4,2)-1 135 end if 136 gmax_a=gbound(4,2) ; gmin_a=gbound(3,2) 137 end if 138 139 igb=5 140 141 ! Consider first every g_a in range 0 ... gmax_a, then gmin_a ... -1 142 gmin(1)=0 ; gmax(1)=gmax_a 143 gmin(2)=gmin_a ; gmax(2)=-1 144 145 do iloop=1,2 146 147 if( gmin(iloop) <= gmax(iloop) )then 148 149 do g_a=gmin(iloop),gmax(iloop) 150 151 if(istwf_k==1 .or. fftalgc/=1)then 152 ! Select the minimal and maximal values, in the selected plane 153 gmin_b=mgfft+1 ! Initialized with a value larger than all possible ones 154 gmax_b=-mgfft-1 ! Initialized with a value smaller than all possible ones 155 do ipw=1,npw 156 if(kg_k(dim_a,ipw)==g_a)then 157 kgk=kg_k(dim_b,ipw) 158 if(kgk<=gmin_b)gmin_b=kgk 159 if(kgk>=gmax_b)gmax_b=kgk 160 end if 161 end do 162 163 else if(istwf_k>=2 .and. fftalgc==1)then 164 165 ! Here, must take into account time-reversal symmetry explicitely 166 167 ! Determine the boundaries for the plane g_a 168 testp=0 169 if(g_a<=gmax_a)then 170 ! Select the minimal and maximal values, in the selected plane 171 gmin_b1=mgfft+1 ! Initialized with a value larger than all possible ones 172 gmax_b1=-mgfft-1 ! Initialized with a value smaller than all possible ones 173 do ipw=1,npw 174 if(kg_k(dim_a,ipw)==g_a)then 175 kgk=kg_k(dim_b,ipw) 176 if(kgk<=gmin_b1)gmin_b1=kgk 177 if(kgk>=gmax_b1)gmax_b1=kgk 178 end if 179 end do 180 181 182 testp=1 183 end if 184 185 ! Determine the boundaries for the plane -g_a or -g_a-1 186 testm=0 187 if( istwf_k==2 .or. istwf_k==3 .or. istwf_k==6 .or. istwf_k==7 )then 188 189 if(-g_a>=gmin_a)then 190 ! Select the minimal and maximal values, in the selected plane 191 ! Warning : there is an inversion of search (might be confusing) 192 gmax_b2=mgfft+1 ! Initialized with a value larger than all possible ones 193 gmin_b2=-mgfft-1 ! Initialized with a value smaller than all possible ones 194 do ipw=1,npw 195 if(kg_k(dim_a,ipw)==-g_a)then 196 kgk=kg_k(dim_b,ipw) 197 if(kgk<=gmax_b2)gmax_b2=kgk 198 if(kgk>=gmin_b2)gmin_b2=kgk 199 end if 200 end do 201 testm=1 202 end if 203 204 else if( istwf_k==4 .or. istwf_k==5 .or. istwf_k==8 .or. istwf_k==9 )then 205 206 if(-g_a-1>=gmin_a)then 207 ! Select the minimal and maximal values, in the selected plane 208 ! Warning : there is an inversion of search (might be confusing) 209 gmax_b2=mgfft+1 ! Initialized with a value larger than all possible ones 210 gmin_b2=-mgfft-1 ! Initialized with a value smaller than all possible ones 211 do ipw=1,npw 212 if(kg_k(dim_a,ipw)==-g_a-1)then 213 kgk=kg_k(dim_b,ipw) 214 if(kgk<=gmax_b2)gmax_b2=kgk 215 if(kgk>=gmin_b2)gmin_b2=kgk 216 end if 217 end do 218 testm=1 219 end if 220 221 end if 222 223 ! Must invert the boundaries, to use them for plane g_a 224 if(testm==1)then 225 ! This is needed to avoid border effect 226 ! if the search did not lead to any element 227 gmin_b2=max(gmin_b2,-mgfft) ; gmax_b2=min(gmax_b2,mgfft) 228 if(istwf_k<=5)then 229 gmax_b2=-gmax_b2 ; gmin_b2=-gmin_b2 230 else 231 gmax_b2=-gmax_b2-1 ; gmin_b2=-gmin_b2-1 232 end if 233 end if 234 235 if( testp==1 .and. testm==1)then 236 gmin_b=min(gmin_b1,gmin_b2) ; gmax_b=max(gmax_b1,gmax_b2) 237 else if( testp==1 )then 238 gmin_b=gmin_b1 ; gmax_b=gmax_b1 239 else if( testm==1 )then 240 gmin_b=gmin_b2 ; gmax_b=gmax_b2 241 end if 242 243 end if ! Endif take into account time-reversal symmetry 244 245 if (igb+1>2*mgfft+4) then 246 write(message, '(2a, 4(a,3(i0,1x),a))' )& 247 "About to overwrite gbound array (FFT mesh too small) ",ch10, & 248 " iloop, igb, mgb = ",iloop,igb,2*mgfft+4, ch10, & 249 " istwfk, mgfft, npw = ",istwf_k, mgfft, npw, ch10, & 250 " minval(kg_k) = ",minval(kg_k, dim=2), ch10, & 251 " maxval(kg_k) = ",maxval(kg_k, dim=2), ch10 252 MSG_BUG(message) 253 end if 254 255 gbound(igb,ii)=gmin_b 256 gbound(igb+1,ii)=gmax_b 257 258 if( iloop==1 .and. istwf_k>=2 .and. istwf_k<=5 .and. fftalgc==2 .and. g_a==0)then 259 ! If k_y=0 , for fftalgc==2, the g_a==0 plane must be completed 260 if(istwf_k==2 .or. istwf_k==4)then 261 gbound(igb+1,ii)=max(gmax_b,-gmin_b) 262 gbound(igb,ii)=-gbound(igb+1,ii) 263 else if(istwf_k==3 .or. istwf_k==5)then 264 gbound(igb+1,ii)=max(gmax_b,-gmin_b-1) 265 gbound(igb,ii)=-gbound(igb+1,ii)-1 266 end if 267 268 end if 269 270 igb=igb+2 271 272 end do ! g_a 273 end if 274 end do ! iloop 275 end do ! ii (fftalgc) 276 277 !DEBUG 278 !write(std_out,'(a)')' sphereoundary : list of plane waves coordinates for 1st k point ' 279 !write(std_out,'(a)')' ipw kg_k(1:3,ipw) ' 280 !do ipw=1,npw 281 !write(std_out, '(i10,a,3i6)' )ipw,' ',kg_k(1:3,ipw) 282 !end do 283 !write(std_out, '(a)' )' sphereboundary : list of boundaries ' 284 !do igb=1,2*mgfft+8 285 !write(std_out, '(i10,a,2i6)' )igb,' ',gbound(igb,1),gbound(igb,2) 286 !end do 287 !write(std_out,*)' sphereboundary : exit ' 288 !ENDDEBUG 289 290 end subroutine sphereboundary