TABLE OF CONTENTS
ABINIT/partial_dos_fractions [ Functions ]
NAME
partial_dos_fractions
FUNCTION
calculate partial DOS fractions to feed to the tetrahedron method 1: project states on angular momenta 2: should be able to choose certain atoms or atom types, slabs of space...
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (MVer,MB,MG) 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
crystal<crystal_t>= data type gathering info on symmetries and unit cell dtset<type(dataset_type)>=all input variables for this dataset npwarr(nkpt)=number of planewaves in basis at this k point kg(3,mpw*mkmem)=reduced planewave coordinates. cg(2,mcg)=planewave coefficients of wavefunctions mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol collect=1 if fractions should be MPI collected at the end, 0 otherwise. mpi_enreg=information about MPI parallelization
SIDE EFFECTS
dos%fractions(ikpt,iband,isppol,ndosfraction) = percentage of s, p, d..
character on each atom for the wavefunction # ikpt,iband, isppol
== if prtdosm /= 0
dos%fractions_m(ikpt,iband,isppol,ndosfraction*mbesslang) = percentage of s, p, d..
character on each atom for the wavefunction # ikpt,iband, isppol (m-resolved)
NOTES
psi(r) = (4pi/sqrt(ucvol)) \sum_{LMG} i**l u(G) e^{i(k+G).Ra} x Y_{LM}^*(k+G) Y_{LM}(r-Ra) j_L(|k+G||r-Ra|)
int_(ratsph) |psi(r)|**2 = \sum_LM rho(LM)
where
rho_{LM} = 4pi \int_o^{rc} dr r**2 ||\sum_G u(G) Y_LM^*(k+G) e^{i(k+G).Ra} j_L(|k+G| r)||**2
where S is a RSH. The final expression is:
When k = G0/2, we have u_{G0/2}(G) = u_{G0/2}(-G-G0)^* and P can be rewritten as
P = (4pi i^L}/sqrt(ucvol) \sum^'_G w(G) S_{LM}(k+G) \int_0^ratsph dr r^2 j_L(|k+G|r) x
2 Re[u_k(G) e^{i(k+G).R_atom}] if L = 2n
2 Im[u_k(G) e^{i(k+G).R_atom}] if L = 2n + 1
where the sum over G is done on the reduced G-sphere and w(G) = 1/2 if G=G0 else 1.
PARENTS
outscfcv
CHILDREN
cg_getspin,cwtime,destroy_mpi_enreg,getkpgnorm,getph,initmpi_seq initylmg,jlspline_free,ph1d3d,recip_ylm,sort_dp,splint,wrtout,xmpi_sum
SOURCE
64 #if defined HAVE_CONFIG_H 65 #include "config.h" 66 #endif 67 68 #include "abi_common.h" 69 70 subroutine partial_dos_fractions(dos,crystal,dtset,npwarr,kg,cg,mcg,collect,mpi_enreg) 71 72 use defs_basis 73 use defs_abitypes 74 use m_profiling_abi 75 use m_errors 76 use m_splines 77 use m_xmpi 78 use m_mpinfo 79 use m_crystal 80 use m_sort 81 82 use m_time, only : cwtime 83 use m_numeric_tools, only : simpson 84 use m_special_funcs, only : jlspline_t, jlspline_new, jlspline_free, jlspline_integral 85 use m_cgtools, only : cg_getspin 86 use m_epjdos, only : recip_ylm, epjdos_t 87 use m_gsphere, only : getkpgnorm 88 use m_kg, only : ph1d3d, getph 89 90 !This section has been created automatically by the script Abilint (TD). 91 !Do not modify the following lines by hand. 92 #undef ABI_FUNC 93 #define ABI_FUNC 'partial_dos_fractions' 94 use interfaces_14_hidewrite 95 use interfaces_32_util 96 use interfaces_51_manage_mpi 97 use interfaces_56_recipspace 98 !End of the abilint section 99 100 implicit none 101 102 !Arguments ------------------------------------ 103 !scalars 104 integer,intent(in) :: mcg,collect 105 type(epjdos_t),intent(inout) :: dos 106 type(MPI_type),intent(inout) :: mpi_enreg 107 type(dataset_type),intent(in) :: dtset 108 type(crystal_t),intent(in) :: crystal 109 !arrays 110 integer,intent(in) :: kg(3,dtset%mpw*dtset%mkmem),npwarr(dtset%nkpt) 111 real(dp),intent(in) :: cg(2,mcg) 112 113 !Local variables------------------------------- 114 !scalars 115 integer,parameter :: prtsphere0=0 ! do not output all the band by band details for projections. 116 integer :: shift_b,shift_sk,iat,iatom,iband,ierr,ikpt,ilang,ioffkg,is1, is2, isoff 117 integer :: ipw,ispinor,isppol,ixint,mbess,mcg_disk,me_kpt,shift_cg 118 integer :: me_g0,mgfft,my_nspinor,n1,n2,n3,natsph_tot,npw_k,nradintmax 119 integer :: comm_pw,comm_kpt,rc_ylm,itypat,nband_k 120 real(dp),parameter :: bessint_delta = 0.1_dp 121 real(dp) :: arg,bessarg,bessargmax,kpgmax,rmax 122 real(dp) :: cpu,wall,gflops 123 character(len=500) :: msg 124 type(jlspline_t) :: jlspl 125 type(MPI_type) :: mpi_enreg_seq 126 !arrays 127 integer :: iindex(dtset%mpw),nband_tmp(1),npwarr_tmp(1) 128 integer,allocatable :: iatsph(:),nradint(:),atindx(:),typat_extra(:),kg_k(:,:) 129 real(dp) :: kpoint(3),spin(3),ylmgr_dum(1) 130 real(dp) :: xfit(dtset%mpw),yfit(dtset%mpw) 131 real(dp),allocatable :: ylm_k(:,:) 132 real(dp),allocatable :: bess_fit(:,:,:) 133 real(dp),allocatable :: cg_1band(:,:),cg_1kpt(:,:),kpgnorm(:),ph1d(:,:) 134 real(dp),allocatable :: ph3d(:,:,:),ratsph(:),rint(:),sum_1atom_1ll(:,:) 135 real(dp),allocatable :: sum_1atom_1lm(:,:) 136 real(dp),allocatable :: xred_sph(:,:),znucl_sph(:),phkxred(:,:) 137 complex(dpc) :: cgcmat(2,2) 138 139 !************************************************************************* 140 141 ! for the moment, only support projection on angular momenta 142 if (dos%partial_dos_flag /= 1 .and. dos%partial_dos_flag /= 2) then 143 write(std_out,*) 'Error: partial_dos_fractions only supports angular ' 144 write(std_out,*) ' momentum projection and spinor components for the moment. return to outscfcv' 145 write(std_out,*) ' partial_dos = ', dos%partial_dos_flag 146 return 147 end if 148 149 ! impose all kpoints have same number of bands 150 do isppol=1,dtset%nsppol 151 do ikpt=1,dtset%nkpt 152 if (dtset%nband((isppol-1)*dtset%nkpt + ikpt) /= dtset%mband) then 153 write(std_out,*) 'Error: partial_dos_fractions wants same number of bands at each kpoint' 154 write(std_out,*) ' isppol, ikpt = ', isppol,ikpt, dtset%nband((isppol-1)*dtset%nkpt + ikpt), dtset%mband 155 write(std_out,*) ' all nband = ', dtset%nband 156 return 157 end if 158 end do 159 end do 160 161 ! Real or complex spherical harmonics? 162 rc_ylm = 2; if (dos%prtdosm == 2) rc_ylm = 1 163 164 comm_kpt = mpi_enreg%comm_kpt 165 me_kpt = mpi_enreg%me_kpt 166 comm_pw = mpi_enreg%comm_bandfft ; me_g0 = mpi_enreg%me_g0 167 my_nspinor = max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 168 169 n1 = dtset%ngfft(1); n2 = dtset%ngfft(2); n3 = dtset%ngfft(3) 170 mgfft = maxval(dtset%ngfft(1:3)) 171 172 call cwtime(cpu, wall, gflops, "start") 173 174 !############################################################## 175 !FIRST CASE: project on angular momenta to get dos parts 176 !############################################################## 177 178 if (dos%partial_dos_flag == 1) then 179 natsph_tot = dtset%natsph + dtset%natsph_extra 180 181 ABI_ALLOCATE(iatsph, (natsph_tot)) 182 ABI_ALLOCATE(typat_extra, (natsph_tot)) 183 ABI_ALLOCATE(ratsph, (natsph_tot)) 184 ABI_ALLOCATE(znucl_sph, (natsph_tot)) 185 ABI_ALLOCATE(nradint, (natsph_tot)) 186 ABI_ALLOCATE(atindx, (natsph_tot)) 187 ABI_ALLOCATE(phkxred, (2,natsph_tot)) 188 189 ! initialize atindx 190 do iatom=1,natsph_tot 191 atindx(iatom) = iatom 192 end do 193 194 iatsph(1:dtset%natsph) = dtset%iatsph(1:dtset%natsph) 195 do iatom=1,dtset%natsph 196 itypat = dtset%typat(iatsph(iatom)) 197 typat_extra(iatom) = itypat 198 ratsph(iatom) = dtset%ratsph(itypat) 199 znucl_sph(iatom) = dtset%znucl(itypat) 200 end do 201 202 ! fictitious atoms are declared with 203 ! %natsph_extra, %ratsph_extra and %xredsph_extra(3, dtset%natsph_extra) 204 ! they have atom index (natom + ii) and itype = ntypat + 1 205 do iatom=1,dtset%natsph_extra 206 typat_extra(iatom+dtset%natsph) = dtset%ntypat + 1 207 ratsph(iatom+dtset%natsph) = dtset%ratsph_extra 208 znucl_sph(iatom+dtset%natsph) = zero 209 iatsph(iatom+dtset%natsph) = dtset%natom + iatom 210 end do 211 212 ! init bessel function integral for recip_ylm max ang mom + 1 213 ABI_ALLOCATE(sum_1atom_1ll,(dos%mbesslang,natsph_tot)) 214 ABI_ALLOCATE(sum_1atom_1lm,(dos%mbesslang**2,natsph_tot)) 215 216 ! Note ecuteff instead of ecut. 217 kpgmax = sqrt(dtset%ecut * dtset%dilatmx**2) 218 rmax = zero; bessargmax = zero; nradintmax = 0 219 do iatom=1,natsph_tot 220 rmax = max(rmax, ratsph(iatom)) 221 bessarg = ratsph(iatom)*two_pi*kpgmax 222 bessargmax = max(bessargmax, bessarg) 223 nradint(iatom) = int (bessarg / bessint_delta) + 1 224 nradintmax = max(nradintmax,nradint(iatom)) 225 end do 226 !write(std_out,*)' partial_dos_fractions: rmax=', rmax,' nradintmax: ", nradintmax 227 ! use same number of grid points to calculate Bessel function and to do the integration later on r 228 ! and make sure bessargmax is a multiple of bessint_delta 229 mbess = nradintmax 230 bessargmax = bessint_delta*mbess 231 232 ABI_ALLOCATE(rint,(nradintmax)) 233 ABI_ALLOCATE(bess_fit,(dtset%mpw,nradintmax,dos%mbesslang)) 234 235 ! initialize general Bessel function array on uniform grid xx, from 0 to (2 \pi |k+G|_{max} |r_{max}|) 236 jlspl = jlspline_new(mbess, bessint_delta, dos%mbesslang) 237 238 ABI_ALLOCATE(xred_sph, (3, natsph_tot)) 239 do iatom=1,dtset%natsph 240 xred_sph(:,iatom) = crystal%xred(:,iatsph(iatom)) 241 end do 242 do iatom=1,dtset%natsph_extra 243 xred_sph(:,dtset%natsph+iatom) = dtset%xredsph_extra(:, iatom) 244 end do 245 246 ABI_ALLOCATE(ph1d,(2,(2*n1+1 + 2*n2+1 + 2*n3+1)*natsph_tot)) 247 call getph(atindx,natsph_tot,n1,n2,n3,ph1d,xred_sph) 248 249 ! Fake MPI data to be used for sequential call to initylmg. 250 call initmpi_seq(mpi_enreg_seq) 251 mpi_enreg_seq%my_natom = dtset%natom 252 253 shift_sk = 0 254 do isppol=1,dtset%nsppol 255 ioffkg = 0 256 do ikpt=1,dtset%nkpt 257 nband_k = dtset%nband((isppol-1)*dtset%nkpt + ikpt) 258 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_kpt)) cycle 259 npw_k = npwarr(ikpt) 260 kpoint(:) = dtset%kpt(:,ikpt) 261 262 ! make phkred for all atoms 263 do iat=1,natsph_tot 264 arg=two_pi*( kpoint(1)*xred_sph(1,iat) + kpoint(2)*xred_sph(2,iat) + kpoint(3)*xred_sph(3,iat) ) 265 phkxred(1,iat)=cos(arg) 266 phkxred(2,iat)=sin(arg) 267 end do 268 269 ABI_MALLOC(kg_k, (3, npw_k)) 270 kg_k = kg(:,ioffkg+1:ioffkg+npw_k) 271 272 ! kpgnorm contains norms only for kpoints used by this processor 273 ABI_ALLOCATE(kpgnorm, (npw_k)) 274 call getkpgnorm(crystal%gprimd, kpoint, kg_k, kpgnorm, npw_k) 275 276 ! Now get Ylm(k, G) factors: returns "real Ylms", which are real (+m) and 277 ! imaginary (-m) parts of actual complex Ylm. Yl-m = Ylm* 278 ! Call initylmg for a single k-point (mind mpi_enreg_seq). 279 ABI_MALLOC(ylm_k, (npw_k, dos%mbesslang**2)) 280 npwarr_tmp(1) = npw_k; nband_tmp(1) = nband_k 281 call initylmg(crystal%gprimd,kg_k,kpoint,1,mpi_enreg_seq,dos%mbesslang,& 282 npw_k,nband_tmp,1,npwarr_tmp,1,0,crystal%rprimd,ylm_k,ylmgr_dum) 283 284 ! get phases exp (2 pi i (k+G).x_tau) in ph3d 285 ABI_ALLOCATE(ph3d,(2,npw_k,natsph_tot)) 286 call ph1d3d(1,natsph_tot,kg_k,natsph_tot,natsph_tot,npw_k,n1,n2,n3,phkxred,ph1d,ph3d) 287 288 ! get Bessel function factors on array of |k+G|*r distances 289 ! since we need many r distances and have a large number of different 290 ! |k+G|, get j_l on uniform grid (above, in array gen_besj), 291 ! and spline it for each kpt Gvector set. 292 ! Note that we use the same step based on rmax, this can lead to (hopefully small) 293 ! inaccuracies when we integrate from 0 up to rmax(iatom) 294 do ixint=1,nradintmax 295 rint(ixint) = (ixint-1)*rmax / (nradintmax-1) 296 do ipw=1,npw_k 297 xfit(ipw) = two_pi * kpgnorm(ipw) * rint(ixint) 298 iindex(ipw) = ipw 299 end do 300 301 call sort_dp(npw_k,xfit,iindex,tol14) 302 do ilang=1,dos%mbesslang 303 call splint(mbess, jlspl%xx, jlspl%bess_spl(:,ilang), jlspl%bess_spl_der(:,ilang), npw_k, xfit, yfit) 304 ! re-order results for different G vectors 305 do ipw=1,npw_k 306 bess_fit(iindex(ipw),ixint,ilang) = yfit(ipw) 307 end do 308 end do 309 end do ! ixint 310 311 shift_b = 0 312 do iband=1,nband_k 313 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband,isppol,me_kpt)) cycle 314 !write(std_out,*)"in band:",iband 315 316 do ispinor=1,my_nspinor 317 ! Select wavefunction in cg array 318 shift_cg = shift_sk + shift_b 319 320 call recip_ylm(bess_fit, cg(:,shift_cg+1:shift_cg+npw_k), comm_pw, dtset%istwfk(ikpt),& 321 & nradint, nradintmax, me_g0, dos%mbesslang , dtset%mpw, natsph_tot, typat_extra, dos%mlang_type,& 322 & npw_k,ph3d, prtsphere0, rint, ratsph, rc_ylm, sum_1atom_1ll, sum_1atom_1lm,& 323 & crystal%ucvol, ylm_k, znucl_sph) 324 325 ! Accumulate 326 do iatom=1,natsph_tot 327 do ilang=1,dos%mbesslang 328 dos%fractions(ikpt,iband,isppol,dos%mbesslang*(iatom-1) + ilang) & 329 & = dos%fractions(ikpt,iband,isppol,dos%mbesslang*(iatom-1) + ilang) & 330 & + sum_1atom_1ll(ilang,iatom) 331 end do 332 end do 333 334 if (dos%prtdosm /= 0) then 335 do iatom=1,natsph_tot 336 do ilang=1,dos%mbesslang**2 337 dos%fractions_m(ikpt,iband,isppol,dos%mbesslang**2*(iatom-1) + ilang) & 338 & = dos%fractions_m(ikpt,iband,isppol,dos%mbesslang**2*(iatom-1) + ilang) & 339 & + sum_1atom_1lm(ilang,iatom) 340 end do 341 end do 342 end if 343 344 ! Increment band, spinor shift 345 shift_b = shift_b + npw_k 346 end do ! spinor 347 end do ! band 348 349 ! Increment kpt and (spin, kpt) shifts 350 ioffkg = ioffkg + npw_k 351 shift_sk = shift_sk + nband_k*my_nspinor*npw_k 352 353 ABI_FREE(kg_k) 354 ABI_FREE(kpgnorm) 355 ABI_FREE(ylm_k) 356 ABI_DEALLOCATE(ph3d) 357 end do ! ikpt 358 end do ! isppol 359 360 ! collect = 1 ==> gather all contributions from different processors 361 if (collect == 1) then 362 call xmpi_sum(dos%fractions,comm_kpt,ierr) 363 if (dos%prtdosm /= 0) call xmpi_sum(dos%fractions_m,comm_kpt,ierr) 364 365 if (mpi_enreg%paral_spinor == 1)then 366 call xmpi_sum(dos%fractions,mpi_enreg%comm_spinor,ierr) 367 if (dos%prtdosm /= 0) call xmpi_sum(dos%fractions_m,mpi_enreg%comm_spinor,ierr) 368 end if 369 end if 370 371 ABI_DEALLOCATE(atindx) 372 ABI_DEALLOCATE(bess_fit) 373 ABI_DEALLOCATE(iatsph) 374 ABI_DEALLOCATE(typat_extra) 375 ABI_DEALLOCATE(nradint) 376 ABI_DEALLOCATE(ph1d) 377 ABI_DEALLOCATE(phkxred) 378 ABI_DEALLOCATE(ratsph) 379 ABI_DEALLOCATE(rint) 380 ABI_DEALLOCATE(sum_1atom_1ll) 381 ABI_DEALLOCATE(sum_1atom_1lm) 382 ABI_DEALLOCATE(xred_sph) 383 ABI_DEALLOCATE(znucl_sph) 384 385 call jlspline_free(jlspl) 386 call destroy_mpi_enreg(mpi_enreg_seq) 387 388 !############################################################## 389 !2ND CASE: project on spinors 390 !############################################################## 391 392 else if (dos%partial_dos_flag == 2) then 393 394 if (dtset%nsppol /= 1 .or. dtset%nspinor /= 2) then 395 MSG_WARNING("spinor projection is meaningless if nsppol==2 or nspinor/=2. Not calculating projections.") 396 return 397 end if 398 if (my_nspinor /= 2) then 399 MSG_WARNING("spinor projection with spinor parallelization is not coded. Not calculating projections.") 400 return 401 end if 402 ABI_CHECK(mpi_enreg%paral_spinor == 0, "prtdos 5 does not support spinor parallelism") 403 404 ! FIXME: We should not allocate such a large chunk of memory! 405 mcg_disk = dtset%mpw*my_nspinor*dtset%mband 406 ABI_ALLOCATE(cg_1kpt,(2,mcg_disk)) 407 shift_sk = 0 408 isppol = 1 409 410 do ikpt=1,dtset%nkpt 411 nband_k = dtset%nband((isppol-1)*dtset%nkpt + ikpt) 412 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_kpt)) cycle 413 npw_k = npwarr(ikpt) 414 415 cg_1kpt(:,:) = cg(:,shift_sk+1:shift_sk+mcg_disk) 416 ABI_ALLOCATE(cg_1band,(2,2*npw_k)) 417 shift_b=0 418 do iband=1,nband_k 419 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband,isppol,me_kpt)) cycle 420 421 ! Select wavefunction in cg array 422 !shift_cg = shift_sk + shift_b 423 cg_1band(:,:) = cg_1kpt(:,shift_b+1:shift_b+2*npw_k) 424 call cg_getspin(cg_1band, npw_k, spin, cgcmat=cgcmat) 425 426 ! MG: TODO: imag part of off-diagonal terms is missing. 427 ! I will add them later on. 428 do is1=1,2 429 do is2=1,2 430 isoff = is2 + (is1-1)*2 431 dos%fractions(ikpt,iband,isppol,isoff) = dos%fractions(ikpt,iband,isppol,isoff) & 432 & + real(cgcmat(is1,is2)) 433 end do 434 end do 435 436 dos%fractions(ikpt,iband,isppol,5) = dos%fractions(ikpt,iband,isppol,5) + spin(1) 437 dos%fractions(ikpt,iband,isppol,6) = dos%fractions(ikpt,iband,isppol,6) + spin(2) 438 dos%fractions(ikpt,iband,isppol,7) = dos%fractions(ikpt,iband,isppol,7) + spin(3) 439 440 shift_b = shift_b + 2*npw_k 441 end do 442 ABI_DEALLOCATE(cg_1band) 443 shift_sk = shift_sk + nband_k*2*npw_k 444 end do 445 ABI_DEALLOCATE(cg_1kpt) 446 447 ! Gather all contributions from different processors 448 if (collect == 1) then 449 call xmpi_sum(dos%fractions,comm_kpt,ierr) 450 call xmpi_sum(dos%fractions,comm_pw,ierr) 451 !below for future use - spinors should not be parallelized for the moment 452 !if (mpi_enreg%paral_spinor == 1)then 453 ! call xmpi_sum(dos%fractions,mpi_enreg%comm_spinor,ierr) 454 !end if 455 end if 456 457 else 458 MSG_WARNING('only partial_dos==1 or 2 is coded') 459 end if 460 461 call cwtime(cpu,wall,gflops,"stop") 462 write(msg,'(2(a,f8.2),a)')" partial_dos_fractions: cpu_time: ",cpu,"[s], walltime: ",wall," [s]" 463 call wrtout(std_out,msg,"PERS") 464 465 end subroutine partial_dos_fractions