TABLE OF CONTENTS
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...
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 cg_set(2,mcg,mset)= plane wave wavefunction coefficients for several sets of wavefunctions (all k points and spinpol) 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
PARENTS
wf_mixing
CHILDREN
dotprod_g,pawcprj_alloc,pawcprj_free,pawcprj_get
SOURCE
69 #if defined HAVE_CONFIG_H 70 #include "config.h" 71 #endif 72 73 #include "abi_common.h" 74 75 subroutine dotprodm_sumdiag_cgcprj(atindx1,cg_set,cprj_set,dimcprj,& 76 & ibg,icg,ikpt,isppol,istwf,mband,mcg,mcprj,mkmem,& 77 & mpi_enreg,mset,natom,nattyp,nbd,npw,nset1,nset2,nspinor,nsppol,ntypat,& 78 & shift_set1,shift_set2,pawtab,smn,usepaw) 79 80 use defs_basis 81 use defs_abitypes 82 use m_cgtools 83 use m_errors 84 use m_xmpi 85 use m_pawtab, only : pawtab_type 86 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_free 87 88 !This section has been created automatically by the script Abilint (TD). 89 !Do not modify the following lines by hand. 90 #undef ABI_FUNC 91 #define ABI_FUNC 'dotprodm_sumdiag_cgcprj' 92 !End of the abilint section 93 94 implicit none 95 96 !Arguments ------------------------------------ 97 !scalars 98 integer, intent(in) :: ibg,icg,ikpt,isppol,istwf 99 integer, intent(in) :: mband,mcg,mcprj,mkmem,mset 100 integer, intent(in) :: natom,nbd,npw,nset1,nset2,nspinor,nsppol,ntypat 101 integer, intent(in) :: shift_set1,shift_set2,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) :: cg_set(2,mcg,mset) 106 real(dp), intent(out) :: smn(2,nset1,nset2) 107 type(pawcprj_type),intent(in) :: cprj_set(natom,mcprj,mset) 108 type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw) 109 110 !Local variables------------------------------- 111 !scalars 112 integer :: ia,iat,itypat,ibd,icgb,ig 113 integer :: ilmn1,ilmn2,ind_set1,ind_set2,iset1,iset2,klmn 114 real(dp) :: dotr,doti 115 !arrays 116 real(dp),allocatable :: cwavef1(:,:),cwavef2(:,:) 117 type(pawcprj_type),allocatable :: cprj1_k(:,:),cprj2_k(:,:) 118 119 ! ************************************************************************* 120 121 !DEBUG 122 !write(std_out,*)' dotprodm_sumdiag_cgcprj : enter ' 123 !call flush(std_out) 124 !ENDDEBUG 125 126 ABI_ALLOCATE(cwavef1,(2,npw*nspinor)) 127 ABI_ALLOCATE(cwavef2,(2,npw*nspinor)) 128 if(usepaw==1) then 129 ABI_DATATYPE_ALLOCATE(cprj1_k,(natom,nspinor*nbd)) 130 ABI_DATATYPE_ALLOCATE(cprj2_k,(natom,nspinor*nbd)) 131 end if 132 if(usepaw==1) then 133 call pawcprj_alloc(cprj1_k,cprj_set(1,1,1)%ncpgr,dimcprj) 134 end if 135 136 smn(:,:,:)=zero 137 138 icgb=icg 139 do ibd=1,nbd 140 141 do iset1=1,nset1 142 143 ind_set1=iset1+shift_set1 144 145 ! Extract wavefunction information 146 do ig=1,npw*nspinor 147 cwavef1(1,ig)=cg_set(1,ig+icgb,ind_set1) 148 cwavef1(2,ig)=cg_set(2,ig+icgb,ind_set1) 149 end do 150 if(usepaw==1) then 151 call pawcprj_get(atindx1,cprj1_k,cprj_set(:,:,ind_set1),natom,1,ibg,ikpt,1,isppol,mband,& 152 & mkmem,natom,nbd,nbd,nspinor,nsppol,0,& 153 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 154 end if 155 156 do iset2=1,nset2 157 158 ind_set2=iset2+shift_set2 159 if(ind_set2<ind_set1 .and. ind_set2>shift_set1)then 160 continue ! These matrix elements have already been computed, the smn matrix will be completed later. 161 162 else if(ind_set1/=ind_set2)then 163 164 ! Extract wavefunction information 165 do ig=1,npw*nspinor 166 cwavef2(1,ig)=cg_set(1,ig+icgb,ind_set2) 167 cwavef2(2,ig)=cg_set(2,ig+icgb,ind_set2) 168 end do 169 170 if(usepaw==1) then 171 call pawcprj_get(atindx1,cprj2_k,cprj_set(:,:,ind_set2),natom,1,ibg,ikpt,1,isppol,mband,& 172 & mkmem,natom,nbd,nbd,nspinor,nsppol,0,& 173 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 174 end if 175 176 ! Calculate Smn=<cg1|cg2> 177 call dotprod_g(dotr,doti,istwf,npw*nspinor,2,cwavef1,cwavef2,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 178 179 if(usepaw==1) then 180 ia =0 181 do itypat=1,ntypat 182 do iat=1+ia,nattyp(itypat)+ia 183 do ilmn1=1,pawtab(itypat)%lmn_size 184 do ilmn2=1,ilmn1 185 klmn=((ilmn1-1)*ilmn1)/2+ilmn2 186 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd)%cp(1,ilmn1)*cprj2_k(iat,ibd)%cp(1,ilmn2)+& 187 & cprj1_k(iat,ibd)%cp(2,ilmn1)*cprj2_k(iat,ibd)%cp(2,ilmn2)) 188 doti=doti+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd)%cp(1,ilmn1)*cprj2_k(iat,ibd)%cp(2,ilmn2)-& 189 & cprj1_k(iat,ibd)%cp(2,ilmn1)*cprj2_k(iat,ibd)%cp(1,ilmn2)) 190 end do 191 do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size 192 klmn=((ilmn2-1)*ilmn2)/2+ilmn1 193 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd)%cp(1,ilmn1)*cprj2_k(iat,ibd)%cp(1,ilmn2)+& 194 & cprj1_k(iat,ibd)%cp(2,ilmn1)*cprj2_k(iat,ibd)%cp(2,ilmn2)) 195 doti=doti+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd)%cp(1,ilmn1)*cprj2_k(iat,ibd)%cp(2,ilmn2)-& 196 & cprj1_k(iat,ibd)%cp(2,ilmn1)*cprj2_k(iat,ibd)%cp(1,ilmn2)) 197 end do 198 end do 199 end do 200 ia=ia+nattyp(itypat) 201 end do 202 end if ! usepaw 203 204 ! if(.false.)then 205 else 206 ! Diagonal part : no need to extract another wavefunction, and the scalar product must be real 207 208 ! Calculate Smn=<cg1|cg1> 209 call dotprod_g(dotr,doti,istwf,npw*nspinor,1,cwavef1,cwavef1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 210 211 if(usepaw==1) then 212 ia =0 213 do itypat=1,ntypat 214 do iat=1+ia,nattyp(itypat)+ia 215 do ilmn1=1,pawtab(itypat)%lmn_size 216 do ilmn2=1,ilmn1 217 klmn=((ilmn1-1)*ilmn1)/2+ilmn2 218 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd)%cp(1,ilmn1)*cprj1_k(iat,ibd)%cp(1,ilmn2)+& 219 & cprj1_k(iat,ibd)%cp(2,ilmn1)*cprj1_k(iat,ibd)%cp(2,ilmn2)) 220 end do 221 do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size 222 klmn=((ilmn2-1)*ilmn2)/2+ilmn1 223 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd)%cp(1,ilmn1)*cprj1_k(iat,ibd)%cp(1,ilmn2)+& 224 & cprj1_k(iat,ibd)%cp(2,ilmn1)*cprj1_k(iat,ibd)%cp(2,ilmn2)) 225 end do 226 end do 227 end do 228 ia=ia+nattyp(itypat) 229 end do 230 end if ! usepaw 231 doti=zero 232 233 end if ! Compare ind_set1 and ind_set2 234 235 smn(1,iset1,iset2)=smn(1,iset1,iset2)+dotr 236 smn(2,iset1,iset2)=smn(2,iset1,iset2)+doti 237 238 end do ! iset2 239 end do ! iset1 240 241 ! End loop over bands ibd 242 icgb=icgb+npw*nspinor 243 end do 244 245 !Complete the matrix, using its hermitian property. 246 do iset1=1,nset1 247 ind_set1=iset1+shift_set1 248 do iset2=1,nset2 249 ind_set2=iset2+shift_set2 250 if(ind_set2<ind_set1 .and. ind_set2>shift_set1)then 251 smn(1,iset1,iset2)= smn(1,iset2,iset1) 252 smn(2,iset1,iset2)=-smn(2,iset2,iset1) 253 end if 254 end do 255 end do 256 257 ABI_DEALLOCATE(cwavef1) 258 ABI_DEALLOCATE(cwavef2) 259 if(usepaw==1) then 260 call pawcprj_free(cprj1_k) 261 call pawcprj_free(cprj2_k) 262 ABI_DATATYPE_DEALLOCATE(cprj1_k) 263 ABI_DATATYPE_DEALLOCATE(cprj2_k) 264 end if 265 266 end subroutine dotprodm_sumdiag_cgcprj