TABLE OF CONTENTS
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...
COPYRIGHT
Copyright (C) 2017-2018 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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
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
PARENTS
cgcprj_cholesky,wf_mixing
CHILDREN
dotprod_g,pawcprj_alloc,pawcprj_free,pawcprj_get,zhpev
SOURCE
71 #if defined HAVE_CONFIG_H 72 #include "config.h" 73 #endif 74 75 #include "abi_common.h" 76 77 subroutine dotprod_set_cgcprj(atindx1,cg1,cg2,cprj1,cprj2,dimcprj,hermitian,& 78 & ibg1,ibg2,icg1,icg2,ikpt,isppol,istwf,mband,mcg1,mcg2,mcprj1,mcprj2,mkmem,& 79 & mpi_enreg,natom,nattyp,nbd1,nbd2,npw,nspinor,nsppol,ntypat,pawtab,smn,usepaw) 80 81 use defs_basis 82 use defs_abitypes 83 use m_cgtools 84 use m_errors 85 use m_xmpi 86 use m_pawtab, only : pawtab_type 87 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_free 88 89 !This section has been created automatically by the script Abilint (TD). 90 !Do not modify the following lines by hand. 91 #undef ABI_FUNC 92 #define ABI_FUNC 'dotprod_set_cgcprj' 93 !End of the abilint section 94 95 implicit none 96 97 !Arguments ------------------------------------ 98 !scalars 99 integer, intent(in) :: hermitian,ibg1,ibg2,icg1,icg2,ikpt,isppol,istwf 100 integer, intent(in) :: mkmem,mband,mcg1,mcg2,mcprj1,mcprj2 101 integer, intent(in) :: natom,nbd1,nbd2,npw,nspinor,nsppol,ntypat,usepaw 102 type(MPI_type),intent(in) :: mpi_enreg 103 !arrays 104 integer, intent(in) :: atindx1(natom),dimcprj(natom),nattyp(ntypat) 105 real(dp), intent(in) :: cg1(2,mcg1),cg2(2,mcg2) 106 real(dp), intent(out) :: smn(2,nbd1,nbd2) 107 type(pawcprj_type),intent(in) :: cprj1(natom,mcprj1),cprj2(natom,mcprj2) 108 type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw) 109 110 !Local variables------------------------------- 111 !scalars 112 integer :: ia,iat,itypat,ibd1,ibd2,icgb1,icgb2,ier,ig,ii,i1,i2 113 integer :: ilmn1,ilmn2,klmn,max_nbd2,nbd 114 real(dp) :: dotr,doti 115 !arrays 116 real(dp),allocatable :: cwavef1(:,:),cwavef2(:,:),proj(:,:,:) 117 real(dp),allocatable :: eigval(:),eigvec(:,:,:),matrx(:,:),zhpev1(:,:),zhpev2(:) 118 type(pawcprj_type),allocatable :: cprj1_k(:,:),cprj2_k(:,:) 119 120 ! ************************************************************************* 121 122 !DEBUG 123 !write(std_out,*)' dotprod_set_cgcprj : enter ' 124 !write(std_out,*)' dotprod_set_cgcprj : npw, nspinor=',npw,nspinor 125 !write(std_out,*)' dotprod_set_cgcprj : usepaw,nbd1,nbd2=',usepaw,nbd1,nbd2 126 !call flush(std_out) 127 !ENDDEBUG 128 129 if(hermitian==1)then 130 if(nbd1/=nbd2)then 131 MSG_ERROR(' With hermitian==1, nb1 and nb2 must be equal ') 132 end if 133 end if 134 135 ABI_ALLOCATE(cwavef1,(2,npw*nspinor)) 136 ABI_ALLOCATE(cwavef2,(2,npw*nspinor)) 137 if(usepaw==1) then 138 ABI_DATATYPE_ALLOCATE(cprj1_k,(natom,nspinor*nbd1)) 139 ABI_DATATYPE_ALLOCATE(cprj2_k,(natom,nspinor*nbd2)) 140 end if 141 if(usepaw==1 .and. ibg1/=0) then 142 call pawcprj_alloc(cprj1_k,cprj1(1,1)%ncpgr,dimcprj) 143 end if 144 if(usepaw==1 .and. ibg2/=0) then 145 call pawcprj_alloc(cprj2_k,cprj1(1,1)%ncpgr,dimcprj) 146 end if 147 148 icgb1=icg1 149 do ibd1=1,nbd1 150 151 ! Extract wavefunction information 152 do ig=1,npw*nspinor 153 cwavef1(1,ig)=cg1(1,ig+icgb1) 154 cwavef1(2,ig)=cg1(2,ig+icgb1) 155 end do 156 if(usepaw==1 .and. ibg1/=0) then 157 call pawcprj_get(atindx1,cprj1_k,cprj1,natom,1,ibg1,ikpt,1,isppol,mband,& 158 & mkmem,natom,nbd1,nbd1,nspinor,nsppol,0,& 159 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 160 end if 161 162 icgb2=icg2 163 max_nbd2=nbd2 164 if(hermitian==1)max_nbd2=ibd1 165 do ibd2=1,max_nbd2 166 167 ! XG171222 Note that this copy step, being inside the ibd1 loop, is quite detrimental. 168 ! It might be reduced by copying several cwavef2, and use a ZGEMM type of approach. 169 170 ! Extract wavefunction information 171 do ig=1,npw*nspinor 172 cwavef2(1,ig)=cg2(1,ig+icgb2) 173 cwavef2(2,ig)=cg2(2,ig+icgb2) 174 end do 175 176 if(usepaw==1 .and. ibg2/=0) then 177 call pawcprj_get(atindx1,cprj2_k,cprj2,natom,1,ibg2,ikpt,1,isppol,mband,& 178 & mkmem,natom,nbd2,nbd2,nspinor,nsppol,0,& 179 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 180 end if 181 182 ! Calculate Smn=<cg1|cg2> 183 call dotprod_g(dotr,doti,istwf,npw*nspinor,2,cwavef1,cwavef2,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 184 185 if(usepaw==1) then 186 ia =0 187 do itypat=1,ntypat 188 if(ibg1/=0 .and. ibg2/=0)then 189 do iat=1+ia,nattyp(itypat)+ia 190 do ilmn1=1,pawtab(itypat)%lmn_size 191 do ilmn2=1,ilmn1 192 klmn=((ilmn1-1)*ilmn1)/2+ilmn2 193 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2)+& 194 & cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2)) 195 doti=doti+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2)-& 196 & cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2)) 197 end do 198 do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size 199 klmn=((ilmn2-1)*ilmn2)/2+ilmn1 200 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2)+& 201 & cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2)) 202 doti=doti+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2)-& 203 & cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2)) 204 end do 205 end do 206 end do 207 else if(ibg1/=0 .and. ibg2==0)then 208 do iat=1+ia,nattyp(itypat)+ia 209 do ilmn1=1,pawtab(itypat)%lmn_size 210 do ilmn2=1,ilmn1 211 klmn=((ilmn1-1)*ilmn1)/2+ilmn2 212 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2)+& 213 & cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2)) 214 doti=doti+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2)-& 215 & cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2)) 216 end do 217 do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size 218 klmn=((ilmn2-1)*ilmn2)/2+ilmn1 219 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2)+& 220 & cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2)) 221 doti=doti+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2)-& 222 & cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2)) 223 end do 224 end do 225 end do 226 else if(ibg1==0 .and. ibg2/=0)then 227 do iat=1+ia,nattyp(itypat)+ia 228 do ilmn1=1,pawtab(itypat)%lmn_size 229 do ilmn2=1,ilmn1 230 klmn=((ilmn1-1)*ilmn1)/2+ilmn2 231 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2)+& 232 & cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2)) 233 doti=doti+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2)-& 234 & cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2)) 235 end do 236 do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size 237 klmn=((ilmn2-1)*ilmn2)/2+ilmn1 238 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2)+& 239 & cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2)) 240 doti=doti+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2)-& 241 & cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2)) 242 end do 243 end do 244 end do 245 else if(ibg1==0 .and. ibg2==0)then 246 do iat=1+ia,nattyp(itypat)+ia 247 do ilmn1=1,pawtab(itypat)%lmn_size 248 do ilmn2=1,ilmn1 249 klmn=((ilmn1-1)*ilmn1)/2+ilmn2 250 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2)+& 251 & cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2)) 252 doti=doti+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2)-& 253 & cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2)) 254 end do 255 do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size 256 klmn=((ilmn2-1)*ilmn2)/2+ilmn1 257 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2)+& 258 & cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2)) 259 doti=doti+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2)-& 260 & cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2)) 261 end do 262 end do 263 end do 264 end if 265 ia=ia+nattyp(itypat) 266 end do 267 end if 268 smn(1,ibd1,ibd2)=dotr 269 smn(2,ibd1,ibd2)=doti 270 ! End loop over bands ibd2 271 icgb2=icgb2+npw*nspinor 272 273 end do 274 275 ! End loop over bands ibd1 276 icgb1=icgb1+npw*nspinor 277 end do 278 279 !Complete the matrix if hermitian 280 if(hermitian==1)then 281 do ibd1=1,nbd1-1 282 do ibd2=ibd1+1,nbd2 283 smn(1,ibd1,ibd2)= smn(1,ibd2,ibd1) 284 smn(2,ibd1,ibd2)=-smn(2,ibd2,ibd1) 285 end do 286 end do 287 end if 288 289 !DEBUG 290 !write(std_out,*)' smn=',smn 291 !ENDDEBUG 292 293 !====== Debugging section ========== 294 if(.false.)then 295 !DEBUG 296 !Compute the eigenvalues of the projector S herm(S) or herm(S) S, depending on which has lowest dimension. 297 !write(std_out,*)' dotprod_set_cgcprj : compute the projector matrix ' 298 nbd=min(nbd1,nbd2) 299 ABI_ALLOCATE(proj,(2,nbd,nbd)) 300 proj(:,:,:)=zero 301 if(nbd1<=nbd2)then 302 do ibd1=1,nbd1 303 do ibd2=1,nbd2 304 proj(1,:,ibd1)=proj(1,:,ibd1)+smn(1,:,ibd2)*smn(1,ibd1,ibd2)+smn(2,:,ibd2)*smn(2,ibd1,ibd2) 305 proj(2,:,ibd1)=proj(2,:,ibd1)-smn(1,:,ibd2)*smn(2,ibd1,ibd2)+smn(2,:,ibd2)*smn(1,ibd1,ibd2) 306 end do 307 end do 308 else 309 do ibd2=1,nbd2 310 do ibd1=1,nbd1 311 proj(1,:,ibd2)=proj(1,:,ibd2)+smn(1,ibd1,:)*smn(1,ibd1,ibd2)+smn(2,ibd1,:)*smn(2,ibd1,ibd2) 312 proj(2,:,ibd2)=proj(2,:,ibd2)+smn(1,ibd1,:)*smn(2,ibd1,ibd2)-smn(2,ibd1,:)*smn(1,ibd1,ibd2) 313 end do 314 end do 315 end if 316 317 !write(std_out,*)' proj=',proj 318 319 !write(std_out,*)' dotprod_set_cgcprj : compute the eigenvalues of the projector ' 320 ABI_ALLOCATE(matrx,(2,(nbd*(nbd+1))/2)) 321 ii=1 322 do i2=1,nbd 323 do i1=1,i2 324 matrx(1,ii)=proj(1,i1,i2) 325 matrx(2,ii)=proj(2,i1,i2) 326 ii=ii+1 327 end do 328 end do 329 330 ABI_ALLOCATE(zhpev1,(2,2*nbd-1)) 331 ABI_ALLOCATE(zhpev2,(3*nbd-2)) 332 ABI_ALLOCATE(eigval,(nbd)) 333 ABI_ALLOCATE(eigvec,(2,nbd,nbd)) 334 335 call ZHPEV ('V','U',nbd,matrx,eigval,eigvec,nbd,zhpev1,zhpev2,ier) 336 337 !write(std_out,*)' eigval=',eigval 338 339 ABI_DEALLOCATE(matrx) 340 ABI_DEALLOCATE(zhpev1) 341 ABI_DEALLOCATE(zhpev2) 342 ABI_DEALLOCATE(eigval) 343 ABI_DEALLOCATE(eigvec) 344 345 ABI_DEALLOCATE(proj) 346 !stop 347 !ENDDEBUG 348 end if 349 !====== End of debugging section ========== 350 351 ABI_DEALLOCATE(cwavef1) 352 ABI_DEALLOCATE(cwavef2) 353 if(usepaw==1) then 354 if(ibg1/=0)then 355 call pawcprj_free(cprj1_k) 356 end if 357 if(ibg2/=0)then 358 call pawcprj_free(cprj2_k) 359 end if 360 ABI_DATATYPE_DEALLOCATE(cprj1_k) 361 ABI_DATATYPE_DEALLOCATE(cprj2_k) 362 end if 363 364 end subroutine dotprod_set_cgcprj