TABLE OF CONTENTS
ABINIT/mlwfovlp_pw [ Functions ]
NAME
mlwfovlp_pw
FUNCTION
Routine which computes PW part of overlap M_{mn}(k,b) for Wannier code (www.wannier.org f90 version).
COPYRIGHT
Copyright (C) 2005-2018 ABINIT group (BAmadon,FJollet,T Rangel) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions. g1(3,nkpt,nntot) = G vector shift which is necessary to obtain k1+b iwav(mband,nkpt,nsppol): shift for pw components in cg. kg(3,mpw*mkmem)=reduced planewave coordinates. mband=maximum number of bands mgfft=maximum size of 1D FFTs mkmem =number of k points treated by this node. mpi_enreg=informations about MPI parallelization mpw=maximum dimensioned size of npw. nfft=(effective) number of FFT grid points (for this processor) (see NOTES at beginning of scfcv) ngfft(18)=contain all needed information about 3D FFT (see NOTES at beginning of scfcv) nkpt=number of k points. npwarr(nkpt)=number of planewaves in basis at this k point nspinor=number of spinorial components of the wavefunctions nsppol=1 for unpolarized, 2 for spin-polarized ovikp(nkpt,nntot)= gives nntot value of k2 (in the BZ) for each k1 (k2=k1+b mod(G)) seed_name= seed_name of files containing cg for all k-points to be used with MPI spin = just used for nsppol>1 ; 0 both, 1 just spin up, 2 just spin down
OUTPUT
cm1(2,mband,mband,nntot,nkpt,nsppol): overlap <u_(nk1)|u_(mk1+b)>.
SIDE EFFECTS
(only writing, printing)
NOTES
PARENTS
mlwfovlp
CHILDREN
wrtout,xmpi_barrier,xmpi_sum
SOURCE
53 #if defined HAVE_CONFIG_H 54 #include "config.h" 55 #endif 56 57 #include "abi_common.h" 58 59 subroutine mlwfovlp_pw(cg,cm1,g1,iwav,kg,mband,mkmem,mpi_enreg,mpw,nfft,ngfft,nkpt,nntot,& 60 & npwarr,nspinor,nsppol,ovikp,seed_name,spin) 61 62 use defs_basis 63 use defs_abitypes 64 use m_profiling_abi 65 use m_errors 66 use m_xmpi 67 68 !This section has been created automatically by the script Abilint (TD). 69 !Do not modify the following lines by hand. 70 #undef ABI_FUNC 71 #define ABI_FUNC 'mlwfovlp_pw' 72 use interfaces_14_hidewrite 73 !End of the abilint section 74 75 implicit none 76 77 !Arguments ------------------------------------ 78 !scalars 79 integer,intent(in) :: mband,mkmem,mpw,nfft,nkpt,nntot 80 integer,intent(in) :: nspinor,nsppol,spin 81 character(len=fnlen) :: seed_name !seed names of files containing cg info used in case of MPI 82 type(MPI_type),intent(in) :: mpi_enreg 83 !arrays 84 integer,intent(in) :: g1(3,nkpt,nntot),kg(3,mpw*mkmem),ngfft(18),npwarr(nkpt) 85 integer,intent(in) :: iwav(mband,nkpt,nsppol) 86 integer,intent(in) :: ovikp(nkpt,nntot) 87 real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol) 88 real(dp),intent(out) :: cm1(2,mband,mband,nntot,nkpt,nsppol) 89 90 !Local variables------------------------------- 91 !scalars 92 integer :: iband1,iband2,ierr,ig,ig1,ig1b,ig2,ig2b 93 integer :: ig3,ig3b,igk1,igk2,igks1,igks2,ii,ikg,ikpt,ikpt1,ikpt2,imntot,index,intot,ios,ipw 94 integer :: ispinor,isppol,iunit,me,n1,n2,n3,npoint,npoint2,npw_k,npw_k2 95 integer :: nprocs,spaceComm 96 integer,allocatable::indpwk(:,:),kg_k(:,:) 97 integer,allocatable :: invpwk(:,:) 98 character(len=500) :: message 99 character(len=fnlen) :: cg_file !file containing cg info used in case of MPI 100 logical::lfile 101 real(dp),allocatable :: cg_read(:,:) !to be used in case of MPI 102 103 !************************************************************************ 104 105 write(message, '(a,a)' ) ch10,& 106 & '** mlwfovlp_pw : compute pw part of overlap' 107 call wrtout(std_out, message,'COLL') 108 109 !initialize flags 110 lfile=.false. 111 !mpi initialization 112 spaceComm=MPI_enreg%comm_cell 113 nprocs=xmpi_comm_size(spaceComm) 114 me=MPI_enreg%me_kpt 115 116 if(nprocs>1) then 117 ABI_ALLOCATE(cg_read,(2,nspinor*mpw*mband)) 118 end if 119 120 121 !****************compute intermediate quantities (index, shifts) ****** 122 !------------compute index for g points-------------------------------- 123 !ig is a plane waves which belongs to the sphere ecut for ikpt (they 124 !are npwarr(ikpt)) 125 !npoint is the position in the grid of planes waves 126 !(they are nfft) 127 !indpwk is a application ig-> npoint 128 !invpwk is not an application (some npoint have no ig corresponding) 129 !cg are ordered with npw_k ! 130 !---------------------------------------------------------------------- 131 !------------compute index for g points-------------------------------- 132 !---------------------------------------------------------------------- 133 write(message, '(a,a)' ) ch10,& 134 & ' first compute index for g-points' 135 call wrtout(std_out, message,'COLL') 136 ! 137 !Allocations 138 ABI_ALLOCATE(kg_k,(3,mpw)) 139 ABI_ALLOCATE(indpwk,(nkpt,mpw)) 140 ABI_ALLOCATE(invpwk,(nkpt,nfft)) 141 ! 142 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 143 invpwk=0 144 indpwk=0 145 kg_k=0 146 do isppol=1,1 !invpwk is not spin dependent 147 ! so we just do it once 148 ikg=0 149 do ikpt=1,nkpt 150 ! 151 ! MPI:cycle over k-points not treated by this node 152 ! 153 if ( ABS(MPI_enreg%proc_distrb(ikpt,1,isppol)-me) /=0) CYCLE 154 155 ! 156 ! write(std_out,*)'me',me,'ikpt',ikpt,'isppol',isppol 157 do npoint=1,nfft 158 if(invpwk(ikpt,npoint)/=0 )then 159 write(std_out,*) "error0 , invpwk is overwritten" 160 write(std_out,*) ikpt,npoint 161 MSG_ERROR("Aborting now") 162 end if 163 end do 164 npw_k=npwarr(ikpt) 165 ! write(std_out,*) ikpt,npw_k,nfft 166 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 167 do ig=1,npw_k 168 if(ig.gt.mpw) then 169 write(std_out,*)"error ig",ig,"greater than mpw ",mpw 170 MSG_ERROR("Aborting now") 171 end if 172 if(indpwk(ikpt,ig)/=0) then 173 write(std_out,*) "error, indpwk is overwritten" 174 write(std_out,*) ikpt,ig,indpwk(ikpt,ig) 175 MSG_ERROR("Aborting now") 176 end if 177 ig1=modulo(kg_k(1,ig),n1) 178 ig2=modulo(kg_k(2,ig),n2) 179 ig3=modulo(kg_k(3,ig),n3) 180 indpwk(ikpt,ig)=ig1+1+n1*(ig2+n2*ig3) 181 npoint=indpwk(ikpt,ig) 182 if(npoint.gt.nfft) then 183 MSG_ERROR("error npoint") 184 end if 185 ! write(std_out,*) ikpt,ig,npoint,invpwk(ikpt,npoint) 186 if(invpwk(ikpt,npoint)/=0) then 187 write(std_out,*) "error, invpwk is overwritten" 188 write(std_out,*) ikpt,ig,npoint,invpwk(ikpt,npoint) 189 MSG_ERROR("Aborting now") 190 end if 191 invpwk(ikpt,npoint)=ig 192 ! write(std_out,*)'ikpt,npoint,invpwk',ikpt,npoint,invpwk(ikpt,npoint) 193 ! if(ikpt.eq.1) write(std_out,*) "ig npoint",ig, npoint 194 ! write(std_out,*) "ikpt ig npoint",ikpt,ig, npoint 195 end do 196 ikg=ikg+npw_k 197 198 end do !ikpt 199 end do !isppol 200 !write(std_out,*) "index for g points has been computed" 201 202 call xmpi_barrier(spaceComm) 203 call xmpi_sum(invpwk,spaceComm,ierr) 204 205 !---------------------------------------------------------------------- 206 !------------test invpwk----------------------------------------------- 207 !---------------------------------------------------------------------- 208 !write(std_out,*) "TEST INVPWK" 209 !ikpt=3 210 !isppol=1 211 !do ig=1,npwarr(ikpt) 212 !npoint=indpwk(ikpt,ig) 213 !write(std_out,*) "ig npoint ",ig, npoint 214 !write(std_out,*) "ig npoint inv",invpwk(ikpt,npoint),npoint 215 !end do 216 !do ig3=1,n3 217 !do ig2=1,n2 218 !do ig1=1,n1 219 !npoint=ig1+(ig2-1)*n1+(ig3-1)*n2*n1 220 !ig=invpwk(ikpt,npoint) 221 !! if(ig/=0) write(std_out,*) "ig npoint",ig, npoint 222 !end do 223 !end do 224 !end do 225 226 227 228 229 ! 230 !Deallocate unused variables 231 ! 232 ABI_DEALLOCATE(kg_k) 233 ABI_DEALLOCATE(indpwk) 234 235 236 !*********************************************************************** 237 !**calculate overlap M_{mn}(k,b)=<\Psi_{k,m}|e^{-ibr}|\Psi_{k+b,n}>***** 238 !*********************************************************************** 239 write(message, '(a,a)' ) ch10,& 240 & ' mlwfovlp_pw : compute overlaps ' 241 call wrtout(std_out, message,'COLL') 242 write(message, '(a,a)' ) ch10,& 243 & " nkpt nntot mband " 244 call wrtout(std_out, message,'COLL') 245 write(message, '(i6,2x,i6,2x,i6,2x,i6)' ) & 246 & nkpt,nntot,mband 247 call wrtout(std_out, message,'COLL') 248 cm1=zero 249 write(message, '(a)' ) ' ' 250 call wrtout(std_out, message,'COLL') 251 do isppol=1,nsppol 252 if(spin.ne.0 .and. spin.ne.isppol) cycle 253 imntot=0 254 do ikpt1=1,nkpt 255 ! 256 ! MPI:cycle over k-points not treated by this node 257 ! 258 if ( ABS(MPI_enreg%proc_distrb(ikpt1,1,isppol)-me) /=0) CYCLE 259 ! 260 write(message, '(a,i6,a,i6,a,i6)' ) & 261 & ' Processor',me,' computes k-point',ikpt1,' and spin=',isppol 262 call wrtout(std_out, message,'COLL') 263 ! write(std_out,*)trim(message) 264 265 do intot=1,nntot 266 lfile=.false. !flag to know if this kpt will be read from a file, see below 267 imntot=imntot+1 268 ikpt2= ovikp(ikpt1,intot) 269 ! write(std_out,*)'me',me,'ikpt1',ikpt1,'ikpt2',ikpt2,'intot',intot,'isppol',isppol 270 271 ! 272 ! MPI: if ikpt2 not found in this processor then 273 ! read info from an unformatted file 274 ! 275 if ( ABS(MPI_enreg%proc_distrb(ikpt2,1,isppol)-me) /=0) then 276 lfile=.true. 277 write(cg_file,'(a,I5.5,".",I1)') trim(seed_name),ikpt2,isppol 278 iunit=1000+ikpt2+ikpt2*(isppol-1) 279 npw_k2=npwarr(ikpt2) 280 open (unit=iunit, file=cg_file,form='unformatted',status='old',iostat=ios) 281 if(ios /= 0) then 282 write(message,*) " mlwfovlp_pw: file",trim(cg_file), "not found" 283 MSG_ERROR(message) 284 end if 285 ! 286 do iband2=1,mband 287 do ipw=1,npw_k2*nspinor 288 index=ipw+(iband2-1)*npw_k2*nspinor 289 read(iunit) (cg_read(ii,index),ii=1,2) 290 ! if(me==0 .and. ikpt2==4)write(300,*)'ipw,iband2,index',ipw,iband2,index,cg_read(:,index) 291 ! if(me==1 .and. ikpt2==4)write(301,*)'ipw,iband2,index',ipw,iband2,index,cg_read(:,index) 292 end do 293 end do 294 close(iunit) 295 end if 296 ! 297 npw_k=npwarr(ikpt1) 298 npw_k2=npwarr(ikpt2) 299 do ig3=1,n3 300 do ig2=1,n2 301 do ig1=1,n1 302 ! write(std_out,*) isppol,ikpt1,iband1,iband2,intot 303 npoint=ig1+(ig2-1)*n1+(ig3-1)*n2*n1 304 if(npoint.gt.nfft) then 305 write(std_out,*) "error npoint b" 306 MSG_ERROR("Aborting now") 307 end if 308 ig1b=ig1+g1(1,ikpt1,intot) 309 ig2b=ig2+g1(2,ikpt1,intot) 310 ig3b=ig3+g1(3,ikpt1,intot) 311 ! write(std_out,*) ig1,ig2,ig3 312 ! write(std_out,*) ig1b,ig2b,ig3b 313 if(ig1b.lt.1) ig1b=ig1b+n1 314 if(ig2b.lt.1) ig2b=ig2b+n2 315 if(ig3b.lt.1) ig3b=ig3b+n3 316 if(ig1b.gt.n1) ig1b=ig1b-n1 317 if(ig2b.gt.n2) ig2b=ig2b-n2 318 if(ig3b.gt.n3) ig3b=ig3b-n3 319 npoint2=ig1b+(ig2b-1)*n1+(ig3b-1)*n2*n1 320 if(npoint2.gt.nfft) then 321 write(std_out,*)"error npoint c" 322 MSG_ERROR("Aborting now") 323 end if 324 igk1=invpwk(ikpt1,npoint) 325 igk2=invpwk(ikpt2,npoint2) 326 327 ! if(intot==10) write(std_out,*)'Before igk1 and igk2',ikpt1,ikpt2,isppol 328 329 if(igk1/=0.and.igk2/=0) then 330 do iband2=1,mband 331 do iband1=1,mband 332 do ispinor=1,nspinor 333 ! igks1= (igk1*nspinor)-(nspinor-ispinor) 334 ! igks2= (igk2*nspinor)-(nspinor-ispinor) 335 igks1= igk1+ (ispinor-1)*npw_k 336 igks2= igk2+ (ispinor-1)*npw_k2 337 338 ! Here the igks is to include the spinor component missing in igk 339 if(lfile) index=igks2+npw_k2*nspinor*(iband2-1) !In case of MPI, see below 340 ! 341 ! If MPI sometimes the info was read from an unformatted file 342 ! If that is the case lfile==.true. 343 ! 344 if(lfile) then 345 cm1(1,iband1,iband2,intot,ikpt1,isppol)=cm1(1,iband1,iband2,intot,ikpt1,isppol)+ & 346 & cg(1,igks1+iwav(iband1,ikpt1,isppol))*cg_read(1,index)& 347 & + cg(2,igks1+iwav(iband1,ikpt1,isppol))*cg_read(2,index) 348 cm1(2,iband1,iband2,intot,ikpt1,isppol)=cm1(2,iband1,iband2,intot,ikpt1,isppol)+ & 349 & cg(1,igks1+iwav(iband1,ikpt1,isppol))*cg_read(2,index)& 350 & - cg(2,igks1+iwav(iband1,ikpt1,isppol))*cg_read(1,index) 351 ! 352 else 353 cm1(1,iband1,iband2,intot,ikpt1,isppol)=cm1(1,iband1,iband2,intot,ikpt1,isppol)+ & 354 & cg(1,igks1+iwav(iband1,ikpt1,isppol))*cg(1,igks2+iwav(iband2,ikpt2,isppol))& 355 & + cg(2,igks1+iwav(iband1,ikpt1,isppol))*cg(2,igks2+iwav(iband2,ikpt2,isppol)) 356 cm1(2,iband1,iband2,intot,ikpt1,isppol)=cm1(2,iband1,iband2,intot,ikpt1,isppol)+ & 357 & cg(1,igks1+iwav(iband1,ikpt1,isppol))*cg(2,igks2+iwav(iband2,ikpt2,isppol))& 358 & - cg(2,igks1+iwav(iband1,ikpt1,isppol))*cg(1,igks2+iwav(iband2,ikpt2,isppol)) 359 end if 360 end do !ispinor 361 end do ! iband1 362 end do ! iband2 363 end if 364 end do 365 end do 366 end do 367 end do ! intot 368 end do ! ikpt1 369 end do ! isppol 370 ! 371 !Deallocations 372 ! 373 ABI_DEALLOCATE(invpwk) 374 if(nprocs>1) then 375 ABI_DEALLOCATE(cg_read) 376 end if 377 378 end subroutine mlwfovlp_pw