TABLE OF CONTENTS
ABINIT/pawmkrhoij [ Functions ]
NAME
pawmkrhoij
FUNCTION
Calculate the PAW quantities rhoij (augmentation occupancies)
Remember:for each atom, rho_ij=Sum_{n,k} {occ(n,k)*<Cnk|p_i><p_j|Cnk>}
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (FJ, 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.txt.
INPUTS
atindx1(natom)=index table for atoms, inverse of atindx
cprj(natom,mcprj)= wave functions projected with non-local projectors:
cprj_nk(i)=<p_i|Cnk> where p_i is a non-local projector.
dimcprj(natom)=array of dimensions of array cprj (ordered by atom-type)
istwfk(nkpt)=parameter that describes the storage of wfs
kptopt=option for the generation of k points
mband=maximum number of bands
mband_cprj=maximum number of bands used in the dimensioning of cprj array (usually mband/nproc_band)
mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
mkmem =number of k points treated by this node.
mpi_enreg=information about MPI parallelization
natom=number of atoms in cell
nband=number of bands for all k points
nkpt=number of k points
nspinor=number of spinorial components of the wavefunctions
nsppol=1 for unpolarized, 2 for spin-polarized
occ(mband*nkpt*nsppol)=occupation number for each band for each k
paral_kgb=Flag related to the kpoint-band-fft parallelism
paw_dmft <type(paw_dmft_type)>= paw+dmft related data
pawprtvol=control print volume and debugging output for PAW
unpaw=unit number for cprj PAW data (if used)
wtk(nkpt)=weight assigned to each k point
SIDE EFFECTS
pawrhoij(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data On input: arrays dimensions On output: pawrhoij(:)%rhoij_(lmn2_size,nspden)= Sum_{n,k} {occ(n,k)*conjugate[cprj_nk(ii)].cprj_nk(jj)} (non symetrized)
PARENTS
afterscfloop,scfcv,vtorho
CHILDREN
pawaccrhoij,pawcprj_alloc,pawcprj_free,pawcprj_gather_spin,pawcprj_get pawio_print_ij,pawrhoij_free,pawrhoij_init_unpacked pawrhoij_mpisum_unpacked,wrtout
NOTES
The cprj are distributed over band processors. Only the mod((iband-1)/mpi_enreg%bandpp,mpi_enreg%nproc_band) projectors are stored on each proc.
SOURCE
64 #if defined HAVE_CONFIG_H 65 #include "config.h" 66 #endif 67 68 #include "abi_common.h" 69 70 subroutine pawmkrhoij(atindx,atindx1,cprj,dimcprj,istwfk,kptopt,mband,mband_cprj,mcprj,mkmem,mpi_enreg,& 71 & natom,nband,nkpt,nspinor,nsppol,occ,paral_kgb,paw_dmft,& 72 & pawprtvol,pawrhoij,unpaw,usewvl,wtk) 73 74 75 use defs_basis 76 use defs_abitypes 77 use m_profiling_abi 78 use m_errors 79 use m_xmpi 80 81 use m_pawrhoij, only : pawrhoij_type, pawrhoij_init_unpacked, pawrhoij_mpisum_unpacked, pawrhoij_free 82 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, & 83 & pawcprj_gather_spin, pawcprj_free 84 use m_paw_io, only : pawio_print_ij 85 use m_paw_dmft, only : paw_dmft_type 86 87 !This section has been created automatically by the script Abilint (TD). 88 !Do not modify the following lines by hand. 89 #undef ABI_FUNC 90 #define ABI_FUNC 'pawmkrhoij' 91 use interfaces_14_hidewrite 92 use interfaces_32_util 93 use interfaces_65_paw, except_this_one => pawmkrhoij 94 !End of the abilint section 95 96 implicit none 97 98 !Arguments --------------------------------------------- 99 !scalars 100 integer,intent(in) :: kptopt,mband,mband_cprj,mcprj,mkmem,natom,nkpt,nspinor,nsppol 101 integer,intent(in) :: paral_kgb,pawprtvol,unpaw,usewvl 102 type(MPI_type),intent(in) :: mpi_enreg 103 !arrays 104 integer,intent(in) :: atindx(natom),atindx1(natom),dimcprj(natom),istwfk(nkpt) 105 integer,intent(in) :: nband(nkpt*nsppol) 106 real(dp),intent(in) :: occ(mband*nkpt*nsppol),wtk(nkpt) 107 type(pawcprj_type),target,intent(in) :: cprj(natom,mcprj) 108 type(paw_dmft_type),intent(in) :: paw_dmft 109 type(pawrhoij_type),intent(inout),target:: pawrhoij(:) 110 111 !Local variables --------------------------------------- 112 !scalars 113 integer,parameter :: max_nband_cprj=100 114 integer :: bdtot_index,cplex 115 integer :: iatom,iatom_tot,ib,ib1,iband,iband1,ibc1,ibg,ib_this_proc,ierr 116 integer :: ikpt,iorder_cprj,isppol,jb_this_proc,jbg,me,my_nspinor,nband_k,nband_k_cprj 117 integer :: nbandc1,nband_k_cprj_read,nband_k_cprj_used,nprocband,nrhoij,nsp2 118 integer :: option,spaceComm,use_nondiag_occup_dmft 119 logical :: locc_test,paral_atom,usetimerev 120 real(dp) :: wtk_k 121 character(len=4) :: wrt_mode 122 character(len=500) :: msg 123 124 !arrays 125 integer :: idum(0) 126 real(dp) :: occup(2) 127 character(len=8),parameter :: dspin(6)=(/"up ","down ","dens (n)","magn (x)","magn (y)","magn (z)"/) 128 type(pawcprj_type),allocatable :: cprj_tmp(:,:),cwaveprj(:,:),cwaveprjb(:,:) 129 type(pawcprj_type),pointer :: cprj_ptr(:,:) 130 type(pawrhoij_type),pointer :: pawrhoij_all(:) 131 132 !************************************************************************ 133 134 DBG_ENTER("COLL") 135 136 ABI_CHECK(mkmem/=0,"mkmem==0 not supported anymore!") 137 138 !Init MPI data 139 ! spaceComm=mpi_enreg%comm_cell 140 ! if(mpi_enreg%paral_kgb==1) spaceComm=mpi_enreg%comm_kpt 141 spaceComm=mpi_enreg%comm_kpt 142 me=mpi_enreg%me_kpt 143 144 !Check size of cprj 145 my_nspinor=max(1,nspinor/mpi_enreg%nproc_spinor) 146 if (mcprj/=my_nspinor*mband_cprj*mkmem*nsppol) then 147 msg=' wrong size for cprj !' 148 MSG_BUG(msg) 149 end if 150 151 !Check if cprj is distributed over bands 152 nprocband=(mband/mband_cprj) 153 if (paral_kgb==1.and.nprocband/=mpi_enreg%nproc_band) then 154 msg=' mband/mband_cprj must be equal to nproc_band!' 155 MSG_BUG(msg) 156 end if 157 if (paw_dmft%use_sc_dmft/=0.and.nprocband/=1) then 158 write(msg,'(4a,e14.3,a)') ch10,& 159 & ' Parallelization over bands is not yet compatible with self-consistency in DMFT !',ch10,& 160 & ' Calculation is thus restricted to nstep =1.' 161 MSG_WARNING(msg) 162 end if 163 164 if( usewvl==1 .and. (nprocband/=1)) then 165 write(msg,'(2a)') ch10,& 166 & ' ERROR: parallelization over bands is not compatible with WAVELETS' 167 MSG_ERROR(msg) 168 end if 169 170 !Initialise and check dmft variables 171 if(paw_dmft%use_sc_dmft/=0) then 172 nbandc1=(paw_dmft%mbandc-1)*paw_dmft%use_sc_dmft+1 173 else 174 nbandc1=1 175 end if 176 177 !Size of pawrhoij datastructure 178 nrhoij=size(pawrhoij) 179 180 !Check if pawrhoij is distributed over atomic sites 181 paral_atom=(nrhoij/=natom.and.mpi_enreg%nproc_atom>1) 182 if (paral_atom.and.nrhoij/=mpi_enreg%my_natom) then 183 msg=' Size of pawrhoij should be natom or my_natom !' 184 MSG_BUG(msg) 185 end if 186 187 !Allocate temporary cwaveprj storage 188 ABI_DATATYPE_ALLOCATE(cwaveprj,(natom,nspinor)) 189 call pawcprj_alloc(cwaveprj,0,dimcprj) 190 if(paw_dmft%use_sc_dmft/=0) then 191 ABI_DATATYPE_ALLOCATE(cwaveprjb,(natom,nspinor)) 192 call pawcprj_alloc(cwaveprjb,0,dimcprj) 193 end if 194 195 !Initialize temporary file (if used) 196 iorder_cprj=0 197 198 !Build and initialize unpacked rhoij (to be computed here) 199 call pawrhoij_init_unpacked(pawrhoij) 200 201 !If pawrhoij is MPI-distributed over atomic sites, gather it 202 if (paral_atom) then 203 ABI_DATATYPE_ALLOCATE(pawrhoij_all,(natom)) 204 else 205 pawrhoij_all => pawrhoij 206 end if 207 208 !LOOP OVER SPINS 209 option=1 210 usetimerev=(kptopt>0.and.kptopt<3) 211 bdtot_index=0;ibg=0;jbg=0 212 do isppol=1,nsppol 213 214 ! LOOP OVER k POINTS 215 do ikpt=1,nkpt 216 217 nband_k=nband(ikpt+(isppol-1)*nkpt) 218 nband_k_cprj=nband_k/nprocband 219 wtk_k=wtk(ikpt) 220 221 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then 222 bdtot_index=bdtot_index+nband_k 223 cycle 224 end if 225 226 cplex=2;if (istwfk(ikpt)>1) cplex=1 227 228 ! In case of spinors parallelism, need some extra storage 229 if (mpi_enreg%paral_spinor==1) then 230 nband_k_cprj_used=min(max_nband_cprj,nband_k_cprj) 231 ABI_DATATYPE_ALLOCATE(cprj_tmp,(natom,my_nspinor*nband_k_cprj_used)) 232 ABI_DATATYPE_ALLOCATE(cprj_ptr,(natom, nspinor*nband_k_cprj_used)) 233 call pawcprj_alloc(cprj_tmp,0,dimcprj) 234 call pawcprj_alloc(cprj_ptr,0,dimcprj) 235 else 236 cprj_ptr => cprj 237 end if 238 239 ! LOOP OVER BANDS 240 ib_this_proc=0;jb_this_proc=0 241 do ib=1,nband_k 242 iband=bdtot_index+ib 243 244 ! Parallelization: treat only some bands 245 if(xmpi_paral==1)then 246 if (paral_kgb==1) then 247 if (mod((ib-1)/mpi_enreg%bandpp,mpi_enreg%nproc_band)/=mpi_enreg%me_band) cycle 248 else 249 if (mpi_enreg%proc_distrb(ikpt,ib,isppol)/=me) cycle 250 end if 251 end if 252 ib_this_proc=ib_this_proc+1 253 254 ! In case of spinors parallelism, gather cprj because we need both components together 255 ! We do that nband_k_cprj_used by nband_k_cprj_used bands 256 if (mpi_enreg%paral_spinor==1) then 257 jb_this_proc=jb_this_proc+1 258 if (mod(jb_this_proc,nband_k_cprj_used)==1) then 259 ib_this_proc=1 260 nband_k_cprj_read=nband_k_cprj_used 261 if (nband_k_cprj<jb_this_proc+nband_k_cprj_used-1) nband_k_cprj_read=nband_k_cprj-jb_this_proc+1 262 call pawcprj_get(atindx1,cprj_tmp,cprj,natom,jb_this_proc,jbg,ikpt,iorder_cprj,isppol,& 263 & mband_cprj,mkmem,natom,nband_k_cprj_read,nband_k_cprj,my_nspinor,nsppol,unpaw,& 264 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 265 call pawcprj_gather_spin(cprj_tmp,cprj_ptr,natom,nband_k_cprj_read,my_nspinor,nspinor,& 266 & mpi_enreg%comm_spinor,ierr) 267 end if 268 end if 269 270 ! DMFT: LOOP ON ADDITIONAL BANDS 271 do ibc1=1,nbandc1 272 ! check if dmft and occupations 273 ! write(std_out,*) 'ib,ibc1 ',ib,ibc1 274 275 ! DMFT stuff: extract cprj and occupations for additional band 276 if(paw_dmft%use_sc_dmft /= 0) then 277 ib1 = paw_dmft%include_bands(ibc1) 278 ! write(std_out,*) 'use_sc_dmft=1 ib,ib1',ib,ib1 279 iband1 = bdtot_index+ib1 280 ! write(std_out,*) 'ib, ib1 ',paw_dmft%band_in(ib),paw_dmft%band_in(ib1) 281 if(paw_dmft%band_in(ib)) then 282 if(.not.paw_dmft%band_in(ib1)) stop 283 use_nondiag_occup_dmft = 1 284 occup(1) = paw_dmft%occnd(1,ib,ib1,ikpt,isppol) 285 if(nspinor==2) occup(2) = paw_dmft%occnd(2,ib,ib1,ikpt,isppol) 286 if(nspinor==1) occup(2) = zero 287 locc_test = abs(paw_dmft%occnd(1,ib,ib1,ikpt,isppol))+abs(paw_dmft%occnd(2,ib,ib1,ikpt,isppol))>tol8 288 ! write(std_out,*) 'use_sc_dmft=1,band_in(ib)=1, ib,ibc1',ib,ib1,locc_test 289 if (locc_test .or. mkmem == 0) then 290 call pawcprj_get(atindx1,cwaveprjb,cprj_ptr,natom,ib1,ibg,ikpt,iorder_cprj,isppol,& 291 & mband_cprj,mkmem,natom,1,nband_k_cprj,nspinor,nsppol,unpaw,& 292 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 293 end if 294 else 295 use_nondiag_occup_dmft = 0 296 locc_test = (abs(occ(iband))>tol8) 297 occup(1) = occ(iband) 298 if(ibc1 /= 1 .and. .not.(paw_dmft%band_in(ib))) cycle 299 end if 300 else ! nbandc1=1 301 use_nondiag_occup_dmft=0 302 locc_test = (abs(occ(iband))>tol8) 303 occup(1) = occ(iband) 304 end if 305 306 ! Extract cprj for current band 307 ! Must read cprj when mkmem=0 (even if unused) to have right pointer inside _PAW file 308 if (locc_test.or.mkmem==0) then 309 call pawcprj_get(atindx1,cwaveprj,cprj_ptr,natom,ib_this_proc,ibg,ikpt,iorder_cprj,isppol,& 310 & mband_cprj,mkmem,natom,1,nband_k_cprj,nspinor,nsppol,unpaw,& 311 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 312 end if 313 314 ! Accumulate contribution from (occupied) current band 315 if (locc_test) then 316 if(use_nondiag_occup_dmft == 1) then 317 call pawaccrhoij(atindx,cplex,cwaveprj,cwaveprjb,0,isppol,nrhoij,natom,& 318 & nspinor,occup(1),option,pawrhoij_all,usetimerev,wtk_k,& 319 & occ_k_2=occup(2)) 320 else 321 call pawaccrhoij(atindx,cplex,cwaveprj,cwaveprj ,0,isppol,nrhoij,natom,& 322 & nspinor,occup(1),option,pawrhoij_all,usetimerev,wtk_k) 323 end if 324 end if 325 end do ! ib1c 326 end do ! ib 327 328 if (mpi_enreg%paral_spinor==1) then 329 call pawcprj_free(cprj_tmp) 330 call pawcprj_free(cprj_ptr) 331 ABI_DATATYPE_DEALLOCATE(cprj_tmp) 332 ABI_DATATYPE_DEALLOCATE(cprj_ptr) 333 else 334 nullify(cprj_ptr) 335 end if 336 337 bdtot_index=bdtot_index+nband_k 338 if (mkmem/=0) then 339 if (mpi_enreg%paral_spinor==0) then 340 ibg=ibg+ nspinor*nband_k_cprj 341 else 342 jbg=jbg+my_nspinor*nband_k_cprj 343 end if 344 end if 345 346 end do ! ikpt 347 end do ! isppol 348 349 !deallocate temporary cwaveprj/cprj storage 350 call pawcprj_free(cwaveprj) 351 ABI_DATATYPE_DEALLOCATE(cwaveprj) 352 if(paw_dmft%use_sc_dmft/=0) then 353 call pawcprj_free(cwaveprjb) 354 ABI_DATATYPE_DEALLOCATE(cwaveprjb) 355 end if 356 357 !MPI: need to exchange rhoij_ between procs 358 if (paral_kgb==1.and.nprocband>1) then 359 call pawrhoij_mpisum_unpacked(pawrhoij_all,spaceComm,comm2=mpi_enreg%comm_band) 360 else 361 call pawrhoij_mpisum_unpacked(pawrhoij_all,spaceComm) 362 end if 363 364 !In case of distribution over atomic sites, dispatch rhoij 365 if (paral_atom) then 366 do iatom=1,nrhoij 367 iatom_tot=mpi_enreg%my_atmtab(iatom) 368 pawrhoij(iatom)%rhoij_(:,:)=pawrhoij_all(iatom_tot)%rhoij_(:,:) 369 end do 370 call pawrhoij_free(pawrhoij_all) 371 ABI_DATATYPE_DEALLOCATE(pawrhoij_all) 372 end if 373 374 !Print info 375 if (abs(pawprtvol)>=1) then 376 wrt_mode='COLL';if (paral_atom) wrt_mode='PERS' 377 do iatom=1,nrhoij 378 iatom_tot=iatom;if (paral_atom) iatom_tot=mpi_enreg%my_atmtab(iatom) 379 if (pawprtvol>=0.and.iatom_tot/=1.and.iatom_tot/=natom) cycle 380 nsp2=pawrhoij(iatom)%nsppol;if (pawrhoij(iatom)%nspden==4) nsp2=4 381 write(msg, '(4a,i3,a)') ch10," PAW TEST:",ch10,& 382 & ' ====== Values of RHOIJ in pawmkrhoij (iatom=',iatom_tot,') ======' 383 if (pawrhoij(iatom)%nspden==2.and.pawrhoij(iatom)%nsppol==1) write(msg,'(3a)') trim(msg),ch10,& 384 & ' (antiferromagnetism case: only one spin component)' 385 call wrtout(std_out,msg,wrt_mode) 386 do isppol=1,nsp2 387 if (pawrhoij(iatom)%nspden/=1) then 388 write(msg,'(3a)') ' Component ',trim(dspin(isppol+2*(pawrhoij(iatom)%nspden/4))),':' 389 call wrtout(std_out,msg,wrt_mode) 390 end if 391 option=2;if (pawrhoij(iatom)%cplex==2.and.pawrhoij(iatom)%nspinor==1) option=1 392 call pawio_print_ij(std_out,pawrhoij(iatom)%rhoij_(:,isppol),pawrhoij(iatom)%lmn2_size,& 393 & pawrhoij(iatom)%cplex,pawrhoij(iatom)%lmn_size,-1,idum,0,pawprtvol,idum,& 394 & -1._dp,1,opt_sym=option,mode_paral=wrt_mode) 395 end do 396 end do 397 end if 398 399 DBG_EXIT("COLL") 400 401 end subroutine pawmkrhoij