TABLE OF CONTENTS
ABINIT/datafordmft [ Functions ]
NAME
datafordmft
FUNCTION
Compute psichi (and print some data for check)
COPYRIGHT
Copyright (C) 2005-2018 ABINIT group (BAmadon) 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
cryst_struc <type(crystal_t)>=crystal structure data -gprimd(3,3)=dimensional reciprocal space primitive translations -indsym(4,nsym,natom)=indirect indexing array for atom labels -symrec(3,3,nsym)=symmetry operations in reciprocal space - nsym= number of symetry operations cprj(natom,nspinor*mband*mkmem*nsppol)= <p_lmn|Cnk> coefficients for each WF |Cnk> and each |p_lmn> non-local projector dimcprj(natom) = dimension for cprj dtset <type(dataset_type)>=all input variables for this dataset eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree) fermie= Fermi energy lda_occup <type(oper_type)> = occupations in the correlated orbitals in LDA mband=maximum number of bands mkmem =number of k points treated by this node mpi_enreg=information about MPI parallelization nkpt=number of k points. my_nspinor=number of spinorial components of the wavefunctions (on current proc) nsppol=1 for unpolarized, 2 for spin-polarized occ(mband*nkpt*nsppol) = occupancies of KS states. paw_dmft <type(paw_dmft_type)>= paw+dmft related data paw_ij(natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels pawang <type(pawang)>=paw angular mesh and related data pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data psps <type(pseudopotential_type)>=variables related to pseudopotentials unpaw = file number for cprj
OUTPUT
paw_dmft%psichi(nsppol,nkpt,mband,my_nspinor,dtset%natom,(2*maxlpawu+1))): projections <Psi|chi> paw_dmft%eigen(paw_dmft%nsppol,paw_dmft%nkpt,paw_dmft%mband)
SIDE EFFECTS
(only writing, printing)
NOTES
PARENTS
outscfcv,vtorho
CHILDREN
SOURCE
59 #if defined HAVE_CONFIG_H 60 #include "config.h" 61 #endif 62 63 64 #include "abi_common.h" 65 66 subroutine datafordmft(cryst_struc,cprj,dimcprj,dtset,eigen,fermie,& 67 & lda_occup,mband,mband_cprj,mkmem,mpi_enreg,nkpt,my_nspinor,nsppol,occ,& 68 & paw_dmft,paw_ij,pawang,pawtab,psps,usecprj,unpaw,nbandkss) 69 70 use defs_basis 71 use defs_datatypes 72 use defs_abitypes 73 use defs_wvltypes 74 use m_profiling_abi 75 use m_errors 76 use m_xmpi 77 78 use m_io_tools, only : open_file 79 use m_matlu, only: matlu_type,init_matlu,sym_matlu,copy_matlu,print_matlu,diff_matlu,destroy_matlu 80 use m_crystal, only : crystal_t 81 use m_oper, only : oper_type 82 83 use m_pawang, only : pawang_type 84 use m_pawtab, only : pawtab_type 85 use m_paw_ij, only : paw_ij_type 86 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_free 87 use m_paw_dmft, only: paw_dmft_type 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 'datafordmft' 93 use interfaces_14_hidewrite 94 use interfaces_32_util 95 use interfaces_68_dmft, except_this_one => datafordmft 96 !End of the abilint section 97 98 implicit none 99 100 !Arguments ------------------------------------ 101 !scalars 102 integer,intent(in) :: mband,mband_cprj,mkmem 103 integer,intent(in) :: nkpt,my_nspinor,nsppol 104 integer,intent(in) :: unpaw,usecprj 105 integer, optional, intent(in) :: nbandkss 106 real(dp),intent(in) :: fermie 107 type(MPI_type),intent(in) :: mpi_enreg 108 type(dataset_type),intent(in) :: dtset 109 type(oper_type), intent(inout) :: lda_occup !vz_i 110 type(pawang_type),intent(in) :: pawang 111 type(pseudopotential_type),intent(in) :: psps 112 type(crystal_t),intent(in) :: cryst_struc 113 !arrays 114 integer, intent(in) :: dimcprj(cryst_struc%natom) 115 real(dp),intent(in) :: eigen(mband*nkpt*nsppol) 116 real(dp),intent(in) :: occ(mband*nkpt*nsppol) 117 type(paw_ij_type),intent(in) :: paw_ij(:) 118 ! type(pawcprj_type) :: cprj(cryst_struc%natom,my_nspinor*mband*mkmem*nsppol) 119 type(pawcprj_type), intent(in) :: cprj(cryst_struc%natom,my_nspinor*mband_cprj*mkmem*nsppol*usecprj) 120 type(paw_dmft_type), intent(inout) :: paw_dmft 121 type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw) 122 !Local variables------------------------------- 123 !scalars 124 integer :: band_index,dimpsichi,facpara 125 integer :: iat,iatom,ib,iband,ibandc,ibg,icat,icount_proj_ilmn,idijeff,ierr,ierrr,ikpt 126 integer :: ilmn,im,im1,iorder_cprj,ispinor,ispinor1,isppol,itypat,ilmn1 127 integer :: jj1,ldim,lmn_size,lpawu 128 integer :: m1,m1_t2g,m1_t2g_mod,maxnproju,me,natom,nband_k,nband_k_cprj 129 integer :: nbandi,nbandf,nnn,nprocband,nsploop,option,opt_renorm,spaceComm,unt 130 real(dp) :: ph0phiint_used 131 character(len=500) :: message 132 !arrays 133 real(dp) :: chinorm 134 complex(dpc), allocatable :: buffer1(:) 135 type(matlu_type), allocatable :: loc_occ_check(:) 136 type(matlu_type), allocatable :: loc_norm_check(:) 137 type(matlu_type), allocatable :: xocc_check(:) 138 type(matlu_type), allocatable :: xnorm_check(:) 139 type(matlu_type), allocatable :: matlu_temp(:) 140 logical :: lprojchi,t2g 141 integer,parameter :: spinor_idxs(2,4)=RESHAPE((/1,1,2,2,1,2,2,1/),(/2,4/)) 142 type(pawcprj_type),allocatable :: cwaveprj(:,:) 143 144 !************************************************************************ 145 146 !DBG_ENTER("COLL") 147 !Fake test to keep fermie as argument. REMOVE IT AS SOON AS POSSIBLE ... 148 if(fermie>huge(zero))chinorm=zero 149 150 facpara=1 !mpi_enreg%nproc 151 if(abs(dtset%pawprtvol)>=3) then 152 write(message,*) " number of k-points used is nkpt=nkpt ", nkpt 153 call wrtout(std_out, message,'COLL') 154 write(message,*) " warning: parallelised version ", nkpt 155 call wrtout(std_out, message,'COLL') 156 write(message,*) " weights k-points used is wtk=wtk" 157 call wrtout(std_out, message,'COLL') 158 end if 159 160 if(usecprj==0) then 161 write(message,*) " usecprj=0 : BUG in datafordmft",usecprj 162 MSG_BUG(message) 163 end if 164 165 if(my_nspinor/=dtset%nspinor) then 166 write(message,*) " my_nspinor=/dtset%nspinor, datafordmft not working in this case",& 167 & my_nspinor,dtset%nspinor 168 MSG_ERROR(message) 169 end if 170 171 172 !do ib=1,my_nspinor*mband_cprj*mkmem*nsppol*usecprj 173 !write(std_out,'(a,i6,3e16.7)') "cprj",ib,cprj(1,ib)%cp(1,19),cprj(1,ib)%cp(2,19),cprj(1,ib)%cp(1,19)**2+cprj(1,ib)%cp(2,19)**2 174 !enddo 175 176 !----------------------------------- MPI------------------------------------- 177 178 !Init parallelism 179 spaceComm=mpi_enreg%comm_cell 180 if(mpi_enreg%paral_kgb==1) spaceComm=mpi_enreg%comm_kpt 181 me=mpi_enreg%me_kpt 182 iorder_cprj=0 183 ABI_CHECK(dtset%mkmem/=0,"mkmem==0 not supported anymore!") 184 !todo_ab: extract cprj from file unpaw in the following.. 185 !call leave_new('COLL') 186 187 !----------------------------------- MPI------------------------------------- 188 189 nbandi=paw_dmft%dmftbandi 190 nbandf=paw_dmft%dmftbandf 191 lprojchi=.false. 192 lprojchi=.true. 193 t2g=(paw_dmft%dmftqmc_t2g==1) 194 natom=cryst_struc%natom 195 196 !if(mpi_enreg%me==0) write(7886,*) "in datafordmft", mpi_enreg%me, mpi_enreg%nproc 197 !if(mpi_enreg%me==1) write(7887,*) "in datafordmft", mpi_enreg%me, mpi_enreg%nproc 198 !if(mpi_enreg%me==2) write(7888,*) "in datafordmft", mpi_enreg%me, mpi_enreg%nproc 199 write(message,'(2a)') ch10,& 200 & ' == Prepare data for DMFT calculation ' 201 call wrtout(std_out,message,'COLL') 202 if(abs(dtset%pawprtvol)>=3) then 203 write(message, '(a,a)' ) ch10,& 204 & '---------------------------------------------------------------' 205 ! call wrtout(ab_out,message,'COLL') 206 ! call wrtout(std_out, message,'COLL') 207 call wrtout(std_out, message,'COLL') 208 write(message, '(a,a,a,a,a,a,a,a,a,a,a,a)' ) ch10,& 209 & ' Print useful data (as a check)',ch10,& 210 & ' - Overlap of KS wfc with atomic orbital inside sphere',ch10,& 211 & ' - Eigenvalues',ch10,& 212 & ' - Weights of k-points',ch10,& 213 & ' - Number of spins ',ch10,& 214 & ' - Number of states' 215 ! call wrtout(ab_out,message,'COLL') 216 ! call wrtout(std_out, message,'COLL') 217 call wrtout(std_out, message,'COLL') 218 write(message, '(a,a)' ) ch10,& 219 & '---------------------------------------------------------------' 220 end if 221 if(dtset%nstep==0.and.dtset%nbandkss==0) then 222 message = 'nstep should be greater than 1' 223 MSG_BUG(message) 224 end if 225 226 !********************* Max Values for U terms. 227 !maxlpawu=0 228 maxnproju=0 229 do iatom=1,natom 230 if(pawtab(dtset%typat(iatom))%lpawu.ne.-1 .and. pawtab(dtset%typat(iatom))%nproju.gt.maxnproju)& 231 & maxnproju=pawtab(dtset%typat(iatom))%nproju 232 end do 233 !***************** in forlb.eig 234 if(me.eq.0.and.abs(dtset%pawprtvol)>=3) then 235 if (open_file('forlb.eig',message,newunit=unt,form='formatted',status='unknown') /= 0) then 236 MSG_ERROR(message) 237 end if 238 rewind(unt) 239 write(unt,*) "Number of bands, spins, and k-point; and spin-orbit flag" 240 write(unt,*) mband,nsppol,nkpt,my_nspinor,nbandi,nbandf 241 write(unt,*) " For each k-point, eigenvalues for each band" 242 write(unt,*) (dtset%wtk(ikpt)*facpara,ikpt=1,nkpt) 243 band_index=0 244 do isppol=1,nsppol 245 write(unt,*) " For spin" 246 write(unt,*) isppol 247 do ikpt=1,nkpt 248 nband_k=dtset%nband(ikpt+(isppol-1)*nkpt) 249 ibandc=0 250 write(unt,*) " For k-point" 251 write(unt,*) ikpt 252 do iband=1,mband 253 if(paw_dmft%band_in(iband)) then 254 ibandc=ibandc+1 255 write(unt, '(2i6,4x,f20.15)' ) ibandc,ikpt,eigen(iband+band_index)*2.d0 256 end if 257 end do 258 band_index=band_index+nband_k 259 end do 260 end do 261 close(unt) 262 end if ! proc=me 263 264 !== put eigen into eigen_lda 265 band_index=0 266 do isppol=1,nsppol 267 do ikpt=1,nkpt 268 nband_k=dtset%nband(ikpt+(isppol-1)*nkpt) 269 ibandc=0 270 do iband=1,mband 271 if(paw_dmft%band_in(iband)) then 272 ibandc=ibandc+1 273 paw_dmft%eigen_lda(isppol,ikpt,ibandc)=eigen(iband+band_index) ! in Ha 274 ! paw_dmft%eigen_lda(isppol,ikpt,ibandc)=fermie 275 end if 276 end do 277 band_index=band_index+nband_k 278 end do 279 end do 280 281 if(abs(dtset%pawprtvol)>=3) then 282 write(message, '(a,a)' ) ch10,& 283 & ' datafordmft : eigenvalues written' 284 call wrtout(std_out, message,'COLL') 285 end if 286 !========================================================================== 287 !***************** Compute <Psi|Chi>=\sum_{proja} <Psi|P_a><phi_a|Chi> 288 !========================================================================== 289 !write(std_out,*) "size(cprj,dim=1)",size(cprj,dim=1),size(cprj,dim=2),dtset%mband,dtset%mkmem,dtset%nkpt 290 291 !Allocate temporary cwaveprj storage 292 ABI_DATATYPE_ALLOCATE(cwaveprj,(natom,my_nspinor)) 293 !write(std_out,*) "before alloc cprj" 294 !write(std_out,*) size(cwaveprj,dim=1),size(cwaveprj,dim=2),size(dimcprj,dim=1) 295 296 call pawcprj_alloc(cwaveprj,0,dimcprj) 297 !write(std_out,*) "after alloc cprj" 298 299 nprocband=(mband/mband_cprj) 300 ibg=0 301 do isppol=1,nsppol 302 do ikpt=1,nkpt 303 nband_k=dtset%nband(ikpt+(isppol-1)*nkpt) 304 nband_k_cprj=nband_k/nprocband 305 306 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) cycle 307 308 ! write(2011,*) ikpt 309 ! ibsp=ibg 310 ibandc=0 311 ! LOOP OVER BANDS 312 ib=0 313 do iband=1,nband_k 314 315 if(paw_dmft%band_in(iband)) ibandc=ibandc+1 316 317 ! Parallelization: treat only some bands 318 if (dtset%paral_kgb==1) then 319 if (mod((iband-1)/mpi_enreg%bandpp,mpi_enreg%nproc_band)/=mpi_enreg%me_band) cycle 320 else 321 if (mpi_enreg%proc_distrb(ikpt,iband,isppol)/=me) cycle 322 end if 323 ib=ib+1 324 325 do ispinor=1,my_nspinor 326 ! ibsp =~ (ikpt-1)*nband*my_nspinor+iband 327 ! ibsp=ibsp+1 328 icat=1 329 ! write(std_out,*) isppol,ikpt,iband,ispinor 330 iat=0 ! to declare 331 do itypat=1,dtset%ntypat 332 lmn_size=pawtab(itypat)%lmn_size 333 ! write(std_out,*) isppol,ikpt,iband,ispinor 334 do iatom=icat,icat+cryst_struc%nattyp(itypat)-1 335 lpawu=pawtab(dtset%typat(iatom))%lpawu 336 ! ---------- t2g case 337 if(paw_dmft%dmftqmc_t2g==1.and.lpawu/=-1) then 338 if(lpawu==2) then 339 ! lpawu==2 must be chosen in input and thus in 340 ! pawtab. On the contrary, paw_dmft now has 341 ! lpawu=1 342 m1_t2g=0 ! index for psichi which has a dimension 3 343 else 344 write(message,'(a,a,i4,i4,2a)') ch10,& 345 & ' Wrong use of dmftqmc_t2g',paw_dmft%dmftqmc_t2g,lpawu,ch10,& 346 & ' Action: desactivate qmftqmc_t2g or use lpawu=1' 347 MSG_ERROR(message) 348 end if 349 end if 350 ! ---------- t2g case 351 ! if(isppol==2) write(std_out,*) "ee",size(cprj(iatom,ibsp)%cp(:,:)) 352 353 iat=iat+1 354 jj1=0 355 if(lpawu.ne.-1) then 356 357 call pawcprj_get(cryst_struc%atindx1,cwaveprj,cprj,natom,ib,ibg,ikpt,& 358 & iorder_cprj,isppol,mband_cprj,mkmem,natom,1,nband_k_cprj,& 359 & my_nspinor,nsppol,unpaw,mpicomm=mpi_enreg%comm_kpt,& 360 & proc_distrb=mpi_enreg%proc_distrb) 361 ! write(std_out,*) "M-2",isppol,ikpt,iband,iatom,& 362 ! & (cwaveprj(iatom,ispinor)%cp(1,13)**2+cwaveprj(iatom,ispinor)%cp(2,13)**2) ,ibsp 363 364 ! chinorm=(pawtab(itypat)%phiphjint(1)) 365 chinorm=1.d0 366 ! write(std_out,*) isppol,ikpt,iband,ispinor,iat 367 do ilmn=1,lmn_size 368 ! write(std_out,*) ilmn 369 ! ------------ Select l=lpawu. 370 if (psps%indlmn(1,ilmn,itypat)==lpawu) then 371 ! ------------ Check that the band is choosen (within nbandi and nbandf) 372 if(paw_dmft%band_in(iband)) then 373 ! if(ilmn==13) then 374 ! write(std_out,*) "M-2c",isppol,ikpt,iband,iatom 375 ! write(std_out,*) "M-2b",isppol,ikpt,ibandc,& 376 ! & (cwaveprj(iatom,ispinor)%cp(1,13)**2+cwaveprj(iatom,ispinor)%cp(2,13)**2) 377 ! endif 378 ! write(std_out,*) "inside paw_dmft%band_in",iband 379 jj1=jj1+1 380 if(jj1>pawtab(dtset%typat(iatom))%nproju*(2*lpawu+1)) then 381 write(message,'(a,a,a,a)') ch10,& 382 & ' jj1 is not correct in datafordmft',ch10,& 383 & ' Action: CONTACT Abinit group' 384 MSG_BUG(message) 385 end if ! jj1 386 icount_proj_ilmn=psps%indlmn(3,ilmn,itypat) ! n: nb of projector 387 ! write(std_out,*) "icount_proj_ilmn",icount_proj_ilmn 388 m1=psps%indlmn(2,ilmn,itypat)+psps%indlmn(1,ilmn,itypat)+1 389 ! ---- if lprochi=true do the sum over every projectors 390 ! ---- if lprochi=false: . do the sum over only ground state projector 391 ! ---- . and take ph0phiint(icount_proj_ilmn)=1 392 393 ! call pawcprj_get(cryst_struc%atindx1,cwaveprj,cprj,natom,& 394 ! & ib,ibg,ikpt,iorder_cprj,isppol,mband_cprj,mkmem,& 395 ! & natom,1,nband_k_cprj,my_nspinor,nsppol,unpaw,& 396 ! & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 397 398 if(lprojchi.or.icount_proj_ilmn==1) then 399 if(.not.lprojchi) ph0phiint_used=one 400 if(lprojchi) ph0phiint_used=pawtab(itypat)%ph0phiint(icount_proj_ilmn) 401 if(paw_dmft%dmftqmc_t2g==1) then ! t2g case 402 403 if(m1==1.or.m1==2.or.m1==4) then ! t2g case 404 m1_t2g=m1_t2g+1 ! t2g case1 405 m1_t2g_mod=mod(m1_t2g-1,3)+1 406 ! write(std_out,*) "M0",isppol,ikpt,iband,ilmn,cprj(iatom,ibsp)%cp(1,ilmn) 407 ! write(std_out,*) "M1",m1,m1_t2g,iband,ilmn,icount_proj_ilmn 408 paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1_t2g_mod)=& ! t2g case 409 & paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1_t2g_mod)+& ! t2g case 410 ! & cmplx(cprj(iatom,ibsp)%cp(1,ilmn)*ph0phiint_used,& ! t2g case 411 ! & cprj(iatom,ibsp)%cp(2,ilmn)*ph0phiint_used,kind=dp) ! t2g case 412 & cmplx(cwaveprj(iatom,ispinor)%cp(1,ilmn)*ph0phiint_used,& ! t2g case 413 & cwaveprj(iatom,ispinor)%cp(2,ilmn)*ph0phiint_used,kind=dp) ! t2g case 414 end if !t2g case 415 else 416 paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1)=& 417 & paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1)+& 418 ! & cmplx(cprj(iatom,ibsp)%cp(1,ilmn)*ph0phiint_used,& 419 ! & cprj(iatom,ibsp)%cp(2,ilmn)*ph0phiint_used,kind=dp) 420 & cmplx(cwaveprj(iatom,ispinor)%cp(1,ilmn)*ph0phiint_used,& 421 & cwaveprj(iatom,ispinor)%cp(2,ilmn)*ph0phiint_used,kind=dp) 422 423 ! if(ibandc==3.and.iat==1.and.m1==1) then 424 ! write(std_out,'(a,3i5)') "psichi integers",isppol,ikpt,ispinor 425 ! write(std_out,'(a,2i5,2e16.7)') "psichi IB3 iAT1 IM1",ilmn,icount_proj_ilmn,& 426 ! & real(paw_dmft%psichi(isppol,ikpt,3,ispinor,1,1)), imag(paw_dmft%psichi(isppol,ikpt,3,ispinor,1,1)) 427 ! write(std_out,'(a,2i5,2e16.7)') "cwaveprj IB3 iAT1 IM1",ilmn,icount_proj_ilmn,cwaveprj(iatom,ispinor)%cp(1,ilmn) & 428 ! & , cwaveprj(iatom,ispinor)%cp(2,ilmn) 429 ! endif 430 end if 431 end if ! lprojchi=.true. (always) 432 end if ! nband belongs to paw_dmft%band_in 433 end if ! L=lpawu 434 end do !ilmn : loop over L,M,N 435 end if ! If lpawu/=-1 436 end do ! iatom 437 icat=icat+cryst_struc%nattyp(itypat) 438 end do ! itypat 439 end do ! ispinor 440 end do !iband 441 ibg=ibg+nband_k_cprj*my_nspinor 442 ! bdtot_index=bdtot_index+nband_k ! useless (only for occ) 443 end do !ikpt 444 end do ! isppol 445 !do isppol=1,nsppol 446 !do ikpt=1,nkpt 447 !do ispinor=1,my_nspinor 448 !write(std_out,*) "psichi integers",isppol,ikpt,ispinor 449 !write(std_out,*) "psichi IB3 iAT1 IM1",& 450 !& real(paw_dmft%psichi(isppol,ikpt,3,ispinor,1,1)), imag(paw_dmft%psichi(isppol,ikpt,3,ispinor,1,1)) 451 ! 452 !enddo 453 !enddo 454 !enddo 455 !call leave_new('COLL') 456 if(abs(dtset%pawprtvol)>=3) then 457 write(message,*) "chinorm used here =",chinorm 458 call wrtout(std_out, message,'COLL') 459 end if 460 461 !deallocate temporary cwaveprj/cprj storage 462 call pawcprj_free(cwaveprj) 463 ABI_DATATYPE_DEALLOCATE(cwaveprj) 464 465 !========================================================================== 466 !********************* Gather information for MPI before printing 467 !========================================================================== 468 469 dimpsichi=2*nsppol*nkpt*mband*my_nspinor*natom*(2*paw_dmft%maxlpawu+1) 470 ABI_ALLOCATE(buffer1,(dimpsichi)) 471 nnn=0 472 !write(176,*) "beg",psichi 473 do isppol=1,nsppol 474 do ikpt=1,nkpt 475 do ibandc=1,paw_dmft%mbandc 476 do ispinor=1,my_nspinor 477 do iat=1,natom 478 do m1=1,2*paw_dmft%maxlpawu+1 479 ! do m=1,2 480 nnn=nnn+1 481 buffer1(nnn)=paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1) 482 ! enddo 483 end do 484 end do 485 end do 486 end do 487 end do 488 end do 489 call xmpi_barrier(spaceComm) 490 call xmpi_sum(buffer1,spaceComm ,ierr) 491 if (dtset%paral_kgb==1.and.nprocband>1) then 492 call xmpi_sum(buffer1,mpi_enreg%comm_band,ierr) !Build sum over band processors 493 end if 494 call xmpi_barrier(spaceComm) 495 nnn=0 496 do isppol=1,nsppol 497 do ikpt=1,nkpt 498 do ibandc=1,paw_dmft%mbandc 499 do ispinor=1,my_nspinor 500 do iat=1,natom 501 do m1=1,2*paw_dmft%maxlpawu+1 502 ! do m=1,2 503 nnn=nnn+1 504 paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1)=buffer1(nnn) 505 ! enddo 506 end do 507 end do 508 end do 509 end do 510 end do 511 end do 512 ABI_DEALLOCATE(buffer1) 513 !write(177,*) "end",psichi 514 515 !do isppol=1,nsppol 516 !do ikpt=1,nkpt 517 !do ibandc=1,paw_dmft%mbandc 518 !do ispinor=1,my_nspinor 519 !write(std_out,*) "psichigather",isppol,ikpt,ibandc,& 520 !& real(paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,1,1))**2+& 521 !& imag(paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,1,1))**2 522 ! 523 !enddo 524 !enddo 525 !enddo 526 !enddo 527 528 call xmpi_barrier(spaceComm) 529 !if(mpi_enreg%me.eq.0) write(177,*) "end",psichi 530 !if(mpi_enreg%me.eq.1) write(178,*) "end",psichi 531 !if(mpi_enreg%me.eq.2) write(179,*) "end",psichi 532 call xmpi_barrier(spaceComm) 533 !========================================================================== 534 !********* WRITE psichi in file for reference 535 !========================================================================== 536 if(me.eq.0) then 537 call psichi_print(dtset,cryst_struc%nattyp,cryst_struc%ntypat,nkpt,my_nspinor,& 538 & nsppol,paw_dmft,pawtab,psps,t2g) 539 end if ! proc=me 540 !========================================================================== 541 !********************* Check normalization and occupations *************** 542 !========================================================================== 543 ABI_DATATYPE_ALLOCATE(xocc_check,(natom)) 544 ABI_DATATYPE_ALLOCATE(xnorm_check,(natom)) 545 call init_matlu(natom,my_nspinor,nsppol,paw_dmft%lpawu,xocc_check) 546 call init_matlu(natom,my_nspinor,nsppol,paw_dmft%lpawu,xnorm_check) 547 call psichi_check(dtset,cryst_struc%nattyp,nkpt,my_nspinor,& 548 & nsppol,cryst_struc%ntypat,paw_dmft,pawtab,psps,xocc_check,xnorm_check) 549 !========================================================================== 550 !*************** write checks ******************************************* 551 !========================================================================== 552 if(abs(dtset%pawprtvol)>=3) then 553 write(message,*) "normalization computed" 554 call wrtout(std_out, message,'COLL') 555 end if 556 557 ABI_DATATYPE_ALLOCATE(loc_occ_check,(natom)) 558 ABI_DATATYPE_ALLOCATE(loc_norm_check,(natom)) 559 call init_matlu(natom,my_nspinor,nsppol,paw_dmft%lpawu,loc_occ_check) 560 call init_matlu(natom,my_nspinor,nsppol,paw_dmft%lpawu,loc_norm_check) 561 call copy_matlu(xocc_check,loc_occ_check,natom) 562 call copy_matlu(xnorm_check,loc_norm_check,natom) 563 564 write(message,'(2a,i4)') ch10," == Check: Occupations and Norm from psichi are" 565 call wrtout(std_out, message,'COLL') 566 567 if(paw_dmft%dmftcheck>=1) then 568 ! print occupations 569 write(message,'(2a,i4)') ch10,' ------ Unsymetrised Occupation' 570 call wrtout(std_out, message,'COLL') 571 572 call print_matlu(xocc_check,natom,dtset%pawprtvol) 573 574 ! print norms 575 write(message,'(2a,i4)') ch10,' ------ Unsymetrised Norm' 576 call wrtout(std_out, message,'COLL') 577 578 call print_matlu(xnorm_check,natom,dtset%pawprtvol) 579 end if 580 581 !symetrise and print occupations 582 call sym_matlu(cryst_struc,loc_occ_check,pawang) 583 584 write(message,'(2a,i4)') ch10,' ------ Symetrised Occupation' 585 call wrtout(std_out, message,'COLL') 586 587 call print_matlu(loc_occ_check,natom,dtset%pawprtvol) 588 589 !symetrise and print norms 590 call sym_matlu(cryst_struc,loc_norm_check,pawang) 591 592 write(message,'(2a,i4)') ch10,' ------ Symetrised Norm' 593 call wrtout(std_out, message,'COLL') 594 595 call print_matlu(loc_norm_check,natom,dtset%pawprtvol) 596 597 !deallocations 598 do iatom=1,natom 599 lda_occup%matlu(iatom)%mat=loc_occ_check(iatom)%mat 600 end do 601 602 !Tests density matrix LDA+U and density matrix computed here. 603 if(paw_dmft%dmftcheck==2.or.(paw_dmft%dmftbandi==1)) then 604 ABI_DATATYPE_ALLOCATE(matlu_temp,(natom)) 605 call init_matlu(natom,paw_dmft%nspinor,paw_dmft%nsppol,paw_dmft%lpawu,matlu_temp) 606 do iatom=1,natom 607 if(paw_dmft%lpawu(iatom).ne.-1) then 608 ldim=2*paw_dmft%lpawu(iatom)+1 609 nsploop=max(paw_dmft%nsppol,paw_dmft%nspinor**2) 610 do idijeff=1,nsploop 611 ispinor=0 612 ispinor1=0 613 if(nsploop==2) then 614 isppol=spinor_idxs(1,idijeff) 615 ispinor=1 616 ispinor1=1 617 else if(nsploop==4) then 618 isppol=1 619 ispinor=spinor_idxs(1,idijeff) 620 ispinor1=spinor_idxs(2,idijeff) 621 else if(nsploop==1) then 622 isppol=1 623 ispinor=1 624 ispinor1=1 625 else 626 write(message,'(2a)') " BUG in datafordmft: nsploop should be equal to 2 or 4" 627 call wrtout(std_out,message,'COLL') 628 end if 629 do im1 = 1 , ldim 630 do im = 1 , ldim 631 if(my_nspinor==2) matlu_temp(iatom)%mat(im,im1,isppol,ispinor,ispinor1)=& 632 & cmplx(paw_ij(iatom)%noccmmp(1,im,im1,idijeff),paw_ij(iatom)%noccmmp(2,im,im1,idijeff),kind=dp) 633 if(my_nspinor==1) matlu_temp(iatom)%mat(im,im1,isppol,ispinor,ispinor1)=& 634 & cmplx(paw_ij(iatom)%noccmmp(1,im,im1,idijeff),zero,kind=dp) 635 end do 636 end do 637 end do 638 end if 639 end do 640 if(paw_dmft%dmftcheck==2) option=1 641 if(paw_dmft%dmftcheck<=1) option=0 642 call diff_matlu("LDA+U density matrix from INPUT wfk",& 643 & "Direct calculation of density matrix with psichi from DIAGONALIZED wfk",& 644 & matlu_temp,lda_occup%matlu,natom,option,tol3,ierrr) !tol1 tol2 tol3 645 if(ierrr==-1) then 646 write(message,'(10a)') ch10,& 647 & ' -> These two quantities should agree if three conditions are fulfilled',ch10,& 648 & ' - input wavefunctions come from the same Hamiltonien (e.g LDA/GGA)',ch10,& 649 & ' - dmatpuopt is equal to 1',ch10,& 650 & ' - all valence states are in the valence',ch10,& 651 & ' (for experts users: it is not compulsary that these conditions are fulfilled)' 652 call wrtout(std_out,message,'COLL') 653 end if 654 ! write(message,'(2a)') ch10,& 655 ! & ' ***** => Calculations of density matrices with projections and in LDA+U are coherent****' 656 ! call wrtout(std_out,message,'COLL') 657 658 call destroy_matlu(matlu_temp,natom) 659 ABI_DATATYPE_DEALLOCATE(matlu_temp) 660 else 661 write(message,'(2a)') ch10,& 662 & ' Warning: Consistency of density matrices computed from projection has not been checked: use dmftcheck>=2 ' 663 call wrtout(std_out,message,'COLL') 664 end if 665 666 call destroy_matlu(loc_norm_check,natom) 667 ABI_DATATYPE_DEALLOCATE(loc_norm_check) 668 call destroy_matlu(loc_occ_check,natom) 669 ABI_DATATYPE_DEALLOCATE(loc_occ_check) 670 671 call destroy_matlu(xocc_check,natom) 672 call destroy_matlu(xnorm_check,natom) 673 ABI_DATATYPE_DEALLOCATE(xocc_check) 674 ABI_DATATYPE_DEALLOCATE(xnorm_check) 675 676 if(present(nbandkss)) then 677 if(me.eq.0.and.nbandkss/=0) then 678 ! opt_renorm=1 ! if ucrpa==1, no need for individual orthonormalization 679 opt_renorm=3 680 if(dtset%ucrpa>=1) opt_renorm=2 681 call psichi_renormalization(cryst_struc,paw_dmft,pawang,opt=opt_renorm) 682 call psichi_print(dtset,cryst_struc%nattyp,cryst_struc%ntypat,nkpt,my_nspinor,& 683 & nsppol,paw_dmft,pawtab,psps,t2g) 684 end if ! proc=me 685 end if 686 687 !********************* 688 !deallocate(psichi) 689 !DEBUG 690 !write(std_out,*)' datafordmft : exit' 691 !stop 692 !ENDDEBUG 693 CONTAINS 694 !===========================================================
datafordmft/psichi_check [ Functions ]
[ Top ] [ datafordmft ] [ Functions ]
NAME
psichi_check
FUNCTION
Check psichi: compute occupations
INPUTS
dtset <type(dataset_type)>=all input variables for this dataset maxnproju = maximum number of projector for LDA+U species nattyp(ntypat)= # atoms of each type mband= number of bands nkpt=number of k points my_nspinor=number of spinorial components of the wavefunctions (on current proc) nsppol=1 for unpolarized, 2 for spin-polarized ntypat= number of species paw_dmft <type(paw_dmft)>=paw data for the self-consistency pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data psps <type(pseudopotential_type)>=variables related to pseudopotentials OUTPUTS: xocc_check(nsppol,my_nspinor,my_nspinor,natom,2*maxlpawu+1,2*maxlpawu+1): density matrix xnorm_check(nsppol,my_nspinor,my_nspinor,natom,2*maxlpawu+1,2*maxlpawu+1): matrix of norms
SIDE EFFECTS
check psichi: compute norm and occupations
PARENTS
datafordmft
CHILDREN
SOURCE
896 subroutine psichi_check(dtset,nattyp,nkpt,my_nspinor,& 897 & nsppol,ntypat,paw_dmft,pawtab,psps,xocc_check,xnorm_check) 898 899 use m_profiling_abi 900 901 use m_matlu, only: matlu_type,init_matlu,sym_matlu 902 903 !This section has been created automatically by the script Abilint (TD). 904 !Do not modify the following lines by hand. 905 #undef ABI_FUNC 906 #define ABI_FUNC 'psichi_check' 907 !End of the abilint section 908 909 implicit none 910 911 !Arguments ------------------------------------ 912 !scalars 913 integer,intent(in) :: nkpt,my_nspinor,nsppol,ntypat 914 !arrays 915 integer, intent(in) :: nattyp(ntypat) 916 type(dataset_type),intent(in) :: dtset 917 type(pseudopotential_type),intent(in) :: psps 918 type(paw_dmft_type), intent(in) :: paw_dmft 919 type(matlu_type), intent(inout) :: xocc_check(dtset%natom) !vz_i 920 type(matlu_type), intent(inout) :: xnorm_check(dtset%natom) !vz_i 921 type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw) 922 !Local variables ------------------------------------ 923 integer :: band_index,facnsppol,ibg,isppol,ikpt,iband,ibandc,ispinor,icat,itypat 924 integer :: iat,iatom,ilmn,lmn_size,m,m1,nband_k 925 complex(dpc) :: psichic,psichic1 926 real(dp) :: chinorm 927 928 ! ********************************************************************* 929 facnsppol=1 930 if(my_nspinor==1.and.nsppol==1) then 931 facnsppol=2 932 end if 933 934 ibg=0 935 band_index=0 936 do isppol=1,nsppol 937 do ikpt=1,nkpt 938 nband_k=dtset%nband(ikpt+(isppol-1)*nkpt) 939 ibandc=0 940 do iband=1,nband_k 941 if(paw_dmft%band_in(iband)) ibandc=ibandc+1 942 do ispinor=1,my_nspinor 943 icat=1 944 iat=0 945 do itypat=1,dtset%ntypat 946 lmn_size=pawtab(itypat)%lmn_size 947 do iatom=icat,icat+nattyp(itypat)-1 948 iat=iat+1 949 ! ------------ Select correlated atoms 950 if(paw_dmft%lpawu(iatom).ne.-1) then 951 chinorm=1.d0 952 do ilmn=1,lmn_size 953 ! ------------ Select l=lpawu. 954 if (psps%indlmn(1,ilmn,itypat)==paw_dmft%lpawu(iatom).and.& 955 & psps%indlmn(3,ilmn,itypat)==1) then 956 do ilmn1=1,lmn_size 957 ! ------------ Select l=lpawu and do not sum over projectors 958 ! (this is already done in paw_dmft%psichi) 959 if (psps%indlmn(1,ilmn1,itypat)==paw_dmft%lpawu(iatom).and.& 960 & psps%indlmn(3,ilmn1,itypat)==1) then 961 ! ------------ Check that the band is choosen (within nbandi and nbandf) 962 if(paw_dmft%band_in(iband)) then 963 m=psps%indlmn(2,ilmn,itypat)+psps%indlmn(1,ilmn,itypat)+1 964 m1=psps%indlmn(2,ilmn1,itypat)+psps%indlmn(1,ilmn,itypat)+1 965 if(psps%indlmn(3,ilmn,itypat)==1) then 966 do ispinor1=1,my_nspinor 967 psichic=paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m) 968 psichic1=paw_dmft%psichi(isppol,ikpt,ibandc,ispinor1,iat,m1) 969 ! ------------ Compute occupation matrix 970 xocc_check(iatom)%mat(m,m1,isppol,ispinor,ispinor1)=& 971 & xocc_check(iatom)%mat(m,m1,isppol,ispinor,ispinor1)& 972 ! & +conjg(psichic)*psichic1*dtset%wtk(ikpt)*facpara*occ(iband+band_index) 973 & +conjg(psichic1)*psichic*dtset%wtk(ikpt)*facpara*occ(iband+band_index)/facnsppol 974 ! ------------ Compute norm (should be equal to noccmmp 975 ! (dmatpuopt=1) if all bands are taken into account) 976 xnorm_check(iatom)%mat(m,m1,isppol,ispinor,ispinor1)=& 977 & xnorm_check(iatom)%mat(m,m1,isppol,ispinor,ispinor1)& 978 ! & +conjg(psichic)*psichic1*dtset%wtk(ikpt)*facpara 979 & +conjg(psichic1)*psichic*dtset%wtk(ikpt)*facpara 980 end do ! ispinor1 981 end if 982 end if ! paw_dmft%band_in 983 end if 984 end do !ilmn1 985 end if 986 end do !ilmn 987 end if ! lpawu.ne.-1 988 end do ! iatom 989 icat=icat+nattyp(itypat) 990 end do ! itypat 991 end do ! ispinor 992 end do !iband 993 band_index=band_index+nband_k 994 ibg=ibg+nband_k*my_nspinor 995 end do !ikpt 996 end do ! isppol 997 998 end subroutine psichi_check 999 !DBG_EXIT("COLL") 1000 1001 end subroutine datafordmft
datafordmft/psichi_print [ Functions ]
[ Top ] [ datafordmft ] [ Functions ]
NAME
psichi_print
FUNCTION
Print psichi for reference
INPUTS
dtset <type(dataset_type)>=all input variables for this dataset maxnproju = maximum number of projector for LDA+U species nattyp(ntypat)= # atoms of each type mband= number of bands nkpt=number of k points my_nspinor=number of spinorial components of the wavefunctions nsppol=1 for unpolarized, 2 for spin-polarized paw_dmft <type(paw_dmft)>=paw data for the self-consistency pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data psichi(2,nsppol,nkpt,mband,my_nspinor,dtset%natom,(2*paw_dmft%maxlpawu+1))) projections for DMFT psps <type(pseudopotential_type)>=variables related to pseudopotentials
SIDE EFFECTS
print psichi in forlb.ovlp
PARENTS
datafordmft
CHILDREN
SOURCE
728 subroutine psichi_print(dtset,nattyp,ntypat,nkpt,my_nspinor,& 729 &nsppol,paw_dmft,pawtab,psps,t2g) 730 731 use m_profiling_abi 732 use m_io_tools, only : open_file 733 734 !This section has been created automatically by the script Abilint (TD). 735 !Do not modify the following lines by hand. 736 #undef ABI_FUNC 737 #define ABI_FUNC 'psichi_print' 738 use interfaces_14_hidewrite 739 !End of the abilint section 740 741 implicit none 742 !Arguments ------------------------------------ 743 !scalars 744 integer,intent(in) :: nkpt,my_nspinor,nsppol,ntypat 745 !arrays 746 integer, intent(in) :: nattyp(ntypat) 747 type(dataset_type),intent(in) :: dtset 748 type(pseudopotential_type),intent(in) :: psps 749 logical t2g 750 type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw) 751 type(paw_dmft_type), intent(in) :: paw_dmft 752 !Local variables ------------------------------------ 753 integer :: ibg,isppol,ikpt,iband,ibandc,ispinor,icat,itypat,lmn_size 754 integer :: iat,iatom,jj1,ilmn,m1,nband_k,unt 755 integer :: m1_t2g,ll 756 real(dp) :: chinorm 757 character(len=500) :: msg 758 759 ! ********************************************************************* 760 ll=1 761 762 if (open_file('forlb.ovlp',msg,newunit=unt,form='formatted',status='unknown') /= 0) then 763 MSG_ERROR(msg) 764 end if 765 rewind(unt) 766 767 ! Header for calc_uCRPA.F90 768 if (COUNT(pawtab(:)%lpawu.NE.-1).EQ.1) then 769 do itypat=1,ntypat 770 if(t2g) then 771 if(pawtab(itypat)%lpawu.ne.-1) write(unt,*) "l= ",ll,itypat 772 else 773 if(pawtab(itypat)%lpawu.ne.-1) write(unt,*) "l= ",pawtab(itypat)%lpawu,itypat 774 end if 775 end do 776 else 777 write(unt,*) "More than one correlated species" 778 end if 779 780 write(unt,*) "Bands ",paw_dmft%dmftbandi,paw_dmft%dmftbandf 781 782 ibg=0 783 do isppol=1,nsppol 784 do ikpt=1,nkpt 785 ! rewind(1023) 786 write(unt,'(a6,2x,i6)') "ikpt =",ikpt 787 nband_k=dtset%nband(ikpt+(isppol-1)*nkpt) 788 ibandc=0 789 do iband=1,nband_k 790 if(paw_dmft%band_in(iband)) then 791 ibandc=ibandc+1 792 write(unt,'(a8,2x,i6)') " iband =",iband 793 end if 794 do ispinor=1,my_nspinor 795 icat=1 796 ! write(std_out,*) isppol,ikpt,iband,ispinor 797 iat=0 ! to declare 798 do itypat=1,dtset%ntypat 799 lmn_size=pawtab(itypat)%lmn_size 800 ! write(std_out,*) isppol,ikpt,iband,ispinor 801 do iatom=icat,icat+nattyp(itypat)-1 802 iat=iat+1 803 jj1=0 804 if(pawtab(dtset%typat(iatom))%lpawu.ne.-1) then 805 ! chinorm=(pawtab(itypat)%phiphjint(1)) 806 chinorm=1.d0 807 ! write(std_out,*) isppol,ikpt,iband,ispinor,iat 808 m1_t2g=0 809 do ilmn=1,lmn_size 810 ! write(std_out,*) ilmn 811 ! ------------ Select l=lpawu. --------------------------------------- 812 if (psps%indlmn(1,ilmn,itypat)==pawtab(dtset%typat(iatom))%lpawu.and.psps%indlmn(3,ilmn,itypat)==1) then 813 ! ------------ Check that the band is choosen (within nbandi and nbandf) 814 if(paw_dmft%band_in(iband)) then 815 jj1=jj1+1 816 if(jj1>(2*pawtab(dtset%typat(iatom))%lpawu+1)) then 817 write(message,'(a,a,i4,i5,i4)') ch10," Error 2 in datafordmft",jj1,pawtab(dtset%typat(iatom))%lpawu 818 call wrtout(std_out, message,'COLL') 819 MSG_ERROR("Aborting now") 820 end if ! jj1 821 ! if(jj1>pawtab(dtset%typat(iatom))%nproju*(2*lpawu+1)) then 822 ! write(message,'(a,a,a,a)') ch10,& 823 ! & 'BUG: jj1 is not correct in datafordmft psichi_print',ch10,& 824 ! & 'Action: CONTACT Abinit group' 825 ! call wrtout(std_out, message,'COLL') 826 ! call leave_new('COLL') 827 ! end if ! jj1 828 m1=psps%indlmn(2,ilmn,itypat)+psps%indlmn(1,ilmn,itypat)+1 829 ! ----- Print only when the sum over projectors is done 830 ! write(std_out,*) ilmn,m1 831 if(t2g) then 832 if(m1==1.or.m1==2.or.m1==4) then 833 m1_t2g=m1_t2g+1 834 write(unt,'(3i6,3x,2f23.15)') isppol, iat, m1,& 835 & real(paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1_t2g))/chinorm,& 836 & aimag(paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1_t2g))/chinorm 837 end if 838 else 839 write(unt,'(3i6,3x,2f23.15)') isppol, iat, m1,& 840 & real(paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1))/chinorm,& 841 & aimag(paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1))/chinorm 842 end if 843 ! if natom=1 then jj1 maximal value should be 2*lpawu+1 844 end if ! paw_dmft%band_in 845 end if 846 end do !ilmn 847 end if ! lpawu.ne.-1 848 end do ! iatom 849 icat=icat+nattyp(itypat) 850 end do ! itypat 851 end do ! ispinor 852 end do !iband 853 ibg=ibg+nband_k*my_nspinor 854 end do !ikpt 855 end do ! isppol 856 ! write(unt,*) "Fermi level (in Ryd)=" 857 ! write(unt,*) fermie*two 858 close(unt) 859 end subroutine psichi_print