TABLE OF CONTENTS
- ABINIT/m_epjdos
- m_epjdos/dens_in_sph
- m_epjdos/dos_calcnwrite
- m_epjdos/epjdos_free
- m_epjdos/epjdos_new
- m_epjdos/epjdos_t
- m_epjdos/fatbands_ncwrite
- m_epjdos/partial_dos_fractions
- m_epjdos/partial_dos_fractions_paw
- m_epjdos/prtfatbands
- m_epjdos/recip_ylm
- m_epjdos/sphericaldens
ABINIT/m_epjdos [ Modules ]
NAME
m_epjdos
FUNCTION
Tools for the computiation of electronic PJDOSes
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group (MVer, XG, SM, MT, BAmadon, MG, MB) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_epjdos 23 24 use defs_basis 25 use m_abicore 26 use m_xmpi 27 use m_errors 28 use m_htetra 29 use m_splines 30 use m_cgtools 31 use m_atomdata 32 use m_crystal 33 use m_ebands 34 use m_nctk 35 #ifdef HAVE_NETCDF 36 use netcdf 37 #endif 38 use m_hdr 39 use m_mpinfo 40 use m_sort 41 use m_dtset 42 43 use defs_abitypes, only : MPI_type 44 use defs_datatypes, only : ebands_t, pseudopotential_type 45 use m_occ, only : dos_hdr_write 46 use m_time, only : cwtime, timab 47 use m_io_tools, only : open_file 48 use m_numeric_tools, only : simpson, simpson_int 49 use m_fstrings, only : int2char4, strcat 50 use m_special_funcs, only : jlspline_t, jlspline_new, jlspline_free, jlspline_integral 51 use m_kpts, only : tetra_from_kptrlatt 52 use m_kg, only : ph1d3d, getph 53 use m_gsphere, only : getkpgnorm 54 use m_fftcore, only : sphereboundary 55 use m_fft, only : fftpac, fourwf, fourdp 56 use m_pawrad, only : pawrad_type, simp_gen 57 use m_pawtab, only : pawtab_type 58 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_free 59 use m_initylmg, only : initylmg 60 61 implicit none 62 63 private
m_epjdos/dens_in_sph [ Functions ]
[ Top ] [ m_epjdos ] [ Functions ]
NAME
dens_in_sph
FUNCTION
Calculate integrated density in sphere around each atom
INPUTS
cg = wavefunction coefficitents in recip space gmet = metric in recip space istwfk = storage mode for cg coefficients kg_k = G vector indices natom = number of atoms mpi_enreg=information about MPI parallelization ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft npw_k = number of plane waves for this kpoint ph1d = phase factors for different atoms for all G vectors rmax(natom) = max radius to integrate to (in bohr)
OUTPUT
cmax = integrated density for each atom for a rmax-radius sphere WARNING cg should not be modified by fourwf.
SOURCE
1238 subroutine dens_in_sph(cmax,cg,gmet,istwfk,kg_k,natom,ngfft,mpi_enreg,npw_k,& 1239 & ph1d,rmax,ucvol) 1240 1241 !Arguments ------------------------------------ 1242 !scalars 1243 integer,intent(in) :: istwfk,natom,npw_k 1244 real(dp),intent(in) :: ucvol 1245 type(MPI_type),intent(in) :: mpi_enreg 1246 !arrays 1247 integer,intent(in) :: kg_k(3,npw_k),ngfft(18) 1248 real(dp),intent(in) :: gmet(3,3) 1249 real(dp),intent(in) :: ph1d(2,(2*ngfft(1)+1+2*ngfft(2)+1+2*ngfft(3)+1)*natom) 1250 real(dp),intent(in) :: rmax(natom) 1251 real(dp),intent(inout) :: cg(2,npw_k) 1252 real(dp),intent(out) :: cmax(natom) 1253 1254 !Local variables ------------------------- 1255 !scalars 1256 integer,parameter :: tim_fourwf=0 1257 integer :: cplex,i1,i2,i3,iatom,id1,id2,id3,ifft,mgfft,n1,n2,n3,n4,n5,n6,nfft,nfftot 1258 real(dp) :: cmaxr,g1,g2,g3,norm,weight 1259 !arrays 1260 integer :: ngfft_here(18) 1261 integer,allocatable :: garr(:,:),gbound(:,:) 1262 real(dp),allocatable :: denpot(:,:,:),fofgout(:,:),fofr(:,:,:,:),gnorm(:) 1263 real(dp),allocatable :: ph3d(:,:,:),phkxred(:,:),rhog(:,:),rhor(:) 1264 real(dp),allocatable :: sphrhog(:,:) 1265 1266 ! ********************************************************************* 1267 1268 n1=ngfft(1) 1269 n2=ngfft(2) 1270 n3=ngfft(3) 1271 n4=ngfft(4) 1272 n5=ngfft(5) 1273 n6=ngfft(6) 1274 nfftot = n1*n2*n3 1275 nfft=n1*n2*n3 1276 ngfft_here(:) = ngfft(:) 1277 !fourwf doesnt work with other options for mode 0 (fft G -> r) 1278 ngfft_here(7)=111 1279 ngfft_here(8)=256 1280 mgfft=maxval(ngfft_here(1:3)) 1281 1282 call sqnorm_g(norm,istwfk,npw_k,cg,mpi_enreg%me_g0,mpi_enreg%comm_fft) 1283 1284 if (abs(one-norm) > tol6) then 1285 write(std_out,'(a,f8.5)' ) ' dens_in_sph : this state is not normalized : norm=',norm 1286 end if 1287 1288 !----------------------------------------------------------------- 1289 !inverse FFT of wavefunction to real space => density in real space 1290 !----------------------------------------------------------------- 1291 ABI_MALLOC(gbound,(2*mgfft+8,2)) 1292 call sphereboundary(gbound,istwfk,kg_k,mgfft,npw_k) 1293 1294 weight = one 1295 cplex=1 1296 ABI_MALLOC(denpot,(cplex*n4,n5,n6)) 1297 denpot(:,:,:)=zero 1298 ABI_MALLOC(fofgout,(2,npw_k)) 1299 ABI_MALLOC(fofr,(2,n4,n5,n6)) 1300 call fourwf(cplex,denpot,cg,fofgout,fofr,gbound,gbound, & 1301 & istwfk,kg_k,kg_k,mgfft,mpi_enreg,1,ngfft_here,npw_k,& 1302 & npw_k,n4,n5,n6,1,tim_fourwf,weight,weight) 1303 ABI_FREE(fofgout) 1304 ABI_FREE(fofr) 1305 ABI_FREE(gbound) 1306 1307 norm = sum(denpot(:,:,:))/nfftot 1308 if (abs(one-norm) > tol6) then 1309 write(std_out,'(a,f8.5)') ' dens_in_sph : this state is not normalized in real space : norm=',norm 1310 end if 1311 1312 !----------------------------------------------------------------- 1313 !FFT of new density: we obtain n(G) in rhog(1,:) 1314 !----------------------------------------------------------------- 1315 1316 !Change the packing of the reciprocal space density 1317 ABI_MALLOC(rhor,(nfft)) 1318 call fftpac(1,mpi_enreg,1,n1,n2,n3,n4,n5,n6,ngfft,rhor,denpot,1) 1319 1320 ABI_MALLOC(rhog,(2,nfft)) 1321 call fourdp(1,rhog,rhor,-1,mpi_enreg,nfft,1,ngfft,0) 1322 1323 ABI_FREE(rhor) 1324 ABI_FREE(denpot) 1325 1326 do ifft=1,nfft 1327 rhog(:,ifft) = rhog(:,ifft) / ucvol 1328 end do 1329 1330 !----------------------------------------------------------------- 1331 !calculate norms of G vectors 1332 !----------------------------------------------------------------- 1333 1334 ABI_MALLOC(garr,(3,nfft)) 1335 ABI_MALLOC(gnorm,(nfft)) 1336 id3=ngfft(3)/2+2 ; id2=ngfft(2)/2+2 ; id1=ngfft(1)/2+2 1337 do i3=1,n3 1338 g3=i3-(i3/(id3))*ngfft(3)-one 1339 do i2=1,n2 1340 g2=i2-(i2/(id2))*ngfft(2)-one 1341 do i1=1,n1 1342 g1=i1-(i1/(id1))*ngfft(1)-one 1343 ifft=i1+(i2-1)*n1+(i3-1)*n1*n2 1344 garr(1,ifft)=nint(g1) 1345 garr(2,ifft)=nint(g2) 1346 garr(3,ifft)=nint(g3) 1347 gnorm(ifft)=sqrt(gmet(1,1)*g1*g1 + & 1348 & two*gmet(2,1)*g2*g1 + & 1349 & two*gmet(3,1)*g3*g1 + & 1350 & gmet(2,2)*g2*g2 + & 1351 & gmet(3,2)*g3*g2 + & 1352 & gmet(3,3)*g3*g3) 1353 end do 1354 end do 1355 end do 1356 1357 !----------------------------------------------------------------- 1358 !For each atom call sphericaldens to calculate 1359 !n(G) * 1/|G|^3 * int_0^2*\pi*r_{max}*|G| 4 \pi y^2 j_0 (y) dy 1360 !for all G vectors put into array sphrhog 1361 !scalar product of phase factors with spherically convoluted density 1362 !----------------------------------------------------------------- 1363 1364 !largest mem occupation = nfft * (2(sphrog) +2*1(ph3d) +3(garr) +2(rhog) +1(gnorm)) = nfft * 10 1365 ABI_MALLOC(sphrhog,(2,nfft)) 1366 ABI_MALLOC(phkxred,(2,natom)) 1367 phkxred(1,:)=one 1368 phkxred(2,:)=zero 1369 ABI_MALLOC(ph3d,(2,nfft,1)) 1370 1371 do iatom=1,natom 1372 1373 call sphericaldens(rhog,gnorm,nfft,rmax(iatom),sphrhog) 1374 ! ----------------------------------------------------------------- 1375 ! Compute the phases for the whole set of fft vectors 1376 ! ----------------------------------------------------------------- 1377 1378 call ph1d3d(iatom,iatom,garr,natom,natom,nfft,ngfft(1),ngfft(2),ngfft(3),phkxred,ph1d,ph3d) 1379 1380 ! For the phase factors, take the compex conjugate, before evaluating the scalar product 1381 do ifft=1,nfft 1382 ph3d(2,ifft,1)=-ph3d(2,ifft,1) 1383 end do 1384 cplex=2 1385 call dotprod_v(cplex,cmaxr,nfft,1,0,ph3d,sphrhog,mpi_enreg%comm_fft) 1386 cmax(iatom) = cmaxr 1387 ! write(std_out,'(a,i4,a,es14.6,a,es12.6)' )' dens_in_sph : At ', iatom, ' has ',cmaxr, ' el.s in a sphere of rad ', rmax 1388 end do 1389 1390 ABI_FREE(rhog) 1391 ABI_FREE(gnorm) 1392 ABI_FREE(garr) 1393 ABI_FREE(sphrhog) 1394 ABI_FREE(ph3d) 1395 ABI_FREE(phkxred) 1396 1397 end subroutine dens_in_sph
m_epjdos/dos_calcnwrite [ Functions ]
[ Top ] [ m_epjdos ] [ Functions ]
NAME
dos_calcnwrite
FUNCTION
calculate DOS and write results to file(s)
INPUTS
dos_fractions= projections of wavefunctions on each angular momentum Ylm
which is the weight going into the DOS for an l-decomposed dos
dos_fractions_m= same as dos_fractions, but m-decomposed not just l-
dos_fractions_paw1= contribution to dos fractions from the PAW partial waves (phi)
dos_fractions_pawt1= contribution to dos fractions from the PAW pseudo partial waves (phi_tild)
dtset structured datatype, in particular one uses :
kptrlatt(3,3)=lattice vectors for full kpoint grid
nshiftk =number of kpoint grid shifts
pawprtdos =option to output the individual contributions to the partial DOS (0, 1 or 2)
shiftk(3,nshiftk)=kpoint shifts
usepaw =option for PAW
crystal<crystal_t>=Object defining the unit cell and its symmetries.
ebands<ebands_t>=Band structure data.
fermie=Fermi energy
fildata=name of the DOS output file
mbesslang=maximum angular momentum for Bessel function expansion
prtdosm=option for the m-contributions to the partial DOS
ndosfraction= number of types of DOS we are calculating, e.g. the number
of l channels. Could be much more general, for other types of partial DOS
paw_dos_flag= option for partial dos in PAW
comm=MPI communicator.
OUTPUT
(no explicit output)
SOURCE
334 subroutine dos_calcnwrite(dos,dtset,crystal,ebands,fildata,comm) 335 336 !Arguments ------------------------------------ 337 !scalars 338 integer,intent(in) :: comm 339 character(len=*),intent(in) :: fildata 340 type(dataset_type),intent(in) :: dtset 341 type(crystal_t),intent(in) :: crystal 342 type(ebands_t),intent(in) :: ebands 343 type(epjdos_t),intent(in) :: dos 344 345 !Local variables------------------------------- 346 !scalars 347 integer,parameter :: bcorr0=0,master=0 348 integer :: iat,iband,iene,ikpt,isppol,natsph,natsph_extra,nkpt,nsppol,i1,i2 349 integer :: nene,prtdos,unitdos,ierr,prtdosm,paw_dos_flag,mbesslang,ndosfraction 350 integer :: my_rank,nprocs,cnt,ifrac,ii 351 real(dp),parameter :: dos_max=9999.9999_dp 352 real(dp) :: buffer,deltaene,enemax,enemin,integral_DOS,max_occ 353 real(dp) :: cpu,wall,gflops 354 logical :: bigDOS,iam_master 355 character(len=10) :: tag 356 character(len=500) :: frmt,frmt_extra,msg 357 type(htetra_t) :: tetra 358 !arrays 359 integer,allocatable :: unt_atsph(:) 360 real(dp) :: list_dp(3) 361 real(dp),allocatable :: tmp_eigen(:),total_dos(:,:,:),eig_dos(:,:) 362 real(dp),allocatable :: dos_m(:,:,:),dos_paw1(:,:,:),dos_pawt1(:,:,:) 363 real(dp),allocatable :: wdt(:,:) 364 365 ! ********************************************************************* 366 367 !DEBUG 368 !write(std_out,*)' m_epjdos%dos_calcncwrite, enter ' 369 !ENDDEBUG 370 371 my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm); iam_master = (my_rank == master) 372 373 prtdosm = dos%prtdosm; paw_dos_flag = dos%paw_dos_flag 374 mbesslang = dos%mbesslang; ndosfraction = dos%ndosfraction 375 376 nkpt = dtset%nkpt; nsppol = dtset%nsppol 377 378 !m-decomposed DOS not compatible with PAW-decomposed DOS 379 if (prtdosm>=1.and.paw_dos_flag==1) then 380 msg = 'm-decomposed DOS (prtdosm>=1) not compatible with PAW-decomposed DOS (pawprtdos=1) !' 381 ABI_ERROR(msg) 382 end if 383 384 !Refuse nband different for different kpoints 385 !Note: This means we can pass ebands%eig(:,:,:) instead of eigen(mband*nkpt*nsppol) in packed form 386 do isppol=1,nsppol 387 do ikpt=1,nkpt 388 if ( dtset%nband(nkpt*(isppol-1) + ikpt) /= dtset%nband(1) ) then 389 write(std_out,*) 'tetrahedron: skip subroutine.' 390 write(std_out,*) 'nband must be the same for all kpoints' 391 write(std_out,*) 'nband=', dtset%nband 392 ABI_WARNING('tetrahedron: skip subroutine. See message above') 393 return 394 end if 395 end do 396 end do 397 398 call cwtime(cpu, wall, gflops, "start") 399 400 tetra = tetra_from_kptrlatt(crystal, dtset%kptopt, dtset%kptrlatt, dtset%nshiftk, & 401 dtset%shiftk, dtset%nkpt, dtset%kpt, comm, msg, ierr) 402 if (ierr /= 0) then 403 call tetra%free() 404 ABI_WARNING(msg) 405 return 406 end if 407 408 natsph=dtset%natsph; natsph_extra=dtset%natsph_extra 409 410 ! Master opens the DOS files. 411 if (iam_master) then 412 if (any(dtset%prtdos == [2, 5])) then 413 if (open_file(fildata, msg, newunit=unitdos, status='unknown', form='formatted', action="write") /= 0) then 414 ABI_ERROR(msg) 415 end if 416 417 else if (dtset%prtdos == 3) then 418 ! unt_atsph(0) is used for the total DOS. 419 ABI_MALLOC(unt_atsph,(0:natsph+natsph_extra)) 420 421 ! Open file for total DOS as well. 422 if (open_file(strcat(fildata, '_TOTAL'), msg, newunit=unt_atsph(0), & 423 status='unknown', form='formatted', action="write") /= 0) then 424 ABI_ERROR(msg) 425 end if 426 427 do iat=1,natsph 428 call int2char4(dtset%iatsph(iat),tag) 429 ABI_CHECK((tag(1:1)/='#'),'Bug: string length too short!') 430 if (open_file(strcat(fildata, '_AT', tag), msg, newunit=unt_atsph(iat), & 431 status='unknown', form='formatted', action="write") /= 0) then 432 ABI_ERROR(msg) 433 end if 434 end do 435 ! do extra spheres in vacuum too. Use _ATEXTRA[NUM] suffix 436 do iat=1,natsph_extra 437 call int2char4(iat,tag) 438 ABI_CHECK((tag(1:1)/='#'),'Bug: string length too short!') 439 if (open_file(strcat(fildata, '_ATEXTRA', tag), msg, newunit=unt_atsph(natsph+iat), & 440 status='unknown', form='formatted', action="write") /= 0) then 441 ABI_ERROR(msg) 442 end if 443 end do 444 end if 445 end if 446 447 ! Write the header of the DOS file, and determine the energy range and spacing 448 prtdos=dtset%prtdos 449 buffer=0.01_dp ! Size of the buffer around the min and max ranges 450 451 ! A Similar section is present is getnel. Should move all DOS stuff to m_ebands 452 ! Choose the lower and upper energies 453 enemax = maxval(ebands%eig) + buffer 454 enemin = minval(ebands%eig) - buffer 455 456 ! Extend the range to a nicer value 457 enemax=0.1_dp*ceiling(enemax*10._dp) 458 enemin=0.1_dp*floor(enemin*10._dp) 459 460 ! Choose the energy increment 461 if(abs(dtset%dosdeltae)<tol10)then 462 deltaene=0.001_dp 463 if(dtset%prtdos>=2)deltaene=0.0005_dp ! Higher resolution possible (and wanted) for tetrahedron 464 else 465 deltaene=dtset%dosdeltae 466 end if 467 nene=nint((enemax-enemin)/deltaene)+1 468 469 call xmpi_bcast(nene, master, comm, ierr) 470 if (iam_master) list_dp(1:3) = [deltaene, enemin, enemax] 471 call xmpi_bcast(list_dp, master, comm, ierr) 472 deltaene = list_dp(1); enemin = list_dp(2); enemax = list_dp(3) 473 474 if (iam_master) then 475 if (any(dtset%prtdos == [2, 5])) then 476 call dos_hdr_write(deltaene,ebands%eig,enemax,enemin,ebands%fermie,ebands%fermih,& 477 dtset%mband,dtset%nband,nene,nkpt,nsppol,dtset%occopt,prtdos,& 478 dtset%tphysel,dtset%tsmear,unitdos) 479 else if (dtset%prtdos == 3) then 480 do iat=0,natsph+natsph_extra 481 call dos_hdr_write(deltaene,ebands%eig,enemax,enemin,ebands%fermie,ebands%fermih,& 482 dtset%mband,dtset%nband,nene,nkpt,nsppol,dtset%occopt,prtdos,& 483 dtset%tphysel,dtset%tsmear,unt_atsph(iat)) 484 end do 485 end if 486 end if 487 488 ! Tetra weights 489 ABI_MALLOC(wdt, (nene, 2)) 490 491 ! Allocate arrays to store DOSes and fill with zeros. 492 ! 1--> DOS , 2--> IDOS 493 ABI_MALLOC(total_dos,(nene,ndosfraction,2)) 494 ABI_MALLOC(eig_dos, (nene, 2)) 495 496 if (paw_dos_flag==1) then 497 ABI_MALLOC(dos_paw1,(nene,ndosfraction,2)) 498 ABI_MALLOC(dos_pawt1,(nene,ndosfraction,2)) 499 end if 500 if (prtdosm>=1) then 501 ABI_MALLOC(dos_m, (nene,ndosfraction*mbesslang,2)) 502 end if 503 504 !Get maximum occupation value (2 or 1) 505 max_occ = one; if (dtset%nspinor == 1 .and. nsppol == 1) max_occ = two 506 507 !------------------------------------------------------------------- 508 !For each spin polarisation and band, interpolate band over kpoints 509 !calculate integration weights and DOS contib from 510 !------------------------------------------------------------------- 511 512 ! Workspace arrays. 513 ABI_MALLOC(tmp_eigen,(nkpt)) 514 515 cnt = 0 516 do isppol=1,nsppol 517 518 total_dos = zero; eig_dos = zero 519 if (prtdosm>=1) dos_m = zero 520 if (paw_dos_flag==1) then 521 dos_paw1 = zero; dos_pawt1 = zero 522 end if 523 524 do ikpt=1,nkpt 525 do iband=1,ebands%nband(ikpt+(isppol-1)*ebands%nkpt) 526 cnt = cnt + 1; if (mod(cnt, nprocs) /= my_rank) cycle ! Mpi parallelism. 527 528 ! Accumulate total DOS from eigenvalues (this is the **exact** total DOS) 529 tmp_eigen(:) = ebands%eig(iband, :, isppol) 530 call tetra%get_onewk(ikpt,bcorr0,nene,nkpt,tmp_eigen,enemin,enemax,max_occ,wdt) 531 wdt = wdt*ebands%wtk(ikpt) 532 eig_dos = eig_dos + wdt 533 534 ! Accumulate L-DOS. 535 do ii=1,2 536 do ifrac=1,ndosfraction 537 total_dos(:,ifrac,ii) = total_dos(:,ifrac,ii) + wdt(:,ii) * dos%fractions(ikpt,iband,isppol,ifrac) 538 end do 539 end do 540 541 if (paw_dos_flag==1) then 542 ! Accumulate L-DOS (on-site terms). 543 do ii=1,2 544 do ifrac=1,ndosfraction 545 dos_paw1(:,ifrac,ii) = dos_paw1(:,ifrac,ii) + wdt(:,ii) * dos%fractions_paw1(ikpt,iband,isppol,ifrac) 546 dos_pawt1(:,ifrac,ii) = dos_pawt1(:,ifrac,ii) + wdt(:,ii) * dos%fractions_pawt1(ikpt,iband,isppol,ifrac) 547 end do 548 end do 549 end if 550 551 if (prtdosm>=1) then 552 ! Accumulate LM-DOS. 553 do ii=1,2 554 do ifrac=1,ndosfraction*mbesslang 555 dos_m(:,ifrac,ii) = dos_m(:,ifrac,ii) + wdt(:,ii) * dos%fractions_m(ikpt, iband, isppol, ifrac) 556 end do 557 end do 558 end if 559 560 end do ! ikpt 561 end do ! iband 562 563 ! Collect results on master 564 call xmpi_sum_master(eig_dos, master, comm, ierr) 565 call xmpi_sum_master(total_dos, master, comm, ierr) 566 bigDOS=(maxval(total_dos(:,:,1))>999._dp) 567 568 if (paw_dos_flag == 1) then 569 call xmpi_sum_master(dos_paw1, master, comm, ierr) 570 call xmpi_sum_master(dos_pawt1, master, comm, ierr) 571 end if 572 if (prtdosm >= 1) call xmpi_sum_master(dos_m, master, comm, ierr) 573 574 ! Write the DOS value in the DOS file 575 ! Print the data for this energy. Note the upper limit (dos_max), to be consistent with the format. 576 ! The use of "E" format is not adequate, for portability of the self-testing procedure. 577 ! header lines depend on the type of DOS (projected etc...) which is output 578 579 if (.not. iam_master) goto 10 580 call write_extra_headers() 581 582 if (prtdos==2) then 583 ! E, DOS, IDOS 584 do iene=1,nene 585 write(unitdos, '(f11.5,1x,2(f10.4,1x))') & 586 enemin + (iene-1)*deltaene, min(total_dos(iene,:,1), dos_max), total_dos(iene,:,2) 587 end do 588 589 else if (prtdos==3) then 590 591 ! Write E, DOS, IDOS 592 do iene=1,nene 593 write(unt_atsph(0), '(f11.5,1x,2(f10.4,1x))') & 594 enemin + (iene-1)*deltaene, min(eig_dos(iene,1), dos_max), eig_dos(iene,2) 595 end do 596 597 ! E, DOS(L=1,LMAX), IDOS(L=1,LMAX) 598 ! Here we assume mpsang = 5 in the format. 599 if (paw_dos_flag/=1.or.dtset%pawprtdos==2) then 600 frmt = '(f11.5,1x,5(f9.4,1x),10x,5(f8.2,1x),10x,25(f8.2,1x))' 601 if (bigDOS) frmt = '(f11.5,1x,5(f10.4,1x),10x,5(f8.2,1x),10x,25(f8.2,1x))' 602 ! for extra atoms in vacuum need more precision 603 frmt_extra = '(f11.5,1x,5(f20.16,1x),10x,5(f20.16,1x),10x,25(f20.16,1x))' 604 605 do iat=1,natsph 606 i1 = (iat-1)*mbesslang+1; i2 = iat*mbesslang 607 if (prtdosm==0) then 608 do iene=1,nene 609 write(unt_atsph(iat), fmt=frmt) enemin + (iene-1)*deltaene, & 610 min(total_dos(iene, i1:i2, 1), dos_max), total_dos(iene, i1:i2,2) 611 end do 612 else 613 do iene=1,nene 614 write(unt_atsph(iat), fmt=frmt) enemin + (iene-1)*deltaene, & 615 min(total_dos(iene, i1:i2, 1), dos_max),& 616 total_dos(iene, i1:i2, 2),& 617 min(dos_m(iene,(iat-1)*mbesslang**2+1:iat*mbesslang**2,1), dos_max) 618 end do 619 end if 620 end do 621 622 ! Extra spheres. 623 do iat=natsph+1,natsph+natsph_extra 624 i1 = (iat-1)*mbesslang+1; i2 = iat*mbesslang 625 if (prtdosm==0) then 626 do iene=1,nene 627 write(unt_atsph(iat), fmt=frmt_extra) enemin + (iene-1)*deltaene, & 628 total_dos(iene, i1:i2, 1), & 629 total_dos(iene, i1:i2, 2) 630 end do 631 else 632 do iene=1,nene 633 write(unt_atsph(iat), fmt=frmt_extra) enemin + (iene-1)*deltaene, & 634 total_dos(iene, i1:i2, 1),& 635 total_dos(iene, i1:i2, 2),& 636 dos_m(iene,(iat-1)*mbesslang**2+1:iat*mbesslang**2, 1) 637 end do 638 end if 639 end do 640 641 else 642 frmt = '(f11.5,1x,5(f9.4,1x),3(6x,5f9.4))' 643 if (bigDOS) frmt = '(f11.5,1x,5(f10.4,1x),3(6x,5f10.4))' 644 ! for extra atom spheres in vacuum need more precision 645 frmt_extra = '(f11.5,1x,5(f20.16,1x),3(6x,5f20.16))' 646 647 do iat=1,natsph 648 i1 = iat*5-4; i2 = iat*5 649 do iene=1,nene 650 write(unt_atsph(iat), fmt=frmt) enemin + (iene-1)*deltaene, & 651 min(total_dos(iene,i1:i2,1), dos_max),& 652 min(total_dos(iene,i1:i2,1) - dos_paw1(iene,i1:i2,1) + dos_pawt1(iene,i1:i2,1), dos_max),& 653 min(dos_paw1(iene,i1:i2,1), dos_max),& 654 min(dos_pawt1(iene,i1:i2,1), dos_max) 655 end do 656 end do 657 658 ! Extra spheres. 659 do iat=natsph+1,natsph+natsph_extra 660 i1 = iat*5-4; i2 = iat*5 661 do iene=1,nene 662 write(unt_atsph(iat), fmt=frmt_extra) enemin + (iene-1)*deltaene, & 663 min(total_dos(iene,i1:i2,1), dos_max),& 664 min(total_dos(iene,i1:i2,1) - dos_paw1(iene,i1:i2,1) + dos_pawt1(iene,i1:i2,1), dos_max),& 665 min(dos_paw1(iene,i1:i2,1), dos_max),& 666 min(dos_pawt1(iene,i1:i2,1), dos_max) 667 end do 668 end do 669 end if 670 671 else if (prtdos==5)then 672 ! E, SPIN-DOS 673 frmt = '(f11.5,1x,7(f9.4,1x),10x,7(f8.2,1x))' 674 if (bigDOS) frmt = '(f11.5,1x,7(f10.4,1x),10x,7(f8.2,1x))' 675 do iene=1,nene 676 write(unitdos, fmt=frmt) enemin + (iene-1)*deltaene, min(total_dos(iene,1:7,1), dos_max), total_dos(iene,1:7,2) 677 end do 678 end if 679 680 10 continue 681 integral_DOS=sum(total_dos(nene,:,2)) 682 write(msg, '(a,es16.8)' ) ' tetrahedron : integrate to',integral_DOS 683 call wrtout(std_out,msg,'COLL') 684 end do ! isppol 685 686 ! Close files. 687 if (iam_master) then 688 if (any(prtdos == [2, 5])) then 689 close(unitdos) 690 else if (prtdos == 3) then 691 do iat=0,natsph+natsph_extra 692 close(unt_atsph(iat)) 693 end do 694 ABI_FREE(unt_atsph) 695 end if 696 end if 697 698 ABI_FREE(tmp_eigen) 699 ABI_FREE(total_dos) 700 ABI_FREE(wdt) 701 ABI_FREE(eig_dos) 702 703 if (prtdosm>=1) then 704 ABI_FREE(dos_m) 705 end if 706 707 if (paw_dos_flag==1) then 708 ABI_FREE(dos_paw1) 709 ABI_FREE(dos_pawt1) 710 end if 711 712 call tetra%free() 713 714 call cwtime(cpu,wall,gflops,"stop") 715 write(msg,'(2(a,f8.2),a)')" tetrahedron: cpu_time: ",cpu,"[s], walltime: ",wall," [s]" 716 call wrtout(std_out,msg,"PERS") 717 718 !DEBUG 719 !write(std_out,*)' m_epjdos%dos_calcncwrite, exit ' 720 !ENDDEBUG 721 722 contains 723 724 subroutine write_extra_headers() 725 726 if (nsppol==2) then 727 if(isppol==1) write(msg,'(a,16x,a)') '#','Spin-up DOS' 728 if(isppol==2) write(msg,'(2a,16x,a)') ch10,'#','Spin-dn DOS' 729 ! NB: dtset%prtdos == 5 should not happen for nsppol==2 730 731 if (any(dtset%prtdos == [2, 5])) then 732 write(unitdos, "(a)")trim(msg) 733 734 else if (dtset%prtdos == 3) then 735 do iat=0,natsph+natsph_extra 736 write(unt_atsph(iat), "(a)")trim(msg) 737 end do 738 end if 739 740 end if 741 742 if (prtdos==2) then 743 write(unitdos, '(a)' )'# energy(Ha) DOS integrated DOS' 744 745 else if (prtdos==3) then 746 747 write(unt_atsph(0), '(a)' )'# energy(Ha) DOS integrated DOS' 748 749 if (paw_dos_flag/=1.or.dtset%pawprtdos==2) then 750 do iat=1,natsph 751 write(unt_atsph(iat), '(3a,i5,a,i5,a,a,es16.6,3a)' ) & 752 '# Local DOS (columns 2-6) and integrated local DOS (columns 7-11),',ch10,& 753 '# for atom number iat=',iat,' iatom=',dtset%iatsph(iat),ch10,& 754 '# inside sphere of radius ratsph=',dtset%ratsph(dtset%typat(dtset%iatsph(iat))),' Bohr.',ch10,"#" 755 756 if (dtset%usepaw==1.and.dtset%pawprtdos==2) then 757 write(unt_atsph(iat), '(3a)' ) & 758 '# PAW: note that only all-electron on-site part has been used to compute DOS !',ch10,"#" 759 end if 760 if (bigDOS) then 761 write(msg, '(a,a)' ) & 762 '# energy(Ha) l=0 l=1 l=2 l=3 l=4',& 763 ' (integral=>) l=0 l=1 l=2 l=3 l=4' 764 else 765 write(msg, '(a,a)' ) & 766 '# energy(Ha) l=0 l=1 l=2 l=3 l=4',& 767 ' (integral=>) l=0 l=1 l=2 l=3 l=4' 768 end if 769 if (prtdosm>=1) then 770 write(msg, '(7a)' ) trim(msg),' ',& 771 ' lm=0 0',& 772 ' lm=1-1 lm=1 0 lm=1 1',& 773 ' lm=2-2 lm=2-1 lm=2 0 lm=2 1 lm=2 2',& 774 ' lm=3-3 lm=3-2 lm=3-1 lm=3 0 lm=3 1 lm=3 2 lm=3 3',& 775 ' lm=4-4 lm=4-3 lm=4-2 lm=4-1 lm=4 0 lm=4 1 lm=4 2 lm=4 3 lm=4 4' 776 end if 777 write(unt_atsph(iat), "(a)")trim(msg) 778 end do 779 else 780 do iat=1,natsph 781 write(unt_atsph(iat), '(9a,i5,a,i5,a,a,es16.6,3a)' ) & 782 '# Local DOS (columns 2-6),',ch10,& 783 '# plane-waves contrib. to DOS (columns 7-11),',ch10,& 784 '# AE on-site contrib. to DOS (columns 12-16),',ch10,& 785 '# -PS on-site contrib. to DOS (columns 17-21),',ch10,& 786 '# for atom number iat=',iat,' iatom=',dtset%iatsph(iat),ch10,& 787 '# inside sphere of radius ratsph=',dtset%ratsph(dtset%typat(dtset%iatsph(iat))),' Bohr.',ch10,"#" 788 if (bigDOS) then 789 write(msg, '(4a)' ) & 790 '#energy(Ha) l=0 l=1 l=2 l=3 l=4',& 791 ' (PW) l=0 l=1 l=2 l=3 l=4',& 792 ' (Phi) l=0 l=1 l=2 l=3 l=4',& 793 ' (tPhi) l=0 l=1 l=2 l=3 l=4' 794 else 795 write(msg, '(4a)' ) & 796 '#energy(Ha) l=0 l=1 l=2 l=3 l=4',& 797 ' (PW) l=0 l=1 l=2 l=3 l=4',& 798 ' (Phi) l=0 l=1 l=2 l=3 l=4',& 799 ' (tPhi) l=0 l=1 l=2 l=3 l=4' 800 end if 801 write(unt_atsph(iat), "(a)")trim(msg) 802 end do 803 end if 804 do iat=1,natsph_extra 805 write(unt_atsph(natsph+iat), '(3a,i5,2a,es16.6,3a)' ) & 806 '# Local DOS (columns 2-6) and integrated local DOS (columns 7-11),',ch10,& 807 '# for non-atomic sphere number iat=',iat,ch10,& 808 '# of radius ratsph=',dtset%ratsph_extra,' Bohr.',ch10,"#" 809 if (bigDOS) then 810 write(msg, '(a,a)' ) & 811 '# energy(Ha) l=0 l=1 l=2 l=3 l=4',& 812 ' (integral=>) l=0 l=1 l=2 l=3 l=4' 813 else 814 write(msg, '(a,a)' ) & 815 '# energy(Ha) l=0 l=1 l=2 l=3 l=4',& 816 ' (integral=>) l=0 l=1 l=2 l=3 l=4' 817 end if 818 if (prtdosm>=1) then 819 write(msg, '(7a)' ) trim(msg),' ',& 820 ' lm=0 0',& 821 ' lm=1-1 lm=1 0 lm=1 1',& 822 ' lm=2-2 lm=2-1 lm=2 0 lm=2 1 lm=2 2',& 823 ' lm=3-3 lm=3-2 lm=3-1 lm=3 0 lm=3 1 lm=3 2 lm=3 3',& 824 ' lm=4-4 lm=4-3 lm=4-2 lm=4-1 lm=4 0 lm=4 1 lm=4 2 lm=4 3 lm=4 4' 825 end if 826 write(unt_atsph(natsph+iat), "(a)")trim(msg) 827 end do 828 829 else if (prtdos==5) then 830 write(unitdos, '(a)' )& 831 '# energy(Ha) DOS up,up up,dn dn,up dn,dn sigma_x sigma_y sigma_z and integrated DOS components' 832 end if ! prtdos value 833 834 end subroutine write_extra_headers 835 836 end subroutine dos_calcnwrite
m_epjdos/epjdos_free [ Functions ]
[ Top ] [ m_epjdos ] [ Functions ]
NAME
epjdos_free
FUNCTION
deallocate memory
SOURCE
280 subroutine epjdos_free(self) 281 282 !Arguments ------------------------------------ 283 class(epjdos_t),intent(inout) :: self 284 285 ! ********************************************************************* 286 287 ! integer 288 ABI_SFREE(self%mlang_type) 289 290 ! real 291 ABI_SFREE(self%fractions) 292 ABI_SFREE(self%fractions_m) 293 ABI_SFREE(self%fractions_paw1) 294 ABI_SFREE(self%fractions_pawt1) 295 296 end subroutine epjdos_free
m_epjdos/epjdos_new [ Functions ]
[ Top ] [ m_epjdos ] [ Functions ]
NAME
epjdos_new
FUNCTION
Create new object from dataset input variables.
INPUTS
dtset <type(dataset_type)>=all input variables for this dataset psps <type(pseudopotential_type)>=variables related to pseudopotentials pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
SOURCE
171 type(epjdos_t) function epjdos_new(dtset, psps, pawtab) result(new) 172 173 !Arguments ------------------------------------ 174 type(dataset_type),intent(in) :: dtset 175 type(pseudopotential_type),intent(in) :: psps 176 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*psps%usepaw) 177 178 !Local variables------------------------------- 179 !scalars 180 integer :: ierr,itypat,iat 181 182 ! ********************************************************************* 183 184 !DEBUG 185 !write(std_out,*)' m_epjdos%epjdos_new, enter ' 186 !ENDDEBUG 187 188 new%nkpt = dtset%nkpt; new%mband = dtset%mband; new%nsppol = dtset%nsppol 189 190 new%prtdos = dtset%prtdos 191 new%partial_dos_flag = 0 192 if (new%prtdos==2) new%partial_dos_flag = 0 ! Standard DOS with tetra. 193 if (new%prtdos==3) new%partial_dos_flag = 1 ! L-DOS with tetra (prtdosm>0 if LM is wanted in Ylm/Slm basis). 194 if (new%prtdos==4) new%partial_dos_flag = 1 ! L-DOS with gaussian (prtdosm if LM is wanted in Ylm/Slm basis). 195 if (new%prtdos==5) new%partial_dos_flag = 2 ! Spin DOS 196 197 new%prtdosm=0 198 if (new%partial_dos_flag==1) new%prtdosm=dtset%prtdosm 199 ! paw_dos_flag= 1 if both PAW contributions are evaluated AND stored 200 new%paw_dos_flag=0 201 if (dtset%usepaw==1 .and. new%partial_dos_flag==1 .and. dtset%pawprtdos==1) new%paw_dos_flag=1 202 203 new%fatbands_flag=0 204 if (dtset%pawfatbnd>0 .and. new%prtdosm==0) new%fatbands_flag=1 205 if (new%prtdosm==1.and.dtset%pawfatbnd>0)then 206 ! because they compute quantities in real and complex harmonics respectively 207 ABI_ERROR('pawfatbnd>0 and prtdosm=1 are not compatible') 208 end if 209 210 ! mjv : initialization is needed as mbesslang is used for allocation below 211 ! NOTE: 10/5/2010 the whole of this could be looped over ndosfraction, 212 ! to store much less in memory. The DOS is accumulated in an array 213 ! and then printed to file at the end. 214 new%mbesslang = 1 215 if (new%partial_dos_flag==1 .or. new%fatbands_flag==1) then 216 217 ABI_MALLOC(new%mlang_type, (dtset%ntypat + dtset%natsph_extra)) 218 new%mlang_type = 0 219 220 ! TODO: Could use mbesslang = 4 or compute it from psps/pawtab 221 ! Increment by one (could underestimate if vloc = vlmax) 222 if (dtset%usepaw == 0) then 223 do iat=1,dtset%natsph 224 itypat = dtset%typat(dtset%iatsph(iat)) 225 new%mlang_type(itypat) = 1 + maxval(psps%indlmn(1, :, itypat)) 226 end do 227 else 228 do iat=1,dtset%natsph 229 itypat= dtset%typat(dtset%iatsph(iat)) 230 new%mlang_type(itypat) = 1 + (pawtab(itypat)%l_size - 1) / 2 231 end do 232 end if 233 234 ! Up to l=g if we have natsph_extra. 235 if (dtset%natsph_extra > 0) new%mlang_type(dtset%ntypat+1:) = 5 236 237 new%mlang_type = 5 ! This is to preserve the old implementation 238 new%mbesslang = maxval(new%mlang_type) 239 new%ndosfraction = (dtset%natsph + dtset%natsph_extra) * new%mbesslang 240 241 else if (new%partial_dos_flag == 2) then 242 new%ndosfraction = 7 243 244 else 245 new%ndosfraction = 1 246 new%mbesslang = 0 247 end if 248 249 ! Check allocations status as these arrays are not distributed and the wavefunctions are still in memory. 250 ABI_MALLOC_OR_DIE(new%fractions, (dtset%nkpt,dtset%mband,dtset%nsppol,new%ndosfraction), ierr) 251 new%fractions = zero 252 253 if (new%prtdosm>=1 .or. new%fatbands_flag==1) then 254 ABI_MALLOC_OR_DIE(new%fractions_m,(dtset%nkpt,dtset%mband,dtset%nsppol,new%ndosfraction*new%mbesslang), ierr) 255 new%fractions_m = zero 256 end if 257 258 if (dtset%usepaw==1 .and. new%partial_dos_flag==1) then 259 ABI_MALLOC_OR_DIE(new%fractions_paw1,(dtset%nkpt,dtset%mband,dtset%nsppol,new%ndosfraction), ierr) 260 ABI_MALLOC_OR_DIE(new%fractions_pawt1,(dtset%nkpt,dtset%mband,dtset%nsppol,new%ndosfraction), ierr) 261 new%fractions_paw1 = zero; new%fractions_pawt1 = zero 262 end if 263 264 !DEBUG 265 !write(std_out,*)' m_epjdos%epjdos_new, exit ' 266 !ENDDEBUG 267 268 end function epjdos_new
m_epjdos/epjdos_t [ Types ]
[ Top ] [ m_epjdos ] [ Types ]
NAME
epjdos_t
FUNCTION
Stores different contributions to the electronic DOS.
NOTES
Please contact gmatteo if you plan to change the internal implementation or add new DOSes. These results are saved in a netcdf file (see fatbands_ncwrite) so that one can read it with python and plot fatbands and PJDOSEs. The python version is able to handle the different cases (L, LM, Spin ...) but any change in the internal Abinit implementation is likely to break the python interface.
SOURCE
91 type,public :: epjdos_t 92 93 integer :: mbesslang 94 ! Max L+1 used in LM-DOS (Bessel function expansion) 95 96 integer :: ndosfraction 97 ! Defines the last dimension of the dos arrays. 98 ! Actual value depends on the other variables. 99 100 integer :: prtdos 101 ! 2 --> Standard DOS with tetra. 102 ! 3 --> L-DOS with tetra (prtdosm>0 if LM is wanted in Ylm/Slm basis). 103 ! 4 --> L-DOS with gaussian (prtdosm if LM is wanted in Ylm/Slm basis). 104 ! 5 --> Spin-DOS 105 106 integer :: prtdosm 107 ! Option for the m-contributions to the partial DOS 108 ! 1 if LM-projection is done onto complex Ylm 109 ! 2 if LM-projection is done onto real Slm 110 111 integer :: partial_dos_flag 112 113 integer :: paw_dos_flag 114 ! 1 if both PAW contributions are evaluated AND stored 115 116 !integer :: pawfatbnd 117 integer :: fatbands_flag 118 119 integer :: nkpt, mband, nsppol 120 ! Used to dimension arrays 121 122 integer,allocatable :: mlang_type(:) 123 ! mlang_type(ntypat + natsph_extra) 124 ! Max L+1 used in LM-DOS for each atom type 125 126 real(dp),allocatable :: fractions(:,:,:,:) 127 ! fractions(nkpt,mband,nsppol,ndosfraction)) 128 ! TODO: replace nsppol with nspden = 1, 2 (nsppol==2) or 4 (nspinor==2) 129 130 real(dp),allocatable :: fractions_m(:,:,:,:) 131 ! fractions_m(nkpt,mband,nsppol,ndosfraction*mbesslang) 132 133 real(dp),allocatable :: fractions_paw1(:,:,:,:) 134 ! fractions_paw1(nkpt,mband,nsppol,ndosfraction) 135 136 real(dp),allocatable :: fractions_pawt1(:,:,:,:) 137 ! fractions_pawt1(nkpt,mband,nsppol,ndosfraction)) 138 139 contains 140 141 procedure :: free => epjdos_free 142 ! Free dynamic memory 143 144 end type epjdos_t 145 146 public :: epjdos_new ! Create new object 147 148 149 public :: prtfatbands ! Print PJDOS contributions in xmgrace format. 150 public :: fatbands_ncwrite ! Write PJDOS contributions to netcdf file. 151 152 !---------------------------------------------------------------------- 153 154 contains !============================================================
m_epjdos/fatbands_ncwrite [ Functions ]
[ Top ] [ m_epjdos ] [ Functions ]
NAME
fatbands_ncwrite
FUNCTION
Write PJDOS contributions to netcdf file.
INPUTS
crystal<crystal_t>=Object defining the unit cell and its symmetries. ebands<ebands_t>=Band structure data. hdr<hdr_t>=Abinit header dtset<dtset_type>=Dataset type psps <type(pseudopotential_type)>=variables related to pseudopotentials pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data ncid=NC file handle.
OUTPUT
Only writing
SOURCE
1715 subroutine fatbands_ncwrite(dos, crystal, ebands, hdr, dtset, psps, pawtab, ncid) 1716 1717 !Arguments ------------------------------------ 1718 !scalars 1719 integer,intent(in) :: ncid 1720 type(epjdos_t),intent(in) :: dos 1721 type(crystal_t),intent(in) :: crystal 1722 type(ebands_t),intent(in) :: ebands 1723 type(hdr_type),intent(in) :: hdr 1724 type(dataset_type),intent(in) :: dtset 1725 type(pseudopotential_type),intent(in) :: psps 1726 !arrays 1727 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw) 1728 1729 #ifdef HAVE_NETCDF 1730 !Local variables------------------------------- 1731 !scalars 1732 integer :: itype,ncerr,fform 1733 real(dp) :: cpu,wall,gflops 1734 character(len=500) :: msg 1735 !arrays 1736 integer :: lmax_type(crystal%ntypat) 1737 1738 !************************************************************************* 1739 1740 ABI_CHECK(dtset%natsph > 0, "natsph <=0") 1741 call cwtime(cpu, wall, gflops, "start") 1742 1743 fform = fform_from_ext("FATBANDS.nc") 1744 ABI_CHECK(fform /= 0, "Cannot find fform associated to FATBANDS.nc") 1745 1746 ! Write header, crystal structure and band energies. 1747 NCF_CHECK(hdr%ncwrite(ncid, fform, nc_define=.True.)) 1748 NCF_CHECK(crystal%ncwrite(ncid)) 1749 NCF_CHECK(ebands_ncwrite(ebands, ncid)) 1750 1751 ! Add fatband-specific quantities 1752 ncerr = nctk_def_dims(ncid, [ & 1753 nctkdim_t("natsph", dtset%natsph), & 1754 nctkdim_t("ndosfraction", dos%ndosfraction)], defmode=.True.) 1755 NCF_CHECK(ncerr) 1756 1757 if (dos%ndosfraction*dos%mbesslang > 0) then 1758 ncerr = nctk_def_dims(ncid, [ & 1759 nctkdim_t("mbesslang", dos%mbesslang), & 1760 nctkdim_t("dos_fractions_m_lastsize", dos%ndosfraction*dos%mbesslang)]) 1761 NCF_CHECK(ncerr) 1762 end if 1763 if (dtset%natsph_extra /= 0) then 1764 NCF_CHECK(nctk_def_dims(ncid, [nctkdim_t("natsph_extra", dtset%natsph_extra)])) 1765 end if 1766 1767 ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "prtdos", "pawprtdos", "prtdosm"]) 1768 NCF_CHECK(ncerr) 1769 ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "ratsph_extra"]) 1770 NCF_CHECK(ncerr) 1771 1772 ncerr = nctk_def_arrays(ncid, [& 1773 nctkarr_t("lmax_type", "int", "number_of_atom_species"), & 1774 nctkarr_t("iatsph", "int", "natsph"), & 1775 nctkarr_t("ratsph", "dp", "number_of_atom_species"), & 1776 nctkarr_t("dos_fractions", "dp", "number_of_kpoints, max_number_of_states, number_of_spins, ndosfraction") & 1777 ]) 1778 NCF_CHECK(ncerr) 1779 1780 if (allocated(dos%fractions_m)) then 1781 ncerr = nctk_def_arrays(ncid, & 1782 nctkarr_t("dos_fractions_m", "dp", & 1783 "number_of_kpoints, max_number_of_states, number_of_spins, dos_fractions_m_lastsize")) 1784 NCF_CHECK(ncerr) 1785 end if 1786 1787 if (allocated(dos%fractions_paw1)) then 1788 ncerr = nctk_def_arrays(ncid, [& 1789 nctkarr_t("dos_fractions_paw1", "dp", "number_of_kpoints, max_number_of_states, number_of_spins, ndosfraction"), & 1790 nctkarr_t("dos_fractions_pawt1", "dp", "number_of_kpoints, max_number_of_states, number_of_spins, ndosfraction") & 1791 ]) 1792 NCF_CHECK(ncerr) 1793 end if 1794 1795 if (dtset%natsph_extra /= 0) then 1796 ncerr = nctk_def_arrays(ncid, [& 1797 nctkarr_t("xredsph_extra", "dp", "number_of_reduced_dimensions, natsph_extra") & 1798 ]) 1799 NCF_CHECK(ncerr) 1800 end if 1801 1802 ! Write variables 1803 NCF_CHECK(nctk_set_datamode(ncid)) 1804 1805 ! scalars 1806 NCF_CHECK(nf90_put_var(ncid, vid("pawprtdos"), dtset%pawprtdos)) 1807 NCF_CHECK(nf90_put_var(ncid, vid("prtdos"), dos%prtdos)) 1808 NCF_CHECK(nf90_put_var(ncid, vid("prtdosm"), dos%prtdosm)) 1809 1810 ! arrays 1811 if (dtset%usepaw == 1) then 1812 lmax_type = (pawtab(:)%l_size - 1) / 2 1813 else 1814 do itype=1,crystal%ntypat 1815 lmax_type(itype) = maxval(psps%indlmn(1, :, itype)) 1816 end do 1817 end if 1818 NCF_CHECK(nf90_put_var(ncid, vid("lmax_type"), lmax_type)) 1819 NCF_CHECK(nf90_put_var(ncid, vid("dos_fractions"), dos%fractions)) 1820 1821 if (dos%prtdos == 3) then 1822 NCF_CHECK(nf90_put_var(ncid, vid("iatsph"), dtset%iatsph(1:dtset%natsph))) 1823 NCF_CHECK(nf90_put_var(ncid, vid("ratsph"), dtset%ratsph(1:dtset%ntypat))) 1824 NCF_CHECK(nf90_put_var(ncid, vid("ratsph_extra"), dtset%ratsph_extra)) 1825 if (dtset%natsph_extra /= 0) then 1826 NCF_CHECK(nf90_put_var(ncid, vid("xredsph_extra"), dtset%xredsph_extra(:, 1:dtset%natsph_extra))) 1827 end if 1828 end if 1829 1830 if (allocated(dos%fractions_m)) then 1831 NCF_CHECK(nf90_put_var(ncid, vid("dos_fractions_m"), dos%fractions_m)) 1832 end if 1833 if (allocated(dos%fractions_paw1)) then 1834 NCF_CHECK(nf90_put_var(ncid, vid("dos_fractions_paw1"), dos%fractions_paw1)) 1835 NCF_CHECK(nf90_put_var(ncid, vid("dos_fractions_pawt1"), dos%fractions_pawt1)) 1836 end if 1837 1838 call cwtime(cpu,wall,gflops,"stop") 1839 write(msg,'(2(a,f8.2),a)')" fatbands_ncwrite: cpu_time: ",cpu,"[s], walltime: ",wall," [s]" 1840 call wrtout(std_out,msg,"PERS") 1841 #endif 1842 1843 contains 1844 integer function vid(vname) 1845 character(len=*),intent(in) :: vname 1846 vid = nctk_idname(ncid, vname) 1847 end function vid 1848 1849 end subroutine fatbands_ncwrite
m_epjdos/partial_dos_fractions [ Functions ]
[ Top ] [ m_epjdos ] [ 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...
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*my_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.
SOURCE
1900 subroutine partial_dos_fractions(dos,crystal,dtset,eigen,occ,npwarr,kg,cg,mcg,collect,mpi_enreg) 1901 1902 !Arguments ------------------------------------ 1903 !scalars 1904 integer,intent(in) :: mcg,collect 1905 type(epjdos_t),intent(inout) :: dos 1906 type(MPI_type),intent(inout) :: mpi_enreg 1907 type(dataset_type),intent(in) :: dtset 1908 type(crystal_t),intent(in) :: crystal 1909 !arrays 1910 integer,intent(in) :: kg(3,dtset%mpw*dtset%mkmem),npwarr(dtset%nkpt) 1911 real(dp),intent(in) :: cg(2,mcg) 1912 real(dp),intent(in) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol) 1913 real(dp),intent(in) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol) 1914 1915 !Local variables------------------------------- 1916 !scalars 1917 integer,parameter :: prtsphere0=0 ! do not output all the band by band details for projections. 1918 integer :: shift_b,shift_sk,iat,iatom,iband,ierr,ikpt,ilang,ioffkg,is1, is2, isoff 1919 integer :: ipw,isppol,ixint,mbess,mcg_disk,me_kpt,shift_cg 1920 integer :: mgfft,my_nspinor,n1,n2,n3,natsph_tot,npw_k,nradintmax 1921 integer :: rc_ylm,itypat,nband_k 1922 integer :: abs_shift_b 1923 integer :: unit_procar, ipauli 1924 real(dp),parameter :: bessint_delta = 0.1_dp 1925 real(dp) :: arg,bessarg,bessargmax,kpgmax,rmax 1926 real(dp) :: cpu,wall,gflops 1927 character(len=500) :: msg 1928 character(len=4) :: ikproc_str 1929 character(len=fnlen) :: filename 1930 type(jlspline_t) :: jlspl 1931 type(MPI_type) :: mpi_enreg_seq 1932 !arrays 1933 integer :: iindex(dtset%mpw),nband_tmp(1),npwarr_tmp(1) 1934 integer,allocatable :: iatsph(:),nradint(:),atindx(:),typat_extra(:),kg_k(:,:) 1935 real(dp) :: kpoint(3),spin(3),ylmgr_dum(1) 1936 real(dp) :: xfit(dtset%mpw),yfit(dtset%mpw) 1937 real(dp),allocatable :: ylm_k(:,:) 1938 real(dp),allocatable :: bess_fit(:,:,:) 1939 real(dp),allocatable :: cg_1band(:,:),cg_1kpt(:,:),kpgnorm(:),ph1d(:,:) 1940 real(dp),allocatable :: ph3d(:,:,:),ratsph(:),rint(:),sum_1ll_1atom(:,:,:) 1941 real(dp),allocatable :: sum_1lm_1atom(:,:,:) 1942 real(dp),allocatable :: cplx_1lm_1atom(:,:,:,:) 1943 real(dp),allocatable :: xred_sph(:,:),znucl_sph(:),phkxred(:,:) 1944 complex(dpc) :: cgcmat(2,2) 1945 1946 !************************************************************************* 1947 1948 !DEBUG 1949 !write(std_out, '(a)') ' m_epjdos%partial_dos_fractions : enter ' 1950 !ENDDEBUG 1951 1952 if(dtset%natsph==dtset%natom)then 1953 write(msg, '(a)') ' Compute the partial DOS fractions. This can be time-consuming. Think using natsph and iatsph.' 1954 else 1955 write(msg, '(a)') ' Compute the partial DOS fractions.' 1956 endif 1957 call wrtout(std_out,msg,'COLL') 1958 1959 ! for the moment, only support projection on angular momenta 1960 if (dos%partial_dos_flag /= 1 .and. dos%partial_dos_flag /= 2) then 1961 write(std_out,*) 'Error: partial_dos_fractions only supports angular ' 1962 write(std_out,*) ' momentum projection and spinor components for the moment. return to outscfcv' 1963 write(std_out,*) ' partial_dos = ', dos%partial_dos_flag 1964 return 1965 end if 1966 1967 ! impose all kpoints have same number of bands 1968 do isppol=1,dtset%nsppol 1969 do ikpt=1,dtset%nkpt 1970 if (dtset%nband((isppol-1)*dtset%nkpt + ikpt) /= dtset%mband) then 1971 write(std_out,*) 'Error: partial_dos_fractions wants same number of bands at each kpoint' 1972 write(std_out,*) ' isppol, ikpt = ', isppol,ikpt, dtset%nband((isppol-1)*dtset%nkpt + ikpt), dtset%mband 1973 write(std_out,*) ' all nband = ', dtset%nband 1974 return 1975 end if 1976 end do 1977 end do 1978 1979 ! Real or complex spherical harmonics? 1980 rc_ylm = 2; if (dos%prtdosm == 2) rc_ylm = 1 1981 1982 me_kpt = mpi_enreg%me_kpt 1983 my_nspinor = max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 1984 1985 n1 = dtset%ngfft(1); n2 = dtset%ngfft(2); n3 = dtset%ngfft(3) 1986 mgfft = maxval(dtset%ngfft(1:3)) 1987 1988 call cwtime(cpu, wall, gflops, "start") 1989 1990 if (dtset%prtprocar /= 0) then 1991 ! open file for each proc, and print header for master node 1992 call int2char4(me_kpt, ikproc_str) 1993 filename = 'PROCAR_'//ikproc_str 1994 if (open_file(filename, msg, newunit=unit_procar, form="formatted", action="write") /= 0) then 1995 ABI_ERROR(msg) 1996 end if 1997 if(mpi_enreg%me==0) then 1998 write (unit_procar,'(a)') 'PROCAR lm decomposed - need to merge files yourself in parallel case!!! Or use pyprocar package' 1999 if (dtset%prtprocar == 2) then 2000 write (unit_procar,'(a)') ' Requested complex output of PROCAR file (prtprocar 2)' 2001 end if 2002 write (unit_procar,'(a,I10,a,I10,a,I10,a)') '# of k-points: ', dtset%nkpt, & 2003 ' # of bands:', dtset%mband, ' # of ions:', dtset%natom, ch10 2004 end if 2005 end if 2006 2007 !############################################################## 2008 !FIRST CASE: project on angular momenta to get dos parts 2009 !############################################################## 2010 2011 if (dos%partial_dos_flag == 1) then 2012 natsph_tot = dtset%natsph + dtset%natsph_extra 2013 2014 ABI_MALLOC(iatsph, (natsph_tot)) 2015 ABI_MALLOC(typat_extra, (natsph_tot)) 2016 ABI_MALLOC(ratsph, (natsph_tot)) 2017 ABI_MALLOC(znucl_sph, (natsph_tot)) 2018 ABI_MALLOC(nradint, (natsph_tot)) 2019 ABI_MALLOC(atindx, (natsph_tot)) 2020 ABI_MALLOC(phkxred, (2,natsph_tot)) 2021 2022 ! initialize atindx 2023 do iatom=1,natsph_tot 2024 atindx(iatom) = iatom 2025 end do 2026 2027 iatsph(1:dtset%natsph) = dtset%iatsph(1:dtset%natsph) 2028 do iatom=1,dtset%natsph 2029 itypat = dtset%typat(iatsph(iatom)) 2030 typat_extra(iatom) = itypat 2031 ratsph(iatom) = dtset%ratsph(itypat) 2032 znucl_sph(iatom) = dtset%znucl(itypat) 2033 end do 2034 2035 ! fictitious atoms are declared with 2036 ! %natsph_extra, %ratsph_extra and %xredsph_extra(3, dtset%natsph_extra) 2037 ! they have atom index (natom + ii) and itype = ntypat + 1 2038 do iatom=1,dtset%natsph_extra 2039 typat_extra(iatom+dtset%natsph) = dtset%ntypat + 1 2040 ratsph(iatom+dtset%natsph) = dtset%ratsph_extra 2041 znucl_sph(iatom+dtset%natsph) = zero 2042 iatsph(iatom+dtset%natsph) = dtset%natom + iatom 2043 end do 2044 2045 ! init bessel function integral for recip_ylm max ang mom + 1 2046 ABI_MALLOC(sum_1ll_1atom,(dtset%nspinor**2,dos%mbesslang,natsph_tot)) 2047 ABI_MALLOC(sum_1lm_1atom,(dtset%nspinor**2,dos%mbesslang**2,natsph_tot)) 2048 ABI_MALLOC(cplx_1lm_1atom,(2,dtset%nspinor,dos%mbesslang**2,natsph_tot)) 2049 2050 ! Note ecuteff instead of ecut. 2051 kpgmax = sqrt(dtset%ecut * dtset%dilatmx**2) 2052 rmax = zero; bessargmax = zero; nradintmax = 0 2053 do iatom=1,natsph_tot 2054 rmax = max(rmax, ratsph(iatom)) 2055 bessarg = ratsph(iatom)*two_pi*kpgmax 2056 bessargmax = max(bessargmax, bessarg) 2057 nradint(iatom) = int (bessarg / bessint_delta) + 1 2058 nradintmax = max(nradintmax,nradint(iatom)) 2059 end do 2060 !write(std_out,*)' partial_dos_fractions: rmax=', rmax,' nradintmax: ", nradintmax 2061 ! use same number of grid points to calculate Bessel function and to do the integration later on r 2062 ! and make sure bessargmax is a multiple of bessint_delta 2063 mbess = nradintmax 2064 bessargmax = bessint_delta*mbess 2065 2066 ABI_MALLOC(rint,(nradintmax)) 2067 ABI_MALLOC(bess_fit,(dtset%mpw,nradintmax,dos%mbesslang)) 2068 2069 ! initialize general Bessel function array on uniform grid xx, from 0 to (2 \pi |k+G|_{max} |r_{max}|) 2070 jlspl = jlspline_new(mbess, bessint_delta, dos%mbesslang) 2071 2072 ABI_MALLOC(xred_sph, (3, natsph_tot)) 2073 do iatom=1,dtset%natsph 2074 xred_sph(:,iatom) = crystal%xred(:,iatsph(iatom)) 2075 end do 2076 do iatom=1,dtset%natsph_extra 2077 xred_sph(:,dtset%natsph+iatom) = dtset%xredsph_extra(:, iatom) 2078 end do 2079 2080 ABI_MALLOC(ph1d,(2,(2*n1+1 + 2*n2+1 + 2*n3+1)*natsph_tot)) 2081 call getph(atindx,natsph_tot,n1,n2,n3,ph1d,xred_sph) 2082 2083 ! Fake MPI data to be used for sequential call to initylmg. 2084 call initmpi_seq(mpi_enreg_seq) 2085 mpi_enreg_seq%my_natom = dtset%natom 2086 2087 shift_sk = 0 2088 abs_shift_b = 0 ! offset to allow for automatic update with +1 below 2089 do isppol=1,dtset%nsppol 2090 ioffkg = 0 2091 2092 do ikpt=1,dtset%nkpt 2093 nband_k = dtset%nband((isppol-1)*dtset%nkpt + ikpt) 2094 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_kpt)) then 2095 abs_shift_b = abs_shift_b + nband_k ! jump the whole kpt in the eig and occ arrays 2096 cycle 2097 end if 2098 npw_k = npwarr(ikpt) 2099 kpoint(:) = dtset%kpt(:,ikpt) 2100 2101 if (dtset%prtprocar /= 0) then 2102 write (unit_procar,'(a,I7,a,3F12.6,a,F12.6,a)') & 2103 ' k-point ', ikpt, ' : ', kpoint(:), ' weight = ', dtset%wtk(ikpt), ch10 2104 end if 2105 2106 ! make phkred for all atoms 2107 do iat=1,natsph_tot 2108 arg=two_pi*( kpoint(1)*xred_sph(1,iat) + kpoint(2)*xred_sph(2,iat) + kpoint(3)*xred_sph(3,iat) ) 2109 phkxred(1,iat)=cos(arg) 2110 phkxred(2,iat)=sin(arg) 2111 end do 2112 2113 ABI_MALLOC(kg_k, (3, npw_k)) 2114 kg_k = kg(:,ioffkg+1:ioffkg+npw_k) 2115 2116 ! kpgnorm contains norms only for kpoints used by this processor 2117 ABI_MALLOC(kpgnorm, (npw_k)) 2118 call getkpgnorm(crystal%gprimd, kpoint, kg_k, kpgnorm, npw_k) 2119 2120 ! Now get Ylm(k, G) factors: returns "real Ylms", which are real (+m) and 2121 ! imaginary (-m) parts of actual complex Ylm. Yl-m = Ylm* 2122 ! Call initylmg for a single k-point (mind mpi_enreg_seq). 2123 ABI_MALLOC(ylm_k, (npw_k, dos%mbesslang**2)) 2124 npwarr_tmp(1) = npw_k; nband_tmp(1) = nband_k 2125 call initylmg(crystal%gprimd,kg_k,kpoint,1,mpi_enreg_seq,dos%mbesslang,& 2126 npw_k,nband_tmp,1,npwarr_tmp,1,0,crystal%rprimd,ylm_k,ylmgr_dum) 2127 2128 ! get phases exp (2 pi i (k+G).x_tau) in ph3d 2129 ABI_MALLOC(ph3d,(2,npw_k,natsph_tot)) 2130 call ph1d3d(1,natsph_tot,kg_k,natsph_tot,natsph_tot,npw_k,n1,n2,n3,phkxred,ph1d,ph3d) 2131 2132 ! get Bessel function factors on array of |k+G|*r distances 2133 ! since we need many r distances and have a large number of different 2134 ! |k+G|, get j_l on uniform grid (above, in array gen_besj), 2135 ! and spline it for each kpt Gvector set. 2136 ! Note that we use the same step based on rmax, this can lead to (hopefully small) 2137 ! inaccuracies when we integrate from 0 up to rmax(iatom) 2138 do ixint=1,nradintmax 2139 rint(ixint) = (ixint-1)*rmax / (nradintmax-1) 2140 do ipw=1,npw_k 2141 xfit(ipw) = two_pi * kpgnorm(ipw) * rint(ixint) 2142 iindex(ipw) = ipw 2143 end do 2144 2145 call sort_dp(npw_k,xfit,iindex,tol14) 2146 do ilang=1,dos%mbesslang 2147 call splint(mbess, jlspl%xx, jlspl%bess_spl(:,ilang), jlspl%bess_spl_der(:,ilang), npw_k, xfit, yfit) 2148 ! re-order results for different G vectors 2149 do ipw=1,npw_k 2150 bess_fit(iindex(ipw),ixint,ilang) = yfit(ipw) 2151 end do 2152 end do 2153 end do ! ixint 2154 2155 shift_b = 0 2156 do iband=1,nband_k 2157 abs_shift_b = abs_shift_b + 1 2158 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband,isppol,me_kpt)) cycle 2159 !write(std_out,*)"in band:",iband 2160 ! TODO: eventually import eig and occ down to here - a pain, but printing outside would imply saving a huge array in memory 2161 if (dtset%prtprocar /= 0) then 2162 write (unit_procar,'(a,I7,a,F12.6,a,F12.6,a)') 'band ', iband, ' # energy ', & 2163 eigen(abs_shift_b), ' # occ. ', occ(abs_shift_b), ch10 2164 end if 2165 2166 ! Select wavefunction in cg array 2167 shift_cg = shift_sk + shift_b 2168 2169 call recip_ylm(bess_fit, cg(:,shift_cg+1:shift_cg+my_nspinor*npw_k), dtset%istwfk(ikpt),& 2170 & mpi_enreg, nradint, nradintmax, dos%mbesslang , dtset%mpw, natsph_tot, typat_extra, dos%mlang_type,& 2171 & npw_k, dtset%nspinor, ph3d, prtsphere0, rint, ratsph, rc_ylm, sum_1ll_1atom, sum_1lm_1atom, cplx_1lm_1atom,& 2172 & crystal%ucvol, ylm_k, znucl_sph) 2173 ! on exit the sum_1atom_* have both spinors counted 2174 2175 ! Accumulate 2176 do iatom=1,natsph_tot 2177 do ilang=1,dos%mbesslang 2178 dos%fractions(ikpt,iband,isppol,dos%mbesslang*(iatom-1) + ilang) & 2179 & = dos%fractions(ikpt,iband,isppol,dos%mbesslang*(iatom-1) + ilang) & 2180 & + sum_1ll_1atom(1,ilang,iatom) 2181 end do 2182 end do 2183 2184 if (dos%prtdosm /= 0) then 2185 do iatom=1,natsph_tot 2186 do ilang=1,dos%mbesslang**2 2187 dos%fractions_m(ikpt,iband,isppol,dos%mbesslang**2*(iatom-1) + ilang) & 2188 & = dos%fractions_m(ikpt,iband,isppol,dos%mbesslang**2*(iatom-1) + ilang) & 2189 & + sum_1lm_1atom(1,ilang,iatom) 2190 end do 2191 end do 2192 end if 2193 2194 ! Increment band, spinor shift 2195 !shift_b = shift_b + npw_k 2196 shift_b = shift_b + my_nspinor*npw_k 2197 2198 ! now we have both spinor components. The first option is for real values projections, eventually decomposed by Pauli spinor components 2199 if (dtset%prtprocar /= 0) then 2200 write (unit_procar,'(a)') 'ion s py pz px dxy dyz dz2 dxz dx2 tot' 2201 do ipauli= 1,dtset%nspinor**2 2202 ! Contract with Pauli matrices to get projections for this k and band, all atoms and ilang 2203 do iatom = 1, natsph_tot 2204 write (unit_procar, '(1x,I5)', advance='no') iatom 2205 do ilang=1,min(dos%mbesslang**2,9) 2206 write (unit_procar, '(F7.3)',advance='no') sum_1lm_1atom(ipauli,ilang,iatom) 2207 end do 2208 write (unit_procar, '(F7.3)',advance='yes') sum(sum_1lm_1atom(ipauli,:,iatom)) 2209 end do 2210 ! final line with sum over atoms 2211 write (unit_procar, '(a)', advance='no') 'tot ' 2212 do ilang=1,min(dos%mbesslang**2,9) 2213 write (unit_procar, '(F7.3)',advance='no') sum(sum_1lm_1atom(ipauli,ilang,:)) 2214 end do 2215 write (unit_procar, '(F7.3)',advance='yes') sum(sum_1lm_1atom(ipauli,:,:)) 2216 end do 2217 2218 ! second option is to also print the complex projection on the atomic like orbital: <psi_nk | Y_lm> in a sphere 2219 ! Two blocks are printed, first real then imaginary part 2220 if (dtset%prtprocar == 2) then 2221 write (unit_procar,'(2a)') 'ion s py pz px',& 2222 & ' dxy dyz dz2 dxz dx2 tot' 2223 do is1= 1,dtset%nspinor 2224 ! Contracted with Pauli matrices to get projections for this k and band, all atoms and ilang 2225 do iatom = 1, natsph_tot 2226 write (unit_procar, '(1x,I5)', advance='no') iatom 2227 do ilang=1,min(dos%mbesslang**2,9) 2228 write (unit_procar, '(2(F7.3,1x))',advance='no') cplx_1lm_1atom(:,is1,ilang,iatom) 2229 end do 2230 write (unit_procar, '(2(F7.3,1x))',advance='yes') sum(cplx_1lm_1atom(1,is1,:,iatom)), & 2231 & sum(cplx_1lm_1atom(2,is1,:,iatom)) 2232 end do 2233 ! final line with sum over atoms 2234 write (unit_procar, '(a)', advance='no') 'charge' 2235 do ilang=1,min(dos%mbesslang**2,9) 2236 write (unit_procar, '(F7.3,9x)',advance='no') sum(cplx_1lm_1atom(:,is1,ilang,:)**2) 2237 end do 2238 write (unit_procar, '(F7.3,9x)',advance='yes') sum(cplx_1lm_1atom(:,is1,:,:)**2) 2239 end do 2240 write (unit_procar,*) 2241 end if 2242 end if 2243 2244 end do ! band 2245 2246 ! Increment kpt and (spin, kpt) shifts 2247 ioffkg = ioffkg + npw_k 2248 shift_sk = shift_sk + nband_k*my_nspinor*npw_k 2249 2250 ABI_FREE(kg_k) 2251 ABI_FREE(kpgnorm) 2252 ABI_FREE(ylm_k) 2253 ABI_FREE(ph3d) 2254 end do ! ikpt 2255 end do ! isppol 2256 2257 ! collect = 1 ==> gather all contributions from different processors 2258 if (collect == 1) then 2259 call xmpi_sum(dos%fractions,mpi_enreg%comm_kpt,ierr) 2260 if (dos%prtdosm /= 0) call xmpi_sum(dos%fractions_m,mpi_enreg%comm_kpt,ierr) 2261 2262 ! this is now done inside recip_ylm 2263 ! if (mpi_enreg%paral_spinor == 1)then 2264 ! call xmpi_sum(dos%fractions,mpi_enreg%comm_spinor,ierr) 2265 ! if (dos%prtdosm /= 0) call xmpi_sum(dos%fractions_m,mpi_enreg%comm_spinor,ierr) 2266 ! end if 2267 end if 2268 2269 ABI_FREE(atindx) 2270 ABI_FREE(bess_fit) 2271 ABI_FREE(iatsph) 2272 ABI_FREE(typat_extra) 2273 ABI_FREE(nradint) 2274 ABI_FREE(ph1d) 2275 ABI_FREE(phkxred) 2276 ABI_FREE(ratsph) 2277 ABI_FREE(rint) 2278 ABI_FREE(sum_1ll_1atom) 2279 ABI_FREE(sum_1lm_1atom) 2280 ABI_FREE(cplx_1lm_1atom) 2281 ABI_FREE(xred_sph) 2282 ABI_FREE(znucl_sph) 2283 2284 call jlspline_free(jlspl) 2285 call destroy_mpi_enreg(mpi_enreg_seq) 2286 2287 !############################################################## 2288 !2ND CASE: project on spinors 2289 !############################################################## 2290 2291 else if (dos%partial_dos_flag == 2) then 2292 2293 if (dtset%nsppol /= 1 .or. dtset%nspinor /= 2) then 2294 ABI_WARNING("spinor projection is meaningless if nsppol==2 or nspinor/=2. Not calculating projections.") 2295 return 2296 end if 2297 if (my_nspinor /= 2) then 2298 ABI_WARNING("spinor projection with spinor parallelization is not coded. Not calculating projections.") 2299 return 2300 end if 2301 ABI_CHECK(mpi_enreg%paral_spinor == 0, "prtdos 5 does not support spinor parallelism") 2302 2303 ! FIXME: We should not allocate such a large chunk of memory! 2304 mcg_disk = dtset%mpw*my_nspinor*dtset%mband 2305 ABI_MALLOC(cg_1kpt,(2,mcg_disk)) 2306 shift_sk = 0 2307 isppol = 1 2308 2309 do ikpt=1,dtset%nkpt 2310 nband_k = dtset%nband((isppol-1)*dtset%nkpt + ikpt) 2311 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_kpt)) cycle 2312 npw_k = npwarr(ikpt) 2313 2314 cg_1kpt(:,:) = cg(:,shift_sk+1:shift_sk+mcg_disk) 2315 ABI_MALLOC(cg_1band,(2,2*npw_k)) 2316 shift_b=0 2317 do iband=1,nband_k 2318 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband,isppol,me_kpt)) cycle 2319 2320 ! Select wavefunction in cg array 2321 !shift_cg = shift_sk + shift_b 2322 cg_1band(:,:) = cg_1kpt(:,shift_b+1:shift_b+2*npw_k) 2323 call cg_getspin(cg_1band, npw_k, spin, cgcmat=cgcmat) 2324 2325 ! MG: TODO: imag part of off-diagonal terms is missing. 2326 ! I will add them later on. 2327 do is1=1,2 2328 do is2=1,2 2329 isoff = is2 + (is1-1)*2 2330 dos%fractions(ikpt,iband,isppol,isoff) = dos%fractions(ikpt,iband,isppol,isoff) & 2331 & + real(cgcmat(is1,is2)) 2332 end do 2333 end do 2334 2335 dos%fractions(ikpt,iband,isppol,5) = dos%fractions(ikpt,iband,isppol,5) + spin(1) 2336 dos%fractions(ikpt,iband,isppol,6) = dos%fractions(ikpt,iband,isppol,6) + spin(2) 2337 dos%fractions(ikpt,iband,isppol,7) = dos%fractions(ikpt,iband,isppol,7) + spin(3) 2338 2339 shift_b = shift_b + 2*npw_k 2340 end do 2341 ABI_FREE(cg_1band) 2342 shift_sk = shift_sk + nband_k*2*npw_k 2343 end do 2344 ABI_FREE(cg_1kpt) 2345 2346 ! Gather all contributions from different processors 2347 if (collect == 1) then 2348 call xmpi_sum(dos%fractions,mpi_enreg%comm_kpt,ierr) 2349 call xmpi_sum(dos%fractions,mpi_enreg%comm_bandfft,ierr) 2350 !below for future use - spinors should not be parallelized for the moment 2351 !if (mpi_enreg%paral_spinor == 1)then 2352 ! call xmpi_sum(dos%fractions,mpi_enreg%comm_spinor,ierr) 2353 !end if 2354 end if 2355 2356 else 2357 ABI_WARNING('only partial_dos==1 or 2 is coded') 2358 end if 2359 2360 !DEBUG 2361 !write(std_out,*) ' m_epjdos%partial_dos_fractions : exit ' 2362 !ENDDEBUG 2363 2364 if (dtset%prtprocar /= 0) close(unit_procar) 2365 2366 call cwtime(cpu,wall,gflops,"stop") 2367 write(msg,'(2(a,f8.2),a)')" partial_dos_fractions: cpu_time: ",cpu,"[s], walltime: ",wall," [s]" 2368 call wrtout(std_out,msg,"PERS") 2369 2370 end subroutine partial_dos_fractions
m_epjdos/partial_dos_fractions_paw [ Functions ]
[ Top ] [ m_epjdos ] [ Functions ]
NAME
partial_dos_fractions_paw
FUNCTION
Calculate PAW contributions to the partial DOS fractions (tetrahedron method)
INPUTS
cprj(natom,mcprj)= <p_lmn|Cnk> coefficients for each WF |Cnk> and each |p_lmn> non-local projector dimcprj(natom)=array of dimensions of array cprj (not ordered) dtset structured datatype, from which one uses : iatsph(nasph)=number of atoms used to project dos kpt(3,nkpt) =irreducible kpoints mband =max number of bands per k-point mkmem =number of kpoints in memory natom =number of atoms in total natsph =number of atoms ofor which the spherical decomposition must be done nband =number of electronic bands for each kpoint nkpt =number of irreducible kpoints nspinor =number of spinor components nsppol =1 or 2 spin polarization channels fatbands_flag =1 if pawfatbnd=1 or 2 mbesslang=maximum angular momentum for Bessel function expansion mpi_enreg=information about MPI parallelization prtdosm=option for the m-contributions to the partial DOS ndosfraction=natsph*mbesslang paw_dos_flag=option for the PAW contributions to the partial DOS pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data: pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
OUTPUT
=== If paw_dos_flag==1: dos%fractions_paw1(ikpt,iband,isppol,natom*mbesslang) = contribution to dos fractions from the PAW partial waves (phi) dos%fractions_pawt1(ikpt,iband,isppol,natom*mbesslang) = contribution to dos fractions from the PAW pseudo partial waves (phi_tild)
SIDE EFFECTS
dos%fractions(ikpt,iband,isppol,ndosfraction) = percentage of s, p, d..
character on each atom for the wavefunction # ikpt,iband, isppol
As input: contains only the pseudo contribution
As output: contains pseudo contribution + PAW corrections
== if prtdosm==1
dos%fractions_m(ikpt,iband,isppol,ndosfraction*mbesslang*prtdosm) =
m discretization of partial DOS fractions
SOURCE
2421 subroutine partial_dos_fractions_paw(dos,cprj,dimcprj,dtset,mcprj,mkmem,mpi_enreg,pawrad,pawtab) 2422 2423 !Arguments ------------------------------------ 2424 !scalars 2425 integer,intent(in) :: mcprj,mkmem 2426 type(MPI_type),intent(in) :: mpi_enreg 2427 type(dataset_type),intent(in) :: dtset 2428 type(epjdos_t),intent(inout) :: dos 2429 !arrays 2430 integer,intent(in) :: dimcprj(dtset%natom) 2431 type(pawcprj_type),intent(in) :: cprj(dtset%natom,mcprj) 2432 type(pawrad_type),intent(in) :: pawrad(dtset%ntypat) 2433 type(pawtab_type),target,intent(in) :: pawtab(dtset%ntypat) 2434 2435 !Local variables------------------------------- 2436 !scalars 2437 integer :: bandpp,basis_size,comm_kptband,cplex,fatbands_flag,iat,iatom,iband,ibg,ibsp 2438 integer :: ierr,ikpt,il,ilang,ilmn,iln,im,iorder_cprj,ispinor,isppol,itypat,j0lmn,j0ln 2439 integer :: jl,jlmn,jln,jm,klmn,kln,lmn_size,mbesslang,me_band,me_kpt,my_nspinor 2440 integer :: nband_cprj_k,nband_k,ndosfraction,nprocband,nproc_spkptband,paw_dos_flag,prtdosm 2441 real(dp) :: cpij,one_over_nproc 2442 !character(len=500) :: msg 2443 !arrays 2444 integer ,allocatable :: dimcprj_atsph(:) 2445 integer,ABI_CONTIGUOUS pointer :: indlmn(:,:) 2446 real(dp) :: tsec(2) 2447 real(dp),allocatable :: int1(:,:),int2(:,:),int1m2(:,:) 2448 type(pawcprj_type),allocatable :: cprj_k(:,:) 2449 2450 !****************************************************************************************** 2451 2452 DBG_ENTER("COLL") 2453 2454 ABI_CHECK(mkmem/=0,"mkmem==0 not supported anymore!") 2455 2456 fatbands_flag = dos%fatbands_flag 2457 mbesslang = dos%mbesslang 2458 prtdosm = dos%prtdosm 2459 ndosfraction = dos%ndosfraction 2460 paw_dos_flag = dos%paw_dos_flag 2461 2462 !m-decomposed DOS not compatible with PAW-decomposed DOS 2463 if(prtdosm>=1.and.paw_dos_flag==1) then 2464 ABI_ERROR('m-decomposed DOS not compatible with PAW-decomposed DOS!') 2465 end if 2466 2467 !Prepare some useful integrals 2468 basis_size=pawtab(1)%basis_size 2469 if (dtset%ntypat>1) then 2470 do itypat=1,dtset%ntypat 2471 basis_size=max(basis_size,pawtab(itypat)%basis_size) 2472 end do 2473 end if 2474 ABI_MALLOC(int1 ,(basis_size*(basis_size+1)/2,dtset%natsph)) 2475 ABI_MALLOC(int2,(basis_size*(basis_size+1)/2,dtset%natsph)) 2476 ABI_MALLOC(int1m2,(basis_size*(basis_size+1)/2,dtset%natsph)) 2477 int1=zero;int2=zero;int1m2=zero 2478 do iat=1,dtset%natsph 2479 iatom=dtset%iatsph(iat) 2480 itypat= dtset%typat(iatom) 2481 do jln=1,pawtab(itypat)%basis_size 2482 j0ln=jln*(jln-1)/2 2483 do iln=1,jln 2484 kln=j0ln+iln 2485 call simp_gen(int1(kln,iat),pawtab(itypat)%phiphj(:,kln),pawrad(itypat)) 2486 if (dtset%pawprtdos<2) then 2487 call simp_gen(int2(kln,iat),pawtab(itypat)%tphitphj(:,kln),pawrad(itypat)) 2488 int1m2(kln,iat)=int1(kln,iat)-int2(kln,iat) 2489 else 2490 int2(kln,iat)=zero;int1m2(kln,iat)=int1(kln,iat) 2491 end if 2492 end do !iln 2493 end do !jln 2494 end do 2495 2496 !Antiferro case 2497 if (dtset%nspden==2.and.dtset%nsppol==1.and.dtset%nspinor==1) then 2498 int1m2(:,:)=half*int1m2(:,:) 2499 if (paw_dos_flag==1.or.fatbands_flag==1.or.prtdosm==2) then 2500 int1(:,:)=half*int1(:,:);int2(:,:)=half*int2(:,:) 2501 end if 2502 end if 2503 2504 !Init parallelism 2505 comm_kptband=mpi_enreg%comm_kptband 2506 nproc_spkptband=xmpi_comm_size(comm_kptband)*mpi_enreg%nproc_spinor 2507 me_kpt=mpi_enreg%me_kpt ; me_band=mpi_enreg%me_band 2508 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 2509 bandpp=1;if (mpi_enreg%paral_kgb==1) bandpp=mpi_enreg%bandpp 2510 !Check if cprj is distributed over bands 2511 nprocband=my_nspinor*dtset%mband*dtset%mkmem*dtset%nsppol/mcprj 2512 if (nprocband/=mpi_enreg%nproc_band) then 2513 ABI_BUG('wrong mcprj/nproc_band!') 2514 end if 2515 2516 !Quick hack: in case of parallelism, dos_fractions have already 2517 ! been reduced over MPI processes; they have to be prepared before 2518 ! the next reduction (at the end of the following loop). 2519 if (nproc_spkptband>1) then 2520 one_over_nproc=one/real(nproc_spkptband,kind=dp) 2521 !$OMP PARALLEL DO COLLAPSE(4) & 2522 !$OMP& DEFAULT(SHARED) PRIVATE(ilang,isppol,iband,ikpt) 2523 do ilang=1,ndosfraction 2524 do isppol=1,dtset%nsppol 2525 do iband=1,dtset%mband 2526 do ikpt=1,dtset%nkpt 2527 dos%fractions(ikpt,iband,isppol,ilang)= & 2528 & one_over_nproc*dos%fractions(ikpt,iband,isppol,ilang) 2529 end do 2530 end do 2531 end do 2532 end do 2533 !$OMP END PARALLEL DO 2534 if (fatbands_flag==1.or.prtdosm==1.or.prtdosm==2) then 2535 !$OMP PARALLEL DO COLLAPSE(4) & 2536 !$OMP& DEFAULT(SHARED) PRIVATE(ilang,isppol,iband,ikpt) 2537 do ilang=1,ndosfraction*mbesslang 2538 do isppol=1,dtset%nsppol 2539 do iband=1,dtset%mband 2540 do ikpt=1,dtset%nkpt 2541 dos%fractions_m(ikpt,iband,isppol,ilang)= & 2542 & one_over_nproc*dos%fractions_m(ikpt,iband,isppol,ilang) 2543 end do 2544 end do 2545 end do 2546 end do 2547 !$OMP END PARALLEL DO 2548 end if 2549 end if 2550 2551 iorder_cprj=0 2552 2553 !LOOPS OVER SPINS,KPTS 2554 ibg=0 2555 do isppol =1,dtset%nsppol 2556 do ikpt=1,dtset%nkpt 2557 2558 nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 2559 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_kpt)) cycle 2560 2561 cplex=2;if (dtset%istwfk(ikpt)>1) cplex=1 2562 nband_cprj_k=nband_k/nprocband 2563 ABI_MALLOC(cprj_k,(dtset%natsph,my_nspinor*nband_cprj_k)) 2564 ABI_MALLOC(dimcprj_atsph,(dtset%natsph)) 2565 do iat=1,dtset%natsph 2566 dimcprj_atsph(iat)=dimcprj(dtset%iatsph(iat)) 2567 end do 2568 call pawcprj_alloc(cprj_k,0,dimcprj_atsph) 2569 ABI_FREE(dimcprj_atsph) 2570 2571 ! Extract cprj for this k-point. 2572 ibsp=0 2573 do iband=1,nband_cprj_k 2574 do ispinor=1,my_nspinor 2575 ibsp=ibsp+1 2576 do iat=1,dtset%natsph 2577 iatom=dtset%iatsph(iat) 2578 cprj_k(iat,ibsp)%cp(:,:)=cprj(iatom,ibsp+ibg)%cp(:,:) 2579 end do 2580 end do 2581 end do 2582 2583 ! LOOP OVER ATOMS (natsph_extra is not included on purpose) 2584 do iat=1,dtset%natsph 2585 iatom=dtset%iatsph(iat) 2586 itypat= dtset%typat(iatom) 2587 lmn_size=pawtab(itypat)%lmn_size 2588 indlmn => pawtab(itypat)%indlmn 2589 2590 ! LOOP OVER BANDS 2591 ibsp=0 2592 do iband=1,nband_k 2593 2594 if (mod((iband-1)/bandpp,nprocband)/=me_band) cycle 2595 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband,isppol,me_kpt)) then 2596 ibsp=ibsp+my_nspinor;cycle 2597 end if 2598 2599 do ispinor=1,my_nspinor 2600 ibsp=ibsp+1 2601 2602 do ilang=1,mbesslang 2603 2604 do jlmn=1,lmn_size 2605 jl=indlmn(1,jlmn) 2606 jm=indlmn(2,jlmn) 2607 j0lmn=jlmn*(jlmn-1)/2 2608 do ilmn=1,jlmn 2609 il=indlmn(1,ilmn) 2610 im=indlmn(2,ilmn) 2611 klmn=j0lmn+ilmn 2612 kln=pawtab(itypat)%indklmn(2,klmn) 2613 2614 if (il==ilang-1.and.jl==ilang-1.and.im==jm) then 2615 2616 cpij=cprj_k(iat,ibsp)%cp(1,ilmn)*cprj_k(iat,ibsp)%cp(1,jlmn) 2617 if (cplex==2) cpij=cpij+cprj_k(iat,ibsp)%cp(2,ilmn)*cprj_k(iat,ibsp)%cp(2,jlmn) 2618 cpij=pawtab(itypat)%dltij(klmn)*cpij 2619 2620 dos%fractions(ikpt,iband,isppol,mbesslang*(iat-1)+ilang)= & 2621 & dos%fractions(ikpt,iband,isppol,mbesslang*(iat-1)+ilang) + & 2622 & cpij*int1m2(kln,iat) 2623 if (prtdosm==1) then 2624 dos%fractions_m(ikpt,iband,isppol,mbesslang**2*(iat-1)+il**2+il+1+im)= & 2625 & dos%fractions_m(ikpt,iband,isppol,mbesslang**2*(iat-1)+il**2+il+1+im) + & 2626 & cpij*int1m2(kln,iat) 2627 end if 2628 if (fatbands_flag==1.or.prtdosm==2) then 2629 dos%fractions_m(ikpt,iband,isppol,mbesslang**2*(iat-1)+il**2+il+1+im)= & 2630 & dos%fractions_m(ikpt,iband,isppol,mbesslang**2*(iat-1)+il**2+il+1+im) + & 2631 & cpij*int1(kln,iat) 2632 end if 2633 if (paw_dos_flag==1) then 2634 dos%fractions_paw1(ikpt,iband,isppol,mbesslang*(iat-1)+ilang)= & 2635 & dos%fractions_paw1(ikpt,iband,isppol,mbesslang*(iat-1)+ilang) + & 2636 & cpij*int1(kln,iat) 2637 dos%fractions_pawt1(ikpt,iband,isppol,mbesslang*(iat-1)+ilang)= & 2638 & dos%fractions_pawt1(ikpt,iband,isppol,mbesslang*(iat-1)+ilang) + & 2639 & cpij*int2(kln,iat) 2640 end if 2641 2642 end if 2643 2644 end do !ilmn 2645 end do !jlmn 2646 2647 end do ! ilang 2648 end do ! ispinor 2649 end do ! iband 2650 2651 end do !iatom 2652 2653 if (mkmem/=0) ibg = ibg + my_nspinor*nband_cprj_k 2654 call pawcprj_free(cprj_k) 2655 ABI_FREE(cprj_k) 2656 end do ! ikpt 2657 end do ! isppol 2658 2659 ABI_FREE(int1) 2660 ABI_FREE(int2) 2661 ABI_FREE(int1m2) 2662 2663 !Reduce data in case of parallelism 2664 call timab(48,1,tsec) 2665 call xmpi_sum(dos%fractions,comm_kptband,ierr) 2666 if (prtdosm>=1.or.fatbands_flag==1) then 2667 call xmpi_sum(dos%fractions_m,comm_kptband,ierr) 2668 end if 2669 if (paw_dos_flag==1) then 2670 call xmpi_sum(dos%fractions_paw1,comm_kptband,ierr) 2671 call xmpi_sum(dos%fractions_pawt1,comm_kptband,ierr) 2672 end if 2673 call timab(48,2,tsec) 2674 if (mpi_enreg%paral_spinor==1) then 2675 call xmpi_sum(dos%fractions,mpi_enreg%comm_spinor,ierr) 2676 if (prtdosm>=1.or.fatbands_flag==1) then 2677 call xmpi_sum(dos%fractions_m,mpi_enreg%comm_spinor,ierr) 2678 end if 2679 if (paw_dos_flag==1) then 2680 call xmpi_sum(dos%fractions_paw1, mpi_enreg%comm_spinor,ierr) 2681 call xmpi_sum(dos%fractions_pawt1, mpi_enreg%comm_spinor,ierr) 2682 end if 2683 end if 2684 2685 !Averaging: A quick hack for m-decomposed LDOS: 2686 !BA: not valid in presence of spin-orbit coupling ! 2687 if (prtdosm==1.and.fatbands_flag==0) then 2688 ! if pawfatbnd is activated, one think in the cubic harmonics basis 2689 ! whereas prtdosm=1 is in the spherical harmonics basis. 2690 ! the following trick is done in order to have everything 2691 ! in the complex spherical basis (not useful for pawfatbnd if we want to 2692 ! have e.g t2g and eg d-orbitals). 2693 do iat=1,dtset%natsph 2694 do il = 0, mbesslang-1 2695 do im = 1, il 2696 dos%fractions_m(:,:,:,mbesslang**2*(iat-1)+il**2+il+1+im) = & 2697 & (dos%fractions_m(:,:,:,mbesslang**2*(iat-1)+il**2+il+1+im) + & 2698 & dos%fractions_m(:,:,:,mbesslang**2*(iat-1)+il**2+il+1-im))/2 2699 dos%fractions_m(:,:,:,mbesslang**2*(iat-1)+il**2+il+1-im) = & 2700 & dos%fractions_m(:,:,:,mbesslang**2*(iat-1)+il**2+il+1+im) 2701 end do 2702 end do 2703 end do !iatom 2704 end if 2705 2706 DBG_EXIT("COLL") 2707 2708 end subroutine partial_dos_fractions_paw
m_epjdos/prtfatbands [ Functions ]
[ Top ] [ m_epjdos ] [ Functions ]
NAME
prtfatbands
FUNCTION
Print dos_fractions_m in order to plot easily fatbands if pawfatbnd=1 1 : fatbands are resolved in L. if pawfatbnd=1 2 : fatbands are resolved in L and M.
INPUTS
dos_fractions_m(nkpt,mband,nsppol,ndosfraction*mbesslang*m_dos_flag)
= m-resolved projected dos inside PAW sphere.
dtset = Input variables
ebands<ebands_t>=Band structure data.
pawfatbnd = keyword for fatbands
mbesslang =maximum angular momentum for Bessel function expansion
m_dos_flag =option for the m-contributions to the partial DOS
ndosfraction =natsph*mbesslang
OUTPUT
(only writing)
NOTES
This routine should be called by master only
SOURCE
1488 subroutine prtfatbands(dos,dtset,ebands,fildata,pawfatbnd,pawtab) 1489 1490 !Arguments ------------------------------------ 1491 !scalars 1492 integer,intent(in) :: pawfatbnd 1493 type(epjdos_t),intent(in) :: dos 1494 type(ebands_t),intent(in) :: ebands 1495 type(dataset_type),intent(in) :: dtset 1496 character(len=fnlen),intent(in) :: fildata 1497 !arrays 1498 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat) 1499 1500 !Local variables------------------------------- 1501 !scalars 1502 integer :: iall,il,iat,natsph,inbfatbands,iband,mband,ixfat,isppol,nkpt,lmax,ll,mm 1503 integer :: ikpt,nband_k,ndosfraction,mbesslang 1504 real(dp) :: xfatband,cpu,wall,gflops 1505 character(len=1) :: tag_l,tag_1m,tag_is 1506 character(len=2) :: tag_2m 1507 character(len=10) :: tag_il,tag_at,tag_grace 1508 character(len=1500) :: msg 1509 character(len=fnlen) :: tmpfil 1510 type(atomdata_t) :: atom 1511 !arrays 1512 integer,allocatable :: unitfatbands_arr(:,:) 1513 real(dp),allocatable :: eigenvalues(:,:,:) 1514 character(len=2) :: symbol 1515 1516 !************************************************************************* 1517 1518 DBG_ENTER("COLL") 1519 1520 ndosfraction = dos%ndosfraction; mbesslang = dos%mbesslang 1521 1522 if(dos%prtdosm.ne.0) then 1523 write(msg,'(3a)')& 1524 'm decomposed dos is activated',ch10, & 1525 'Action: deactivate it with prtdosm=0 !' 1526 ABI_ERROR(msg) 1527 end if 1528 1529 if(dtset%nspinor==2) then 1530 ABI_WARNING("Fatbands are not yet available in the case nspinor==2!") 1531 end if 1532 1533 ABI_CHECK(allocated(dos%fractions_m), "dos%fractions_m is not allocated!") 1534 1535 natsph=dtset%natsph 1536 nkpt=dtset%nkpt 1537 mband=dtset%mband 1538 1539 if(natsph>1000) then 1540 write(msg,'(3a)')& 1541 'Too big number of fat bands!',ch10, & 1542 'Action: decrease natsph in input file !' 1543 ABI_ERROR(msg) 1544 end if 1545 1546 !-------------- PRINTING IN LOG 1547 call cwtime(cpu, wall, gflops, "start") 1548 write(msg,'(a,a,a,a,i5,a,a,1000i5)') ch10," ***** Print of fatbands activated ****** ",ch10,& 1549 " Number of atom: natsph = ",natsph,ch10, & 1550 " atoms are = ",(dtset%iatsph(iat),iat=1,natsph) 1551 call wrtout(std_out,msg,'COLL') 1552 call wrtout(ab_out,msg,'COLL') 1553 iall=0;inbfatbands=0 1554 1555 if(pawfatbnd==1) then 1556 inbfatbands=mbesslang-1 1557 write(msg,'(3a)')" (fatbands are in eV and are given for each value of L)",ch10 1558 else if(pawfatbnd==2) then 1559 write(msg,'(3a)')" (fatbands are in eV and are given for each value of L and M)",ch10 1560 inbfatbands=(mbesslang-1)**2 1561 end if 1562 call wrtout(std_out,msg,'COLL') 1563 call wrtout(ab_out,msg,'COLL') 1564 1565 write(msg,'(a,e12.5,a,e12.5,a)') " Fermi energy is ",ebands%fermie*Ha_eV," eV = ",ebands%fermie," Ha" 1566 call wrtout(std_out,msg,'COLL') 1567 1568 !-------------- OPEN AND NAME FILES FOR FATBANDS 1569 ABI_MALLOC(unitfatbands_arr,(natsph*inbfatbands,dtset%nsppol)) 1570 unitfatbands_arr = -3 1571 1572 do iat=1,natsph 1573 lmax=(pawtab(dtset%typat(dtset%iatsph(iat)))%l_size-1)/2 1574 call int2char4(dtset%iatsph(iat),tag_at) 1575 ABI_CHECK((tag_at(1:1)/='#'),'Bug: string length too short!') 1576 call atomdata_from_znucl(atom,dtset%znucl(dtset%typat(dtset%iatsph(iat)))) 1577 symbol = atom%symbol 1578 do il=1,inbfatbands 1579 iall=iall+1 1580 ll=int(sqrt(float(il-1))) ! compute l 1581 if(ll.le.lmax) then ! print only angular momentum included in the PAW data 1582 do isppol=1,dtset%nsppol 1583 write(tag_is,'(i1)')isppol 1584 if(pawfatbnd==1) then 1585 call int2char4(il-1,tag_il) 1586 ABI_CHECK((tag_il(1:1)/='#'),'Bug: string length too short!') 1587 tmpfil = trim(fildata)//'_at'//trim(tag_at)//'_'//trim(adjustl(symbol))//'_is'//tag_is//'_l'//trim(tag_il) 1588 else if (pawfatbnd==2) then 1589 write(tag_l,'(i1)') ll 1590 mm=il-(ll**2+ll+1) ! compute m 1591 if(mm<0) write(tag_2m,'(i2)') mm 1592 if(mm>=0) write(tag_1m,'(i1)') mm 1593 if(mm<0) tmpfil = trim(fildata)// & 1594 & '_at'//trim(tag_at)//'_'//trim(adjustl(symbol))//'_is'//tag_is//'_l'//tag_l//'_m'//tag_2m 1595 if(mm>=0) tmpfil = trim(fildata)// & 1596 & '_at'//trim(tag_at)//'_'//trim(adjustl(symbol))//'_is'//tag_is//'_l'//tag_l//'_m+'//tag_1m 1597 end if 1598 !unitfatbands_arr(iall,isppol)=tmp_unit+100+iall-1+(natsph*inbfatbands)*(isppol-1) 1599 !open (unit=unitfatbands_arr(iall,isppol),file=trim(tmpfil),status='unknown',form='formatted') 1600 if (open_file(tmpfil, msg, newunit=unitfatbands_arr(iall,isppol), status='unknown',form='formatted') /= 0) then 1601 ABI_ERROR(msg) 1602 end if 1603 1604 write(msg,'(a,a,a,i4)') 'opened file : ', trim(tmpfil), ' unit', unitfatbands_arr(iall,isppol) 1605 call wrtout(std_out,msg,'COLL') 1606 write(msg,'(9a)') "# ",ch10,"# ABINIT package : FATBAND file ", ch10,& 1607 "# It contains, for each band: the eigenvalues in eV (and the character of the band) as a function of the k-point",& 1608 ch10,"# This file can be read with xmgrace (http://plasma-gate.weizmann.ac.il/Grace/) ",ch10,"# " 1609 write(unitfatbands_arr(iall,isppol), "(a)")trim(msg) 1610 do iband=1,mband 1611 call int2char4(iband-1,tag_grace) 1612 ABI_CHECK((tag_grace(1:1)/='#'),'Bug: string length too short!') 1613 write(msg,'(16a)') ch10,"@ s",trim(tag_grace)," line color 1",& 1614 ch10,"@ s",trim(tag_grace)," errorbar color 2",& 1615 ch10,"@ s",trim(tag_grace)," errorbar riser linewidth 5.0", & 1616 ch10,"@ s",trim(tag_grace)," errorbar linestyle 0" 1617 write(unitfatbands_arr(iall,isppol), "(a)")trim(msg) 1618 end do !iband 1619 write(unitfatbands_arr(iall,isppol), '(a,a)') ch10,'@type xydy' 1620 end do ! isppol 1621 end if ! ll=<lmax 1622 end do ! il 1623 end do ! iat 1624 1625 if(iall.ne.(natsph*inbfatbands)) then 1626 ABI_ERROR("error1 ") 1627 end if 1628 1629 1630 1631 !-------------- WRITE FATBANDS IN FILES 1632 if (pawfatbnd>0) then 1633 ! Store eigenvalues with nkpt as first dimension for efficiency reasons 1634 ABI_MALLOC(eigenvalues,(nkpt,mband,dtset%nsppol)) 1635 do isppol=1,dtset%nsppol 1636 do ikpt=1,nkpt 1637 nband_k=dtset%nband(ikpt+(isppol-1)*nkpt) 1638 do iband=1,mband 1639 eigenvalues(ikpt,iband,isppol)= ebands%eig(iband, ikpt, isppol) - ebands%fermie 1640 end do 1641 end do 1642 end do 1643 iall=0 1644 do iat=1,natsph 1645 lmax=(pawtab(dtset%typat(dtset%iatsph(iat)))%l_size-1)/2 1646 do il=1,inbfatbands 1647 iall=iall+1 1648 ll=int(sqrt(float(il-1))) 1649 if(ll.le.lmax) then 1650 do isppol=1,dtset%nsppol 1651 do iband=1,mband 1652 write(unitfatbands_arr(iall,isppol),'(a,a,i8)') ch10,"# BAND number :",iband 1653 do ikpt=1,nkpt 1654 if(pawfatbnd==1) then 1655 xfatband=0.d0 1656 do ixfat=(il-1)**2+1,il**2 1657 xfatband=xfatband+dos%fractions_m(ikpt,iband,isppol,(iat-1)*mbesslang**2+ixfat) 1658 end do ! ixfat 1659 else if (pawfatbnd==2) then 1660 xfatband=dos%fractions_m(ikpt,iband,isppol,(iat-1)*mbesslang**2+il) 1661 end if 1662 write(unitfatbands_arr(iall,isppol),'(i5,e20.5,e20.5)')& 1663 ikpt-1,eigenvalues(ikpt,iband,isppol)*Ha_eV,xfatband 1664 end do ! ikpt 1665 end do !iband 1666 write(unitfatbands_arr(iall,isppol),'(a)') '&' 1667 !close(unitfatbands_arr(iall,isppol)) 1668 end do !isppol 1669 end if 1670 end do ! il 1671 end do ! iat 1672 ABI_FREE(eigenvalues) 1673 end if 1674 1675 do isppol=1,size(unitfatbands_arr, dim=2) 1676 do iat=1,size(unitfatbands_arr, dim=1) 1677 if (unitfatbands_arr(iat, isppol) /= -3) close (unitfatbands_arr(iat, isppol)) 1678 end do 1679 end do 1680 1681 ABI_FREE(unitfatbands_arr) 1682 1683 call cwtime(cpu,wall,gflops,"stop") 1684 write(msg,'(2(a,f8.2),a)')" prtfatbands: cpu_time: ",cpu,"[s], walltime: ",wall," [s]" 1685 call wrtout(std_out,msg,"PERS") 1686 1687 DBG_EXIT("COLL") 1688 1689 end subroutine prtfatbands
m_epjdos/recip_ylm [ Functions ]
[ Top ] [ m_epjdos ] [ Functions ]
NAME
recip_ylm
FUNCTION
Project input wavefunctions in reciprocal space on to Ylm (real or complex harmonics depending on rc_ylm).
INPUTS
bess_fit(mpw,nradintmax,ll) = Bessel functions for L, splined with arguments $2 \pi |k+G| \Delta r$, for all G vectors in sphere and all points on radial grid. cg_1band(2,npw_k)=wavefunction in recip space (note that nspinor is missing, see Notes). istwfk= storage mode of cg_1band nradint(natsph)=number of points on radial real-space grid for a given atom. nradintmax=dimension of rint array. me_g0=1 if this processor has G=0, 0 otherwise mlang=maximum angular momentum in Bessel functions. mpw=Maximum number of planewaves. Used to dimension bess_fit natsph=number of atoms around which ang mom projection has to be done typat_extra(natsph)=Type of each atom. ntypat + 1 if empty sphere mlang_type(ntypat + natsph_extra)=Max L+1 for each atom type npw_k=number of plane waves for this kpt nspinor=number of spinor components ph3d(2,npw_k,natsph)=3-dim structure factors, for each atom and plane wave. prtsphere= if 1, print a complete analysis of the angular momenta in atomic spheres rint(nradintmax) = points on radial real-space grid for integration rmax(natsph)=maximum radius for real space integration sphere rc_ylm= 1 for real spherical harmonics. 2 for complex spherical harmonics, ucvol=unit cell volume in bohr**3. ylm_k(npw_k,mlang**2)=real spherical harmonics for each G and LM. znucl_sph(natsph)=gives the nuclear number for each type of atom
OUTPUT
sum_1ll_1atom(mlang,natsph)= projected scalars for each atom and ang. mom. sum_1lm_1atom(mlang*mlang,natsph)= projected scalars for each atom and LM component. cplx_1lm_1atom(2,dtset%nspinor**2,dos%mbesslang**2,natsph_tot) = complex projection of wave function on atomic like orbital
NOTES
* ph3d atoms are ordered with natsph and must be provided by the caller in the correct order!
* spinor components are not treated here. This facilitates the implementation of spinor parallelism
because the caller can easily call the routine inside a loop over spinors and then sum the
different contributions outside the loop thus reducing the number of MPI calls.
SOURCE
886 subroutine recip_ylm (bess_fit, cg_1band, istwfk, mpi_enreg, nradint, nradintmax, mlang,& 887 & mpw, natsph, typat_extra, mlang_type, npw_k, nspinor, ph3d, prtsphere, rint, rmax,& 888 & rc_ylm, sum_1ll_1atom, sum_1lm_1atom, cplx_1lm_1atom, ucvol, ylm_k, znucl_sph) 889 890 !Arguments ------------------------------------ 891 !scalars 892 integer,intent(in) :: istwfk,mlang,mpw,natsph,npw_k,nradintmax 893 integer,intent(in) :: nspinor 894 integer,intent(in) :: prtsphere,rc_ylm 895 real(dp),intent(in) :: ucvol 896 !arrays 897 integer,intent(in) :: nradint(natsph),typat_extra(natsph),mlang_type(:) 898 real(dp),intent(in) :: bess_fit(mpw,nradintmax,mlang),cg_1band(:,:) !(2,my_nspinor*npw_k) 899 real(dp),intent(in) :: ph3d(2,npw_k,natsph),rint(nradintmax) 900 real(dp),intent(in) :: rmax(natsph),ylm_k(npw_k,mlang*mlang) 901 real(dp),intent(in) :: znucl_sph(natsph) 902 type(MPI_type),intent(in) :: mpi_enreg 903 real(dp),intent(out) :: sum_1ll_1atom(nspinor**2, mlang, natsph) 904 real(dp),intent(out) :: sum_1lm_1atom(nspinor**2, mlang*mlang,natsph) 905 real(dp),intent(out) :: cplx_1lm_1atom(2,nspinor, mlang*mlang,natsph) 906 907 !Local variables------------------------------- 908 !scalars 909 integer :: ilm,iat,ipw,ixint,ll,mm,il,jlm,ierr,lm_size,itypat 910 integer :: ispinor 911 integer :: ipauli, is, isp 912 integer :: my_nspinor 913 real(dp),parameter :: invsqrt2=one/sqrt2 914 real(dp) :: sum_all, dr, fact 915 type(atomdata_t) :: atom 916 character(len=500) :: msg 917 !arrays 918 integer :: ilang(mlang**2) 919 real(dp) :: c1(2),c2(2) 920 real(dp) :: sum_1atom(natsph),sum_1ll(mlang),sum_1lm(mlang**2) 921 real(dp) :: func(nradintmax) 922 real(dp) :: func_cplx(nradintmax,2) 923 complex(dpc) :: vect(npw_k) 924 complex(dpc),allocatable :: tmppsia(:,:),tmppsim(:,:),dotc(:) 925 integer, allocatable :: ispinors(:) 926 complex(dpc),allocatable :: values(:,:,:,:) 927 928 ! ************************************************************************* 929 930 ! Workspace array (used to reduce the number of MPI communications) 931 ! One could reduce a bit the memory requirement by using non-blocking operations ... 932 ABI_MALLOC_OR_DIE(values, (nradintmax, nspinor, mlang**2, natsph), ierr) 933 values = czero 934 935 my_nspinor = max(1,nspinor/mpi_enreg%nproc_spinor) 936 ABI_MALLOC(tmppsia, (npw_k,my_nspinor)) 937 ABI_MALLOC(tmppsim, (npw_k,my_nspinor)) 938 ABI_MALLOC(dotc, (my_nspinor)) 939 ABI_MALLOC(ispinors, (my_nspinor)) 940 if (my_nspinor == 2) then 941 ispinors(1) = 1 942 ispinors(2) = 2 943 else 944 ispinors(1) = mpi_enreg%me_spinor+1 945 end if 946 947 sum_1lm_1atom = zero 948 cplx_1lm_1atom = zero 949 950 do ll=0,mlang-1 951 do mm=-ll,ll 952 ilm = (ll+1)**2-ll+mm 953 ilang(ilm) = ll+1 954 end do 955 end do 956 957 ! Big loop on all atoms 958 do iat=1,natsph 959 itypat = typat_extra(iat) 960 lm_size = mlang_type(itypat) ** 2 961 dr = rmax(iat) / (nradint(iat)-1) 962 963 ! u(G) e^{i(k+G).Ra} 964 ! tmppsia = Temporary array for part which depends only on iat 965 do ispinor=1,my_nspinor 966 do ipw=1,npw_k 967 tmppsia(ipw,ispinor) = dcmplx(cg_1band(1,ipw+(ispinor-1)*npw_k),cg_1band(2,ipw+(ispinor-1)*npw_k)) & 968 & * dcmplx(ph3d(1,ipw,iat), ph3d(2,ipw,iat)) 969 end do 970 end do 971 972 ! tmppsim = temporary arrays for part of psi which does not depend on ixint = tmppsia * ylm. 973 ! could remove this intermediate array to save memory... 974 ! u(G) Y_LM^*(k+G) e^{i(k+G).Ra} 975 ! Take into account the fact that ylm_k are REAL spherical harmonics, see initylmg.f 976 ! For time-reversal states, detailed treatment show that only the real or imaginary 977 ! part of tmppsia is needed here, depending on l being even or odd: only one of the coef is 1, the other 0 978 do ilm=1,lm_size 979 il = ilang(ilm) 980 ll = ilang(ilm) - 1 981 mm = ilm - (ll+1)**2 + ll 982 983 select case (rc_ylm) 984 case (1) 985 ! to get PDOS for real spherical harmonics, simply multiply here by ylm_k 986 do ispinor=1,my_nspinor 987 do ipw=1,npw_k 988 tmppsim(ipw,ispinor) = tmppsia(ipw,ispinor) * ylm_k(ipw,ilm) 989 end do 990 end do 991 992 ! Handle time-reversal 993 ! TODO: check if time reversal with spinors is special and may need some treatment. 994 ! normally SOC will simply impose istwfk 1 in appropriate cases (kptopt 4). 995 if (istwfk /= 1) then 996 if (mod(ll, 2) == 0) then 997 tmppsim(:,:) = dcmplx(real(tmppsim(:,:)),zero) 998 else 999 tmppsim(:,:) = dcmplx(aimag(tmppsim(:,:)),zero) 1000 end if 1001 ! f2008 version: 1002 ! if (mod(ll, 2) == 0) then 1003 ! tmppsim(:,:)%im = zero 1004 ! else 1005 ! tmppsim(:,:)%re = tmppsim(:,:)%im 1006 ! tmppsim(:,:)%im = zero 1007 ! end if 1008 end if 1009 1010 case (2) 1011 ! to get PDOS for complex spherical harmonics, build linear combination of real ylm_k 1012 jlm = (ll+1)**2-ll-mm ! index of (l, -m) 1013 if (mm == 0) then 1014 vect(:) = dcmplx(ylm_k(1:npw_k,ilm),zero) 1015 else if (mm > 0) then 1016 !vect(1,:) = invsqrt2 * ylm_k(1:npw_k,ilm) * (-1)**mm 1017 !vect(2,:) = +invsqrt2 * ylm_k(1:npw_k,jlm) * (-1)**mm 1018 c1 = sy(ll, mm, mm) 1019 c2 = sy(ll,-mm, mm) 1020 vect(:) = dcmplx(c1(1) * ylm_k(1:npw_k,ilm) + c2(1) * ylm_k(1:npw_k,jlm), & 1021 c1(2) * ylm_k(1:npw_k,ilm) + c2(2) * ylm_k(1:npw_k,jlm)) 1022 1023 else if (mm < 0) then 1024 !vect(1,:) = invsqrt2 * ylm_k(1:npw_k,jlm) !* (-1)**mm 1025 !vect(2,:) = -invsqrt2 * ylm_k(1:npw_k,ilm) !* (-1)**mm 1026 c1 = sy(ll, mm, mm) 1027 c2 = sy(ll,-mm, mm) 1028 vect(:) = dcmplx(c1(1) * ylm_k(1:npw_k,ilm) + c2(1) * ylm_k(1:npw_k,jlm),& 1029 c1(2) * ylm_k(1:npw_k,ilm) + c2(2) * ylm_k(1:npw_k,jlm)) 1030 end if 1031 vect(:) = dcmplx(real(vect(:)), -aimag(vect(:))) 1032 !vect(:)%im = -vect(:)%im 1033 1034 if (istwfk == 1) then 1035 do ispinor=1,my_nspinor 1036 do ipw=1,npw_k 1037 tmppsim(ipw,ispinor) = tmppsia(ipw,ispinor) * vect(ipw) 1038 end do 1039 end do 1040 else 1041 ! Handle time-reversal 1042 if (mod(ll, 2) == 0) then 1043 do ispinor=1,my_nspinor 1044 do ipw=1,npw_k 1045 tmppsim(ipw,ispinor) = real(tmppsia(ipw,ispinor)) * vect(ipw) 1046 end do 1047 end do 1048 else 1049 do ispinor=1,my_nspinor 1050 do ipw=1,npw_k 1051 tmppsim(ipw,ispinor) = aimag(tmppsia(ipw,ispinor)) * vect(ipw) 1052 end do 1053 end do 1054 end if 1055 end if 1056 1057 case default 1058 ABI_ERROR("Wrong value for rc_ylm") 1059 end select 1060 1061 ! Compute integral $ \int_0^{rc} dr r**2 ||\sum_G u(G) Y_LM^*(k+G) e^{i(k+G).Ra} j_L(|k+G| r)||**2 $ 1062 ! or more general spinor case integral 1063 ! $ \int_0^{rc} dr r**2 dotc^*_s \sigma^x_{ss'} dotc_{s'} 1064 ! where dotc_s = \sum_G u_s (G) Y_LM^*(k+G) e^{i(k+G).Ra} j_L(|k+G| r) 1065 do ixint=1,nradint(iat) 1066 dotc = czero 1067 do ispinor=1, my_nspinor 1068 do ipw=1,npw_k 1069 dotc(ispinor) = dotc(ispinor) + bess_fit(ipw, ixint, il) * tmppsim(ipw, ispinor) 1070 end do 1071 end do 1072 if (istwfk /= 1) then 1073 dotc = two * dotc 1074 if (istwfk == 2 .and. mpi_enreg%me_g0 == 1) then 1075 dotc(:) = dotc(:) - bess_fit(1, ixint, il) * tmppsim(1, :) 1076 end if 1077 end if 1078 1079 ! Store results to reduce number of xmpi_sum calls if MPI 1080 do ispinor=1, my_nspinor 1081 values(ixint, ispinors(ispinor), ilm, iat) = dotc(ispinor) 1082 end do 1083 end do ! ixint 1084 1085 end do ! ilm 1086 end do ! iat 1087 1088 ABI_FREE(tmppsia) 1089 ABI_FREE(tmppsim) 1090 ABI_FREE(dotc) 1091 ABI_FREE(ispinors) 1092 1093 ! Collect results in comm_pw (data are distributed over plane waves) 1094 call xmpi_sum(values, mpi_enreg%comm_bandfft, ierr) 1095 ! ! Collect results in mpi_enreg%comm_spinor (data are distributed over spinor components) 1096 call xmpi_sum(values, mpi_enreg%comm_spinor, ierr) 1097 1098 ! Multiply by r**2 and take norm, integrate 1099 do iat=1,natsph 1100 itypat = typat_extra(iat) 1101 lm_size = mlang_type(itypat) ** 2 1102 do ilm=1,lm_size 1103 1104 do ipauli=0,nspinor**2-1 1105 do ixint=1,nradint(iat) 1106 func(ixint) = zero 1107 do is=1,nspinor 1108 do isp=1,nspinor 1109 func(ixint) = func(ixint) + real(conjg(values(ixint, is, ilm, iat)) * pauli_mat(is,isp,ipauli)*& 1110 & values(ixint, isp, ilm, iat)) 1111 end do 1112 end do 1113 func(ixint) = rint(ixint)**2 * func(ixint) 1114 end do 1115 ! Here I should treat the case in which the last point /= rcut 1116 ! NB: indexing is from 1 not 0 for spin matrix components 1117 sum_1lm_1atom (ipauli+1, ilm, iat) = simpson(dr, func(1:nradint(iat))) 1118 end do ! ipauli 1119 1120 do is = 1, nspinor 1121 func_cplx(:,1) = real(values(:, is, ilm, iat)) 1122 func_cplx(:,2) = aimag(values(:, is, ilm, iat)) 1123 do ixint=1,nradint(iat) 1124 func_cplx(ixint,:) = rint(ixint)**2 * func_cplx(ixint,:) 1125 end do 1126 cplx_1lm_1atom(1, is, ilm, iat) = simpson(dr, func_cplx(1:nradint(iat),1)) 1127 cplx_1lm_1atom(2, is, ilm, iat) = simpson(dr, func_cplx(1:nradint(iat),2)) 1128 end do 1129 1130 end do ! ilm 1131 end do ! iat 1132 1133 ! Normalize with unit cell volume and include 4pi term coming from Rayleigh expansion. 1134 fact = four_pi**2 / ucvol 1135 sum_1lm_1atom = fact * sum_1lm_1atom 1136 cplx_1lm_1atom = fact * cplx_1lm_1atom 1137 1138 ! sum up the m-independent fractions 1139 sum_1ll_1atom = zero 1140 do iat=1,natsph 1141 itypat = typat_extra(iat) 1142 lm_size = mlang_type(itypat) ** 2 1143 do ilm=1,lm_size 1144 il = ilang(ilm) 1145 sum_1ll_1atom(:,il, iat) = sum_1ll_1atom(:,il, iat) + sum_1lm_1atom(:,ilm, iat) 1146 end do 1147 end do 1148 1149 ABI_FREE(values) 1150 1151 ! Output 1152 if (prtsphere == 1) then 1153 sum_1ll = zero 1154 sum_1lm = zero 1155 sum_1atom = zero 1156 do iat=1,natsph 1157 sum_1atom(iat) = sum(sum_1lm_1atom(1,:,iat)) 1158 sum_1ll(:)=sum_1ll(:)+sum_1ll_1atom(1,:,iat) 1159 sum_1lm(:)=sum_1lm(:)+sum_1lm_1atom(1,:,iat) 1160 end do 1161 sum_all = sum(sum_1atom) 1162 1163 if (rc_ylm == 1) msg = " Angular analysis (real spherical harmonics)" 1164 if (rc_ylm == 2) msg = " Angular analysis (complex spherical harmonics)" 1165 call wrtout(std_out, msg) 1166 do iat=1,natsph 1167 call atomdata_from_znucl(atom, znucl_sph(iat)) 1168 call wrtout(std_out, " ") 1169 write(msg,'(a,i3,a,a,a,f10.6)' )' Atom # ',iat, ' is ', atom%symbol,', in-sphere charge =',sum_1atom(iat) 1170 call wrtout(std_out, msg) 1171 do ll=0,mlang-1 1172 write(msg,'(a,i1,a,f9.6,a,9f6.3)' )& 1173 ' l=',ll,', charge=',sum_1ll_1atom(1,ll+1,iat),& 1174 ', m=-l,l splitting:',sum_1lm_1atom(1,1+ll**2:(ll+1)**2,iat) 1175 call wrtout(std_out, msg) 1176 end do ! ll 1177 end do ! iat 1178 write(msg,'(a,a)') ch10,' Sum of angular contributions for all atomic spheres ' 1179 call wrtout(std_out, msg) 1180 do ll=0,mlang-1 1181 write(msg,'(a,i1,a,f9.6,a,f9.6)' )& 1182 ' l=',ll,', charge =',sum_1ll(ll+1),' proportion =',sum_1ll(ll+1)/sum_all 1183 call wrtout(std_out, msg) 1184 end do 1185 write(msg,'(a,a,f10.6)' ) ch10,' Total over all atoms and l=0 to 4 :',sum_all 1186 call wrtout(std_out, msg) 1187 call wrtout(std_out, " ") 1188 end if 1189 1190 contains 1191 1192 function sy(ll, mm, mp) 1193 use m_paw_sphharm, only : ys 1194 ! Computes the matrix element <Slm|Ylm'> 1195 integer,intent(in) :: ll,mm, mp 1196 1197 real(dp) :: sy(2) 1198 complex(dpc) :: ys_val 1199 1200 ! Computes the matrix element <Yl'm'|Slm> 1201 call ys(ll,mp,ll,mm,ys_val) 1202 !call ys(ll,mm,ll,mp,ys_val) 1203 sy(1) = real(ys_val) 1204 sy(2) = -aimag(ys_val) 1205 1206 end function sy 1207 1208 end subroutine recip_ylm
m_epjdos/sphericaldens [ Functions ]
[ Top ] [ m_epjdos ] [ Functions ]
NAME
sphericaldens
FUNCTION
Compute the convolution of a function with the unity constant function over a sphere of radius rmax . The function is to be given in reciprocal space, the resulting function is also given in reciprocal space. The routine needs the norm of the reciprocal space vectors. The resulting function in reciprocal space can give the integral of the density in any sphere of that radius, centered on any point, by a simple scalar product.
INPUTS
fofg(2,nfft)=initial function, in reciprocal space gnorm(nfft)=norm of the reciprocal space vectors nfft=(effective) number of FFT grid points (for this processor) rmax=radius of the sphere
OUTPUT
sphfofg(2,nfft)=convoluted function, in reciprocal space
SOURCE
1426 subroutine sphericaldens(fofg,gnorm,nfft,rmax,sphfofg) 1427 1428 !Arguments ------------------------------------ 1429 !scalars 1430 integer,intent(in) :: nfft 1431 real(dp),intent(in) :: rmax 1432 !arrays 1433 real(dp),intent(in) :: fofg(2,nfft),gnorm(nfft) 1434 real(dp),intent(out) :: sphfofg(2,nfft) 1435 1436 !Local variables------------------------------- 1437 !scalars 1438 integer :: ifft 1439 real(dp) :: factor,int0yy,rmax_2pi,yy 1440 1441 ! ************************************************************************* 1442 1443 rmax_2pi=two_pi*rmax 1444 factor=four_pi/(two_pi)**3 1445 1446 do ifft=1,nfft 1447 if(abs(gnorm(ifft)) < tol12)then 1448 sphfofg(1,ifft)=fofg(1,ifft)*four_pi*third*rmax**3 1449 sphfofg(2,ifft)=fofg(2,ifft)*four_pi*third*rmax**3 1450 else 1451 yy=gnorm(ifft)*rmax_2pi 1452 int0yy=factor*(sin(yy)-yy*cos(yy))/(gnorm(ifft)**3) 1453 sphfofg(1,ifft)=fofg(1,ifft)*int0yy 1454 sphfofg(2,ifft)=fofg(2,ifft)*int0yy 1455 end if 1456 end do 1457 1458 end subroutine sphericaldens