TABLE OF CONTENTS
ABINIT/posratecore [ Functions ]
NAME
posratecore
FUNCTION
Calculate the annihilataion rate of a given core state
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (JW,GJ,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
dtset <type(dataset_type)>=all input variables for this dataset | nspden=number of spin-density components | ntypat=number of atom types | paral_kgb=flag controlling (k,g,bands) parallelization | pawxcdev=Choice of XC development (0=no dev. (use of angular mesh) ; 1 or 2=dev. on moments) | usepaw=flag for PAW iatom= index of the current atom in posdoppler mesh_sizej= size of the radial mesh for the current atom in posdoppler mpi_enreg= information about MPI parallelization my_natom=number of atoms treated by current processor option= if 1, use gamma if 2, use IPM (gamma=1) pawang <type(pawang)>=paw angular mesh and related data pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= -PAW only- atomic occupancies pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
OUTPUT
rate= annihilation rate of a given core state needed for state dependent scheme for doppler broadening
SIDE EFFECTS
PARENTS
posdoppler
CHILDREN
gammapositron,mkdenpos,nderiv_gen,pawdensities,pawxcsum,simp_gen xmpi_sum
SOURCE
48 #if defined HAVE_CONFIG_H 49 #include "config.h" 50 #endif 51 52 #include "abi_common.h" 53 54 subroutine posratecore(dtset,electronpositron,iatom,my_natom,mesh_sizej,mpi_enreg,& 55 & option,pawang,pawrad,pawrhoij,pawrhoij_ep,& 56 & pawtab,rate,rhocorej) 57 58 use defs_basis 59 use defs_abitypes 60 use m_profiling_abi 61 use m_errors 62 use m_xmpi 63 64 use m_electronpositron 65 66 use m_pawang, only : pawang_type 67 use m_pawrad, only : pawrad_type, simp_gen, nderiv_gen 68 use m_pawtab, only : pawtab_type 69 use m_pawrhoij, only : pawrhoij_type 70 use m_pawxc, only: pawxcsum 71 72 !This section has been created automatically by the script Abilint (TD). 73 !Do not modify the following lines by hand. 74 #undef ABI_FUNC 75 #define ABI_FUNC 'posratecore' 76 use interfaces_41_xc_lowlevel 77 use interfaces_56_xc 78 use interfaces_65_paw 79 !End of the abilint section 80 81 implicit none 82 83 !Arguments ------------------------------------ 84 !scalars 85 integer,intent(in) :: iatom,my_natom,option,mesh_sizej 86 real(dp),intent(out) :: rate 87 type(dataset_type), intent(in) :: dtset 88 type(electronpositron_type),pointer :: electronpositron 89 type(MPI_type),intent(in) :: mpi_enreg 90 type(pawang_type), intent(in) :: pawang 91 !arrays 92 real(dp),intent(in) :: rhocorej(mesh_sizej) 93 type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*dtset%usepaw) 94 type(pawrhoij_type),intent(in) :: pawrhoij(my_natom*dtset%usepaw) 95 type(pawrhoij_type),intent(in),target :: pawrhoij_ep(my_natom*dtset%usepaw) 96 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw) 97 98 !Local variables------------------------------- 99 !scalars 100 integer :: cplex,ierr,igamma,ii,ilm,ipt,ir 101 integer :: itypat,iwarn,iwarnj,iwarnp,lm_size,lmn2_size,mesh_size 102 integer :: ngr,ngrad,nspden_ep,opt_dens 103 logical,parameter :: include_nhat_in_gamma=.false. 104 real(dp),parameter :: delta=1.d-4 105 real(dp) :: fact,fact2,intg 106 real(dp) :: mpibuf,rdum,sqfpi 107 character(len=500) :: msg 108 !arrays 109 logical,allocatable :: lmselect(:),lmselect_ep(:),lmselect_dum(:) 110 real(dp),parameter :: qphon(3)=(/zero,zero,zero/),lsign(2)=(/one,-one/) 111 real(dp),allocatable :: d1gam(:,:),d2gam(:,:),ff(:),gam_(:,:,:),gamma(:,:),gammam(:,:),gg(:,:) 112 real(dp),allocatable :: grhocore2(:),grhocor2_(:),grhoe2(:),grho2_(:) 113 real(dp),allocatable :: nhat1(:,:,:),nhat1_ep(:,:,:) 114 real(dp),allocatable :: rho_(:),rho_ep_(:),rho1(:,:,:),rho1_ep(:,:,:) 115 real(dp),allocatable :: rhoarr1(:),rhoarr1_ep(:),rhoarr2(:) 116 real(dp),allocatable :: rhocore(:),rhocor_(:) 117 real(dp),allocatable :: rhosph(:),rhosph_ep(:),rhotot(:,:),rhotot_ep(:,:) 118 real(dp),allocatable :: trho1(:,:,:),trho1_ep(:,:,:) 119 real(dp),allocatable :: v1sum(:,:),v2sum(:,:,:) 120 type(pawrhoij_type),pointer :: pawrhoij_ep_(:) 121 122 ! ************************************************************************* 123 124 DBG_ENTER("COLL") 125 126 !Tests for developers 127 if (.not.associated(electronpositron)) then 128 msg='electronpositron variable must be associated!' 129 MSG_BUG(msg) 130 end if 131 !Constants 132 fact=0.0 133 cplex=1;nspden_ep=1 134 ngrad=1;if (electronpositron%ixcpositron==3.or.electronpositron%ixcpositron==31) ngrad=2 135 iwarn=0;iwarnj=0;iwarnp=1 136 sqfpi=sqrt(four_pi) 137 138 !Compatibility tests 139 if (electronpositron%particle==EP_NOTHING) then 140 msg='Not valid for electronpositron%particle=NOTHING!' 141 MSG_BUG(msg) 142 end if 143 144 if (dtset%usepaw==1) then 145 if(dtset%pawxcdev==0.and.ngrad==2) then 146 msg='GGA is not implemented for pawxcdev=0 (use dtset%pawxcdev/=0)!' 147 MSG_BUG(msg) 148 end if 149 end if 150 151 !Select type(s) of enhancement factor 152 if (electronpositron%ixcpositron==-1) igamma=0 153 if (electronpositron%ixcpositron== 2) igamma=4 154 if (electronpositron%ixcpositron==11.or.electronpositron%ixcpositron==31) igamma=3 155 if (electronpositron%ixcpositron==1.or.electronpositron%ixcpositron==3) igamma=2 156 if (option==2) igamma=0 157 158 pawrhoij_ep_ => pawrhoij_ep 159 160 rate=zero 161 162 itypat=pawrhoij(iatom)%itypat 163 lmn2_size=pawtab(itypat)%lmn2_size 164 mesh_size=pawtab(itypat)%mesh_size 165 lm_size=pawtab(itypat)%lcut_size**2 166 cplex=1 167 ngr=0;if (ngrad==2) ngr=mesh_size 168 169 !Allocations of "on-site" densities 170 ABI_ALLOCATE(rho1 ,(cplex*mesh_size,lm_size,nspden_ep)) 171 ABI_ALLOCATE(trho1,(cplex*mesh_size,lm_size,nspden_ep)) 172 ABI_ALLOCATE(rho1_ep ,(cplex*mesh_size,lm_size,nspden_ep)) 173 ABI_ALLOCATE(trho1_ep,(cplex*mesh_size,lm_size,nspden_ep)) 174 ABI_ALLOCATE(lmselect,(lm_size)) 175 ABI_ALLOCATE(lmselect_ep,(lm_size)) 176 ABI_ALLOCATE(lmselect_dum,(lm_size)) 177 if (include_nhat_in_gamma) then 178 ABI_ALLOCATE(nhat1,(cplex*mesh_size,lm_size,nspden_ep)) 179 ABI_ALLOCATE(nhat1_ep,(cplex*mesh_size,lm_size,nspden_ep)) 180 else 181 ABI_ALLOCATE(nhat1,(0,0,0)) 182 ABI_ALLOCATE(nhat1_ep,(0,0,0)) 183 end if 184 185 !Compute "on-site" densities (n1, ntild1, nhat1) for electron and positron ===== 186 lmselect(:)=.true. 187 opt_dens=1;if (include_nhat_in_gamma) opt_dens=0 188 call pawdensities(rdum,cplex,iatom,lmselect,lmselect_dum,lm_size,nhat1,nspden_ep,1,& 189 & 0,opt_dens,-1,0,pawang,0,pawrad(itypat),pawrhoij(iatom),& 190 & pawtab(itypat),rho1,trho1) 191 lmselect_ep(:)=.true. 192 call pawdensities(rdum,cplex,iatom,lmselect_ep,lmselect_dum,lm_size,nhat1_ep,nspden_ep,1,& 193 & 0,opt_dens,-1,0,pawang,0,pawrad(itypat),pawrhoij_ep_(iatom),& 194 & pawtab(itypat),rho1_ep,trho1_ep) 195 !Compute contribution to annihilation rate 196 197 ABI_ALLOCATE(rhocore,(mesh_size)) 198 199 !First formalism: use densities on r,theta,phi 200 if (dtset%pawxcdev==0) then 201 202 ABI_ALLOCATE(gamma,(mesh_size,2)) 203 ABI_ALLOCATE(rhoarr1,(mesh_size)) 204 ABI_ALLOCATE(rhoarr1_ep,(mesh_size)) 205 206 ! Loop on the angular part 207 do ipt=1,pawang%angl_size 208 ! Build densities 209 rhoarr1=zero;rhoarr1_ep=zero;rhocore=zero 210 do ilm=1,lm_size 211 if (lmselect(ilm)) rhoarr1(:)=rhoarr1(:)+rho1(:,ilm,1)*pawang%ylmr(ilm,ipt) 212 end do 213 do ilm=1,lm_size 214 if (lmselect_ep(ilm)) rhoarr1_ep(:)=rhoarr1_ep(:)+rho1_ep(:,ilm,1)*pawang%ylmr(ilm,ipt) 215 end do 216 rhocore(:)=pawtab(itypat)%coredens(:) 217 ! Make the densities positive 218 if (electronpositron%particle==EP_ELECTRON) then 219 call mkdenpos(iwarnp,mesh_size,1,1,rhoarr1 ,dtset%xc_denpos) 220 call mkdenpos(iwarn ,mesh_size,1,1,rhoarr1_ep,dtset%xc_denpos) 221 else if (electronpositron%particle==EP_POSITRON) then 222 call mkdenpos(iwarn ,mesh_size,1,1,rhoarr1 ,dtset%xc_denpos) 223 call mkdenpos(iwarnp,mesh_size,1,1,rhoarr1_ep,dtset%xc_denpos) 224 end if 225 ! Compute Gamma 226 ABI_ALLOCATE(grhoe2,(ngr)) 227 ABI_ALLOCATE(grhocore2,(ngr)) 228 if (electronpositron%particle==EP_ELECTRON) then 229 call gammapositron(gamma,grhocore2,grhoe2,igamma,ngr,mesh_size,& 230 & rhocore,rhoarr1_ep,rhoarr1,1) 231 else if (electronpositron%particle==EP_POSITRON) then 232 call gammapositron(gamma,grhocore2,grhoe2,igamma,ngr,mesh_size,& 233 & rhocore,rhoarr1,rhoarr1_ep,1) 234 end if 235 ABI_DEALLOCATE(grhoe2) 236 ABI_DEALLOCATE(grhocore2) 237 ! Compute contribution to annihilation rates 238 239 ABI_ALLOCATE(ff,(mesh_size)) 240 ff(1:mesh_size)=rhoarr1_ep(1:mesh_size)*rhocorej(1:mesh_size) & 241 & *gamma(1:mesh_size,1)*pawrad(itypat)%rad(1:mesh_size)**2 242 call simp_gen(intg,ff,pawrad(itypat)) 243 intg=intg*pawang%angwgth(ipt)*four_pi 244 rate =rate +intg 245 ABI_DEALLOCATE(ff) 246 end do ! ipt 247 ABI_DEALLOCATE(gamma) 248 ABI_DEALLOCATE(rhoarr1) 249 ABI_DEALLOCATE(rhoarr1_ep) 250 251 !Second formalism: use (l,m) moments for densities 252 else if (dtset%pawxcdev/=0) then 253 254 ! Build densities 255 ABI_ALLOCATE(gammam,(mesh_size,lm_size)) 256 ABI_ALLOCATE(rhotot,(mesh_size,lm_size)) 257 ABI_ALLOCATE(rhotot_ep,(mesh_size,lm_size)) 258 ABI_ALLOCATE(rhosph,(mesh_size)) 259 ABI_ALLOCATE(rhosph_ep,(mesh_size)) 260 261 rhotot (:,:)=rho1 (:,:,1) 262 rhotot_ep(:,:)=rho1_ep(:,:,1) 263 rhocore(:)=pawtab(itypat)%coredens(:) 264 rhosph (:)=rhotot (:,1)/sqfpi 265 rhosph_ep(:)=rhotot_ep(:,1)/sqfpi 266 ! Make spherical densities positive 267 if (electronpositron%particle==EP_ELECTRON) then 268 call mkdenpos(iwarnp,mesh_size,1,1,rhosph ,dtset%xc_denpos) 269 call mkdenpos(iwarn ,mesh_size,1,1,rhosph_ep,dtset%xc_denpos) 270 else if (electronpositron%particle==EP_POSITRON) then 271 call mkdenpos(iwarn ,mesh_size,1,1,rhosph ,dtset%xc_denpos) 272 call mkdenpos(iwarnp,mesh_size,1,1,rhosph_ep,dtset%xc_denpos) 273 end if 274 275 ! Need gradients of electronic densities for GGA 276 ABI_ALLOCATE(grhoe2,(ngr)) 277 ABI_ALLOCATE(grhocore2,(ngr)) 278 if (ngr>0) then 279 if (electronpositron%particle==EP_ELECTRON) then 280 call nderiv_gen(grhoe2,rhosph_ep,pawrad(itypat)) 281 else if (electronpositron%particle==EP_POSITRON) then 282 call nderiv_gen(grhoe2,rhosph,pawrad(itypat)) 283 end if 284 grhoe2(:)=grhoe2(:)**2 285 call nderiv_gen(grhocore2,rhocore,pawrad(itypat)) 286 grhocore2(:)=grhocore2(:)**2 287 end if 288 ! Compute Gamma for (rho-,rho+), 289 ! (rho- +drho-,rho+), (rho- -drho-,rho+), 290 ! (rho-,rho+ +drho+), (rho-,rho+ -drho+), 291 ! (rho- +drho-,rho+ +drho+), (rho- -drho-,rho+ -drho+) 292 ! Do a seven steps loop 293 ABI_ALLOCATE(gam_,(mesh_size,2,7)) 294 ABI_ALLOCATE(rho_,(mesh_size)) 295 ABI_ALLOCATE(rho_ep_,(mesh_size)) 296 ABI_ALLOCATE(rhocor_,(mesh_size)) 297 ABI_ALLOCATE(grho2_,(ngr)) 298 ABI_ALLOCATE(grhocor2_,(ngr)) 299 300 do ii=1,7 301 ! Apply delta to get perturbed densities 302 rho_(:)=rhosph(:);rho_ep_(:)=rhosph_ep(:);rhocor_(:)=rhocore(:) 303 if (ngr>0) grho2_(:)=grhoe2(:) 304 if (ngr>0) grhocor2_(:)=grhocore2(:) 305 if (ii==2.or.ii==4.or.ii==6) fact=(one+delta) 306 if (ii==3.or.ii==5.or.ii==7) fact=(one-delta) 307 fact2=fact**2 308 if (ii==2.or.ii==3.or.ii==6.or.ii==7) then 309 rho_(:)=fact*rho_(:) 310 if (electronpositron%particle==EP_POSITRON) then 311 if (ngr>0) grho2_(:)=fact2*grho2_(:) 312 rhocor_(:)=fact*rhocor_(:) 313 if (ngr>0) grhocor2_(:)=fact2*grhocor2_(:) 314 end if 315 end if 316 317 if (ii==4.or.ii==5.or.ii==6.or.ii==7) then 318 rho_ep_(:)=fact*rho_ep_(:) 319 if (electronpositron%particle==EP_ELECTRON) then 320 if (ngr>0) grho2_(:)=fact2*grho2_(:) 321 rhocor_(:)=fact*rhocor_(:) 322 if (ngr>0) grhocor2_(:)=fact2*grhocor2_(:) 323 end if 324 end if 325 ! Compute gamma for these perturbed densities 326 if (electronpositron%particle==EP_ELECTRON) then 327 call gammapositron(gam_(:,:,ii),grhocor2_,grho2_,igamma,ngr,mesh_size,rhocor_,rho_ep_,rho_,1) 328 else if (electronpositron%particle==EP_POSITRON) then 329 call gammapositron(gam_(:,:,ii),grhocor2_,grho2_,igamma,ngr,mesh_size,rhocor_,rho_,rho_ep_,1) 330 end if 331 332 end do ! end loop ii=1,7 333 334 ABI_DEALLOCATE(rhocor_) 335 ABI_DEALLOCATE(grho2_) 336 ABI_DEALLOCATE(grhocor2_) 337 ABI_DEALLOCATE(grhoe2) 338 ABI_DEALLOCATE(grhocore2) 339 rho_ (:)=rhosph (:);if (electronpositron%particle==EP_POSITRON) rho_ (:)=rho_ (:)+rhocore(:) 340 rho_ep_(:)=rhosph_ep(:);if (electronpositron%particle==EP_ELECTRON) rho_ep_(:)=rho_ep_(:)+rhocore(:) 341 ! Compute numerical first and second derivatives of Gamma 342 ! d1gam(1) = dgam/drho+ (particle=ELECTRON), dgam/drho- (particle=POSITRON) 343 ! d1gam(2) = dgam/drho- (particle=ELECTRON), dgam/drho+ (particle=POSITRON) 344 ABI_ALLOCATE(d1gam,(mesh_size,2)) 345 d1gam(:,:)=zero 346 do ir=1,mesh_size 347 if (rho_ (ir)>tol14) d1gam(ir,1)=(gam_(ir,1,2)-gam_(ir,1,3))*half/(delta*rho_ (ir)) 348 if (rho_ep_ (ir)>tol14) d1gam(ir,2)=(gam_(ir,1,4)-gam_(ir,1,5))*half/(delta*rho_ep_ (ir)) 349 end do 350 351 ! d2gam(1) = d2gam/drho+_drho+ (particle=ELECTRON), dgam/drho-_drho- (particle=POSITRON) 352 ! d2gam(2) = d2gam/drho-_drho+ (particle=ELECTRON), dgam/drho+_drho- (particle=POSITRON) 353 ! d2gam(3) = d2gam/drho-_drho- (particle=ELECTRON), dgam/drho+_drho+ (particle=POSITRON) 354 ABI_ALLOCATE(d2gam,(mesh_size,3)) 355 d2gam(:,:)=zero 356 do ir=1,mesh_size 357 if (rho_ (ir)>tol14) d2gam(ir,1)=(gam_(ir,1,2)+gam_(ir,1,3)-two*gam_(ir,1,1))/(delta*rho_ (ir))**2 358 if (rho_ep_(ir)>tol14) then 359 d2gam(ir,3)=(gam_(ir,1,4)+gam_(ir,1,5)-two*gam_(ir,1,1))/(delta*rho_ep_(ir))**2 360 if (rho_(ir)>tol14) then 361 d2gam(ir,2)=(gam_(ir,1,6)+gam_(ir,1,7)+two*gam_(ir,1,1) & 362 & -gam_(ir,1,2)-gam_(ir,1,3)-gam_(ir,1,4)-gam_(ir,1,5)) & 363 & *half/(delta*rho_(ir))/(delta*rho_ep_(ir)) 364 end if 365 end if 366 end do 367 368 ABI_DEALLOCATE(rho_) 369 ABI_DEALLOCATE(rho_ep_) 370 ! Compute useful sums of densities 371 ABI_ALLOCATE(v1sum,(mesh_size,3)) 372 if ( dtset%pawxcdev>=2) then 373 ABI_ALLOCATE(v2sum,(mesh_size,lm_size,3)) 374 else 375 ABI_ALLOCATE(v2sum,(0,0,0)) 376 end if 377 rhotot(:,1)=sqfpi*rhosph(:);rhotot_ep(:,1)=sqfpi*rhosph_ep(:) 378 call pawxcsum(1,1,1,lmselect,lmselect_ep,lm_size,mesh_size,3,dtset%pawxcdev,& 379 & pawang,rhotot,rhotot_ep,v1sum,v2sum) 380 ! Compute final developpment of gamma moments 381 gammam(:,:)=zero 382 gammam(:,1)=gam_(:,1,1)*sqfpi 383 gammam(:,1)=gammam(:,1)+(d2gam(:,2)*v1sum(:,2) & 384 & +half*(d2gam(:,1)*v1sum(:,1)+d2gam(:,3)*v1sum(:,3)))/sqfpi 385 do ilm=2,lm_size 386 if (lmselect(ilm)) then 387 gammam(:,ilm)=gammam(:,ilm)+d1gam(:,1)*rhotot(:,ilm) 388 end if 389 if (lmselect_ep(ilm)) then 390 gammam(:,ilm)=gammam(:,ilm)+d1gam(:,2)*rhotot_ep(:,ilm) 391 end if 392 end do 393 if (dtset%pawxcdev>1) then 394 do ilm=2,lm_size 395 gammam(:,ilm)=gammam(:,ilm)+d2gam(:,2)*v2sum(:,ilm,2) & 396 & +half*(d2gam(:,1)*v2sum(:,ilm,1)+d2gam(:,3)*v2sum(:,ilm,3)) 397 end do 398 end if 399 400 ABI_DEALLOCATE(gam_) 401 ABI_DEALLOCATE(d1gam) 402 ABI_DEALLOCATE(d2gam) 403 ABI_DEALLOCATE(v1sum) 404 ABI_DEALLOCATE(v2sum) 405 ! Compute contribution to annihilation rate 406 ABI_ALLOCATE(gg,(mesh_size,4)) 407 gg=zero 408 ABI_ALLOCATE(rhoarr2,(mesh_size)) 409 do ilm=1,lm_size 410 if (lmselect_ep(ilm)) gg(:,1)=gg(:,1)+rhotot_ep(:,ilm)*rhocorej(:)*gammam(:,ilm) 411 end do 412 ABI_DEALLOCATE(rhoarr2) 413 gg(1:mesh_size,1)=gg(1:mesh_size,1)*pawrad(itypat)%rad(1:mesh_size)**2 414 call simp_gen(intg,gg(:,1),pawrad(itypat)) 415 rate =rate +intg 416 ABI_DEALLOCATE(gg) 417 ABI_DEALLOCATE(gammam) 418 ABI_DEALLOCATE(rhotot) 419 ABI_DEALLOCATE(rhotot_ep) 420 ABI_DEALLOCATE(rhosph) 421 ABI_DEALLOCATE(rhosph_ep) 422 423 end if ! dtset%pawxcdev 424 ABI_DEALLOCATE(rhocore) 425 426 ABI_DEALLOCATE(rho1) 427 ABI_DEALLOCATE(trho1) 428 ABI_DEALLOCATE(rho1_ep) 429 ABI_DEALLOCATE(trho1_ep) 430 ABI_DEALLOCATE(lmselect) 431 ABI_DEALLOCATE(lmselect_ep) 432 ABI_DEALLOCATE(lmselect_dum) 433 ABI_DEALLOCATE(nhat1) 434 ABI_DEALLOCATE(nhat1_ep) 435 436 !Reduction in case of distribution over atomic sites 437 if (mpi_enreg%nproc_atom>1) then 438 mpibuf=rate 439 call xmpi_sum(mpibuf,mpi_enreg%comm_atom,ierr) 440 rate=mpibuf 441 end if 442 443 DBG_EXIT("COLL") 444 445 end subroutine posratecore