TABLE OF CONTENTS
ABINIT/setrhoijpbe0 [ Functions ]
NAME
setrhoijpbe0
FUNCTION
PAW local exact exchange only: Impose value of rhoij for f electrons using an auxiliairy file
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (FJ,MT,MD) 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
dtset <type(dataset_type)>=all input variables for this dataset initialized= if 0, the initialization of the gstate run is not yet finished istep=index of the number of steps in the routine scfcv istep_mix=index of the number of steps for the SCF mixing (can be <istep) mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc comm_atom=--optional-- MPI communicator over atoms mpi_comm_read=MPI communicator containing all the processes reading the PBE0 file my_natom=number of atoms treated by current processor natom=number of atoms in cell ntypat=number of types of atoms in unit cell pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data typat(natom)=type integer for each atom in cell
SIDE EFFECTS
istep_mix=index of the number of steps for the SCF mixing (can be <istep) pawrhoij(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
NOTES
Only valid for f electrons !!!
PARENTS
scfcv
CHILDREN
free_my_atmtab,get_my_atmtab,pawio_print_ij,wrtout,xmpi_sum
SOURCE
46 #if defined HAVE_CONFIG_H 47 #include "config.h" 48 #endif 49 50 #include "abi_common.h" 51 52 subroutine setrhoijpbe0(dtset,initialized,istep,istep_mix,& 53 & mpi_comm_read,my_natom,natom,ntypat,pawrhoij,pawtab,typat,& 54 & mpi_atmtab,comm_atom) ! optional arguments (parallelism) 55 56 use defs_basis 57 use defs_abitypes 58 use m_profiling_abi 59 use m_errors 60 use m_xmpi 61 62 use m_io_tools, only : open_file 63 use m_pawtab, only : pawtab_type 64 use m_pawrhoij, only : pawrhoij_type 65 use m_paw_io, only : pawio_print_ij 66 use m_paral_atom, only : get_my_atmtab, free_my_atmtab 67 68 !This section has been created automatically by the script Abilint (TD). 69 !Do not modify the following lines by hand. 70 #undef ABI_FUNC 71 #define ABI_FUNC 'setrhoijpbe0' 72 use interfaces_14_hidewrite 73 !End of the abilint section 74 75 implicit none 76 77 !Arguments --------------------------------------------- 78 !scalars 79 integer,intent(in) :: initialized,istep,mpi_comm_read,my_natom,natom,ntypat 80 integer,intent(inout) :: istep_mix 81 integer,optional,intent(in) :: comm_atom 82 type(dataset_type),intent(in) :: dtset 83 !arrays 84 integer,intent(in) :: typat(natom) 85 integer,optional,target,intent(in) :: mpi_atmtab(:) 86 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom) 87 type(pawtab_type),intent(in) :: pawtab(ntypat) 88 89 !Local variables --------------------------------------- 90 !scalars 91 integer,parameter :: ll=3 92 integer :: iatom,iatom_tot,ierr,ii,ios,iread,irhoij,ispden,itypat,jj,klmn,my_comm_atom,my_rank,nselect 93 integer :: nstep1,nstep1_abs,rhoijshft,rhoijsz 94 logical :: my_atmtab_allocated,paral_atom,test0 95 character(len=9),parameter :: filnam='rhoijpbe0' 96 character(len=9),parameter :: dspin(6)=(/"up ","down ","up-up ","down-down","Re[up-dn]","Im[up-dn]"/) 97 character(len=500) :: strg, message 98 !arrays 99 integer, allocatable :: nspden_tmp(:) 100 integer,pointer :: my_atmtab(:) 101 real(dp),allocatable :: rhoijtmp(:,:),rhoijtmp1(:,:),rhoijtmp2(:,:,:,:) 102 103 ! ********************************************************************* 104 105 DBG_ENTER("COLL") 106 107 !Test existence of file and open it 108 inquire(file=filnam,iostat=ios,exist=test0) 109 if(.not.test0) return 110 111 !Look for parallelisation over atomic sites 112 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 113 114 !Test if exact-exch. is on f electrons 115 test0=.false. 116 do itypat=1,ntypat 117 if (pawtab(itypat)%useexexch>0.and.pawtab(itypat)%lexexch/=ll) test0=.true. 118 end do 119 if (test0) then 120 write(message, '(3a,i1,a)' ) & 121 & ' Local exact exchange: occ. matrix can only be imposed for l=',ll,' !' 122 MSG_ERROR(message) 123 end if 124 125 !============================================================ 126 !===== First case: no parallelisation over atomic sites ===== 127 !============================================================ 128 129 if (.not.paral_atom) then 130 131 ! Open file 132 if (open_file(filnam,message,unit=77,form='formatted') /= 0) then 133 MSG_ERROR(message) 134 end if 135 136 ! Read step number and eventually exit 137 nstep1=0;test0=.false. 138 do while (.not.test0) 139 read(77,'(A)') strg 140 test0=(strg(1:1)/="#") 141 if (test0) read(unit=strg,fmt=*) nstep1 142 end do 143 nstep1_abs=abs(nstep1) 144 if (nstep1_abs==0.or.istep>nstep1_abs.or.(nstep1>0.and.initialized/=0)) then 145 close(77) 146 ! Reinitalize mixing when rhoij is allowed to change; for experimental purpose... 147 if (dtset%userib==1234.and.istep==1+nstep1_abs.and.(nstep1<0.or.initialized==0)) istep_mix=1 148 return 149 end if 150 151 ! Loop on atoms 152 do iatom=1,natom 153 itypat=typat(iatom) 154 if (pawtab(itypat)%useexexch>0) then 155 156 ! Set sizes depending on ll 157 rhoijsz=4*ll+2 158 rhoijshft=2*ll*ll 159 160 ! Uncompress rhoij 161 ABI_ALLOCATE(rhoijtmp,(pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden)) 162 do ispden=1,pawrhoij(iatom)%nspden 163 rhoijtmp=zero 164 do irhoij=1,pawrhoij(iatom)%nrhoijsel 165 klmn=pawrhoij(iatom)%rhoijselect(irhoij) 166 rhoijtmp(klmn,ispden)=pawrhoij(iatom)%rhoijp(irhoij,ispden) 167 end do 168 end do 169 ! Read rhoij from file 170 ABI_ALLOCATE(rhoijtmp1,(rhoijsz,rhoijsz)) 171 do ispden=1,pawrhoij(iatom)%nspden 172 do ii=1,rhoijsz 173 test0=.false. 174 do while (.not.test0) 175 read(77,'(A)') strg 176 test0=(strg(1:1)/="#") 177 if (test0) read(unit=strg,fmt=*) (rhoijtmp1(ii,jj), jj=1,rhoijsz) 178 end do 179 end do 180 181 ! Impose rhoij 182 do jj=1,rhoijsz 183 do ii=1,jj 184 rhoijtmp((jj+rhoijshft)*((jj+rhoijshft)-1)/2+ii+rhoijshft,ispden)=rhoijtmp1(ii,jj) 185 end do 186 end do 187 188 end do 189 ABI_DEALLOCATE(rhoijtmp1) 190 191 ! Compress rhoij 192 nselect=0 193 do klmn=1,pawrhoij(iatom)%lmn2_size 194 if (any(abs(rhoijtmp(klmn,:))>tol10)) then 195 nselect=nselect+1 196 do ispden=1,pawrhoij(iatom)%nspden 197 pawrhoij(iatom)%rhoijp(nselect,ispden)=rhoijtmp(klmn,ispden) 198 end do 199 pawrhoij(iatom)%rhoijselect(nselect)=klmn 200 end if 201 end do 202 pawrhoij(iatom)%nrhoijsel=nselect 203 ABI_DEALLOCATE(rhoijtmp) 204 205 ! Print new rhoij 206 do ispden=1,pawrhoij(iatom)%nspden 207 write(message,'(2a,i3,a)') ch10,'== Atom ',iatom,& 208 & ' == Imposed occupation matrix' 209 if (pawrhoij(iatom)%nspden==1) write(message,fmt='(2a)') trim(message)," for spin up ==" 210 if (pawrhoij(iatom)%nspden==2) write(message,fmt='(2a,i3,a)')trim(message)," for spin ",ispden," ==" 211 if (pawrhoij(iatom)%nspden==4) write(message,fmt='(4a)') trim(message)," for component ", & 212 & trim(dspin(ispden+2*(pawrhoij(iatom)%nspden/4)))," ==" 213 call wrtout(std_out,message,'COLL') 214 call pawio_print_ij(std_out,pawrhoij(iatom)%rhoijp(:,ispden),pawrhoij(iatom)%nrhoijsel,& 215 & pawrhoij(iatom)%cplex,pawrhoij(iatom)%lmn_size,ll,& 216 & pawtab(itypat)%indlmn(1,1:pawtab(itypat)%lmn_size),& 217 & 1,-1,pawrhoij(iatom)%rhoijselect(:),-1.d0,1,mode_paral='COLL') 218 end do 219 220 ! End loop on atoms 221 end if 222 end do 223 224 ! Close file 225 close (77) 226 227 else 228 229 ! ============================================================ 230 ! ====== 2nd case: no parallelisation over atomic sites ===== 231 ! ============================================================ 232 233 my_rank=xmpi_comm_rank(mpi_comm_read) 234 235 ! Read step number and eventually exit 236 iread=0 237 if (my_rank==0) then 238 if (open_file(filnam,message,unit=77,form='formatted') /=0 ) then 239 MSG_ERROR(message) 240 end if 241 nstep1=0;test0=.false. 242 do while (.not.test0) 243 read(77,'(A)') strg 244 test0=(strg(1:1)/="#") 245 if (test0) read(unit=strg,fmt=*) nstep1 246 end do 247 nstep1_abs=abs(nstep1) 248 if (nstep1_abs==0.or.istep>nstep1_abs.or.(nstep1>0.and.initialized/=0)) then 249 close(77) 250 ! Reinitalize mixing when rhoij is allowed to change; for experimental purpose... 251 if (dtset%userib==1234.and.istep==1+nstep1_abs.and.(nstep1<0.or.initialized==0)) istep_mix=1 252 iread=1 253 end if 254 end if 255 call xmpi_sum(iread,mpi_comm_read,ierr) 256 if (iread/=0) return 257 258 ! Set up parallelism over atoms 259 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 260 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 261 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom) 262 263 ! Store number of component for rhoij 264 ABI_ALLOCATE(nspden_tmp,(natom)) 265 nspden_tmp(:)=zero 266 do iatom=1,my_natom 267 iatom_tot=my_atmtab(iatom) 268 nspden_tmp(iatom_tot)=pawrhoij(iatom)%nspden 269 end do 270 call xmpi_sum(nspden_tmp,mpi_comm_read,ierr) 271 272 ! To be improve if too much memory 273 ABI_ALLOCATE(rhoijtmp2,(natom,4,rhoijsz,rhoijsz)) 274 rhoijtmp2=zero 275 276 ! Read rhoij from file 277 if (my_rank==0) then 278 do iatom=1,natom 279 itypat=typat(iatom) 280 if (pawtab(itypat)%useexexch>0) then 281 rhoijsz=4*ll+2 282 do ispden=1,nspden_tmp(iatom) 283 do ii=1,rhoijsz 284 test0=.false. 285 do while (.not.test0) 286 read(77,'(A)') strg 287 test0=(strg(1:1)/="#") 288 if (test0) read(unit=strg,fmt=*) (rhoijtmp2(iatom,ispden,ii,jj),jj=1,rhoijsz) 289 end do 290 end do 291 end do 292 end if 293 end do 294 end if 295 call xmpi_sum(rhoijtmp2,mpi_comm_read,ierr) 296 297 ! Now, distribute rhoij 298 do iatom=1,my_natom 299 iatom_tot=my_atmtab(iatom) 300 itypat=pawrhoij(iatom)%itypat 301 302 if (pawtab(itypat)%useexexch>0) then 303 304 ! Set sizes depending on ll 305 rhoijsz=4*ll+2 306 rhoijshft=2*ll*ll 307 308 ! Uncompress rhoij 309 ABI_ALLOCATE(rhoijtmp,(pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden)) 310 do ispden=1,pawrhoij(iatom)%nspden 311 rhoijtmp=zero 312 do irhoij=1,pawrhoij(iatom)%nrhoijsel 313 klmn=pawrhoij(iatom)%rhoijselect(irhoij) 314 rhoijtmp(klmn,ispden)=pawrhoij(iatom)%rhoijp(irhoij,ispden) 315 end do 316 317 ! Impose rhoij 318 do jj=1,rhoijsz 319 do ii=1,jj 320 rhoijtmp((jj+rhoijshft)*((jj+rhoijshft)-1)/2+ii+rhoijshft,ispden)=rhoijtmp2(iatom_tot,ispden,ii,jj) 321 end do 322 end do 323 324 end do 325 326 ! Compress rhoij 327 nselect=0 328 do klmn=1,pawrhoij(iatom)%lmn2_size 329 if (any(abs(rhoijtmp(klmn,:))>tol10)) then 330 nselect=nselect+1 331 do ispden=1,pawrhoij(iatom)%nspden 332 pawrhoij(iatom)%rhoijp(nselect,ispden)=rhoijtmp(klmn,ispden) 333 end do 334 pawrhoij(iatom)%rhoijselect(nselect)=klmn 335 end if 336 end do 337 pawrhoij(iatom)%nrhoijsel=nselect 338 ABI_DEALLOCATE(rhoijtmp) 339 340 end if ! useexexch>0 341 342 ! Print new rhoij 343 do ispden=1,pawrhoij(iatom)%nspden 344 write(message,'(2a,i3,a)') ch10,'== Atom ',iatom,' == Imposed occupation matrix' 345 if (pawrhoij(iatom)%nspden==1) write(message,fmt='(2a)') trim(message)," for spin up ==" 346 if (pawrhoij(iatom)%nspden==2) write(message,fmt='(2a,i3,a)')trim(message)," for spin ",ispden," ==" 347 if (pawrhoij(iatom)%nspden==4) write(message,fmt='(4a)') trim(message)," for component ", & 348 & trim(dspin(ispden+2*(pawrhoij(iatom)%nspden/4)))," ==" 349 call wrtout(std_out,message,'PERS') 350 call pawio_print_ij(std_out,pawrhoij(iatom)%rhoijp(:,ispden),pawrhoij(iatom)%nrhoijsel,& 351 & pawrhoij(iatom)%cplex,pawrhoij(iatom)%lmn_size,ll,& 352 & pawtab(itypat)%indlmn(1,1:pawtab(itypat)%lmn_size),& 353 & 1,-1,pawrhoij(iatom)%rhoijselect(:),-1.d0,1,mode_paral='PERS') 354 end do 355 356 ! end loop on atoms 357 end do 358 359 ABI_DEALLOCATE(nspden_tmp) 360 ABI_DEALLOCATE(rhoijtmp2) 361 362 ! Destroy atom table used for parallelism 363 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 364 365 ! ============================================================ 366 end if ! paral_atom 367 368 DBG_EXIT("COLL") 369 370 end subroutine setrhoijpbe0