TABLE OF CONTENTS
- ABINIT/cgcprj_cholesky
- ABINIT/dotprod_set_cgcprj
- ABINIT/dotprodm_sumdiag_cgcprj
- ABINIT/lincom_cgcprj
- ABINIT/m_cgcprj
ABINIT/cgcprj_cholesky [ Modules ]
NAME
cgcprj_cholesky
FUNCTION
Cholesky orthonormalization of the vectors stored in cg+cprj mode. This implementation is NOT band-parallelized Also, it is far of being optimal at the level of linear algebra
INPUTS
atindx1(natom)=index table for atoms, inverse of atindx dimcprj(natom)=number of lmn components in the <p_{lmn}^i|\psi> for the i-th atom icg=shift in cg array to locate current k-point and spinpol ikpt=current k point index isppol=current spin polarization index istwf=input option parameter that describes the storage of wfs mcg=second dimension of cg array (mpw*nspinor*mband*mkmem*nsppol) mcprj=second dimension of cprj_k array mkmem=number of k points which can fit in memory mpi_enreg=information about MPI parallelization natom=number of atoms nattyp(ntypat)=number of atoms of each type in cell. nband=number of bands npw=number of planewaves in basis at this k point nspinor=number of spinor components nsppol=number of spin polarizations ntypat=number of types of atoms pawtab(dtset%ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data usepaw=1 if PAW is activated
SIDE EFFECTS
cg(2,mcg)= plane wave wavefunction coefficients for the set of input wavefunctions (all k points and spinpol) cprj_k(natom,mcprj) <type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk> with NL projectors for the specific k point and spinpol
SOURCE
786 subroutine cgcprj_cholesky(atindx1,cg,cprj_k,dimcprj,icg,ikpt,isppol,istwf,mcg,mcprj,mkmem,& 787 & mpi_enreg,natom,nattyp,nband,npw,nspinor,nsppol,ntypat,pawtab,usepaw) 788 789 !Arguments ------------------------------------ 790 !scalars 791 integer,intent(in) :: icg,ikpt,isppol,istwf,mcg,mcprj,mkmem 792 integer,intent(in) :: natom,nband,npw,nspinor,nsppol,ntypat,usepaw 793 !arrays 794 integer, intent(in) :: atindx1(natom),dimcprj(natom),nattyp(ntypat) 795 real(dp), intent(inout) :: cg(2,mcg) 796 type(pawcprj_type),intent(inout) :: cprj_k(natom,mcprj) 797 type(MPI_type),intent(in) :: mpi_enreg 798 type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw) 799 800 !Local variables ------------------------------ 801 !scalars 802 integer :: hermitian,ierr,ii,inplace 803 !arrays 804 real(dp), allocatable :: dmn(:,:,:),smn(:,:,:) 805 806 ! ************************************************************************* 807 808 ABI_MALLOC(smn,(2,nband,nband)) 809 ABI_MALLOC(dmn,(2,nband,nband)) 810 811 hermitian=1 812 call dotprod_set_cgcprj(atindx1,cg,cg,cprj_k,cprj_k,dimcprj,hermitian,& 813 & 0,0,icg,icg,ikpt,isppol,istwf,nband,mcg,mcg,mcprj,mcprj,mkmem,& 814 & mpi_enreg,natom,nattyp,nband,nband,npw,nspinor,nsppol,ntypat,pawtab,smn,usepaw) 815 816 !Cholesky factorization: O = U^H U with U upper triangle matrix. 817 call ZPOTRF('U',nband,smn,nband,ierr) 818 819 !Solve X U = 1. 820 dmn=zero 821 do ii=1,nband 822 dmn(1,ii,ii)=one 823 end do 824 call ZTRSM('Right','Upper','Normal','Normal',nband,nband,cone,smn,nband,dmn,nband) 825 826 inplace=1 827 !This call does not take into account the fact that X=dmn is an upper triangular matrix... 828 !The number of operations might be divided by two. 829 call lincom_cgcprj(dmn,cg,cprj_k,dimcprj,& 830 & icg,inplace,mcg,mcprj,natom,nband,nband,npw,nspinor,usepaw) 831 832 ABI_FREE(smn) 833 ABI_FREE(dmn) 834 835 end subroutine cgcprj_cholesky
ABINIT/dotprod_set_cgcprj [ Functions ]
NAME
dotprod_set_cgcprj
FUNCTION
For one k point and spinpol, compute the matrix of scalar products between two sets of nband wavefunctions that are known in the cg+cprj representation Smn=<wf1m|wf2n> Can also treat the case of the computation of scalar products within one set of nband wavefunctions without recomputing already computed matrix elements (define hermitian=1) This implementation is NOT band-parallelized Also, it is far of being optimal at the level of linear algebra, and involves extra copying that are detrimental for performance...
INPUTS
atindx1(natom)=index table for atoms, inverse of atindx cg1(2,mcg1)= plane wave wavefunction coefficients for the first set of wavefunctions (all k points and spinpol) cg2(2,mcg2)= plane wave wavefunction coefficients for the second set of wavefunctions (all k points and spinpol) cprj1(natom,mcprj1) <type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk> with NL projectors in the first set, cprj2(natom,mcprj2) <type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk> with NL projectors in the second set, dimcprj(natom)=number of lmn components in the <p_{lmn}^i|\psi> for the i-th atom hermitian= if 1, consider that the Smn matrix is hermitian, and do not recompute already computed matrix elements. ibg1=shift in cprj1 array to locate current k-point. Might be 0, in which case cprj1 is not copied internally, which saves some time/space ibg2=shift in cprj2 array to locate current k-point. Might be 0, in which case cprj2 is not copied internally, which saves some time/space icg1=shift in cg1 array to locate current k-point icg2=shift in cg2 array to locate current k-point ikpt=current k point index isppol=current spin polarization index istwf=input option parameter that describes the storage of wfs mcg1=second dimension of cg1 array (mpw*nspinor*mband1*mkmem*nsppol) mcg2=second dimension of cg2 array (mpw*nspinor*mband2*mkmem*nsppol) mcprj1=second dimension of cprj1 array mcprj2=second dimension of cprj2 array mkmem=number of k points which can fit in memory mpi_enreg=information about MPI parallelization natom=number of atoms nattyp(ntypat)=number of atoms of each type in cell. nbd1=number of bands for the first set of wavefunctions nbd2=number of bands for the second set of wavefunctions npw=number of planewaves in basis at this k point nspinor=number of spinor components nsppol=number of spin polarizations ntypat=number of types of atoms pawtab(dtset%ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data usepaw=1 if PAW is activated
OUTPUT
smn(2,nbd1,nbd2)=matrix of scalar products between the first set of wavefunctions and the second set of wavefunctions
SIDE EFFECTS
SOURCE
104 subroutine dotprod_set_cgcprj(atindx1,cg1,cg2,cprj1,cprj2,dimcprj,hermitian,& 105 & ibg1,ibg2,icg1,icg2,ikpt,isppol,istwf,mband,mcg1,mcg2,mcprj1,mcprj2,mkmem,& 106 & mpi_enreg,natom,nattyp,nbd1,nbd2,npw,nspinor,nsppol,ntypat,pawtab,smn,usepaw) 107 108 !Arguments ------------------------------------ 109 !scalars 110 integer, intent(in) :: hermitian,ibg1,ibg2,icg1,icg2,ikpt,isppol,istwf 111 integer, intent(in) :: mkmem,mband,mcg1,mcg2,mcprj1,mcprj2 112 integer, intent(in) :: natom,nbd1,nbd2,npw,nspinor,nsppol,ntypat,usepaw 113 type(MPI_type),intent(in) :: mpi_enreg 114 !arrays 115 integer, intent(in) :: atindx1(natom),dimcprj(natom),nattyp(ntypat) 116 real(dp), intent(in) :: cg1(2,mcg1),cg2(2,mcg2) 117 real(dp), intent(out) :: smn(2,nbd1,nbd2) 118 type(pawcprj_type),intent(in) :: cprj1(natom,mcprj1),cprj2(natom,mcprj2) 119 type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw) 120 121 !Local variables------------------------------- 122 !scalars 123 integer :: ia,iat,itypat,ibd1,ibd2,icgb1,icgb2,ier,ig,ii,i1,i2,iorder 124 integer :: ilmn1,ilmn2,klmn,max_nbd2,nbd 125 real(dp) :: dotr,doti 126 !arrays 127 real(dp),allocatable :: cwavef1(:,:),cwavef2(:,:),proj(:,:,:) 128 real(dp),allocatable :: eigval(:),eigvec(:,:,:),matrx(:,:),zhpev1(:,:),zhpev2(:) 129 type(pawcprj_type),allocatable :: cprj1_k(:,:),cprj2_k(:,:) 130 131 ! ************************************************************************* 132 133 !DEBUG 134 !write(std_out,*)' dotprod_set_cgcprj : enter ' 135 !write(std_out,*)' dotprod_set_cgcprj : npw, nspinor=',npw,nspinor 136 !write(std_out,*)' dotprod_set_cgcprj : usepaw,nbd1,nbd2=',usepaw,nbd1,nbd2 137 !call flush(std_out) 138 !ENDDEBUG 139 140 if(hermitian==1)then 141 if(nbd1/=nbd2)then 142 ABI_ERROR(' With hermitian==1, nb1 and nb2 must be equal ') 143 end if 144 end if 145 146 ABI_MALLOC(cwavef1,(2,npw*nspinor)) 147 ABI_MALLOC(cwavef2,(2,npw*nspinor)) 148 if(usepaw==1) then 149 ABI_MALLOC(cprj1_k,(natom,nspinor*nbd1)) 150 ABI_MALLOC(cprj2_k,(natom,nspinor*nbd2)) 151 iorder=0 ! There is no change of ordering of cprj when copying wavefunctions 152 end if 153 if(usepaw==1 .and. ibg1/=0) then 154 call pawcprj_alloc(cprj1_k,cprj1(1,1)%ncpgr,dimcprj) 155 end if 156 if(usepaw==1 .and. ibg2/=0) then 157 call pawcprj_alloc(cprj2_k,cprj1(1,1)%ncpgr,dimcprj) 158 end if 159 160 icgb1=icg1 161 do ibd1=1,nbd1 162 163 ! Extract wavefunction information 164 do ig=1,npw*nspinor 165 cwavef1(1,ig)=cg1(1,ig+icgb1) 166 cwavef1(2,ig)=cg1(2,ig+icgb1) 167 end do 168 if(usepaw==1 .and. ibg1/=0) then 169 call pawcprj_get(atindx1,cprj1_k,cprj1,natom,1,ibg1,ikpt,iorder,isppol,mband,& 170 & mkmem,natom,nbd1,nbd1,nspinor,nsppol,0,& 171 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 172 end if 173 174 icgb2=icg2 175 max_nbd2=nbd2 176 if(hermitian==1)max_nbd2=ibd1 177 do ibd2=1,max_nbd2 178 179 ! XG171222 Note that this copy step, being inside the ibd1 loop, is quite detrimental. 180 ! It might be reduced by copying several cwavef2, and use a ZGEMM type of approach. 181 182 ! Extract wavefunction information 183 do ig=1,npw*nspinor 184 cwavef2(1,ig)=cg2(1,ig+icgb2) 185 cwavef2(2,ig)=cg2(2,ig+icgb2) 186 end do 187 188 if(usepaw==1 .and. ibg2/=0) then 189 call pawcprj_get(atindx1,cprj2_k,cprj2,natom,1,ibg2,ikpt,iorder,isppol,mband,& 190 & mkmem,natom,nbd2,nbd2,nspinor,nsppol,0,& 191 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 192 end if 193 194 ! Calculate Smn=<cg1|cg2> 195 call dotprod_g(dotr,doti,istwf,npw*nspinor,2,cwavef1,cwavef2,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 196 197 if(usepaw==1) then 198 ia =0 199 do itypat=1,ntypat 200 if(ibg1/=0 .and. ibg2/=0)then 201 do iat=1+ia,nattyp(itypat)+ia 202 do ilmn1=1,pawtab(itypat)%lmn_size 203 do ilmn2=1,ilmn1 204 klmn=((ilmn1-1)*ilmn1)/2+ilmn2 205 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2)+& 206 & cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2)) 207 doti=doti+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2)-& 208 & cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2)) 209 end do 210 do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size 211 klmn=((ilmn2-1)*ilmn2)/2+ilmn1 212 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2)+& 213 & cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2)) 214 doti=doti+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2)-& 215 & cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2)) 216 end do 217 end do 218 end do 219 else if(ibg1/=0 .and. ibg2==0)then 220 do iat=1+ia,nattyp(itypat)+ia 221 do ilmn1=1,pawtab(itypat)%lmn_size 222 do ilmn2=1,ilmn1 223 klmn=((ilmn1-1)*ilmn1)/2+ilmn2 224 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2)+& 225 & cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2)) 226 doti=doti+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2)-& 227 & cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2)) 228 end do 229 do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size 230 klmn=((ilmn2-1)*ilmn2)/2+ilmn1 231 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2)+& 232 & cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2)) 233 doti=doti+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2)-& 234 & cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2)) 235 end do 236 end do 237 end do 238 else if(ibg1==0 .and. ibg2/=0)then 239 do iat=1+ia,nattyp(itypat)+ia 240 do ilmn1=1,pawtab(itypat)%lmn_size 241 do ilmn2=1,ilmn1 242 klmn=((ilmn1-1)*ilmn1)/2+ilmn2 243 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2)+& 244 & cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2)) 245 doti=doti+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2)-& 246 & cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2)) 247 end do 248 do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size 249 klmn=((ilmn2-1)*ilmn2)/2+ilmn1 250 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2)+& 251 & cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2)) 252 doti=doti+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2)-& 253 & cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2)) 254 end do 255 end do 256 end do 257 else if(ibg1==0 .and. ibg2==0)then 258 do iat=1+ia,nattyp(itypat)+ia 259 do ilmn1=1,pawtab(itypat)%lmn_size 260 do ilmn2=1,ilmn1 261 klmn=((ilmn1-1)*ilmn1)/2+ilmn2 262 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2)+& 263 & cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2)) 264 doti=doti+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2)-& 265 & cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2)) 266 end do 267 do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size 268 klmn=((ilmn2-1)*ilmn2)/2+ilmn1 269 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2)+& 270 & cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2)) 271 doti=doti+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2)-& 272 & cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2)) 273 end do 274 end do 275 end do 276 end if 277 ia=ia+nattyp(itypat) 278 end do 279 end if 280 smn(1,ibd1,ibd2)=dotr 281 smn(2,ibd1,ibd2)=doti 282 ! End loop over bands ibd2 283 icgb2=icgb2+npw*nspinor 284 285 end do 286 287 ! End loop over bands ibd1 288 icgb1=icgb1+npw*nspinor 289 end do 290 291 !Complete the matrix if hermitian 292 if(hermitian==1)then 293 do ibd1=1,nbd1-1 294 do ibd2=ibd1+1,nbd2 295 smn(1,ibd1,ibd2)= smn(1,ibd2,ibd1) 296 smn(2,ibd1,ibd2)=-smn(2,ibd2,ibd1) 297 end do 298 end do 299 end if 300 301 !DEBUG 302 !write(std_out,*)' smn=',smn 303 !ENDDEBUG 304 305 !====== Debugging section ========== 306 if(.false.)then 307 !DEBUG 308 !Compute the eigenvalues of the projector S herm(S) or herm(S) S, depending on which has lowest dimension. 309 !write(std_out,*)' dotprod_set_cgcprj : compute the projector matrix ' 310 nbd=min(nbd1,nbd2) 311 ABI_MALLOC(proj,(2,nbd,nbd)) 312 proj(:,:,:)=zero 313 if(nbd1<=nbd2)then 314 do ibd1=1,nbd1 315 do ibd2=1,nbd2 316 proj(1,:,ibd1)=proj(1,:,ibd1)+smn(1,:,ibd2)*smn(1,ibd1,ibd2)+smn(2,:,ibd2)*smn(2,ibd1,ibd2) 317 proj(2,:,ibd1)=proj(2,:,ibd1)-smn(1,:,ibd2)*smn(2,ibd1,ibd2)+smn(2,:,ibd2)*smn(1,ibd1,ibd2) 318 end do 319 end do 320 else 321 do ibd2=1,nbd2 322 do ibd1=1,nbd1 323 proj(1,:,ibd2)=proj(1,:,ibd2)+smn(1,ibd1,:)*smn(1,ibd1,ibd2)+smn(2,ibd1,:)*smn(2,ibd1,ibd2) 324 proj(2,:,ibd2)=proj(2,:,ibd2)+smn(1,ibd1,:)*smn(2,ibd1,ibd2)-smn(2,ibd1,:)*smn(1,ibd1,ibd2) 325 end do 326 end do 327 end if 328 329 !write(std_out,*)' proj=',proj 330 331 !write(std_out,*)' dotprod_set_cgcprj : compute the eigenvalues of the projector ' 332 ABI_MALLOC(matrx,(2,(nbd*(nbd+1))/2)) 333 ii=1 334 do i2=1,nbd 335 do i1=1,i2 336 matrx(1,ii)=proj(1,i1,i2) 337 matrx(2,ii)=proj(2,i1,i2) 338 ii=ii+1 339 end do 340 end do 341 342 ABI_MALLOC(zhpev1,(2,2*nbd-1)) 343 ABI_MALLOC(zhpev2,(3*nbd-2)) 344 ABI_MALLOC(eigval,(nbd)) 345 ABI_MALLOC(eigvec,(2,nbd,nbd)) 346 347 call ZHPEV ('V','U',nbd,matrx,eigval,eigvec,nbd,zhpev1,zhpev2,ier) 348 349 !write(std_out,*)' eigval=',eigval 350 351 ABI_FREE(matrx) 352 ABI_FREE(zhpev1) 353 ABI_FREE(zhpev2) 354 ABI_FREE(eigval) 355 ABI_FREE(eigvec) 356 357 ABI_FREE(proj) 358 !stop 359 !ENDDEBUG 360 end if 361 !====== End of debugging section ========== 362 363 ABI_FREE(cwavef1) 364 ABI_FREE(cwavef2) 365 if(usepaw==1) then 366 if(ibg1/=0)then 367 call pawcprj_free(cprj1_k) 368 end if 369 if(ibg2/=0)then 370 call pawcprj_free(cprj2_k) 371 end if 372 ABI_FREE(cprj1_k) 373 ABI_FREE(cprj2_k) 374 end if 375 376 end subroutine dotprod_set_cgcprj
ABINIT/dotprodm_sumdiag_cgcprj [ Functions ]
NAME
dotprodm_sumdiag_cgcprj
FUNCTION
For one k point and spinpol, compute the matrix of sum of diagonal scalar products between sets of nband wavefunctions that are known in the cg+cprj representation The sets of wavefunctions must be contained in a big array of sets of wavefunctions. Sij=Sum_m <wfm(set i)|wfm(set j)> Can also treat the case of the computation of scalar products within one set of wavefunctions This implementation is NOT band-parallelized Also, it is far of being optimal at the level of linear algebra, and involves extra copying that are detrimental for performance...
INPUTS
atindx1(natom)=index table for atoms, inverse of atindx cg_set(2,mcg,mset)= plane wave wavefunction coefficients for several sets of wavefunctions (all k points and spins) cprj_set(natom,mcprj,mset) <type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk> with NL projectors in the different sets dimcprj(natom)=number of lmn components in the <p_{lmn}^i|\psi> for the i-th atom ibg=shift in cprj_set array to locate current k-point icg=shift in cg_set array to locate current k-point ikpt=current k point index isppol=current spin polarization index istwf=input option parameter that describes the storage of wfs mband=maximum number of bands (used in the dimensioning of cprj_set) mcg=second dimension of cg array (mpw*nspinor*mband*mkmem*nsppol) mcprj=second dimension of cprj array mkmem=number of k points which can fit in memory mpi_enreg=information about MPI parallelization mset=third dimension of cg_set and cprj_set, maximum number of sets natom=number of atoms nattyp(ntypat)=number of atoms of each type in cell. nbd=number of bands for each set of wavefunctions npw=number of planewaves in basis at this k point nset1=number of sets of wavefunctions to be considered in the left side of the scalar products nset2=number of sets of wavefunctions to be considered in the right side of the scalar products nspinor=number of spinor components nsppol=number of spin polarizations ntypat=number of types of atoms pawtab(dtset%ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data shift_set1=shift that defines the first set of wavefunctions to be considered in the left side of the scalar products shift_set2=shift that defines the first set of wavefunctions to be considered in the right side of the scalar products usepaw=1 if PAW is activated
OUTPUT
smn(2,nset1,nset2)=matrix of sum of diagonal scalar products between the first set of wavefunctions and the second set of wavefunctions
SIDE EFFECTS
SOURCE
435 subroutine dotprodm_sumdiag_cgcprj(atindx1,cg_set,cprj_set,dimcprj,& 436 & ibg,icg,ikpt,isppol,istwf,mband,mcg,mcprj,mkmem,& 437 & mpi_enreg,mset,natom,nattyp,nbd,npw,nset1,nset2,nspinor,nsppol,ntypat,& 438 & shift_set1,shift_set2,pawtab,smn,usepaw) 439 440 !Arguments ------------------------------------ 441 !scalars 442 integer, intent(in) :: ibg,icg,ikpt,isppol,istwf 443 integer, intent(in) :: mband,mcg,mcprj,mkmem,mset 444 integer, intent(in) :: natom,nbd,npw,nset1,nset2,nspinor,nsppol,ntypat 445 integer, intent(in) :: shift_set1,shift_set2,usepaw 446 type(MPI_type),intent(in) :: mpi_enreg 447 !arrays 448 integer, intent(in) :: atindx1(natom),dimcprj(natom),nattyp(ntypat) 449 real(dp), intent(in) :: cg_set(2,mcg,mset) 450 real(dp), intent(out) :: smn(2,nset1,nset2) 451 type(pawcprj_type),intent(in) :: cprj_set(natom,mcprj,mset) 452 type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw) 453 454 !Local variables------------------------------- 455 !scalars 456 integer :: ia,iat,itypat,ibd,icgb,ig,iorder 457 integer :: ilmn1,ilmn2,ind_set1,ind_set2,iset1,iset2,klmn 458 real(dp) :: dotr,doti 459 !arrays 460 real(dp),allocatable :: cwavef1(:,:),cwavef2(:,:) 461 type(pawcprj_type),allocatable :: cprj1_k(:,:),cprj2_k(:,:) 462 463 ! ************************************************************************* 464 465 !DEBUG 466 !write(std_out,*)' dotprodm_sumdiag_cgcprj : enter ' 467 !call flush(std_out) 468 !ENDDEBUG 469 470 ABI_MALLOC(cwavef1,(2,npw*nspinor)) 471 ABI_MALLOC(cwavef2,(2,npw*nspinor)) 472 if(usepaw==1) then 473 ABI_MALLOC(cprj1_k,(natom,nspinor*nbd)) 474 ABI_MALLOC(cprj2_k,(natom,nspinor*nbd)) 475 iorder=0 ! There is no change of ordering in the copy of wavefunctions 476 call pawcprj_alloc(cprj1_k,cprj_set(1,1,1)%ncpgr,dimcprj) 477 end if 478 479 smn(:,:,:)=zero 480 481 icgb=icg 482 do ibd=1,nbd 483 484 do iset1=1,nset1 485 486 ind_set1=iset1+shift_set1 487 488 ! Extract wavefunction information 489 do ig=1,npw*nspinor 490 cwavef1(1,ig)=cg_set(1,ig+icgb,ind_set1) 491 cwavef1(2,ig)=cg_set(2,ig+icgb,ind_set1) 492 end do 493 if(usepaw==1) then 494 call pawcprj_get(atindx1,cprj1_k,cprj_set(:,:,ind_set1),natom,1,ibg,ikpt,iorder,isppol,mband,& 495 & mkmem,natom,nbd,nbd,nspinor,nsppol,0,& 496 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 497 end if 498 499 do iset2=1,nset2 500 501 ind_set2=iset2+shift_set2 502 if(ind_set2<ind_set1 .and. ind_set2>shift_set1)then 503 continue ! These matrix elements have already been computed, the smn matrix will be completed later. 504 505 else if(ind_set1/=ind_set2)then 506 507 ! Extract wavefunction information 508 do ig=1,npw*nspinor 509 cwavef2(1,ig)=cg_set(1,ig+icgb,ind_set2) 510 cwavef2(2,ig)=cg_set(2,ig+icgb,ind_set2) 511 end do 512 513 if(usepaw==1) then 514 call pawcprj_get(atindx1,cprj2_k,cprj_set(:,:,ind_set2),natom,1,ibg,ikpt,iorder,isppol,mband,& 515 & mkmem,natom,nbd,nbd,nspinor,nsppol,0,& 516 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 517 end if 518 519 ! Calculate Smn=<cg1|cg2> 520 call dotprod_g(dotr,doti,istwf,npw*nspinor,2,cwavef1,cwavef2,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 521 522 if(usepaw==1) then 523 ia =0 524 do itypat=1,ntypat 525 do iat=1+ia,nattyp(itypat)+ia 526 do ilmn1=1,pawtab(itypat)%lmn_size 527 do ilmn2=1,ilmn1 528 klmn=((ilmn1-1)*ilmn1)/2+ilmn2 529 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd)%cp(1,ilmn1)*cprj2_k(iat,ibd)%cp(1,ilmn2)+& 530 & cprj1_k(iat,ibd)%cp(2,ilmn1)*cprj2_k(iat,ibd)%cp(2,ilmn2)) 531 doti=doti+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd)%cp(1,ilmn1)*cprj2_k(iat,ibd)%cp(2,ilmn2)-& 532 & cprj1_k(iat,ibd)%cp(2,ilmn1)*cprj2_k(iat,ibd)%cp(1,ilmn2)) 533 end do 534 do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size 535 klmn=((ilmn2-1)*ilmn2)/2+ilmn1 536 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd)%cp(1,ilmn1)*cprj2_k(iat,ibd)%cp(1,ilmn2)+& 537 & cprj1_k(iat,ibd)%cp(2,ilmn1)*cprj2_k(iat,ibd)%cp(2,ilmn2)) 538 doti=doti+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd)%cp(1,ilmn1)*cprj2_k(iat,ibd)%cp(2,ilmn2)-& 539 & cprj1_k(iat,ibd)%cp(2,ilmn1)*cprj2_k(iat,ibd)%cp(1,ilmn2)) 540 end do 541 end do 542 end do 543 ia=ia+nattyp(itypat) 544 end do 545 end if ! usepaw 546 547 ! if(.false.)then 548 else 549 ! Diagonal part : no need to extract another wavefunction, and the scalar product must be real 550 551 ! Calculate Smn=<cg1|cg1> 552 call dotprod_g(dotr,doti,istwf,npw*nspinor,1,cwavef1,cwavef1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 553 554 if(usepaw==1) then 555 ia =0 556 do itypat=1,ntypat 557 do iat=1+ia,nattyp(itypat)+ia 558 do ilmn1=1,pawtab(itypat)%lmn_size 559 do ilmn2=1,ilmn1 560 klmn=((ilmn1-1)*ilmn1)/2+ilmn2 561 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd)%cp(1,ilmn1)*cprj1_k(iat,ibd)%cp(1,ilmn2)+& 562 & cprj1_k(iat,ibd)%cp(2,ilmn1)*cprj1_k(iat,ibd)%cp(2,ilmn2)) 563 end do 564 do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size 565 klmn=((ilmn2-1)*ilmn2)/2+ilmn1 566 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd)%cp(1,ilmn1)*cprj1_k(iat,ibd)%cp(1,ilmn2)+& 567 & cprj1_k(iat,ibd)%cp(2,ilmn1)*cprj1_k(iat,ibd)%cp(2,ilmn2)) 568 end do 569 end do 570 end do 571 ia=ia+nattyp(itypat) 572 end do 573 end if ! usepaw 574 doti=zero 575 576 end if ! Compare ind_set1 and ind_set2 577 578 smn(1,iset1,iset2)=smn(1,iset1,iset2)+dotr 579 smn(2,iset1,iset2)=smn(2,iset1,iset2)+doti 580 581 end do ! iset2 582 end do ! iset1 583 584 ! End loop over bands ibd 585 icgb=icgb+npw*nspinor 586 end do 587 588 !Complete the matrix, using its hermitian property. 589 do iset1=1,nset1 590 ind_set1=iset1+shift_set1 591 do iset2=1,nset2 592 ind_set2=iset2+shift_set2 593 if(ind_set2<ind_set1 .and. ind_set2>shift_set1)then 594 smn(1,iset1,iset2)= smn(1,iset2,iset1) 595 smn(2,iset1,iset2)=-smn(2,iset2,iset1) 596 end if 597 end do 598 end do 599 600 ABI_FREE(cwavef1) 601 ABI_FREE(cwavef2) 602 if(usepaw==1) then 603 call pawcprj_free(cprj1_k) 604 call pawcprj_free(cprj2_k) 605 ABI_FREE(cprj1_k) 606 ABI_FREE(cprj2_k) 607 end if 608 609 end subroutine dotprodm_sumdiag_cgcprj
ABINIT/lincom_cgcprj [ Functions ]
NAME
lincom_cgcprj
FUNCTION
For one k point and spin, compute a set (size nband_out) of linear combinations of nband_in wavefunctions, that are known in the cg+cprj representation : cgout_n(:,:) <--- Sum_m [ cg_m(:,:) . alpha_mn ] cprjout_n(:,:) <--- Sum_m [ cprj_m(:,:) . alpha_mn ] If nband_out is smaller or equal to nband_in, the result might be in-place output in cg instead of cgout, and in cprj instead of cprjout). Otherwise, it is contained in the optional cgout+cprjout pair. In the present status, the cg and cgout relates to all the k points and spins, and rely on the icg index, while it is assumed that cprj and cprjout refer to the specific k point and spin. This is not coherent. THIS MIGHT BE CHANGED IN THE FUTURE ! This implementation is NOT band-parallelized Also, it is far of being optimal at the level of linear algebra, and involves extra copying that are detrimental for performance...
INPUTS
alpha_mn(2,nband_in,nband_out)=complex matrix of coefficients of the linear combinations to be computed dimcprj(natom)=number of lmn components in the <p_{lmn}^i|\psi> for the i-th atom icg=shift in cg array to locate current k-point and spinpol (for input, and possibly for in-place output) inplace= if 0, output in cgout and cprjout ; if 1, output in cg and cprj mcg=second dimension of cg array (mpw*nspinor*mband*mkmem*nsppol) mcprj=second dimension of cprj array natom=number of atoms nband_in=number of bands, size of the input set of wavefunctions nband_out=number of bands, size of the output set of wavefunctions (should be equal to nband_in if inplace==1) npw=number of planewaves in basis at this k point nspinor=number of spinor components usepaw=1 if PAW is activated [icgout= shift in cgout array to locate current k-point and spinpol (for output)] [mcgout=second dimension of cgout array (mpw*nspinor*mband*mkmem*nsppol)] [mcprjout=second dimension of cprjout array]
OUTPUT
[cgout(2,mcgout)= plane wave wavefunction coefficients for the set of output wavefunctions] [cprjout(natom,mcprjout) <type(pawcprj_type)>= projected output wave functions <Proj_i|Cnk> with NL projectors]
SIDE EFFECTS
(this quantities are input, and possibly updated output when inplace==1) cg(2,mcg)= plane wave wavefunction coefficients for the set of input wavefunctions (all k points and spinpol) cprj(natom,mcprj) <type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk> with NL projectors
SOURCE
663 subroutine lincom_cgcprj(alpha_mn,cg,cprj,dimcprj,& 664 & icg,inplace,mcg,mcprj,natom,nband_in,nband_out,npw,nspinor,usepaw, & 665 & cgout,cprjout,icgout) ! optional args 666 667 !Arguments ------------------------------------ 668 !scalars 669 integer, intent(in) :: icg,inplace,mcg,mcprj 670 integer, intent(in) :: natom,nband_in,nband_out,npw,nspinor,usepaw 671 integer, intent(in),optional :: icgout 672 !arrays 673 integer, intent(in) :: dimcprj(natom) 674 real(dp), intent(inout) :: cg(2,mcg) 675 real(dp), intent(in) :: alpha_mn(2,nband_in,nband_out) 676 real(dp), intent(out),optional :: cgout(:,:) 677 type(pawcprj_type),intent(inout) :: cprj(natom,mcprj) 678 type(pawcprj_type),intent(inout),optional :: cprjout(:,:) ! ifort and others are buggy for optional intent(out) structured types 679 680 !Local variables------------------------------- 681 !scalars 682 integer :: iband_in,iband_out,ii 683 !arrays 684 real(dp),allocatable :: al(:,:),cgout_(:,:) 685 type(pawcprj_type),allocatable :: cprjout_(:,:) 686 687 ! ************************************************************************* 688 689 !DEBUG 690 !write(std_out,*)' lincom_cgcprj : enter ' 691 !write(std_out,*)' lincom_cgcprj : npw, nspinor=',npw,nspinor 692 !write(std_out,*)' lincom_cgcprj : icgout=',icgout 693 !ENDDEBUG 694 695 if(inplace==0)then 696 if(.not.present(cgout))then 697 ABI_ERROR(' inplace==0 while .not.present(cgout) is not permitted ') 698 end if 699 if(usepaw==1) then 700 if(.not.present(cprjout))then 701 ABI_ERROR(' inplace==0 and usepaw==1 while .not.present(cprjout) is not permitted ') 702 end if 703 end if 704 end if 705 706 !Take care of the plane wave part 707 ABI_MALLOC(cgout_,(2,npw*nspinor*nband_out)) 708 709 call zgemm('N','N',npw*nspinor,nband_out,nband_in,dcmplx(1._dp), & 710 & cg(:,icg+1:icg+npw*nspinor*nband_in),npw*nspinor, & 711 & alpha_mn,nband_in,dcmplx(0._dp),cgout_,npw*nspinor) 712 713 if(inplace==1)then 714 cg(:,icg+1:icg+npw*nspinor*nband_out)=cgout_ 715 else 716 cgout(:,icgout+1:icgout+npw*nspinor*nband_out)=cgout_ 717 end if 718 ABI_FREE(cgout_) 719 720 !Take care of the cprj part 721 if(usepaw==1) then 722 723 ABI_MALLOC(cprjout_,(natom,nspinor*nband_out)) 724 call pawcprj_alloc(cprjout_,cprj(1,1)%ncpgr,dimcprj) 725 ABI_MALLOC(al,(2,nband_in)) 726 do iband_out=1,nband_out 727 ii=(iband_out-1)*nspinor 728 do iband_in=1,nband_in 729 al(1,iband_in)=alpha_mn(1,iband_in,iband_out) 730 al(2,iband_in)=alpha_mn(2,iband_in,iband_out) 731 end do 732 call pawcprj_lincom(al,cprj,cprjout_(:,ii+1:ii+nspinor),nband_in) 733 end do 734 ABI_FREE(al) 735 736 if(inplace==1)then 737 cprj=cprjout_ 738 else 739 cprjout=cprjout_ 740 end if 741 call pawcprj_free(cprjout_) 742 ABI_FREE(cprjout_) 743 744 end if 745 746 end subroutine lincom_cgcprj
ABINIT/m_cgcprj [ Modules ]
NAME
m_cgcprj
FUNCTION
Functions operating on wavefunctions in the cg+cprj representation.
COPYRIGHT
Copyright (C) 2008-2022 ABINIT group (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 .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_cgcprj 23 24 use defs_basis 25 use defs_abitypes 26 use m_cgtools 27 use m_errors 28 use m_xmpi 29 30 use m_pawtab, only : pawtab_type 31 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_free, pawcprj_lincom 32 33 implicit none 34 35 private