TABLE OF CONTENTS
ABINIT/poslifetime [ Functions ]
NAME
poslifetime
FUNCTION
Calculate the positron lifetime
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (GJ,MT,JW) 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 gprimd(3,3)= dimensional reciprocal space primitive translations mpi_enreg= information about MPI parallelization my_natom=number of atoms treated by current processor n3xccc= dimension of the xccc3d array (0 or nfft). nfft= number of FFT grid points ngfft(18)= contain all needed information about 3D FFT nhat(nfft,nspden)=charge compensation density (content depends on electronpositron%particle) option= if 1, calculate positron lifetime for whole density if 2, calculate positron lifetime for given state if 3, calculate positron lifetime for given state with IPM 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 rhor(nfft,nspden)=total electron/positron density (content depends on electronpositron%particle) ucvol=unit cell volume in bohr**3. xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3 ===== Optional arguments, used only if option>1 ===== pawrhoij_dop_el(my_natom*usepaw) <type(pawrhoij_type)>= -PAW only- atomic occupancies of one state rhor_dop_el(nfft)=electron density of given state for the state dependent scheme ===== Optional argument ===== pawrhoij_ep(my_natom*usepaw) <type(pawrhoij_type)>= atomic occupancies to be used in place of electronpositron%pawrhoij_ep
OUTPUT
rate= annihilation rate of a given state needed for state dependent scheme for doppler broadening
SIDE EFFECTS
electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation
PARENTS
outscfcv,posdoppler
CHILDREN
gammapositron,gammapositron_fft,mkdenpos,nderiv_gen,pawdensities pawxcsum,simp_gen,wrtout,xmpi_sum
SOURCE
62 #if defined HAVE_CONFIG_H 63 #include "config.h" 64 #endif 65 66 #include "abi_common.h" 67 68 subroutine poslifetime(dtset,electronpositron,gprimd,my_natom,mpi_enreg,n3xccc,nfft,ngfft,nhat,& 69 & option,pawang,pawrad,pawrhoij,pawtab,rate,rate_paw,rhor,ucvol,xccc3d,& 70 & rhor_dop_el,pawrhoij_dop_el,pawrhoij_ep) ! optional arguments 71 72 73 use defs_basis 74 use defs_abitypes 75 use m_profiling_abi 76 use m_errors 77 use m_xmpi 78 79 use m_electronpositron 80 81 use m_pawang, only : pawang_type 82 use m_pawrad, only : pawrad_type, simp_gen, nderiv_gen 83 use m_pawtab, only : pawtab_type 84 use m_pawrhoij, only : pawrhoij_type 85 use m_pawxc, only: pawxcsum 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 'poslifetime' 91 use interfaces_14_hidewrite 92 use interfaces_41_xc_lowlevel 93 use interfaces_56_xc 94 use interfaces_65_paw 95 !End of the abilint section 96 97 implicit none 98 99 !Arguments ------------------------------------ 100 !scalars 101 integer,intent(in) :: my_natom,n3xccc,nfft,option 102 real(dp),intent(in) :: ucvol 103 real(dp),intent(out) :: rate,rate_paw 104 type(dataset_type), intent(in) :: dtset 105 type(electronpositron_type),pointer :: electronpositron 106 type(MPI_type),intent(in) :: mpi_enreg 107 type(pawang_type), intent(in) :: pawang 108 !arrays 109 integer,intent(in) :: ngfft(18) 110 real(dp),intent(in) :: gprimd(3,3),nhat(nfft,dtset%nspden*dtset%usepaw),xccc3d(n3xccc) 111 real(dp),intent(in),target :: rhor(nfft,dtset%nspden) 112 real(dp),optional,intent(in) :: rhor_dop_el(nfft) 113 type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*dtset%usepaw) 114 type(pawrhoij_type),intent(in) :: pawrhoij(my_natom*dtset%usepaw) 115 type(pawrhoij_type),optional,intent(in) :: pawrhoij_dop_el(my_natom*dtset%usepaw) 116 type(pawrhoij_type),optional,target,intent(in) :: pawrhoij_ep(my_natom*dtset%usepaw) 117 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw) 118 119 !Local variables------------------------------- 120 !scalars 121 integer :: cplex,iatom,ierr,ifft,igam,ii,ilm,ilm1,ilm2,iloop,ipt,ir,isel 122 integer :: itypat,iwarn,iwarnj,iwarnp,lm_size,lmn2_size,mesh_size 123 integer :: nfftot,ngamma,ngr,ngrad,nspden_ep,opt_dens,usecore 124 logical,parameter :: include_nhat_in_gamma=.false. 125 real(dp),parameter :: delta=1.d-4 126 real(dp) :: fact,fact2,intg 127 real(dp) :: lambda_core ,lambda_core_ipm ,lambda ,lambda_ipm 128 real(dp) :: lambda_core_paw,lambda_core_paw_ipm,lambda_paw,lambda_paw_ipm 129 real(dp) :: lifetime,lifetime_ipm,nbec,nbev,nbp,rdum,sqfpi,units 130 character(len=500) :: msg 131 !arrays 132 integer,allocatable :: igamma(:) 133 logical,allocatable :: lmselect(:),lmselect_ep(:),lmselect_dum(:) 134 real(dp) :: mpibuf(4) 135 real(dp),parameter :: qphon(3)=(/zero,zero,zero/),lsign(2)=(/one,-one/) 136 real(dp),allocatable :: d1gam(:,:,:),d2gam(:,:,:),ff(:),gam_(:,:,:),gamma(:,:),gammam(:,:,:),gg(:,:) 137 real(dp),allocatable :: grhocore2(:),grhocor2_(:),grhoe2(:),grho2_(:) 138 real(dp),allocatable :: nhat1(:,:,:),nhat1_ep(:,:,:),nhat1_j(:,:,:) 139 real(dp),allocatable :: rho_(:),rho_ep_(:),rho1(:,:,:),rho1_ep(:,:,:),rho1_j(:,:,:) 140 real(dp),allocatable :: rhoarr1(:),rhoarr1_ep(:),rhoarr1_j(:),rhoarr2(:) 141 real(dp),allocatable :: rhocore(:),rhocor_(:),rhoe(:,:),rhop(:,:),rhor_dop_el_(:) 142 real(dp),allocatable :: rhosph(:),rhosph_ep(:),rhosph_j(:),rhotot(:,:),rhotot_ep(:,:) 143 real(dp),allocatable :: rhotot_j(:,:),trho1(:,:,:),trho1_ep(:,:,:),trho1_j(:,:,:) 144 real(dp),allocatable :: v1sum(:,:),v2sum(:,:,:) 145 real(dp),pointer :: rhor_(:,:),rhor_ep_(:,:) 146 type(pawrhoij_type),pointer :: pawrhoij_ep_(:) 147 148 ! ************************************************************************* 149 150 DBG_ENTER("COLL") 151 152 !Tests for developers 153 if (.not.associated(electronpositron)) then 154 msg='electronpositron variable must be associated!' 155 MSG_BUG(msg) 156 end if 157 if (option/=1) then 158 if ((.not.present(rhor_dop_el)).or.(.not.present(pawrhoij_dop_el))) then 159 msg='when option/=1, rhor_dop_el and pawrhoij_dop_el must be present!' 160 MSG_BUG(msg) 161 end if 162 end if 163 164 !Constants 165 fact=0.0 166 cplex=1;nspden_ep=1 167 usecore=n3xccc/nfft 168 nfftot=ngfft(1)*ngfft(2)*ngfft(3) 169 ngrad=1;if (electronpositron%ixcpositron==3.or.electronpositron%ixcpositron==31) ngrad=2 170 iwarn=0;iwarnj=0;iwarnp=1 171 sqfpi=sqrt(four_pi) 172 173 !Compatibility tests 174 if (electronpositron%particle==EP_NOTHING) then 175 msg='Not valid for electronpositron%particle=NOTHING!' 176 MSG_BUG(msg) 177 end if 178 if (electronpositron%nfft/=nfft) then 179 msg='nfft/=electronpositron%nfft!' 180 MSG_BUG(msg) 181 end if 182 if (dtset%usepaw==1) then 183 if(dtset%pawxcdev==0.and.ngrad==2) then 184 msg='GGA is not implemented for pawxcdev=0 (use dtset%pawxcdev/=0)!' 185 MSG_BUG(msg) 186 end if 187 end if 188 189 !Select type(s) of enhancement factor 190 if ((electronpositron%ixcpositron==1.or.electronpositron%ixcpositron==3).and.option==1) then 191 ngamma=2 192 ABI_ALLOCATE(igamma,(ngamma)) 193 igamma(1)=1;igamma(2)=2 194 else 195 ngamma=1 196 ABI_ALLOCATE(igamma,(ngamma)) 197 if (electronpositron%ixcpositron==-1) igamma(1)=0 198 if (electronpositron%ixcpositron== 2) igamma(1)=4 199 if (electronpositron%ixcpositron==11.or.electronpositron%ixcpositron==31) igamma(1)=3 200 if (electronpositron%ixcpositron==1.or.electronpositron%ixcpositron==3) igamma(1)=2 201 end if 202 203 !Select density according to nhat choice 204 if (dtset%usepaw==0.or.include_nhat_in_gamma) then 205 rhor_ => rhor 206 rhor_ep_ => electronpositron%rhor_ep 207 else 208 ABI_ALLOCATE(rhor_,(nfft,dtset%nspden)) 209 ABI_ALLOCATE(rhor_ep_,(nfft,dtset%nspden)) 210 rhor_=rhor-nhat 211 rhor_ep_=electronpositron%rhor_ep-electronpositron%nhat_ep 212 end if 213 214 !Eventually overwrite electronpositron%pawrhoij_ep 215 if (present(pawrhoij_ep)) then 216 pawrhoij_ep_ => pawrhoij_ep 217 else 218 pawrhoij_ep_ => electronpositron%pawrhoij_ep 219 end if 220 221 !Loop on different enhancement factors 222 do igam=1,ngamma 223 224 ! Compute electron-positron annihilation rate using pseudo densities (plane waves) 225 ! ---------------------------------------------------------------------------------------- 226 227 ! Select the densities and make them positive 228 ABI_ALLOCATE(rhoe,(nfft,nspden_ep)) 229 ABI_ALLOCATE(rhop,(nfft,nspden_ep)) 230 if (electronpositron%particle==EP_ELECTRON) then 231 rhoe(:,1)=rhor_ep_(:,1);rhop(:,1)=rhor_(:,1) 232 else if (electronpositron%particle==EP_POSITRON) then 233 rhoe(:,1)=rhor_(:,1);rhop(:,1)=rhor_ep_(:,1) 234 end if 235 call mkdenpos(iwarn ,nfft,nspden_ep,1,rhoe,dtset%xc_denpos) 236 call mkdenpos(iwarnp,nfft,nspden_ep,1,rhop,dtset%xc_denpos) 237 if (option/=1) then 238 ABI_ALLOCATE(rhor_dop_el_,(nfft)) 239 rhor_dop_el_(:)=rhor_dop_el(:) 240 call mkdenpos(iwarnp,nfft,1,1,rhor_dop_el_,dtset%xc_denpos) 241 end if 242 243 ! Compute enhancement factor at each FFT grid point 244 ! gamma(:,1): using total electronic density 245 ! gamma(:,2): using valence electronic density 246 ABI_ALLOCATE(gamma,(nfft,2)) 247 if (option==1.or.option==2) then 248 call gammapositron_fft(electronpositron,gamma,gprimd,igamma(igam),mpi_enreg,& 249 & n3xccc,nfft,ngfft,rhoe,rhop,xccc3d) 250 else 251 gamma=one 252 end if 253 254 ! Compute positron annihilation rates 255 lambda =zero;lambda_ipm =zero 256 lambda_core=zero;lambda_core_ipm=zero 257 if (option==1) then 258 do ifft=1,nfft 259 lambda =lambda +rhop(ifft,1)*rhoe(ifft,1)*gamma(ifft,1) 260 lambda_ipm=lambda_ipm+rhop(ifft,1)*rhoe(ifft,1)*gamma(ifft,2) 261 end do 262 else 263 do ifft=1,nfft 264 lambda =lambda +rhop(ifft,1)*rhor_dop_el_(ifft)*gamma(ifft,1) 265 lambda_ipm=lambda_ipm+rhop(ifft,1)*rhor_dop_el_(ifft)*gamma(ifft,2) 266 end do 267 end if 268 if (usecore==1) then 269 do ifft=1,nfft 270 lambda_core =lambda_core +rhop(ifft,1)*xccc3d(ifft)*gamma(ifft,1) 271 lambda_core_ipm=lambda_core_ipm+rhop(ifft,1)*xccc3d(ifft) 272 end do 273 end if 274 lambda =lambda *ucvol/dble(nfftot) 275 lambda_ipm =lambda_ipm *ucvol/dble(nfftot) 276 lambda_core =lambda_core *ucvol/dble(nfftot) 277 lambda_core_ipm=lambda_core_ipm*ucvol/dble(nfftot) 278 ABI_DEALLOCATE(gamma) 279 ABI_DEALLOCATE(rhoe) 280 ABI_DEALLOCATE(rhop) 281 if (option/=1) then 282 ABI_DEALLOCATE(rhor_dop_el_) 283 end if 284 ! NC pseudopotential: check electrons/positron number 285 if (dtset%usepaw==0.and.igam==ngamma) then 286 nbec=zero;nbev=zero;nbp=zero 287 if (electronpositron%particle==EP_ELECTRON) then 288 do ifft=1,nfft 289 nbec=nbec+xccc3d(ifft) 290 nbev=nbev+electronpositron%rhor_ep(ifft,1) 291 nbp =nbp +rhor(ifft,1) 292 end do 293 else 294 do ifft=1,nfft 295 nbec=nbec+xccc3d(ifft) 296 nbev=nbev+rhor(ifft,1) 297 nbp =nbp +electronpositron%rhor_ep(ifft,1) 298 end do 299 end if 300 nbec=nbec*ucvol/dble(nfftot) 301 nbev=nbev*ucvol/dble(nfftot) 302 nbp =nbp *ucvol/dble(nfftot) 303 end if 304 305 ! MPI parallelization 306 if(mpi_enreg%nproc_fft>1)then 307 call xmpi_sum(lambda ,mpi_enreg%comm_fft,ierr) 308 call xmpi_sum(lambda_ipm,mpi_enreg%comm_fft,ierr) 309 call xmpi_sum(lambda_core ,mpi_enreg%comm_fft,ierr) 310 call xmpi_sum(lambda_core_ipm,mpi_enreg%comm_fft,ierr) 311 if (dtset%usepaw==0.and.igam==ngamma) then 312 call xmpi_sum(nbec,mpi_enreg%comm_fft,ierr) 313 call xmpi_sum(nbev,mpi_enreg%comm_fft,ierr) 314 call xmpi_sum(nbp ,mpi_enreg%comm_fft,ierr) 315 end if 316 end if 317 318 319 ! PAW: add on-site contributions to electron-positron annihilation rate 320 ! ---------------------------------------------------------------------------------------- 321 if (dtset%usepaw==1) then 322 323 lambda_paw =zero;lambda_paw_ipm =zero 324 lambda_core_paw=zero;lambda_core_paw_ipm=zero 325 326 ! Loop on atoms 327 do iatom=1,my_natom 328 329 itypat=pawrhoij(iatom)%itypat 330 lmn2_size=pawtab(itypat)%lmn2_size 331 mesh_size=pawtab(itypat)%mesh_size 332 lm_size=pawtab(itypat)%lcut_size**2 333 cplex=1 334 ngr=0;if (ngrad==2) ngr=mesh_size 335 336 ! Allocations of "on-site" densities 337 ABI_ALLOCATE(rho1 ,(cplex*mesh_size,lm_size,nspden_ep)) 338 ABI_ALLOCATE(trho1,(cplex*mesh_size,lm_size,nspden_ep)) 339 ABI_ALLOCATE(rho1_ep ,(cplex*mesh_size,lm_size,nspden_ep)) 340 ABI_ALLOCATE(trho1_ep,(cplex*mesh_size,lm_size,nspden_ep)) 341 if (option/=1) then 342 ABI_ALLOCATE(rho1_j ,(cplex*mesh_size,lm_size,nspden_ep)) 343 ABI_ALLOCATE(trho1_j,(cplex*mesh_size,lm_size,nspden_ep)) 344 else 345 ABI_ALLOCATE(rho1_j ,(0,0,0)) 346 ABI_ALLOCATE(trho1_j,(0,0,0)) 347 end if 348 if (include_nhat_in_gamma) then 349 ABI_ALLOCATE(nhat1,(cplex*mesh_size,lm_size,nspden_ep)) 350 ABI_ALLOCATE(nhat1_ep,(cplex*mesh_size,lm_size,nspden_ep)) 351 else 352 ABI_ALLOCATE(nhat1,(0,0,0)) 353 ABI_ALLOCATE(nhat1_ep,(0,0,0)) 354 end if 355 if (include_nhat_in_gamma.and.option/=1) then 356 ABI_ALLOCATE(nhat1_j,(cplex*mesh_size,lm_size,nspden_ep)) 357 else 358 ABI_ALLOCATE(nhat1_j,(0,0,0)) 359 end if 360 ABI_ALLOCATE(lmselect,(lm_size)) 361 ABI_ALLOCATE(lmselect_ep,(lm_size)) 362 ABI_ALLOCATE(lmselect_dum,(lm_size)) 363 364 ! Compute "on-site" densities (n1, ntild1, nhat1) for electron and positron ===== 365 lmselect(:)=.true. 366 opt_dens=1;if (include_nhat_in_gamma) opt_dens=0 367 call pawdensities(rdum,cplex,iatom,lmselect,lmselect_dum,lm_size,nhat1,nspden_ep,1,& 368 & 0,opt_dens,-1,0,pawang,0,pawrad(itypat),pawrhoij(iatom),& 369 & pawtab(itypat),rho1,trho1) 370 lmselect_ep(:)=.true. 371 call pawdensities(rdum,cplex,iatom,lmselect_ep,lmselect_dum,lm_size,nhat1_ep,nspden_ep,1,& 372 & 0,opt_dens,-1,0,pawang,0,pawrad(itypat),pawrhoij_ep_(iatom),& 373 & pawtab(itypat),rho1_ep,trho1_ep) 374 375 ! For state dependent scheme in Doppler ===== 376 ! Compute "on-site" densities (n1, ntild1, nhat1) for a given electron state j===== 377 if (option/=1) then 378 opt_dens=1;if (include_nhat_in_gamma) opt_dens=0 379 call pawdensities(rdum,cplex,iatom,lmselect,lmselect_dum,lm_size,nhat1_j,nspden_ep,1,& 380 & 0,opt_dens,-1,0,pawang,0,pawrad(itypat),pawrhoij_dop_el(iatom),& 381 & pawtab(itypat),rho1_j,trho1_j) 382 end if 383 384 ! Compute contribution to annihilation rate: 385 ! Loop: first step: compute all-electron contribution (from n^1, n_c) 386 ! 2nd step: compute pseudo contribution (from tild_n^1, hat_n^1, tild_n_c) 387 do iloop=1,2 388 if (iloop==1) usecore=1 389 if (iloop==2) usecore=pawtab(itypat)%usetcore 390 ABI_ALLOCATE(rhocore,(mesh_size)) 391 392 ! First formalism: use densities on r,theta,phi 393 if (dtset%pawxcdev==0) then 394 395 ABI_ALLOCATE(gamma,(mesh_size,2)) 396 ABI_ALLOCATE(rhoarr1,(mesh_size)) 397 ABI_ALLOCATE(rhoarr1_ep,(mesh_size)) 398 if (option/=1) then 399 ABI_ALLOCATE(rhoarr1_j,(mesh_size)) 400 end if 401 ! Loop on the angular part 402 do ipt=1,pawang%angl_size 403 ! Build densities 404 rhoarr1=zero;rhoarr1_ep=zero;rhocore=zero 405 if (option/=1) rhoarr1_j=zero 406 if (iloop==1) then 407 do ilm=1,lm_size 408 if (lmselect(ilm)) rhoarr1(:)=rhoarr1(:)+rho1(:,ilm,1)*pawang%ylmr(ilm,ipt) 409 end do 410 if (option/=1) then 411 do ilm=1,lm_size 412 if (lmselect(ilm)) rhoarr1_j(:)=rhoarr1_j(:)+rho1_j(:,ilm,1)*pawang%ylmr(ilm,ipt) 413 end do 414 end if 415 do ilm=1,lm_size 416 if (lmselect_ep(ilm)) rhoarr1_ep(:)=rhoarr1_ep(:)+rho1_ep(:,ilm,1)*pawang%ylmr(ilm,ipt) 417 end do 418 if (usecore==1) rhocore(:)=pawtab(itypat)%coredens(:) 419 else 420 if (include_nhat_in_gamma) then 421 do ilm=1,lm_size 422 if (lmselect(ilm)) rhoarr1(:)=rhoarr1(:)+(trho1(:,ilm,1)+nhat1(:,ilm,1))*pawang%ylmr(ilm,ipt) 423 end do 424 if (option/=1) then 425 do ilm=1,lm_size 426 if (lmselect(ilm)) rhoarr1_j(:)=rhoarr1_j(:)+(trho1_j(:,ilm,1)+nhat1_j(:,ilm,1))*pawang%ylmr(ilm,ipt) 427 end do 428 end if 429 do ilm=1,lm_size 430 if (lmselect_ep(ilm)) rhoarr1_ep(:)=rhoarr1_ep(:)+(trho1_ep(:,ilm,1)+nhat1_ep(:,ilm,1))*pawang%ylmr(ilm,ipt) 431 end do 432 else 433 do ilm=1,lm_size 434 if (lmselect(ilm)) rhoarr1(:)=rhoarr1(:)+trho1(:,ilm,1)*pawang%ylmr(ilm,ipt) 435 end do 436 if (option/=1) then 437 do ilm=1,lm_size 438 if (lmselect(ilm)) rhoarr1_j(:)=rhoarr1(:)+trho1_j(:,ilm,1)*pawang%ylmr(ilm,ipt) 439 end do 440 end if 441 do ilm=1,lm_size 442 if (lmselect_ep(ilm)) rhoarr1_ep(:)=rhoarr1_ep(:)+trho1_ep(:,ilm,1)*pawang%ylmr(ilm,ipt) 443 end do 444 end if 445 if (usecore==1) rhocore(:)=pawtab(itypat)%tcoredens(:,1) 446 end if 447 ! Make the densities positive 448 if (electronpositron%particle==EP_ELECTRON) then 449 call mkdenpos(iwarnp,mesh_size,1,1,rhoarr1 ,dtset%xc_denpos) 450 call mkdenpos(iwarn ,mesh_size,1,1,rhoarr1_ep,dtset%xc_denpos) 451 else if (electronpositron%particle==EP_POSITRON) then 452 call mkdenpos(iwarn ,mesh_size,1,1,rhoarr1 ,dtset%xc_denpos) 453 call mkdenpos(iwarnp,mesh_size,1,1,rhoarr1_ep,dtset%xc_denpos) 454 if (option/=1) then 455 call mkdenpos(iwarnj,mesh_size,1,1,rhoarr1_j,dtset%xc_denpos) 456 end if 457 end if 458 ! Compute Gamma 459 ABI_ALLOCATE(grhoe2,(ngr)) 460 ABI_ALLOCATE(grhocore2,(ngr)) 461 if (option==1.or.option==2) then 462 if (electronpositron%particle==EP_ELECTRON) then 463 call gammapositron(gamma,grhocore2,grhoe2,igamma(igam),ngr,mesh_size,& 464 & rhocore,rhoarr1_ep,rhoarr1,usecore) 465 else if (electronpositron%particle==EP_POSITRON) then 466 call gammapositron(gamma,grhocore2,grhoe2,igamma(igam),ngr,mesh_size,& 467 & rhocore,rhoarr1,rhoarr1_ep,usecore) 468 end if 469 else 470 gamma(:,:)=one 471 end if 472 ABI_DEALLOCATE(grhoe2) 473 ABI_DEALLOCATE(grhocore2) 474 ! Compute contribution to annihilation rates 475 ABI_ALLOCATE(ff,(mesh_size)) 476 if (option/=1) rhoarr1(:)=rhoarr1_j(:) 477 do ii=1,4 478 if (ii==1) ff(1:mesh_size)=rhoarr1(1:mesh_size)*rhoarr1_ep(1:mesh_size) & 479 & *gamma(1:mesh_size,1)*pawrad(itypat)%rad(1:mesh_size)**2 480 if (ii==2) ff(1:mesh_size)=rhoarr1(1:mesh_size)*rhoarr1_ep(1:mesh_size) & 481 & *gamma(1:mesh_size,2)*pawrad(itypat)%rad(1:mesh_size)**2 482 if (electronpositron%particle==EP_ELECTRON) then 483 if (ii==3) ff(1:mesh_size)=rhoarr1(1:mesh_size)*rhocore(1:mesh_size) & 484 & *gamma(1:mesh_size,1)*pawrad(itypat)%rad(1:mesh_size)**2 485 if (ii==4) ff(1:mesh_size)=rhoarr1(1:mesh_size)*rhocore(1:mesh_size) & 486 & *pawrad(itypat)%rad(1:mesh_size)**2 487 else 488 if (ii==3) ff(1:mesh_size)=rhoarr1_ep(1:mesh_size)*rhocore(1:mesh_size) & 489 & *gamma(1:mesh_size,1)*pawrad(itypat)%rad(1:mesh_size)**2 490 if (ii==4) ff(1:mesh_size)=rhoarr1_ep(1:mesh_size)*rhocore(1:mesh_size) & 491 & *pawrad(itypat)%rad(1:mesh_size)**2 492 end if 493 call simp_gen(intg,ff,pawrad(itypat)) 494 intg=intg*pawang%angwgth(ipt)*four_pi 495 if (ii==1) lambda_paw =lambda_paw +lsign(iloop)*intg 496 if (ii==2) lambda_paw_ipm =lambda_paw_ipm +lsign(iloop)*intg 497 if (ii==3) lambda_core_paw =lambda_core_paw +lsign(iloop)*intg 498 if (ii==4) lambda_core_paw_ipm=lambda_core_paw_ipm+lsign(iloop)*intg 499 end do 500 ABI_DEALLOCATE(ff) 501 end do ! ipt 502 ABI_DEALLOCATE(gamma) 503 ABI_DEALLOCATE(rhoarr1) 504 ABI_DEALLOCATE(rhoarr1_ep) 505 if (option/=1) then 506 ABI_DEALLOCATE(rhoarr1_j) 507 end if 508 509 ! Second formalism: use (l,m) moments for densities 510 else if (dtset%pawxcdev/=0) then 511 512 ! Build densities 513 ABI_ALLOCATE(gammam,(mesh_size,2,lm_size)) 514 ABI_ALLOCATE(rhotot,(mesh_size,lm_size)) 515 ABI_ALLOCATE(rhotot_ep,(mesh_size,lm_size)) 516 ABI_ALLOCATE(rhosph,(mesh_size)) 517 ABI_ALLOCATE(rhosph_ep,(mesh_size)) 518 if (option/=1) then 519 ABI_ALLOCATE(rhotot_j,(mesh_size,lm_size)) 520 ABI_ALLOCATE(rhosph_j,(mesh_size)) 521 end if 522 if (usecore==0) rhocore(:)=zero 523 if (iloop==1) then 524 rhotot (:,:)=rho1 (:,:,1) 525 rhotot_ep(:,:)=rho1_ep(:,:,1) 526 if (option/=1) rhotot_j (:,:)=rho1_j (:,:,1) 527 if (usecore==1) rhocore(:)=pawtab(itypat)%coredens(:) 528 else 529 if (include_nhat_in_gamma) then 530 rhotot (:,:)=trho1 (:,:,1)+nhat1 (:,:,1) 531 rhotot_ep(:,:)=trho1_ep(:,:,1)+nhat1_ep(:,:,1) 532 if (option/=1) rhotot_j (:,:)=trho1_j (:,:,1)+nhat1_j (:,:,1) 533 else 534 rhotot (:,:)=trho1 (:,:,1) 535 rhotot_ep(:,:)=trho1_ep(:,:,1) 536 if (option/=1) rhotot_j (:,:)=trho1_j (:,:,1) 537 end if 538 if (usecore==1) rhocore(:)=pawtab(itypat)%tcoredens(:,1) 539 end if 540 rhosph (:)=rhotot (:,1)/sqfpi 541 rhosph_ep(:)=rhotot_ep(:,1)/sqfpi 542 if (option/=1) rhosph_j (:)=rhotot_j(:,1)/sqfpi 543 ! Make spherical densities positive 544 if (electronpositron%particle==EP_ELECTRON) then 545 call mkdenpos(iwarnp,mesh_size,1,1,rhosph ,dtset%xc_denpos) 546 call mkdenpos(iwarn ,mesh_size,1,1,rhosph_ep,dtset%xc_denpos) 547 else if (electronpositron%particle==EP_POSITRON) then 548 call mkdenpos(iwarn ,mesh_size,1,1,rhosph ,dtset%xc_denpos) 549 call mkdenpos(iwarnp,mesh_size,1,1,rhosph_ep,dtset%xc_denpos) 550 if (option/=1) then 551 call mkdenpos(iwarnp,mesh_size,1,1,rhosph_j,dtset%xc_denpos) 552 end if 553 end if 554 ! Need gradients of electronic densities for GGA 555 ABI_ALLOCATE(grhoe2,(ngr)) 556 ABI_ALLOCATE(grhocore2,(ngr)) 557 if (ngr>0) then 558 if (electronpositron%particle==EP_ELECTRON) then 559 call nderiv_gen(grhoe2,rhosph_ep,pawrad(itypat)) 560 else if (electronpositron%particle==EP_POSITRON) then 561 call nderiv_gen(grhoe2,rhosph,pawrad(itypat)) 562 end if 563 grhoe2(:)=grhoe2(:)**2 564 if (usecore==1) then 565 call nderiv_gen(grhocore2,rhocore,pawrad(itypat)) 566 grhocore2(:)=grhocore2(:)**2 567 end if 568 end if 569 ! Compute Gamma for (rho-,rho+), 570 ! (rho- +drho-,rho+), (rho- -drho-,rho+), 571 ! (rho-,rho+ +drho+), (rho-,rho+ -drho+), 572 ! (rho- +drho-,rho+ +drho+), (rho- -drho-,rho+ -drho+) 573 ! Do a seven steps loop 574 ABI_ALLOCATE(gam_,(mesh_size,2,7)) 575 ABI_ALLOCATE(rho_,(mesh_size)) 576 ABI_ALLOCATE(rho_ep_,(mesh_size)) 577 ABI_ALLOCATE(rhocor_,(mesh_size)) 578 ABI_ALLOCATE(grho2_,(ngr)) 579 ABI_ALLOCATE(grhocor2_,(ngr)) 580 do ii=1,7 581 ! Apply delta to get perturbed densities 582 rho_(:)=rhosph(:);rho_ep_(:)=rhosph_ep(:);if (usecore==1) rhocor_(:)=rhocore(:) 583 if (ngr>0) grho2_(:)=grhoe2(:) 584 if (ngr>0) grhocor2_(:)=grhocore2(:) 585 if (ii==2.or.ii==4.or.ii==6) fact=(one+delta) 586 if (ii==3.or.ii==5.or.ii==7) fact=(one-delta) 587 fact2=fact**2 588 if (ii==2.or.ii==3.or.ii==6.or.ii==7) then 589 rho_(:)=fact*rho_(:) 590 if (electronpositron%particle==EP_POSITRON) then 591 if (ngr>0) grho2_(:)=fact2*grho2_(:) 592 if (usecore==1)rhocor_(:)=fact*rhocor_(:) 593 if (ngr>0.and.usecore==1) grhocor2_(:)=fact2*grhocor2_(:) 594 end if 595 end if 596 if (ii==4.or.ii==5.or.ii==6.or.ii==7) then 597 rho_ep_(:)=fact*rho_ep_(:) 598 if (electronpositron%particle==EP_ELECTRON) then 599 if (ngr>0) grho2_(:)=fact2*grho2_(:) 600 if (usecore==1)rhocor_(:)=fact*rhocor_(:) 601 if (ngr>0.and.usecore==1) grhocor2_(:)=fact2*grhocor2_(:) 602 end if 603 end if 604 ! Compute gamma for these perturbed densities 605 if (option==1.or.option==2) then 606 if (electronpositron%particle==EP_ELECTRON) then 607 call gammapositron(gam_(:,:,ii),grhocor2_,grho2_,igamma(igam),ngr,mesh_size,rhocor_,rho_ep_,rho_,usecore) 608 else if (electronpositron%particle==EP_POSITRON) then 609 call gammapositron(gam_(:,:,ii),grhocor2_,grho2_,igamma(igam),ngr,mesh_size,rhocor_,rho_,rho_ep_,usecore) 610 end if 611 else 612 gam_(:,:,:)=one 613 end if 614 end do ! end loop ii=1,7 615 616 ABI_DEALLOCATE(rhocor_) 617 ABI_DEALLOCATE(grho2_) 618 ABI_DEALLOCATE(grhocor2_) 619 ABI_DEALLOCATE(grhoe2) 620 ABI_DEALLOCATE(grhocore2) 621 rho_ (:)=rhosph (:);if (electronpositron%particle==EP_POSITRON.and.usecore==1) rho_ (:)=rho_ (:)+rhocore(:) 622 rho_ep_(:)=rhosph_ep(:);if (electronpositron%particle==EP_ELECTRON.and.usecore==1) rho_ep_(:)=rho_ep_(:)+rhocore(:) 623 ! Compute numerical first and second derivatives of Gamma 624 ! d1gam(1) = dgam/drho+ (particle=ELECTRON), dgam/drho- (particle=POSITRON) 625 ! d1gam(2) = dgam/drho- (particle=ELECTRON), dgam/drho+ (particle=POSITRON) 626 ABI_ALLOCATE(d1gam,(mesh_size,2,2)) 627 d1gam(:,:,:)=zero 628 do ir=1,mesh_size 629 if (rho_ (ir)>tol14) d1gam(ir,1,1)=(gam_(ir,1,2)-gam_(ir,1,3))*half/(delta*rho_ (ir)) 630 if (rhosph (ir)>tol14) d1gam(ir,2,1)=(gam_(ir,2,2)-gam_(ir,2,3))*half/(delta*rhosph (ir)) 631 if (rho_ep_ (ir)>tol14) d1gam(ir,1,2)=(gam_(ir,1,4)-gam_(ir,1,5))*half/(delta*rho_ep_ (ir)) 632 if (rhosph_ep(ir)>tol14) d1gam(ir,2,2)=(gam_(ir,2,4)-gam_(ir,2,5))*half/(delta*rhosph_ep(ir)) 633 end do 634 635 ! d2gam(1) = d2gam/drho+_drho+ (particle=ELECTRON), dgam/drho-_drho- (particle=POSITRON) 636 ! d2gam(2) = d2gam/drho-_drho+ (particle=ELECTRON), dgam/drho+_drho- (particle=POSITRON) 637 ! d2gam(3) = d2gam/drho-_drho- (particle=ELECTRON), dgam/drho+_drho+ (particle=POSITRON) 638 ABI_ALLOCATE(d2gam,(mesh_size,2,3)) 639 d2gam(:,:,:)=zero 640 do ir=1,mesh_size 641 if (rho_ (ir)>tol14) d2gam(ir,1,1)=(gam_(ir,1,2)+gam_(ir,1,3)-two*gam_(ir,1,1))/(delta*rho_ (ir))**2 642 if (rhosph(ir)>tol14) d2gam(ir,2,1)=(gam_(ir,2,2)+gam_(ir,2,3)-two*gam_(ir,2,1))/(delta*rhosph(ir))**2 643 if (rho_ep_(ir)>tol14) then 644 d2gam(ir,1,3)=(gam_(ir,1,4)+gam_(ir,1,5)-two*gam_(ir,1,1))/(delta*rho_ep_(ir))**2 645 if (rho_(ir)>tol14) then 646 d2gam(ir,1,2)=(gam_(ir,1,6)+gam_(ir,1,7)+two*gam_(ir,1,1) & 647 & -gam_(ir,1,2)-gam_(ir,1,3)-gam_(ir,1,4)-gam_(ir,1,5)) & 648 & *half/(delta*rho_(ir))/(delta*rho_ep_(ir)) 649 end if 650 end if 651 if (rhosph_ep(ir)>tol14) then 652 d2gam(ir,2,3)=(gam_(ir,2,4)+gam_(ir,2,5)-two*gam_(ir,2,1))/(delta*rhosph_ep(ir))**2 653 if (rhosph(ir)>tol14) then 654 d2gam(ir,2,2)=(gam_(ir,2,6)+gam_(ir,2,7)+two*gam_(ir,2,1) & 655 & -gam_(ir,2,2)-gam_(ir,2,3)-gam_(ir,2,4)-gam_(ir,2,5)) & 656 & *half/(delta*rhosph(ir))/(delta*rhosph_ep(ir)) 657 end if 658 end if 659 end do 660 ABI_DEALLOCATE(rho_) 661 ABI_DEALLOCATE(rho_ep_) 662 ! Compute useful sums of densities 663 ABI_ALLOCATE(v1sum,(mesh_size,3)) 664 if ( dtset%pawxcdev>=2) then 665 ABI_ALLOCATE(v2sum,(mesh_size,lm_size,3)) 666 else 667 ABI_ALLOCATE(v2sum,(0,0,0)) 668 end if 669 rhotot(:,1)=sqfpi*rhosph(:);rhotot_ep(:,1)=sqfpi*rhosph_ep(:) 670 call pawxcsum(1,1,1,lmselect,lmselect_ep,lm_size,mesh_size,3,dtset%pawxcdev,& 671 & pawang,rhotot,rhotot_ep,v1sum,v2sum) 672 ! Compute final developpment of gamma moments 673 gammam(:,:,:)=zero 674 gammam(:,:,1)=gam_(:,:,1)*sqfpi 675 gammam(:,1,1)=gammam(:,1,1)+(d2gam(:,1,2)*v1sum(:,2) & 676 & +half*(d2gam(:,1,1)*v1sum(:,1)+d2gam(:,1,3)*v1sum(:,3)))/sqfpi 677 gammam(:,2,1)=gammam(:,2,1)+(d2gam(:,2,2)*v1sum(:,2) & 678 & +half*(d2gam(:,2,1)*v1sum(:,1)+d2gam(:,2,3)*v1sum(:,3)))/sqfpi 679 do ilm=2,lm_size 680 if (lmselect(ilm)) then 681 gammam(:,1,ilm)=gammam(:,1,ilm)+d1gam(:,1,1)*rhotot(:,ilm) 682 gammam(:,2,ilm)=gammam(:,2,ilm)+d1gam(:,2,1)*rhotot(:,ilm) 683 end if 684 if (lmselect_ep(ilm)) then 685 gammam(:,1,ilm)=gammam(:,1,ilm)+d1gam(:,1,2)*rhotot_ep(:,ilm) 686 gammam(:,2,ilm)=gammam(:,2,ilm)+d1gam(:,2,2)*rhotot_ep(:,ilm) 687 end if 688 end do 689 if (dtset%pawxcdev>1) then 690 do ilm=2,lm_size 691 gammam(:,1,ilm)=gammam(:,1,ilm)+d2gam(:,1,2)*v2sum(:,ilm,2) & 692 & +half*(d2gam(:,1,1)*v2sum(:,ilm,1)+d2gam(:,1,3)*v2sum(:,ilm,3)) 693 gammam(:,2,ilm)=gammam(:,2,ilm)+d2gam(:,2,2)*v2sum(:,ilm,2) & 694 & +half*(d2gam(:,2,1)*v2sum(:,ilm,1)+d2gam(:,2,3)*v2sum(:,ilm,3)) 695 end do 696 end if 697 ABI_DEALLOCATE(gam_) 698 ABI_DEALLOCATE(d1gam) 699 ABI_DEALLOCATE(d2gam) 700 ABI_DEALLOCATE(v1sum) 701 ABI_DEALLOCATE(v2sum) 702 703 ! Compute contribution to annihilation rate 704 705 ! In state dependent scheme for Doppler, replace electronic density with one state 706 if (option/=1) then 707 rhotot(:,:) = rhotot_j(:,:) 708 rhosph (:) = rhosph_j (:) 709 end if 710 711 ABI_ALLOCATE(gg,(mesh_size,4)) 712 gg=zero 713 ABI_ALLOCATE(rhoarr1,(mesh_size)) 714 ABI_ALLOCATE(rhoarr2,(mesh_size)) 715 do ilm=1,lm_size 716 do ilm1=1,lm_size 717 if (lmselect(ilm1)) then 718 if (ilm1==1) rhoarr1(:)=sqfpi*rhosph(:) 719 if (ilm1/=1) rhoarr1(:)=rhotot(:,ilm1) 720 do ilm2=1,lm_size 721 if (lmselect_ep(ilm2)) then 722 if (ilm2==1) rhoarr2(:)=sqfpi*rhosph_ep(:) 723 if (ilm2/=1) rhoarr2(:)=rhotot_ep(:,ilm2) 724 if (ilm1>=ilm2) then 725 isel=pawang%gntselect(ilm,ilm2+ilm1*(ilm1-1)/2) 726 else 727 isel=pawang%gntselect(ilm,ilm1+ilm2*(ilm2-1)/2) 728 end if 729 if (isel>0) then 730 fact=pawang%realgnt(isel) 731 gg(:,1)=gg(:,1)+fact*rhoarr1(:)*rhoarr2(:)*gammam(:,1,ilm) 732 gg(:,2)=gg(:,2)+fact*rhoarr1(:)*rhoarr2(:)*gammam(:,2,ilm) 733 end if 734 end if 735 end do 736 end if 737 end do 738 end do 739 ABI_DEALLOCATE(rhoarr1) 740 ABI_DEALLOCATE(rhoarr2) 741 if (electronpositron%particle==EP_ELECTRON) then 742 do ilm=1,lm_size 743 if (lmselect(ilm)) gg(:,3)=gg(:,3)+rhotot(:,ilm)*rhocore(:)*gammam(:,1,ilm) 744 end do 745 gg(:,4)=sqfpi*rhotot(:,1)*rhocore(:) 746 else if (electronpositron%particle==EP_POSITRON) then 747 do ilm=1,lm_size 748 if (lmselect_ep(ilm)) gg(:,3)=gg(:,3)+rhotot_ep(:,ilm)*rhocore(:)*gammam(:,1,ilm) 749 end do 750 gg(:,4)=sqfpi*rhotot_ep(:,1)*rhocore(:) 751 end if 752 do ii=1,4 753 gg(1:mesh_size,ii)=gg(1:mesh_size,ii)*pawrad(itypat)%rad(1:mesh_size)**2 754 call simp_gen(intg,gg(:,ii),pawrad(itypat)) 755 if (ii==1) lambda_paw =lambda_paw +lsign(iloop)*intg 756 if (ii==2) lambda_paw_ipm =lambda_paw_ipm +lsign(iloop)*intg 757 if (ii==3) lambda_core_paw =lambda_core_paw +lsign(iloop)*intg 758 if (ii==4) lambda_core_paw_ipm=lambda_core_paw_ipm+lsign(iloop)*intg 759 end do 760 ABI_DEALLOCATE(gg) 761 ABI_DEALLOCATE(gammam) 762 ABI_DEALLOCATE(rhotot) 763 ABI_DEALLOCATE(rhotot_ep) 764 ABI_DEALLOCATE(rhosph) 765 ABI_DEALLOCATE(rhosph_ep) 766 if (option/=1) then 767 ABI_DEALLOCATE(rhotot_j) 768 ABI_DEALLOCATE(rhosph_j) 769 end if 770 771 end if ! dtset%pawxcdev 772 773 ABI_DEALLOCATE(rhocore) 774 775 end do ! iloop 776 777 ABI_DEALLOCATE(rho1) 778 ABI_DEALLOCATE(trho1) 779 ABI_DEALLOCATE(rho1_ep) 780 ABI_DEALLOCATE(trho1_ep) 781 ABI_DEALLOCATE(rho1_j) 782 ABI_DEALLOCATE(trho1_j) 783 ABI_DEALLOCATE(nhat1) 784 ABI_DEALLOCATE(nhat1_ep) 785 ABI_DEALLOCATE(nhat1_j) 786 ABI_DEALLOCATE(lmselect) 787 ABI_DEALLOCATE(lmselect_ep) 788 ABI_DEALLOCATE(lmselect_dum) 789 790 end do ! iatom 791 792 ! Reduction in case of distribution over atomic sites 793 if (mpi_enreg%nproc_atom>1) then 794 mpibuf(1)=lambda_paw ;mpibuf(2)=lambda_paw_ipm 795 mpibuf(3)=lambda_core_paw;mpibuf(4)=lambda_core_paw_ipm 796 call xmpi_sum(mpibuf,mpi_enreg%comm_atom,ierr) 797 lambda_paw=mpibuf(1) ;lambda_paw_ipm=mpibuf(2) 798 lambda_core_paw=mpibuf(3);lambda_core_paw_ipm=mpibuf(4) 799 end if 800 801 ! Add plane-wave and PAW contributions to annihilation rates 802 lambda =lambda +lambda_paw 803 lambda_ipm =lambda_ipm +lambda_paw_ipm 804 lambda_core =lambda_core +lambda_core_paw 805 lambda_core_ipm=lambda_core_ipm+lambda_core_paw_ipm 806 end if ! dtset%usepaw 807 808 809 ! Convert into proper units and print 810 ! --------------------------------------------------------------------------------------- 811 812 ! Sum valence and core contributions to annihilation rates 813 lambda =lambda +lambda_core 814 lambda_ipm =lambda_ipm +lambda_core_ipm 815 if (dtset%usepaw==1) then 816 lambda_paw =lambda_paw +lambda_core_paw 817 lambda_paw_ipm=lambda_paw_ipm+lambda_core_paw_ipm 818 end if 819 820 ! Set annihilation rate in proper unit (picosec.) 821 units=pi*(one/InvFineStruct)**3/Time_Sec/1.e12_dp/electronpositron%posocc 822 823 lambda =lambda *units 824 lambda_ipm =lambda_ipm *units 825 lambda_core =lambda_core *units 826 lambda_core_ipm=lambda_core_ipm*units 827 lifetime =one/lambda 828 lifetime_ipm=one/lambda_ipm 829 electronpositron%lambda=lambda 830 electronpositron%lifetime=lifetime 831 if (dtset%usepaw==1) then 832 lambda_paw =lambda_paw *units 833 lambda_paw_ipm =lambda_paw_ipm *units 834 lambda_core_paw =lambda_core_paw *units 835 lambda_core_paw_ipm=lambda_core_paw_ipm*units 836 end if 837 rate=lambda-lambda_core-lambda_paw+lambda_core_paw 838 rate_paw=lambda_paw-lambda_core_paw 839 ! Print life time and additional information 840 if (option==1) then 841 if (igam==1) then 842 write(msg,'(a,80("-"),2a)') ch10,ch10,' Results for electron-positron annihilation:' 843 call wrtout(ab_out,msg,'COLL') 844 call wrtout(std_out,msg,'COLL') 845 end if 846 if (ngamma>1.and.igam==1) then 847 write(msg,'(a,i1,a)') ch10,ngamma,& 848 & ' computations of positron lifetime have been performed (with different enhancement factors).' 849 call wrtout(ab_out,msg,'COLL') 850 call wrtout(std_out,msg,'COLL') 851 end if 852 if (ngamma>1) then 853 write(msg,'(2a,i1)') ch10,"########## Lifetime computation ",igam 854 call wrtout(ab_out,msg,'COLL') 855 call wrtout(std_out,msg,'COLL') 856 end if 857 if (abs(electronpositron%ixcpositron)==1) then 858 write(msg,'(4a)') ch10,' # Zero-positron density limit of Arponen and Pajanne provided by Boronski & Nieminen',& 859 & ch10,' Ref.: Boronski and R.M. Nieminen, Phys. Rev. B 34, 3820 (1986)' 860 else if (electronpositron%ixcpositron==11) then 861 write(msg,'(4a)') ch10,' # Zero-positron density limit of Arponen and Pajanne fitted by Sterne & Kaiser',& 862 & ch10,' Ref.: P.A. Sterne and J.H. Kaiser, Phys. Rev. B 43, 13892 (1991)' 863 else if (electronpositron%ixcpositron==2) then 864 write(msg,'(4a)') ch10,' # Electron-positron correlation provided by Puska, Seitsonen, and Nieminen',& 865 & ch10,' Ref: M.J. Puska, A.P. Seitsonen and R.M. Nieminen, Phys. Rev. B 52, 10947 (1994)' 866 else if (electronpositron%ixcpositron==3) then 867 write(msg,'(8a)') ch10,' # Zero-positron density limit of Arponen and Pajanne provided by Boronski & Nieminen',& 868 & ch10,' + GGA corrections',& 869 & ch10,' Ref.: Boronski and R.M. Nieminen, Phys. Rev. B 34, 3820 (1986)',& 870 & ch10,' B. Barbiellini, M.J. Puska, T. Torsti and R.M.Nieminen, Phys. Rev. B 51, 7341 (1994)' 871 else if (electronpositron%ixcpositron==31) then 872 write(msg,'(8a)') ch10,' # Zero-positron density limit of Arponen and Pajanne fitted by Sterne & Kaiser',& 873 & ch10,' + GGA corrections',& 874 & ch10,' Ref.: P.A. Sterne and J.H. Kaiser, Phys. Rev. B 43, 13892 (1991)',& 875 & ch10,' B. Barbiellini, M.J. Puska, T. Torsti and R.M. Nieminen, Phys. Rev. B 51, 7341 (1994)' 876 end if 877 call wrtout(ab_out,msg,'COLL') 878 call wrtout(std_out, msg,'COLL') 879 if (igamma(igam)==0) then 880 write(msg,'(a)') ' # Enhancement factor set to one (test)' 881 else if (igamma(igam)==1) then 882 write(msg,'(3a)') ' # Enhancement factor of Boronski & Nieminen',& 883 & ch10,' Ref.: Boronski and R.M. Nieminen, Phys. Rev. B 34, 3820 (1986)' 884 else if (igamma(igam)==2) then 885 write(msg,'(3a)') ' # Enhancement factor of Boronski & Nieminen IN THE RPA LIMIT',& 886 & ch10,' Ref.: Boronski and R.M. Nieminen, Phys. Rev. B 34, 3820 (1986)' 887 else if (igamma(igam)==3) then 888 write(msg,'(3a)') ' # Enhancement factor of Sterne & Kaiser',& 889 & ch10,' Ref.: P.A. Sterne and J.H. Kaiser, Phys. Rev. B 43, 13892 (1991)' 890 else if (igamma(igam)==4) then 891 write(msg,'(3a)') ' # Enhancement factor of Puska, Seitsonen, and Nieminen',& 892 & ch10,' Ref.: M.J. Puska, A.P. Seitsonen and R.M. Nieminen, Phys. Rev. B 52, 10947 (1994)' 893 end if 894 call wrtout(ab_out,msg,'COLL') 895 call wrtout(std_out,msg,'COLL') 896 write(msg, '(4(2a,es16.8))' ) ch10,& 897 & ' Positron lifetime (ps) =',lifetime ,ch10,& 898 & ' Positron lifetime with IPM for core elec. (ps) =',lifetime_ipm,ch10,& 899 & ' Annihilation rate (ns-1) =',lambda *1000._dp,ch10,& 900 & ' Annihilation rate with IPM for core elec. (ns-1) =',lambda_ipm*1000._dp 901 call wrtout(ab_out,msg,'COLL') 902 call wrtout(std_out,msg,'COLL') 903 write(msg,'(2a,5(2a,es16.8))' ) ch10,& 904 & ' Annihilation rate core/valence decomposition:',ch10,& 905 & ' Core contribution to ann.rate (ns-1) =', lambda_core *1000._dp,ch10,& 906 & ' Valence contribution to ann.rate (ns-1) =',(lambda-lambda_core) *1000._dp,ch10,& 907 & ' Core contribution to ann.rate with IPM (ns-1) =', lambda_core_ipm *1000._dp,ch10,& 908 & ' Valence contribution to ann.rate with IPM (ns-1) =',(lambda_ipm-lambda_core_ipm) *1000._dp 909 call wrtout(ab_out,msg,'COLL') 910 call wrtout(std_out,msg,'COLL') 911 if (dtset%usepaw==1) then 912 write(msg, '(2a,6(2a,es16.8))' ) ch10,& 913 & ' Annihilation rate PAW decomposition:',ch10,& 914 & ' Plane-wave contribution to ann.rate (ns-1) =',(lambda-lambda_paw)*1000._dp,ch10,& 915 & ' Plane-wave valence contribution to ann.rate (ns-1) =',(lambda-lambda_paw-lambda_core+lambda_core_paw)*1000._dp,ch10,& 916 & ' On-site core contribution to ann.rate (ns-1) =', lambda_core_paw*1000._dp,ch10,& 917 & ' On-site valence contribution to ann.rate (ns-1) =',(lambda_paw-lambda_core_paw)*1000._dp,ch10,& 918 & ' Plane-wave contribution to ann.rate with IPM (ns-1) =',(lambda_ipm-lambda_paw_ipm)*1000._dp,ch10,& 919 & ' Plane-wave core contrb. to ann.rate with IPM (ns-1) =',(lambda_core_ipm-lambda_core_paw_ipm)*1000._dp 920 call wrtout(ab_out,msg,'COLL') 921 call wrtout(std_out,msg,'COLL') 922 end if 923 if (dtset%usepaw==0.and.igam==ngamma) then ! These tests are not relevant with PAW 924 write(msg, '(2a,3(2a,es16.8))' ) ch10,& 925 & ' ########## Some checks, for testing purpose:',ch10,& 926 & ' Number of core electrons =',nbec,ch10,& 927 & ' Number of valence electrons =',nbev,ch10,& 928 & ' Number of positrons =',nbp 929 call wrtout(ab_out,msg,'COLL') 930 call wrtout(std_out,msg,'COLL') 931 end if 932 end if !end if option 933 end do ! Big loop on igam 934 935 if (option==1) then 936 write(msg, '(3a)' ) ch10,' (*) IPM=Independent particle Model',ch10 937 call wrtout(ab_out,msg,'COLL') 938 call wrtout(std_out,msg,'COLL') 939 end if !end if option 940 941 !Deallocate memory 942 ABI_DEALLOCATE(igamma) 943 if (dtset%usepaw==1.and.(.not.include_nhat_in_gamma)) then 944 ABI_DEALLOCATE(rhor_) 945 ABI_DEALLOCATE(rhor_ep_) 946 end if 947 948 DBG_EXIT("COLL") 949 950 end subroutine poslifetime