TABLE OF CONTENTS
ABINIT/partial_dos_fractions_paw [ Functions ]
NAME
partial_dos_fractions_paw
FUNCTION
Calculate PAW contributions to the partial DOS fractions (tetrahedron method)
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (SM,MT) 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 .
INPUTS
cprj(natom,mcprj)= <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) dtset structured datatype, from which one uses : iatsph(nasph)=number of atoms used to project dos kpt(3,nkpt) =irreducible kpoints mband =max number of bands per k-point mkmem =number of kpoints in memory natom =number of atoms in total natsph =number of atoms ofor which the spherical decomposition must be done nband =number of electronic bands for each kpoint nkpt =number of irreducible kpoints nspinor =number of spinor components nsppol =1 or 2 spin polarization channels fatbands_flag =1 if pawfatbnd=1 or 2 mbesslang=maximum angular momentum for Bessel function expansion mpi_enreg=informations about MPI parallelization prtdosm=option for the m-contributions to the partial DOS ndosfraction=natsph*mbesslang paw_dos_flag=option for the PAW contributions to the partial DOS pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data: pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
OUTPUT
=== If paw_dos_flag==1: dos%fractions_paw1(ikpt,iband,isppol,natom*mbesslang) = contribution to dos fractions from the PAW partial waves (phi) dos%fractions_pawt1(ikpt,iband,isppol,natom*mbesslang) = contribution to dos fractions from the PAW pseudo partial waves (phi_tild)
SIDE EFFECTS
dos%fractions(ikpt,iband,isppol,ndosfraction) = percentage of s, p, d..
character on each atom for the wavefunction # ikpt,iband, isppol
As input: contains only the pseudo contribution
As output: contains pseudo contribution + PAW corrections
== if prtdosm==1
dos%fractions_m(ikpt,iband,isppol,ndosfraction*mbesslang*prtdosm) =
m discretization of partial DOS fractions
PARENTS
outscfcv
CHILDREN
pawcprj_alloc,pawcprj_free,simp_gen,timab,xmpi_sum
SOURCE
63 #if defined HAVE_CONFIG_H 64 #include "config.h" 65 #endif 66 67 #include "abi_common.h" 68 69 subroutine partial_dos_fractions_paw(dos,cprj,dimcprj,dtset,mcprj,mkmem,mpi_enreg,pawrad,pawtab) 70 71 use defs_basis 72 use defs_abitypes 73 use m_xmpi 74 use m_errors 75 use m_profiling_abi 76 77 use m_pawrad, only : pawrad_type, simp_gen 78 use m_pawtab, only : pawtab_type 79 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_free 80 use m_epjdos, only : epjdos_t 81 82 !This section has been created automatically by the script Abilint (TD). 83 !Do not modify the following lines by hand. 84 #undef ABI_FUNC 85 #define ABI_FUNC 'partial_dos_fractions_paw' 86 use interfaces_18_timing 87 use interfaces_32_util 88 !End of the abilint section 89 90 implicit none 91 92 !Arguments ------------------------------------ 93 !scalars 94 integer,intent(in) :: mcprj,mkmem 95 type(MPI_type),intent(in) :: mpi_enreg 96 type(dataset_type),intent(in) :: dtset 97 type(epjdos_t),intent(inout) :: dos 98 !arrays 99 integer,intent(in) :: dimcprj(dtset%natom) 100 type(pawcprj_type),intent(in) :: cprj(dtset%natom,mcprj) 101 type(pawrad_type),intent(in) :: pawrad(dtset%ntypat) 102 type(pawtab_type),target,intent(in) :: pawtab(dtset%ntypat) 103 104 !Local variables------------------------------- 105 !scalars 106 integer :: bandpp,basis_size,comm_kptband,cplex,fatbands_flag,iat,iatom,iband,ibg,ibsp 107 integer :: ierr,ikpt,il,ilang,ilmn,iln,im,iorder_cprj,ispinor,isppol,itypat,j0lmn,j0ln 108 integer :: jl,jlmn,jln,jm,klmn,kln,lmn_size,mbesslang,me_band,me_kpt,my_nspinor 109 integer :: nband_cprj_k,nband_k,ndosfraction,nprocband,nproc_spkptband,paw_dos_flag,prtdosm 110 real(dp) :: cpij,one_over_nproc 111 character(len=500) :: message 112 !arrays 113 integer ,allocatable :: dimcprj_atsph(:) 114 integer,ABI_CONTIGUOUS pointer :: indlmn(:,:) 115 real(dp) :: tsec(2) 116 real(dp),allocatable :: int1(:,:),int2(:,:),int1m2(:,:) 117 type(pawcprj_type),allocatable :: cprj_k(:,:) 118 119 !****************************************************************************************** 120 121 DBG_ENTER("COLL") 122 !return 123 124 ABI_CHECK(mkmem/=0,"mkmem==0 not supported anymore!") 125 126 fatbands_flag = dos%fatbands_flag 127 mbesslang = dos%mbesslang 128 prtdosm = dos%prtdosm 129 ndosfraction = dos%ndosfraction 130 paw_dos_flag = dos%paw_dos_flag 131 132 !m-decomposed DOS not compatible with PAW-decomposed DOS 133 if(prtdosm>=1.and.paw_dos_flag==1) then 134 message = 'm-decomposed DOS not compatible with PAW-decomposed DOS !' 135 MSG_ERROR(message) 136 end if 137 138 !Prepare some useful integrals 139 basis_size=pawtab(1)%basis_size 140 if (dtset%ntypat>1) then 141 do itypat=1,dtset%ntypat 142 basis_size=max(basis_size,pawtab(itypat)%basis_size) 143 end do 144 end if 145 ABI_ALLOCATE(int1 ,(basis_size*(basis_size+1)/2,dtset%natsph)) 146 ABI_ALLOCATE(int2,(basis_size*(basis_size+1)/2,dtset%natsph)) 147 ABI_ALLOCATE(int1m2,(basis_size*(basis_size+1)/2,dtset%natsph)) 148 int1=zero;int2=zero;int1m2=zero 149 do iat=1,dtset%natsph 150 iatom=dtset%iatsph(iat) 151 itypat= dtset%typat(iatom) 152 do jln=1,pawtab(itypat)%basis_size 153 j0ln=jln*(jln-1)/2 154 do iln=1,jln 155 kln=j0ln+iln 156 call simp_gen(int1(kln,iat),pawtab(itypat)%phiphj(:,kln),pawrad(itypat)) 157 if (dtset%pawprtdos<2) then 158 call simp_gen(int2(kln,iat),pawtab(itypat)%tphitphj(:,kln),pawrad(itypat)) 159 int1m2(kln,iat)=int1(kln,iat)-int2(kln,iat) 160 else 161 int2(kln,iat)=zero;int1m2(kln,iat)=int1(kln,iat) 162 end if 163 end do !iln 164 end do !jln 165 end do 166 167 !Antiferro case 168 if (dtset%nspden==2.and.dtset%nsppol==1.and.dtset%nspinor==1) then 169 int1m2(:,:)=half*int1m2(:,:) 170 if (paw_dos_flag==1.or.fatbands_flag==1.or.prtdosm==2) then 171 int1(:,:)=half*int1(:,:);int2(:,:)=half*int2(:,:) 172 end if 173 end if 174 175 !Init parallelism 176 comm_kptband=mpi_enreg%comm_kptband 177 nproc_spkptband=xmpi_comm_size(comm_kptband)*mpi_enreg%nproc_spinor 178 me_kpt=mpi_enreg%me_kpt ; me_band=mpi_enreg%me_band 179 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 180 bandpp=1;if (mpi_enreg%paral_kgb==1) bandpp=mpi_enreg%bandpp 181 !Check if cprj is distributed over bands 182 nprocband=my_nspinor*dtset%mband*dtset%mkmem*dtset%nsppol/mcprj 183 if (nprocband/=mpi_enreg%nproc_band) then 184 message='wrong mcprj/nproc_band!' 185 MSG_BUG(message) 186 end if 187 188 !Quick hack: in case of parallelism, dos_fractions have already 189 ! been reduced over MPI processes; they have to be prepared before 190 ! the next reduction (at the end of the following loop). 191 if (nproc_spkptband>1) then 192 one_over_nproc=one/real(nproc_spkptband,kind=dp) 193 !$OMP PARALLEL DO COLLAPSE(4) & 194 !$OMP& DEFAULT(SHARED) PRIVATE(ilang,isppol,iband,ikpt) 195 do ilang=1,ndosfraction 196 do isppol=1,dtset%nsppol 197 do iband=1,dtset%mband 198 do ikpt=1,dtset%nkpt 199 dos%fractions(ikpt,iband,isppol,ilang)= & 200 & one_over_nproc*dos%fractions(ikpt,iband,isppol,ilang) 201 end do 202 end do 203 end do 204 end do 205 !$OMP END PARALLEL DO 206 if (fatbands_flag==1.or.prtdosm==1.or.prtdosm==2) then 207 !$OMP PARALLEL DO COLLAPSE(4) & 208 !$OMP& DEFAULT(SHARED) PRIVATE(ilang,isppol,iband,ikpt) 209 do ilang=1,ndosfraction*mbesslang 210 do isppol=1,dtset%nsppol 211 do iband=1,dtset%mband 212 do ikpt=1,dtset%nkpt 213 dos%fractions_m(ikpt,iband,isppol,ilang)= & 214 & one_over_nproc*dos%fractions_m(ikpt,iband,isppol,ilang) 215 end do 216 end do 217 end do 218 end do 219 !$OMP END PARALLEL DO 220 end if 221 end if 222 223 iorder_cprj=0 224 225 !LOOPS OVER SPINS,KPTS 226 ibg=0 227 do isppol =1,dtset%nsppol 228 do ikpt=1,dtset%nkpt 229 230 nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 231 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_kpt)) cycle 232 233 cplex=2;if (dtset%istwfk(ikpt)>1) cplex=1 234 nband_cprj_k=nband_k/nprocband 235 ABI_DATATYPE_ALLOCATE(cprj_k,(dtset%natsph,my_nspinor*nband_cprj_k)) 236 ABI_ALLOCATE(dimcprj_atsph,(dtset%natsph)) 237 do iat=1,dtset%natsph 238 dimcprj_atsph(iat)=dimcprj(dtset%iatsph(iat)) 239 end do 240 call pawcprj_alloc(cprj_k,0,dimcprj_atsph) 241 ABI_DEALLOCATE(dimcprj_atsph) 242 243 ! Extract cprj for this k-point. 244 ibsp=0 245 do iband=1,nband_cprj_k 246 do ispinor=1,my_nspinor 247 ibsp=ibsp+1 248 do iat=1,dtset%natsph 249 iatom=dtset%iatsph(iat) 250 cprj_k(iat,ibsp)%cp(:,:)=cprj(iatom,ibsp+ibg)%cp(:,:) 251 end do 252 end do 253 end do 254 255 ! LOOP OVER ATOMS (natsph_extra is not included on purpose) 256 do iat=1,dtset%natsph 257 iatom=dtset%iatsph(iat) 258 itypat= dtset%typat(iatom) 259 lmn_size=pawtab(itypat)%lmn_size 260 indlmn => pawtab(itypat)%indlmn 261 262 ! LOOP OVER BANDS 263 ibsp=0 264 do iband=1,nband_k 265 266 if (mod((iband-1)/bandpp,nprocband)/=me_band) cycle 267 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband,isppol,me_kpt)) then 268 ibsp=ibsp+my_nspinor;cycle 269 end if 270 271 do ispinor=1,my_nspinor 272 ibsp=ibsp+1 273 274 do ilang=1,mbesslang 275 276 do jlmn=1,lmn_size 277 jl=indlmn(1,jlmn) 278 jm=indlmn(2,jlmn) 279 j0lmn=jlmn*(jlmn-1)/2 280 do ilmn=1,jlmn 281 il=indlmn(1,ilmn) 282 im=indlmn(2,ilmn) 283 klmn=j0lmn+ilmn 284 kln=pawtab(itypat)%indklmn(2,klmn) 285 286 if (il==ilang-1.and.jl==ilang-1.and.im==jm) then 287 288 cpij=cprj_k(iat,ibsp)%cp(1,ilmn)*cprj_k(iat,ibsp)%cp(1,jlmn) 289 if (cplex==2) cpij=cpij+cprj_k(iat,ibsp)%cp(2,ilmn)*cprj_k(iat,ibsp)%cp(2,jlmn) 290 cpij=pawtab(itypat)%dltij(klmn)*cpij 291 292 dos%fractions(ikpt,iband,isppol,mbesslang*(iat-1)+ilang)= & 293 & dos%fractions(ikpt,iband,isppol,mbesslang*(iat-1)+ilang) + & 294 & cpij*int1m2(kln,iat) 295 if (prtdosm==1) then 296 dos%fractions_m(ikpt,iband,isppol,mbesslang**2*(iat-1)+il**2+il+1+im)= & 297 & dos%fractions_m(ikpt,iband,isppol,mbesslang**2*(iat-1)+il**2+il+1+im) + & 298 & cpij*int1m2(kln,iat) 299 end if 300 if (fatbands_flag==1.or.prtdosm==2) then 301 dos%fractions_m(ikpt,iband,isppol,mbesslang**2*(iat-1)+il**2+il+1+im)= & 302 & dos%fractions_m(ikpt,iband,isppol,mbesslang**2*(iat-1)+il**2+il+1+im) + & 303 & cpij*int1(kln,iat) 304 end if 305 if (paw_dos_flag==1) then 306 dos%fractions_paw1(ikpt,iband,isppol,mbesslang*(iat-1)+ilang)= & 307 & dos%fractions_paw1(ikpt,iband,isppol,mbesslang*(iat-1)+ilang) + & 308 & cpij*int1(kln,iat) 309 dos%fractions_pawt1(ikpt,iband,isppol,mbesslang*(iat-1)+ilang)= & 310 & dos%fractions_pawt1(ikpt,iband,isppol,mbesslang*(iat-1)+ilang) + & 311 & cpij*int2(kln,iat) 312 end if 313 314 end if 315 316 end do !ilmn 317 end do !jlmn 318 319 end do ! ilang 320 end do ! ispinor 321 end do ! iband 322 323 end do !iatom 324 325 if (mkmem/=0) ibg = ibg + my_nspinor*nband_cprj_k 326 call pawcprj_free(cprj_k) 327 ABI_DATATYPE_DEALLOCATE(cprj_k) 328 end do ! ikpt 329 end do ! isppol 330 331 ABI_DEALLOCATE(int1) 332 ABI_DEALLOCATE(int2) 333 ABI_DEALLOCATE(int1m2) 334 335 !Reduce data in case of parallelism 336 call timab(48,1,tsec) 337 call xmpi_sum(dos%fractions,comm_kptband,ierr) 338 if (prtdosm>=1.or.fatbands_flag==1) then 339 call xmpi_sum(dos%fractions_m,comm_kptband,ierr) 340 end if 341 if (paw_dos_flag==1) then 342 call xmpi_sum(dos%fractions_paw1,comm_kptband,ierr) 343 call xmpi_sum(dos%fractions_pawt1,comm_kptband,ierr) 344 end if 345 call timab(48,2,tsec) 346 if (mpi_enreg%paral_spinor==1) then 347 call xmpi_sum(dos%fractions,mpi_enreg%comm_spinor,ierr) 348 if (prtdosm>=1.or.fatbands_flag==1) then 349 call xmpi_sum(dos%fractions_m,mpi_enreg%comm_spinor,ierr) 350 end if 351 if (paw_dos_flag==1) then 352 call xmpi_sum(dos%fractions_paw1, mpi_enreg%comm_spinor,ierr) 353 call xmpi_sum(dos%fractions_pawt1, mpi_enreg%comm_spinor,ierr) 354 end if 355 end if 356 357 !Averaging: A quick hack for m-decomposed LDOS: 358 !BA: not valid in presence of spin-orbit coupling ! 359 if (prtdosm==1.and.fatbands_flag==0) then 360 ! if pawfatbnd is activated, one think in the cubic harmonics basis 361 ! whereas prtdosm=1 is in the spherical harmonics basis. 362 ! the following trick is done in order to have everything 363 ! in the complex spherical basis (not useful for pawfatbnd if we want to 364 ! have e.g t2g and eg d-orbitals). 365 do iat=1,dtset%natsph 366 do il = 0, mbesslang-1 367 do im = 1, il 368 dos%fractions_m(:,:,:,mbesslang**2*(iat-1)+il**2+il+1+im) = & 369 & (dos%fractions_m(:,:,:,mbesslang**2*(iat-1)+il**2+il+1+im) + & 370 & dos%fractions_m(:,:,:,mbesslang**2*(iat-1)+il**2+il+1-im))/2 371 dos%fractions_m(:,:,:,mbesslang**2*(iat-1)+il**2+il+1-im) = & 372 & dos%fractions_m(:,:,:,mbesslang**2*(iat-1)+il**2+il+1+im) 373 end do 374 end do 375 end do !iatom 376 end if 377 378 DBG_EXIT("COLL") 379 380 end subroutine partial_dos_fractions_paw