TABLE OF CONTENTS
ABINIT/smatrix_pawinit [ Functions ]
NAME
smatrix_pawinit
FUNCTION
Routine which computes paw part of the overlap used to compute LMWF wannier functions and berryphase
COPYRIGHT
Copyright (C) 2005-2018 ABINIT group (T Rangel,BAmadon,FJollet,PHermet) 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
atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f) cprj(natom,nspinor*mband*mkmem*nsppol)= <p_lmn|Cnk> coefficients for each WF |Cnk> and each |p_lmn> non-local projector dimcprj(natom)=array of dimensions of array cprj (not ordered) g1(3)= reciprocal vector to put k1+b inside the BZ. bb=k2-k1=b-G1 ("b" is the true b, so we have to correct bb with G1). gprimd(3,3)=dimensional reciprocal space primitive translations ikpt1(3)=cartesian coordinates of k1 ikpt2(3)=cartesian coordinates of k2 isppol = spin polarization mband=maximum number of bands mkmem =number of k points treated by this node. mpi_enreg=information about MPI parallelization natom=number of atoms in cell. nkpt=number of k points. nspinor=number of spinorial components of the wavefunctions (on current proc) nsppol=1 for unpolarized, 2 for spin-polarized ntypat=number of types of atoms in unit cell. rprimd(3,3)=dimensional primitive translations for real space (bohr) seed_name= seed_name of files containing cg for all k-points to be used with MPI xred(3,natom)=reduced dimensionless atomic coordinates
OUTPUT
cm2: Inside sphere part of the overlap needed for constructing wannier function
SIDE EFFECTS
(only writing, printing)
NOTES
The mpi part will work with mlwfovlp but not for berryphase_new
PARENTS
mlwfovlp
CHILDREN
initylmr,pawcprj_alloc,pawcprj_free,pawcprj_get,pawcprj_getdim,sbf8 simp_gen
SOURCE
57 #if defined HAVE_CONFIG_H 58 #include "config.h" 59 #endif 60 61 #include "abi_common.h" 62 63 subroutine smatrix_pawinit(atindx1,cm2,cprj,ikpt1,ikpt2,isppol,& 64 & g1,gprimd,kpt,mband,mbandw,mkmem,mpi_enreg,& 65 & natom,nband,nkpt,nspinor,nsppol,ntypat,pawang,pawrad,pawtab,rprimd,& 66 & seed_name,typat,xred) 67 68 use defs_basis 69 use defs_abitypes 70 use m_errors 71 use m_profiling_abi 72 use m_xmpi 73 74 use m_special_funcs, only : sbf8 75 use m_pawang, only : pawang_type 76 use m_pawrad, only : pawrad_type, simp_gen 77 use m_pawtab, only : pawtab_type 78 use m_paw_sphharm, only : initylmr 79 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_free, pawcprj_getdim 80 81 !This section has been created automatically by the script Abilint (TD). 82 !Do not modify the following lines by hand. 83 #undef ABI_FUNC 84 #define ABI_FUNC 'smatrix_pawinit' 85 !End of the abilint section 86 87 implicit none 88 89 !Arguments--------------------------- 90 !scalars 91 integer,intent(in) :: ikpt1,ikpt2,isppol,mband,mbandw,mkmem,natom,nkpt,nspinor,nsppol 92 integer,intent(in) :: ntypat 93 character(len=fnlen) :: seed_name !seed names of files containing cg info used in case of MPI 94 type(MPI_type),intent(in) :: mpi_enreg 95 type(pawang_type),intent(in) :: pawang 96 97 !arrays 98 integer,intent(in) :: atindx1(natom),g1(3),nband(nsppol*nkpt),typat(natom) 99 real(dp),intent(in) :: gprimd(3,3),kpt(3,nkpt),rprimd(3,3),xred(3,natom) 100 real(dp),intent(inout) :: cm2(2,mbandw,mbandw) 101 type(pawcprj_type) :: cprj(natom,nspinor*mband*mkmem*nsppol) 102 type(pawrad_type),intent(in) :: pawrad(ntypat) 103 type(pawtab_type),intent(in) :: pawtab(ntypat) 104 105 !Local variables--------------------------- 106 !scalars 107 integer :: dummy 108 integer :: iatom,iband1,iband2,icg1,icg2,idx1,idx2,ii 109 integer :: ilmn,ios,iunit,ir 110 integer :: iorder_cprj,isel,ispinor,itypat,j0lmn,jj,jlmn,klm,klmn,kln,ll,lm0,lmax 111 integer :: lmin,lmn_size,max_lmn,mesh_size,mm,nband_k 112 integer :: nprocs,spaceComm,rank !for mpi 113 real(dp) :: arg,bnorm,delta,intg,ppi,ppr,qijbtemp,qijtot,x1 114 real(dp) :: x2,xsum,xtemp,xx,yy,zz 115 character(len=500) :: message 116 character(len=fnlen) :: cprj_file !file containing cg info used in case of MPI 117 logical::lfile 118 119 !arrays 120 integer,allocatable :: dimcprj(:),nattyp_dum(:) 121 real(dp),parameter :: ili(7)=(/zero,-one,zero,one,zero,-one,zero/) 122 real(dp),parameter :: ilr(7)=(/one,zero,-one,zero,one,zero,-one/) 123 real(dp) :: bb(3),bb1(3),bbn(3),qijb(2),xcart(3,natom) 124 real(dp),allocatable :: ff(:),j_bessel(:,:),ylmb(:),ylmrgr_dum(:,:,:) 125 real(dp),allocatable :: sb_out(:) 126 type(pawcprj_type),allocatable :: cprj_k1(:,:) 127 type(pawcprj_type),allocatable :: cprj_k2(:,:) 128 129 ! ************************************************************************* 130 131 DBG_ENTER("COLL") 132 133 ! 134 !Allocate cprj_k1 and cprj_k2 135 ! 136 ABI_ALLOCATE(dimcprj,(natom)) 137 call pawcprj_getdim(dimcprj,natom,nattyp_dum,ntypat,typat,pawtab,'R') 138 139 nband_k=nband(ikpt1) 140 ABI_DATATYPE_ALLOCATE(cprj_k1,(natom,nband_k*nspinor)) 141 call pawcprj_alloc(cprj_k1,0,dimcprj) 142 143 nband_k=nband(ikpt2) 144 ABI_DATATYPE_ALLOCATE(cprj_k2,(natom,nband_k*nspinor)) 145 call pawcprj_alloc(cprj_k2,0,dimcprj) 146 ABI_DEALLOCATE(dimcprj) 147 148 !mpi initialization 149 spaceComm=MPI_enreg%comm_cell 150 nprocs=xmpi_comm_size(spaceComm) 151 rank=MPI_enreg%me_kpt 152 153 lfile=.false. 154 ! 155 !write(std_out,*) "compute PAW overlap for k-points",ikpt1,ikpt2 156 do iatom=1,natom 157 xcart(:,iatom)=rprimd(:,1)*xred(1,iatom)+& 158 & rprimd(:,2)*xred(2,iatom)+& 159 & rprimd(:,3)*xred(3,iatom) 160 end do 161 162 ! 163 !Calculate indices icg1 and icg2 164 ! 165 icg1=0 166 do ii=1,isppol 167 ll=nkpt 168 if(ii==isppol) ll=ikpt1-1 169 do jj=1,ll 170 ! MPI: cycle over kpts not treated by this node 171 if ( ABS(MPI_enreg%proc_distrb(jj,1,ii)-rank)/=0) CYCLE 172 ! write(std_out,'("kpt loop2: ikpt",i3," rank ",i3)') jj,rank 173 icg1=icg1+nspinor*nband(jj+(ii-1)*nkpt) 174 end do 175 end do 176 icg2=0 177 do ii=1,isppol 178 ll=nkpt 179 if(isppol==ii) ll=ikpt2-1 180 do jj=1,ll 181 ! MPI: cycle over kpts not treated by this node 182 if (ABS(MPI_enreg%proc_distrb(jj,1,ii)-rank)/=0) CYCLE 183 ! write(std_out,'("kpt loop2: ikpt",i3," rank ",i3)') jj,rank 184 icg2=icg2+nspinor*nband(jj+(ii-1)*nkpt) 185 end do 186 end do 187 ! 188 !MPI: if ikpt2 not found in this processor then 189 !read info from an unformatted file 190 ! 191 if (nprocs>1) then 192 if (ABS(MPI_enreg%proc_distrb(ikpt2,1,isppol)-rank)/=0) then 193 lfile=.true. 194 ! 195 ! get maximum of lmn_size 196 max_lmn=0 197 do itypat=1,ntypat 198 lmn_size=pawtab(itypat)%lmn_size 199 if(lmn_size>max_lmn) max_lmn=lmn_size 200 end do 201 ! 202 ! get file name and open it 203 ! 204 write(cprj_file,'(a,I5.5,".",I1)') trim(seed_name),ikpt2,isppol 205 iunit=1000 206 ! write(std_out,*)'reading file',trim(cprj_file) 207 open (unit=iunit, file=cprj_file,form='unformatted',status='old',iostat=ios) 208 if(ios /= 0) then 209 write(message,*) " smatrix_pawinit: file",trim(cprj_file), "not found" 210 MSG_ERROR(message) 211 end if 212 ! 213 ! start reading 214 do ii=1,mband*nspinor 215 do iatom=1,natom 216 itypat=typat(iatom) 217 lmn_size=pawtab(itypat)%lmn_size 218 do ilmn=1,lmn_size 219 read(iunit)(cprj_k2(iatom,ii)%cp(jj,ilmn),jj=1,2) 220 end do !ilmn 221 end do 222 end do 223 ! 224 ! close file 225 ! 226 close (unit=iunit,iostat=ios) 227 if(ios /= 0) then 228 write(message,*) " smatrix_pawinit: error closing file ",trim(cprj_file) 229 MSG_ERROR(message) 230 end if 231 ! 232 end if 233 end if !mpi 234 235 !Extract cprj_k1 and cprj_k2 236 !these contain the projectors cprj for just one k-point (ikpt1 or ikpt2) 237 238 !Extract cprj for k-point 1 239 iorder_cprj=0 !do not change the ordering of cprj 240 nband_k=nband(ikpt1) 241 dummy=1000 !index of file not implemented here, mkmem==0 not implemented 242 call pawcprj_get(atindx1,cprj_k1,cprj,natom,1,icg1,ikpt1,iorder_cprj,isppol,& 243 & mband,mkmem,natom,nband_k,nband_k,nspinor,nsppol,dummy,& 244 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 245 246 !Extract cprj for k-point 2 247 if( lfile .eqv. .false. ) then !if it was not already read above 248 iorder_cprj=0 !do not change the ordering of cprj 249 nband_k=nband(ikpt2) 250 dummy=1000 !index of file not implemented here, mkmem==0 not implemented 251 call pawcprj_get(atindx1,cprj_k2,cprj,natom,1,icg2,ikpt2,iorder_cprj,isppol,& 252 & mband,mkmem,natom,nband_k,nband_k,nspinor,nsppol,dummy,& 253 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 254 end if 255 256 !DEBUG 257 !if(ikpt2==2) then 258 !if(ikpt1==1) then 259 !do iband1=1,mbandw 260 !do iatom=1,natom 261 !itypat=typat(atindx1(iatom)) 262 !lmn_size=pawtab(itypat)%lmn_size 263 !do ilmn=1,1!lmn_size 264 !!write(500,'(a,i3,a,i3,a,i3,a,2f13.7)')'iband ',iband1,' iatom ',iatom,' ilmn ',ilmn,' cprj',cprj(iatom,iband1+icg1)%cp(:,ilmn) 265 !write(500,'(a,i3,a,i3,a,i3,a,2f13.7)')'iband ',iband1,' iatom ',atindx1(iatom),' ilmn ',ilmn,' cprj',cprj(atindx1(iatom),iband1+icg1)%cp(:,ilmn) 266 !write(500,'(a,i3,a,i3,a,i3,a,2f13.7)')'iband ',iband1,' iatom ',iatom,' ilmn ',ilmn,' cprj_k1 ',cprj_k1(iatom,iband1)%cp(:,ilmn) 267 !end do 268 !end do 269 !end do 270 !end if 271 !end if 272 !NDDEBUG 273 274 !!!!!!!!!!!!!!!!! 275 !--- Compute intermediate quantities: "b" vector=k2-k1 and its 276 !normalized value: bbn (and its norm: bnorm) 277 !compute also Ylm(b). 278 ABI_ALLOCATE(ylmb,(pawang%l_size_max*pawang%l_size_max)) 279 ABI_ALLOCATE(ylmrgr_dum,(1,1,0)) 280 bb(:)=kpt(:,ikpt2)-kpt(:,ikpt1)+g1(:) 281 bb1=bb 282 xx=gprimd(1,1)*bb(1)+gprimd(1,2)*bb(2)+gprimd(1,3)*bb(3) 283 yy=gprimd(2,1)*bb(1)+gprimd(2,2)*bb(2)+gprimd(2,3)*bb(3) 284 zz=gprimd(3,1)*bb(1)+gprimd(3,2)*bb(2)+gprimd(3,3)*bb(3) 285 bnorm=two_pi*dsqrt(xx**2+yy**2+zz**2) 286 if(bnorm<tol8) then 287 ! write(std_out,*) "WARNING: bnorm=",bnorm 288 bbn(:)=zero 289 else 290 xx=xx*two_pi 291 yy=yy*two_pi 292 zz=zz*two_pi 293 bb(1)=xx 294 bb(2)=yy 295 bb(3)=zz 296 bbn(1)=xx/bnorm 297 bbn(2)=yy/bnorm 298 bbn(3)=zz/bnorm 299 end if 300 301 !debug bbn=0 302 !debug bnorm=0 303 !bbn has to ne normalized 304 call initylmr(pawang%l_size_max,0,1,(/one/),1,bbn(:),ylmb(:),ylmrgr_dum) 305 !write(std_out,*) "ylmb(:)",ylmb(:) 306 !write(std_out,*) pawang%l_size_max 307 !write(std_out,*) "bbn",bbn(:) 308 !write(std_out,*) "xx,yy,zz",xx,yy,zz 309 !write(std_out,*) "bnorm",bnorm 310 ABI_DEALLOCATE(ylmrgr_dum) 311 312 !------- First Compute Qij(b)- 313 ABI_ALLOCATE(sb_out, (pawang%l_size_max)) 314 cm2=zero 315 316 do iatom=1,natom 317 itypat=typat(iatom) 318 lmn_size=pawtab(itypat)%lmn_size 319 ! --- en coordonnnes reelles cartesiennes (espace reel) 320 ! --- first radial part(see pawinit) 321 mesh_size=pawtab(itypat)%mesh_size 322 ABI_ALLOCATE(j_bessel,(mesh_size,pawang%l_size_max)) 323 324 325 ! --- compute bessel function for (br) for all angular momenta necessary 326 ! --- and for all value of r. 327 ! --- they are needed for radial part 328 ! --- of the integration => j_bessel(ir,:) 329 do ir=1,mesh_size 330 arg=bnorm*pawrad(itypat)%rad(ir) 331 call sbf8(pawang%l_size_max,arg,sb_out) 332 j_bessel(ir,:) = sb_out 333 end do 334 335 ! do jlmn=1,pawang%l_size_max 336 ! write(665,*) "j_bessel",j_bessel(1:mesh_size,jlmn) 337 ! enddo 338 ! write(std_out,*) "bessel function computed" 339 ! --- Compute \Sum b.R=xsum for future use 340 xtemp=zero 341 do mm=1,3 342 xtemp=xtemp+xred(mm,iatom)*bb1(mm) 343 end do 344 xtemp=xtemp*two_pi 345 xsum=zero 346 do mm=1,3 347 xsum=xsum+xcart(mm,iatom)*bbn(mm)*bnorm 348 end do 349 ! write(std_out,*)'xsum',xsum,xtemp,lmn_size 350 351 ! --- Loop on jlmn and ilmn 352 qijtot=zero 353 do jlmn=1,lmn_size 354 j0lmn=jlmn*(jlmn-1)/2 355 do ilmn=1,jlmn 356 357 klmn=j0lmn+ilmn 358 klm=pawtab(itypat)%indklmn(1,klmn);kln=pawtab(itypat)%indklmn(2,klmn) 359 lmin=pawtab(itypat)%indklmn(3,klmn);lmax=pawtab(itypat)%indklmn(4,klmn) 360 ! --- Sum over angular momenta 361 ! --- compute radial part integration for each angular momentum => intg 362 ! --- (3j) symbols follows the rule: l belongs to abs(li-lj), li+lj. 363 qijb=zero 364 do ll=lmin,lmax,2 365 lm0=ll*ll+ll+1 366 ABI_ALLOCATE(ff,(mesh_size)) 367 ff(1:mesh_size)=(pawtab(itypat)%phiphj(1:mesh_size,kln)& 368 & -pawtab(itypat)%tphitphj(1:mesh_size,kln))& 369 & *j_bessel(1:mesh_size,ll+1) 370 call simp_gen(intg,ff,pawrad(itypat)) 371 ABI_DEALLOCATE(ff) 372 qijbtemp=zero 373 do mm=-ll,ll 374 isel=pawang%gntselect(lm0+mm,klm) 375 if (isel>0) qijbtemp=qijbtemp& 376 & +pawang%realgnt(isel)*ylmb(lm0+mm) 377 end do ! mm 378 ! --- compute angular part with a summation 379 ! --- qijb =\sum_{lm} intg(lm)*qijbtemp 380 qijb(1)=qijb(1) +intg*qijbtemp*ilr(ll+1) 381 qijb(2)=qijb(2) +intg*qijbtemp*ili(ll+1) 382 ! if(ilmn==jlmn) write(std_out,*) "intg, qij",intg,qijbtemp 383 end do ! ll 384 385 ! --- Add exp(-i.b*R) for each atom. 386 if(ilmn==jlmn) qijtot=qijtot+qijb(1) 387 ! if(ilmn==jlmn) write(std_out,*) "qijtot",qijtot 388 x1=qijb(1)*dcos(-xsum)-qijb(2)*dsin(-xsum) 389 x2=qijb(1)*dsin(-xsum)+qijb(2)*dcos(-xsum) 390 ! x1 x2 necessary to avoid changing qijb(1) before 391 ! computing qijb(2) 392 qijb(1)=x1 393 qijb(2)=x2 ! 394 ! if(ilmn==jlmn) write(std_out,*) "qij",jlmn,ilmn,qijb(1),qijb(2) 395 396 do iband1=1,mbandw ! limite inferieure a preciser 397 do iband2=1,mbandw 398 ppr=0.d0 399 ppi=0.d0 400 do ispinor=1,nspinor 401 idx1=iband1*nspinor-(nspinor-ispinor) 402 idx2=iband2*nspinor-(nspinor-ispinor) !to take into account spinors 403 ! write(std_out,*) "iband2",iband2 404 ! product of (a1+ia2)*(b1-ib2) (minus sign because conjugated) 405 ppr=ppr+& 406 ! real part a_1*b_1+a_2*b_2 407 & cprj_k1(iatom,idx1)%cp(1,ilmn)*cprj_k2(iatom,idx2)%cp(1,jlmn)+& 408 & cprj_k1(iatom,idx1)%cp(2,ilmn)*cprj_k2(iatom,idx2)%cp(2,jlmn)+& 409 ! & cprj(iatom,idx1+icg1)%cp(1,ilmn)*cprj(iatom,idx2+icg2)%cp(1,jlmn)+& 410 ! & cprj(iatom,idx1+icg1)%cp(2,ilmn)*cprj(iatom,idx2+icg2)%cp(2,jlmn)+& 411 ! add term on the other triangle of the matrix 412 ! qij is the same for this part because phi are real. 413 & cprj_k1(iatom,idx1)%cp(1,jlmn)*cprj_k2(iatom,idx2)%cp(1,ilmn)+& 414 & cprj_k1(iatom,idx1)%cp(2,jlmn)*cprj_k2(iatom,idx2)%cp(2,ilmn) 415 ! & cprj(iatom,idx1+icg1)%cp(1,jlmn)*cprj(iatom,idx2+icg2)%cp(1,ilmn)+& 416 ! & cprj(iatom,idx1+icg1)%cp(2,jlmn)*cprj(iatom,idx2+icg2)%cp(2,ilmn) 417 ppi=ppi+& 418 ! imaginary part a_1*b_2-a_2*b_1 419 & cprj_k1(iatom,idx1)%cp(1,ilmn)*cprj_k2(iatom,idx2)%cp(2,jlmn)-& 420 & cprj_k1(iatom,idx1)%cp(2,ilmn)*cprj_k2(iatom,idx2)%cp(1,jlmn)+& 421 ! & cprj(iatom,idx1+icg1)%cp(1,ilmn)*cprj(iatom,idx2+icg2)%cp(2,jlmn)-& 422 ! & cprj(iatom,idx1+icg1)%cp(2,ilmn)*cprj(iatom,idx2+icg2)%cp(1,jlmn)+& 423 ! add term on the other triangle of the matrix 424 & cprj_k1(iatom,idx1)%cp(1,jlmn)*cprj_k2(iatom,idx2)%cp(2,ilmn)-& 425 & cprj_k1(iatom,idx1)%cp(2,jlmn)*cprj_k2(iatom,idx2)%cp(1,ilmn) 426 ! & cprj(iatom,idx1+icg1)%cp(1,jlmn)*cprj(iatom,idx2+icg2)%cp(2,ilmn)-& 427 ! & cprj(iatom,idx1+icg1)%cp(2,jlmn)*cprj(iatom,idx2+icg2)%cp(1,ilmn) 428 end do !ispinor 429 ! 430 ! delta: diagonal terms are counted twice ! so 431 ! we need a 0.5 factor for diagonal elements. 432 delta=one 433 ! write(std_out,*) "ppr and ppi computed",ikpt1,ikpt2,iband1,iband2 434 if(ilmn==jlmn) delta=half 435 cm2(1,iband1,iband2)= cm2(1,iband1,iband2)+ & 436 & (qijb(1)*ppr-qijb(2)*ppi)*delta 437 cm2(2,iband1,iband2)= cm2(2,iband1,iband2)+ & 438 & (qijb(2)*ppr+qijb(1)*ppi)*delta 439 end do ! iband2 440 end do ! iband1 441 442 end do ! ilmn 443 end do ! jlmn 444 ! write(std_out,*) "final qijtot",qijtot 445 ABI_DEALLOCATE(j_bessel) 446 end do ! iatom 447 448 ABI_DEALLOCATE(sb_out) 449 ABI_DEALLOCATE(ylmb) 450 call pawcprj_free(cprj_k1) 451 call pawcprj_free(cprj_k2) 452 ABI_DATATYPE_DEALLOCATE(cprj_k1) 453 ABI_DATATYPE_DEALLOCATE(cprj_k2) 454 455 DBG_EXIT("COLL") 456 457 end subroutine smatrix_pawinit