TABLE OF CONTENTS
ABINIT/pawdensities [ Functions ]
NAME
pawdensities
FUNCTION
Compute PAW on-site densities (all-electron ,pseudo and compensation) for a given atom
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (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
cplex: if 1, on-site densities are REAL, if 2, COMPLEX (response function only) iatom=index of current atom (note: this is the absolute index, not the index on current proc) lm_size=number of (l,m) moments lmselectin(lm_size)=flags selecting the non-zero LM-moments of on-site densities (value of these flags at input; must be .TRUE. for nzlmopt/=1) nspden=number of spin-density components nzlmopt=if -1, compute all LM-moments of densities (lmselectin=.true. forced) initialize "lmselectout" (index of non-zero LM-moments of densities) if 0, compute all LM-moments of densities (lmselectin=.true. forced) force "lmselectout" to .true. (index of non-zero LM-moments of densities) if 1, compute only non-zero LM-moments of densities (stored before in "lmselectin") one_over_rad2(mesh_size)= contains 1/r**2 for each point of the radial grid -optional argument- opt_compch=flag controlling the accumulation of compensation charge density moments inside PAW spheres (compch_sph) opt_dens=flag controlling which on-site density(ies) is (are) computed 0: all on-site densities (all-electron, pseudo and compensation) 1: all-electron and pseudo densities (no compensation) 2: only all-electron density opt_l=controls which l-moment(s) contribute to the density: <0 : all l contribute >=0: only l=opt_l contributes Note: opt_l>=0 is only compatible with opt_dens=2 opt_print=1 if the densities moments have to be printed out (only if pawprtvol>=2) pawang <type(pawang_type)>=paw angular mesh and related data pawprtvol=control print volume and debugging output for PAW pawrad <type(pawrad_type)>=paw radial mesh and related data (for the current atom type) pawrhoij <type(pawrhoij_type)>= paw rhoij occupancies and related data (for the current atom) pawtab <type(pawtab_type)>=paw tabulated starting data (for the current atom type)
OUTPUT
nhat1(cplex*mesh_size,lm_size,nspden)= compensation charge on-site density for current atom rho1(cplex*mesh_size,lm_size,nspden)= all electron on-site density for current atom trho1(cplex*mesh_size,lm_size,nspden)= pseudo on-site density for current atom ==== if nzlmopt/=1 lmselectout(lm_size)=flags selecting the non-zero LM-moments of on-site densities (value of these flags at output if updated, i.e. if nzlmopt<1)
SIDE EFFECTS
==== if opt_compch==1 compch_sph=compensation charge integral inside spheres computed over spherical meshes updated with the contribution of current atom
PARENTS
make_efg_onsite,pawdenpot,pawdfptenergy,poslifetime,posratecore
CHILDREN
pawrad_deducer0,simp_gen,wrtout
SOURCE
67 #if defined HAVE_CONFIG_H 68 #include "config.h" 69 #endif 70 71 #include "abi_common.h" 72 73 subroutine pawdensities(compch_sph,cplex,iatom,lmselectin,lmselectout,lm_size,nhat1,nspden,nzlmopt,& 74 & opt_compch,opt_dens,opt_l,opt_print,pawang,pawprtvol,pawrad,pawrhoij,pawtab,rho1,trho1,& 75 & one_over_rad2) ! optional 76 77 use m_profiling_abi 78 use defs_basis 79 use m_errors 80 use m_pawang, only : pawang_type 81 use m_pawrad, only : pawrad_type, pawrad_deducer0, simp_gen 82 use m_pawtab, only : pawtab_type 83 use m_pawrhoij, only : pawrhoij_type 84 85 !This section has been created automatically by the script Abilint (TD). 86 !Do not modify the following lines by hand. 87 #undef ABI_FUNC 88 #define ABI_FUNC 'pawdensities' 89 use interfaces_14_hidewrite 90 !End of the abilint section 91 92 implicit none 93 94 !Arguments --------------------------------------------- 95 !scalars 96 integer,intent(in) :: cplex,iatom,lm_size,nspden,nzlmopt,opt_compch,opt_dens,opt_l,opt_print,pawprtvol 97 ! jmb real(dp),intent(out) :: compch_sph 98 real(dp),intent(inout) :: compch_sph 99 type(pawang_type),intent(in) :: pawang 100 type(pawrad_type),intent(in) :: pawrad 101 type(pawrhoij_type),intent(in) :: pawrhoij 102 type(pawtab_type),intent(in) :: pawtab 103 !arrays 104 logical,intent(in) :: lmselectin(lm_size) 105 logical,intent(inout) :: lmselectout(lm_size) 106 real(dp),intent(in),target,optional :: one_over_rad2(pawtab%mesh_size) 107 real(dp),intent(out) :: nhat1(cplex*pawtab%mesh_size,lm_size,nspden*(1-((opt_dens+1)/2))) 108 real(dp),intent(out) :: rho1(cplex*pawtab%mesh_size,lm_size,nspden) 109 real(dp),intent(out) :: trho1(cplex*pawtab%mesh_size,lm_size,nspden*(1-(opt_dens/2))) 110 111 !Local variables --------------------------------------- 112 !scalars 113 integer :: dplex,ii,ilm,iplex,ir,irhoij,isel,ispden,jrhoij 114 integer :: klm,klmn,kln,ll,lmax,lmin,mesh_size 115 real(dp) :: m1,mt1,rdum 116 character(len=500) :: msg 117 !arrays 118 real(dp) :: compchspha(cplex),compchsphb(cplex),ro(cplex),ro_ql(cplex),ro_rg(cplex) 119 real(dp),allocatable :: aa(:),bb(:) 120 real(dp),pointer :: one_over_rad2_(:) 121 122 ! ************************************************************************* 123 124 DBG_ENTER("COLL") 125 126 !Compatibility tests 127 if (opt_dens/=2.and.opt_l>=0) then 128 msg=' opt_dens/=2 incompatible with opt_l>=0 !' 129 MSG_BUG(msg) 130 end if 131 if(nzlmopt/=0.and.nzlmopt/=1.and.nzlmopt/=-1) then 132 msg=' invalid value for variable "nzlmopt".' 133 MSG_BUG(msg) 134 end if 135 if(nspden>pawrhoij%nspden) then 136 msg=' nspden must be <= pawrhoij%nspden !' 137 MSG_BUG(msg) 138 end if 139 if (cplex>pawrhoij%cplex) then 140 msg=' cplex must be <= pawrhoij%cplex !' 141 MSG_BUG(msg) 142 end if 143 if (nzlmopt/=1) then 144 if (any(.not.lmselectin(1:lm_size))) then 145 msg=' With nzlmopt/=1, lmselectin must be true !' 146 MSG_BUG(msg) 147 end if 148 end if 149 if (pawang%gnt_option==0) then 150 msg=' pawang%gnt_option=0 !' 151 MSG_BUG(msg) 152 end if 153 154 !Various inits 155 rho1=zero 156 if (opt_dens <2) trho1=zero 157 if (opt_dens==0) nhat1=zero 158 mesh_size=pawtab%mesh_size;dplex=cplex-1 159 if (nzlmopt<1) lmselectout(1:lm_size)=.true. 160 if (present(one_over_rad2)) then 161 one_over_rad2_ => one_over_rad2 162 else 163 ABI_ALLOCATE(one_over_rad2_,(mesh_size)) 164 one_over_rad2_(1)=zero 165 one_over_rad2_(2:mesh_size)=one/pawrad%rad(2:mesh_size)**2 166 end if 167 168 !===== Compute "on-site" densities (n1, ntild1, nhat1) ===== 169 !========================================================== 170 171 do ispden=1,nspden 172 173 ! -- Loop over ij channels (basis components) 174 jrhoij=1 175 do irhoij=1,pawrhoij%nrhoijsel 176 klmn=pawrhoij%rhoijselect(irhoij) 177 klm =pawtab%indklmn(1,klmn) 178 kln =pawtab%indklmn(2,klmn) 179 lmin=pawtab%indklmn(3,klmn) 180 lmax=pawtab%indklmn(4,klmn) 181 182 183 ! Retrieve rhoij 184 if (pawrhoij%nspden/=2) then 185 ro(1:cplex)=pawrhoij%rhoijp(jrhoij:jrhoij+dplex,ispden) 186 else 187 if (ispden==1) then 188 ro(1:cplex)=pawrhoij%rhoijp(jrhoij:jrhoij+dplex,1)& 189 & +pawrhoij%rhoijp(jrhoij:jrhoij+dplex,2) 190 else if (ispden==2) then 191 ro(1:cplex)=pawrhoij%rhoijp(jrhoij:jrhoij+dplex,1) 192 end if 193 end if 194 ro(1:cplex)=pawtab%dltij(klmn)*ro(1:cplex) 195 196 ! First option: all on-site densities are computed (opt_dens==0) 197 ! -------------------------------------------------------------- 198 if (opt_dens==0) then 199 do ll=lmin,lmax,2 200 do ilm=ll**2+1,(ll+1)**2 201 if (lmselectin(ilm)) then 202 isel=pawang%gntselect(ilm,klm) 203 if (isel>0) then 204 ro_ql(1:cplex)=ro(1:cplex)*pawtab%qijl(ilm,klmn) 205 ro_rg(1:cplex)=ro(1:cplex)*pawang%realgnt(isel) 206 ! == nhat1(r=0) 207 nhat1(1:cplex,ilm,ispden)=nhat1(1:cplex,ilm,ispden) & 208 & +ro_ql(1:cplex)*pawtab%shapefunc(1,ll+1) 209 ! == rho1(r>0), trho1(r>0), nhat1(r>0) 210 do ir=2,mesh_size 211 rho1(cplex*ir-dplex:ir*cplex,ilm,ispden) =rho1(cplex*ir-dplex:ir*cplex,ilm,ispden)& 212 & +ro_rg(1:cplex)*pawtab%phiphj (ir,kln)*one_over_rad2_(ir) 213 trho1(cplex*ir-dplex:ir*cplex,ilm,ispden)=trho1(cplex*ir-dplex:ir*cplex,ilm,ispden)& 214 & +ro_rg(1:cplex)*pawtab%tphitphj(ir,kln)*one_over_rad2_(ir) 215 nhat1(cplex*ir-dplex:ir*cplex,ilm,ispden)=nhat1(cplex*ir-dplex:ir*cplex,ilm,ispden)& 216 & +ro_ql(1:cplex)*pawtab%shapefunc(ir,ll+1) 217 end do 218 end if 219 end if 220 end do ! End loops over ll,lm 221 end do 222 223 ! 2nd option: AE and pseudo densities are computed (opt_dens==1) 224 ! -------------------------------------------------------------- 225 else if (opt_dens==1) then 226 do ll=lmin,lmax,2 227 do ilm=ll**2+1,(ll+1)**2 228 if (lmselectin(ilm)) then 229 isel=pawang%gntselect(ilm,klm) 230 if (isel>0) then 231 ro_rg(1:cplex)=ro(1:cplex)*pawang%realgnt(isel) 232 ! == rho1(r>0), trho1(r>0) 233 do ir=2,mesh_size 234 rho1(cplex*ir-dplex:ir*cplex,ilm,ispden) =rho1(cplex*ir-dplex:ir*cplex,ilm,ispden)& 235 & +ro_rg(1:cplex)*pawtab%phiphj (ir,kln)*one_over_rad2_(ir) 236 trho1(cplex*ir-dplex:ir*cplex,ilm,ispden)=trho1(cplex*ir-dplex:ir*cplex,ilm,ispden)& 237 & +ro_rg(1:cplex)*pawtab%tphitphj(ir,kln)*one_over_rad2_(ir) 238 end do 239 end if 240 end if 241 end do ! End loops over ll,lm 242 end do 243 244 ! 3rd option: only all-electron on-site density is computed (opt_dens==2) 245 ! ----------------------------------------------------------------------- 246 else if (opt_dens==2) then 247 if (opt_l<0.or.(pawtab%indklmn(3,klmn)==0.and.pawtab%indklmn(4,klmn)==2*opt_l)) then 248 do ll=lmin,lmax,2 249 do ilm=ll**2+1,(ll+1)**2 250 if (lmselectin(ilm)) then 251 isel=pawang%gntselect(ilm,klm) 252 if (isel>0) then 253 ro_rg(1:cplex)=ro(1:cplex)*pawang%realgnt(isel) 254 ! == rho1(r>0) 255 do ir=2,mesh_size 256 rho1(cplex*ir-dplex:ir*cplex,ilm,ispden) =rho1(cplex*ir-dplex:ir*cplex,ilm,ispden)& 257 & +ro_rg(1:cplex)*pawtab%phiphj(ir,kln)*one_over_rad2_(ir) 258 end do 259 end if 260 end if 261 end do ! End loops over ll, lm 262 end do 263 end if 264 end if 265 266 ! -- End loop over ij channels 267 jrhoij=jrhoij+pawrhoij%cplex 268 end do 269 270 ! Scale densities with 1/r**2 and compute rho1(r=0) and trho1(r=0) 271 if (cplex==2) then 272 ABI_ALLOCATE(aa,(5)) 273 ABI_ALLOCATE(bb,(5)) 274 end if 275 if (opt_dens==0.or.opt_dens==1) then 276 do ll=0,pawtab%lcut_size-1 277 do ilm=ll**2+1,(ll+1)**2 278 if (lmselectin(ilm)) then 279 if (cplex==1) then 280 call pawrad_deducer0(rho1 (:,ilm,ispden),mesh_size,pawrad) 281 call pawrad_deducer0(trho1(:,ilm,ispden),mesh_size,pawrad) 282 else 283 do ii=0,1 284 do ir=2,5 285 aa(ir)=rho1 (2*ir-ii,ilm,ispden) 286 bb(ir)=trho1(2*ir-ii,ilm,ispden) 287 end do 288 call pawrad_deducer0(aa,5,pawrad) 289 call pawrad_deducer0(bb,5,pawrad) 290 rho1 (2-ii,ilm,ispden)=aa(1) 291 trho1(2-ii,ilm,ispden)=bb(1) 292 end do 293 end if 294 end if 295 end do 296 end do 297 else 298 do ll=0,pawtab%lcut_size-1 299 do ilm=ll**2+1,(ll+1)**2 300 if (lmselectin(ilm)) then 301 if (cplex==1) then 302 call pawrad_deducer0(rho1(:,ilm,ispden),mesh_size,pawrad) 303 else 304 do ii=0,1 305 do ir=2,5 306 aa(ir)=rho1 (2*ir-ii,ilm,ispden) 307 end do 308 call pawrad_deducer0(aa,5,pawrad) 309 rho1(2-ii,ilm,ispden)=aa(1) 310 end do 311 end if 312 end if 313 end do 314 end do 315 end if 316 if (cplex==2) then 317 ABI_DEALLOCATE(aa) 318 ABI_DEALLOCATE(bb) 319 end if 320 321 ! -- Test moments of densities and store non-zero ones 322 if (nzlmopt==-1) then 323 do ll=0,pawtab%lcut_size-1 324 do ilm=ll**2+1,(ll+1)**2 325 m1=zero;mt1=zero 326 if (cplex==1) then 327 m1=maxval(abs(rho1 (1:mesh_size,ilm,ispden))) 328 if (opt_dens<2) mt1=maxval(abs(trho1(1:mesh_size,ilm,ispden))) 329 else 330 do ir=1,mesh_size 331 rdum=sqrt(rho1(2*ir-1,ilm,ispden)**2+rho1(2*ir,ilm,ispden)**2) 332 m1=max(m1,rdum) 333 end do 334 if (opt_dens<2) then 335 do ir=1,mesh_size 336 rdum=sqrt(trho1(2*ir-1,ilm,ispden)**2+trho1(2*ir,ilm,ispden)**2) 337 mt1=max(mt1,rdum) 338 end do 339 end if 340 end if 341 if (ispden==1) then 342 if ((ilm>1).and.(m1<tol16).and.(mt1<tol16)) then 343 lmselectout(ilm)=.false. 344 end if 345 else if (.not.(lmselectout(ilm))) then 346 lmselectout(ilm)=((m1>=tol16).or.(mt1>=tol16)) 347 end if 348 end do 349 end do 350 end if 351 352 ! -- Compute integral of (n1-tn1) inside spheres 353 if (opt_compch==1.and.ispden==1.and.opt_dens<2) then 354 ABI_ALLOCATE(aa,(mesh_size)) 355 aa(1:mesh_size)=(rho1(1:mesh_size,1,1)-trho1(1:mesh_size,1,1)) & 356 & *pawrad%rad(1:mesh_size)**2 357 ! jmb 358 ! compchspha(1)=zero 359 call simp_gen(compchspha(1),aa,pawrad) 360 compch_sph=compch_sph+compchspha(1)*sqrt(four_pi) 361 ABI_DEALLOCATE(aa) 362 end if 363 364 ! -- Print out moments of densities (if requested) 365 if (abs(pawprtvol)>=2.and.opt_print==1.and.opt_dens<2) then 366 ABI_ALLOCATE(aa,(cplex*mesh_size)) 367 ABI_ALLOCATE(bb,(cplex*mesh_size)) 368 if (opt_dens==0) then 369 write(msg,'(2a,i3,a,i1,3a)') ch10, & 370 & ' Atom ',iatom,' (ispden=',ispden,'):',ch10,& 371 & ' ******* Moment of (n1-tn1) ** Moment of (n1-tn1-nhat1)' 372 else 373 write(msg,'(2a,i3,a,i1,3a)') ch10, & 374 & ' Atom ',iatom,' (ispden=',ispden,'):',ch10,& 375 & ' ******* Moment of (n1-tn1)' 376 end if 377 call wrtout(std_out,msg,'PERS') 378 do ll=0,pawtab%lcut_size-1 379 do ilm=ll**2+1,(ll+1)**2 380 if (lmselectin(ilm)) then 381 do iplex=1,cplex 382 if (opt_dens==0) then 383 do ir=1,mesh_size 384 ii=cplex*(ir-1)+iplex 385 ro(1)=pawrad%rad(ir)**(2+ll) 386 aa(ir)=ro(1)*(rho1(ii,ilm,ispden)-trho1(ii,ilm,ispden)) 387 bb(ir)=ro(1)*nhat1(ii,ilm,ispden) 388 end do 389 call simp_gen(compchspha(iplex),aa,pawrad) 390 call simp_gen(compchsphb(iplex),bb,pawrad) 391 else 392 do ir=1,mesh_size 393 ii=cplex*(ir-1)+iplex 394 ro(1)=pawrad%rad(ir)**(2+ll) 395 aa(ir)=ro(1)*(rho1(ii,ilm,ispden)-trho1(ii,ilm,ispden)) 396 end do 397 call simp_gen(compchspha(iplex),aa,pawrad) 398 end if 399 end do 400 if (opt_dens==0) then 401 if (cplex==1) then 402 write(msg,'(3x,a,2i2,2(a,es14.7))') & 403 & 'l,m=',ll,ilm-(ll**2+ll+1),': M=',compchspha(1),& 404 & ' ** M=',compchspha(1)-compchsphb(1) 405 else 406 write(msg,'(3x,a,2i2,2(a,2es14.7))') & 407 & 'l,m=',ll,ilm-(ll**2+ll+1),': M=',compchspha(1:2),& 408 & ' ** M=',compchspha(1:2)-compchsphb(1:2) 409 end if 410 else 411 if (cplex==1) then 412 write(msg,'(3x,a,2i2,a,es14.7)') & 413 & 'l,m=',ll,ilm-(ll**2+ll+1),': M=',compchspha(1) 414 else 415 write(msg,'(3x,a,2i2,a,2es14.7)') & 416 & 'l,m=',ll,ilm-(ll**2+ll+1),': M=',compchspha(1:2) 417 end if 418 end if 419 call wrtout(std_out,msg,'PERS') 420 end if 421 end do 422 end do 423 ABI_DEALLOCATE(aa) 424 ABI_DEALLOCATE(bb) 425 end if 426 427 ! ----- End loop over spin components 428 end do 429 430 if (.not.present(one_over_rad2)) then 431 ABI_DEALLOCATE(one_over_rad2_) 432 end if 433 434 DBG_EXIT("COLL") 435 436 end subroutine pawdensities