TABLE OF CONTENTS


ABINIT/m_epjdos [ Modules ]

[ Top ] [ Modules ]

NAME

  m_epjdos

FUNCTION

  Tools for the computiation of electronic PJDOSes

COPYRIGHT

  Copyright (C) 2008-2018 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 .

PARENTS

SOURCE

19 #if defined HAVE_CONFIG_H
20 #include "config.h"
21 #endif
22 
23 #include "abi_common.h"
24 
25 module m_epjdos
26 
27  use defs_basis
28  use defs_abitypes
29  use m_abicore
30  use m_xmpi
31  use m_errors
32  use m_tetrahedron
33  use m_splines
34  use m_cgtools
35  use m_atomdata
36  use m_crystal
37  use m_crystal_io
38  use m_ebands
39  use m_nctk
40 #ifdef HAVE_NETCDF
41  use netcdf
42 #endif
43  use m_hdr
44  use m_mpinfo
45  use m_sort
46 
47  use defs_datatypes,   only : ebands_t, pseudopotential_type
48  use m_occ,            only : dos_hdr_write
49  use m_time,           only : cwtime, timab
50  use m_io_tools,       only : open_file
51  use m_numeric_tools,  only : simpson, simpson_int
52  use m_fstrings,       only : int2char4, strcat
53  use m_special_funcs,  only : jlspline_t, jlspline_new, jlspline_free, jlspline_integral
54  use m_kpts,           only : tetra_from_kptrlatt
55  use m_kg,             only : ph1d3d, getph
56  use m_gsphere,        only : getkpgnorm
57  use m_fftcore,        only : sphereboundary
58  use m_fft,            only : fftpac, fourwf, fourdp
59  use m_pawrad,         only : pawrad_type, simp_gen
60  use m_pawtab,         only : pawtab_type
61  use m_pawcprj,        only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_free
62  use m_initylmg,       only : initylmg
63 
64  implicit none
65 
66  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.

PARENTS

      m_cut3d

CHILDREN

      cwtime,wrtout

SOURCE

1291 subroutine dens_in_sph(cmax,cg,gmet,istwfk,kg_k,natom,ngfft,mpi_enreg,npw_k,&
1292 &                       paral_kgb,ph1d,rmax,ucvol)
1293 
1294 
1295 !This section has been created automatically by the script Abilint (TD).
1296 !Do not modify the following lines by hand.
1297 #undef ABI_FUNC
1298 #define ABI_FUNC 'dens_in_sph'
1299 !End of the abilint section
1300 
1301  implicit none
1302 
1303 !Arguments ------------------------------------
1304 !scalars
1305  integer,intent(in) :: istwfk,natom,npw_k,paral_kgb
1306  real(dp),intent(in) :: ucvol
1307  type(MPI_type),intent(in) :: mpi_enreg
1308 !arrays
1309  integer,intent(in) :: kg_k(3,npw_k),ngfft(18)
1310  real(dp),intent(in) :: gmet(3,3)
1311  real(dp),intent(in) :: ph1d(2,(2*ngfft(1)+1+2*ngfft(2)+1+2*ngfft(3)+1)*natom)
1312  real(dp),intent(in) :: rmax(natom)
1313  real(dp),intent(inout) :: cg(2,npw_k)
1314  real(dp),intent(out) :: cmax(natom)
1315 
1316 !Local variables -------------------------
1317 !scalars
1318  integer,parameter :: tim_fourwf=0
1319  integer :: cplex,i1,i2,i3,iatom,id1,id2,id3,ifft,mgfft,n1,n2,n3,n4,n5,n6,nfft,nfftot
1320  real(dp) :: cmaxr,g1,g2,g3,norm,weight
1321 !arrays
1322  integer :: ngfft_here(18)
1323  integer,allocatable :: garr(:,:),gbound(:,:)
1324  real(dp),allocatable :: denpot(:,:,:),fofgout(:,:),fofr(:,:,:,:),gnorm(:)
1325  real(dp),allocatable :: ph3d(:,:,:),phkxred(:,:),rhog(:,:),rhor(:)
1326  real(dp),allocatable :: sphrhog(:,:)
1327 
1328 ! *********************************************************************
1329 
1330  n1=ngfft(1)
1331  n2=ngfft(2)
1332  n3=ngfft(3)
1333  n4=ngfft(4)
1334  n5=ngfft(5)
1335  n6=ngfft(6)
1336  nfftot = n1*n2*n3
1337  nfft=n1*n2*n3
1338  ngfft_here(:) = ngfft(:)
1339 !fourwf doesnt work with other options for mode 0 (fft G -> r)
1340  ngfft_here(7)=111
1341  ngfft_here(8)=256
1342  mgfft=maxval(ngfft_here(1:3))
1343 
1344  call sqnorm_g(norm,istwfk,npw_k,cg,mpi_enreg%me_g0,mpi_enreg%comm_fft)
1345 
1346  if (abs(one-norm) > tol6) then
1347    write(std_out,'(a,f8.5)' ) ' dens_in_sph : this state is not normalized : norm=',norm
1348  end if
1349 
1350 !-----------------------------------------------------------------
1351 !inverse FFT of wavefunction to real space => density in real space
1352 !-----------------------------------------------------------------
1353  ABI_ALLOCATE(gbound,(2*mgfft+8,2))
1354  call sphereboundary(gbound,istwfk,kg_k,mgfft,npw_k)
1355 
1356  weight = one
1357  cplex=1
1358  ABI_ALLOCATE(denpot,(cplex*n4,n5,n6))
1359  denpot(:,:,:)=zero
1360  ABI_ALLOCATE(fofgout,(2,npw_k))
1361  ABI_ALLOCATE(fofr,(2,n4,n5,n6))
1362  call fourwf(cplex,denpot,cg,fofgout,fofr,gbound,gbound, &
1363 & istwfk,kg_k,kg_k,mgfft,mpi_enreg,1,ngfft_here,npw_k,&
1364 & npw_k,n4,n5,n6,1,paral_kgb,tim_fourwf,weight,weight)
1365  ABI_DEALLOCATE(fofgout)
1366  ABI_DEALLOCATE(fofr)
1367  ABI_DEALLOCATE(gbound)
1368 
1369  norm = sum(denpot(:,:,:))/nfftot
1370  if (abs(one-norm) > tol6) then
1371    write(std_out,'(a,f8.5)') ' dens_in_sph : this state is not normalized in real space : norm=',norm
1372  end if
1373 
1374 !-----------------------------------------------------------------
1375 !FFT of new density: we obtain n(G) in rhog(1,:)
1376 !-----------------------------------------------------------------
1377 
1378 !Change the packing of the reciprocal space density
1379  ABI_ALLOCATE(rhor,(nfft))
1380  call fftpac(1,mpi_enreg,1,n1,n2,n3,n4,n5,n6,ngfft,rhor,denpot,1)
1381 
1382  ABI_ALLOCATE(rhog,(2,nfft))
1383  call fourdp(1,rhog,rhor,-1,mpi_enreg,nfft,ngfft,paral_kgb,0)
1384 
1385  ABI_DEALLOCATE(rhor)
1386  ABI_DEALLOCATE(denpot)
1387 
1388  do ifft=1,nfft
1389    rhog(:,ifft) = rhog(:,ifft) / ucvol
1390  end do
1391 
1392 !-----------------------------------------------------------------
1393 !calculate norms of G vectors
1394 !-----------------------------------------------------------------
1395 
1396  ABI_ALLOCATE(garr,(3,nfft))
1397  ABI_ALLOCATE(gnorm,(nfft))
1398  id3=ngfft(3)/2+2 ; id2=ngfft(2)/2+2 ; id1=ngfft(1)/2+2
1399  do i3=1,n3
1400    g3=i3-(i3/(id3))*ngfft(3)-one
1401    do i2=1,n2
1402      g2=i2-(i2/(id2))*ngfft(2)-one
1403      do i1=1,n1
1404        g1=i1-(i1/(id1))*ngfft(1)-one
1405        ifft=i1+(i2-1)*n1+(i3-1)*n1*n2
1406        garr(1,ifft)=nint(g1)
1407        garr(2,ifft)=nint(g2)
1408        garr(3,ifft)=nint(g3)
1409        gnorm(ifft)=sqrt(gmet(1,1)*g1*g1 + &
1410 &       two*gmet(2,1)*g2*g1 + &
1411 &       two*gmet(3,1)*g3*g1 + &
1412 &       gmet(2,2)*g2*g2 + &
1413 &       gmet(3,2)*g3*g2 + &
1414 &       gmet(3,3)*g3*g3)
1415      end do
1416    end do
1417  end do
1418 
1419 !-----------------------------------------------------------------
1420 !For each atom call sphericaldens to calculate
1421 !n(G) * 1/|G|^3  *  int_0^2*\pi*r_{max}*|G| 4 \pi y^2 j_0 (y) dy
1422 !for all G vectors put into array sphrhog
1423 !scalar product of phase factors with spherically convoluted density
1424 !-----------------------------------------------------------------
1425 
1426 !largest mem occupation = nfft * (2(sphrog) +2*1(ph3d) +3(garr) +2(rhog) +1(gnorm)) = nfft * 10
1427  ABI_ALLOCATE(sphrhog,(2,nfft))
1428  ABI_ALLOCATE(phkxred,(2,natom))
1429  phkxred(1,:)=one
1430  phkxred(2,:)=zero
1431  ABI_ALLOCATE(ph3d,(2,nfft,1))
1432 
1433  do iatom=1,natom
1434 
1435    call sphericaldens(rhog,gnorm,nfft,rmax(iatom),sphrhog)
1436 !  -----------------------------------------------------------------
1437 !  Compute the phases for the whole set of fft vectors
1438 !  -----------------------------------------------------------------
1439 
1440    call ph1d3d(iatom,iatom,garr,natom,natom,nfft,ngfft(1),ngfft(2),ngfft(3),&
1441 &   phkxred,ph1d,ph3d)
1442 
1443 !  For the phase factors, take the compex conjugate, before evaluating the scalar product
1444    do ifft=1,nfft
1445      ph3d(2,ifft,1)=-ph3d(2,ifft,1)
1446    end do
1447    cplex=2
1448    call dotprod_v(cplex,cmaxr,nfft,1,0,ph3d,sphrhog,mpi_enreg%comm_fft)
1449    cmax(iatom) = cmaxr
1450 !  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
1451  end do
1452 
1453  ABI_DEALLOCATE(rhog)
1454  ABI_DEALLOCATE(gnorm)
1455  ABI_DEALLOCATE(garr)
1456  ABI_DEALLOCATE(sphrhog)
1457  ABI_DEALLOCATE(ph3d)
1458  ABI_DEALLOCATE(phkxred)
1459 
1460 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)

PARENTS

      outscfcv

CHILDREN

      cwtime,wrtout

SOURCE

370 subroutine dos_calcnwrite(dos,dtset,crystal,ebands,fildata,comm)
371 
372 
373 !This section has been created automatically by the script Abilint (TD).
374 !Do not modify the following lines by hand.
375 #undef ABI_FUNC
376 #define ABI_FUNC 'dos_calcnwrite'
377 !End of the abilint section
378 
379  implicit none
380 
381 !Arguments ------------------------------------
382 !scalars
383  integer,intent(in) :: comm
384  character(len=*),intent(in) :: fildata
385  type(dataset_type),intent(in) :: dtset
386  type(crystal_t),intent(in) :: crystal
387  type(ebands_t),intent(in) :: ebands
388  type(epjdos_t),intent(in) :: dos
389 
390 !Local variables-------------------------------
391 !scalars
392  integer,parameter :: bcorr0=0,master=0
393  integer :: iat,iband,iene,ikpt,isppol,natsph,natsph_extra,nkpt,nsppol,i1,i2
394  integer :: nene,prtdos,unitdos,ierr,prtdosm,paw_dos_flag,mbesslang,ndosfraction
395  integer :: my_rank,nprocs,cnt,ifrac,ii
396  real(dp),parameter :: dos_max=9999.9999_dp
397  real(dp) :: buffer,deltaene,enemax,enemin,integral_DOS,max_occ
398  real(dp) :: cpu,wall,gflops
399  logical :: bigDOS,iam_master
400  character(len=10) :: tag
401  character(len=500) :: frmt,frmt_extra,msg
402  type(t_tetrahedron) :: tetra
403 !arrays
404  integer,allocatable :: unt_atsph(:)
405  real(dp) :: list_dp(3)
406  real(dp),allocatable :: tmp_eigen(:),total_dos(:,:,:),eig_dos(:,:)
407  real(dp),allocatable :: dos_m(:,:,:),dos_paw1(:,:,:),dos_pawt1(:,:,:)
408  real(dp),allocatable :: wdt(:,:)
409 
410 ! *********************************************************************
411 
412  my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm); iam_master = (my_rank == master)
413 
414  prtdosm = dos%prtdosm; paw_dos_flag = dos%paw_dos_flag
415  mbesslang = dos%mbesslang; ndosfraction = dos%ndosfraction
416 
417  nkpt = dtset%nkpt; nsppol = dtset%nsppol
418 
419 !m-decomposed DOS not compatible with PAW-decomposed DOS
420  if (prtdosm>=1.and.paw_dos_flag==1) then
421    msg = 'm-decomposed DOS (prtdosm>=1) not compatible with PAW-decomposed DOS (pawprtdos=1) !'
422    MSG_ERROR(msg)
423  end if
424 
425 !Refuse nband different for different kpoints
426 !Note: This means we can pass ebands%eig(:,:,:) instead of eigen(mband*nkpt*nsppol) in packed form
427  do isppol=1,nsppol
428    do ikpt=1,nkpt
429      if ( dtset%nband(nkpt*(isppol-1) + ikpt) /= dtset%nband(1) ) then
430        write(std_out,*) 'tetrahedron: skip subroutine.'
431        write(std_out,*) 'nband must be the same for all kpoints'
432        write(std_out,*) 'nband=', dtset%nband
433        MSG_WARNING('tetrahedron: skip subroutine. See message above')
434        return
435      end if
436    end do
437  end do
438 
439  call cwtime(cpu, wall, gflops, "start")
440 
441  tetra = tetra_from_kptrlatt(crystal, dtset%kptopt, dtset%kptrlatt, dtset%nshiftk, &
442    dtset%shiftk, dtset%nkpt, dtset%kpt, msg, ierr)
443  if (ierr /= 0) then
444    call destroy_tetra(tetra)
445    MSG_WARNING(msg)
446    return
447  end if
448 
449  natsph=dtset%natsph; natsph_extra=dtset%natsph_extra
450 
451  ! Master opens the DOS files.
452  if (iam_master) then
453    if (any(dtset%prtdos == [2, 5])) then
454      if (open_file(fildata, msg, newunit=unitdos, status='unknown', form='formatted', action="write") /= 0) then
455        MSG_ERROR(msg)
456      end if
457 
458    else if (dtset%prtdos == 3) then
459      ! unt_atsph(0) is used for the total DOS.
460      ABI_MALLOC(unt_atsph,(0:natsph+natsph_extra))
461 
462      ! Open file for total DOS as well.
463      if (open_file(strcat(fildata, '_TOTAL'), msg, newunit=unt_atsph(0), &
464          status='unknown', form='formatted', action="write") /= 0) then
465        MSG_ERROR(msg)
466      end if
467 
468      do iat=1,natsph
469        call int2char4(dtset%iatsph(iat),tag)
470        ABI_CHECK((tag(1:1)/='#'),'Bug: string length too short!')
471        if (open_file(strcat(fildata, '_AT', tag), msg, newunit=unt_atsph(iat), &
472            status='unknown', form='formatted', action="write") /= 0) then
473          MSG_ERROR(msg)
474        end if
475      end do
476      ! do extra spheres in vacuum too. Use _ATEXTRA[NUM] suffix
477      do iat=1,natsph_extra
478        call int2char4(iat,tag)
479        ABI_CHECK((tag(1:1)/='#'),'Bug: string length too short!')
480        if (open_file(strcat(fildata, '_ATEXTRA', tag), msg, newunit=unt_atsph(natsph+iat), &
481            status='unknown', form='formatted', action="write") /= 0) then
482          MSG_ERROR(msg)
483        end if
484      end do
485    end if
486  end if
487 
488  ! Write the header of the DOS file, and determine the energy range and spacing
489  prtdos=dtset%prtdos
490  buffer=0.01_dp ! Size of the buffer around the min and max ranges
491 
492  ! A Similar section is present is getnel. Should move all DOS stuff to m_ebands
493  ! Choose the lower and upper energies
494  enemax = maxval(ebands%eig) + buffer
495  enemin = minval(ebands%eig) - buffer
496 
497  ! Extend the range to a nicer value
498  enemax=0.1_dp*ceiling(enemax*10._dp)
499  enemin=0.1_dp*floor(enemin*10._dp)
500 
501  ! Choose the energy increment
502  if(abs(dtset%dosdeltae)<tol10)then
503    deltaene=0.001_dp
504    if(dtset%prtdos>=2)deltaene=0.0005_dp ! Higher resolution possible (and wanted) for tetrahedron
505  else
506    deltaene=dtset%dosdeltae
507  end if
508  nene=nint((enemax-enemin)/deltaene)+1
509 
510  call xmpi_bcast(nene, master, comm, ierr)
511  if (iam_master) list_dp(1:3) = [deltaene, enemin, enemax]
512  call xmpi_bcast(list_dp, master, comm, ierr)
513  deltaene = list_dp(1); enemin = list_dp(2); enemax = list_dp(3)
514 
515  if (iam_master) then
516    if (any(dtset%prtdos == [2, 5])) then
517      call dos_hdr_write(deltaene,ebands%eig,enemax,enemin,ebands%fermie,dtset%mband,&
518      dtset%nband,nene,nkpt,nsppol,dtset%occopt,prtdos,&
519      dtset%tphysel,dtset%tsmear,unitdos)
520    else if (dtset%prtdos == 3) then
521      do iat=0,natsph+natsph_extra
522        call dos_hdr_write(deltaene,ebands%eig,enemax,enemin,ebands%fermie,dtset%mband,&
523        dtset%nband,nene,nkpt,nsppol,dtset%occopt,prtdos,&
524        dtset%tphysel,dtset%tsmear,unt_atsph(iat))
525      end do
526    end if
527  end if
528 
529  ! Tetra weights
530  ABI_MALLOC(wdt, (nene, 2))
531 
532  ! Allocate arrays to store DOSes and fill with zeros.
533  ! 1--> DOS , 2--> IDOS
534  ABI_MALLOC(total_dos,(nene,ndosfraction,2))
535  ABI_MALLOC(eig_dos, (nene, 2))
536 
537  if (paw_dos_flag==1) then
538    ABI_MALLOC(dos_paw1,(nene,ndosfraction,2))
539    ABI_MALLOC(dos_pawt1,(nene,ndosfraction,2))
540  end if
541  if (prtdosm>=1) then
542    ABI_MALLOC(dos_m, (nene,ndosfraction*mbesslang,2))
543  end if
544 
545 !Get maximum occupation value (2 or 1)
546  max_occ = one; if (dtset%nspinor == 1 .and. nsppol == 1) max_occ = two
547 
548 !-------------------------------------------------------------------
549 !For each spin polarisation and band, interpolate band over kpoints
550 !calculate integration weights and DOS contib from
551 !-------------------------------------------------------------------
552 
553  ! Workspace arrays.
554  ABI_MALLOC(tmp_eigen,(nkpt))
555 
556  cnt = 0
557  do isppol=1,nsppol
558 
559    total_dos = zero; eig_dos = zero
560    if (prtdosm>=1) dos_m = zero
561    if (paw_dos_flag==1) then
562      dos_paw1 = zero; dos_pawt1 = zero
563    end if
564 
565    do ikpt=1,nkpt
566       do iband=1,ebands%nband(ikpt+(isppol-1)*ebands%nkpt)
567         cnt = cnt + 1; if (mod(cnt, nprocs) /= my_rank) cycle ! Mpi parallelism.
568 
569         ! Accumulate total DOS from eigenvalues (this is the **exact** total DOS)
570         tmp_eigen(:) = ebands%eig(iband, :, isppol)
571         call tetra_get_onewk(tetra,ikpt,bcorr0,nene,nkpt,tmp_eigen,enemin,enemax,max_occ,wdt)
572         eig_dos = eig_dos + wdt
573 
574         ! Accumulate L-DOS.
575         do ii=1,2
576           do ifrac=1,ndosfraction
577             total_dos(:,ifrac,ii) = total_dos(:,ifrac,ii) + wdt(:,ii) * dos%fractions(ikpt,iband,isppol,ifrac)
578           end do
579         end do
580 
581         if (paw_dos_flag==1) then
582           ! Accumulate L-DOS (on-site terms).
583           do ii=1,2
584             do ifrac=1,ndosfraction
585               dos_paw1(:,ifrac,ii) = dos_paw1(:,ifrac,ii) + wdt(:,ii) * dos%fractions_paw1(ikpt,iband,isppol,ifrac)
586               dos_pawt1(:,ifrac,ii) = dos_pawt1(:,ifrac,ii) + wdt(:,ii) * dos%fractions_pawt1(ikpt,iband,isppol,ifrac)
587             end do
588           end do
589         end if
590 
591         if (prtdosm>=1) then
592          ! Accumulate LM-DOS.
593          do ii=1,2
594            do ifrac=1,ndosfraction*mbesslang
595              dos_m(:,ifrac,ii) = dos_m(:,ifrac,ii) + wdt(:,ii) * dos%fractions_m(ikpt, iband, isppol, ifrac)
596            end do
597          end do
598         end if
599 
600       end do ! ikpt
601    end do ! iband
602 
603    ! Collect results on master
604    call xmpi_sum_master(eig_dos, master, comm, ierr)
605    call xmpi_sum_master(total_dos, master, comm, ierr)
606    bigDOS=(maxval(total_dos(:,:,1))>999._dp)
607 
608    if (paw_dos_flag == 1) then
609      call xmpi_sum_master(dos_paw1, master, comm, ierr)
610      call xmpi_sum_master(dos_pawt1, master, comm, ierr)
611    end if
612    if (prtdosm >= 1) call xmpi_sum_master(dos_m, master, comm, ierr)
613 
614    ! Write the DOS value in the DOS file
615    ! Print the data for this energy. Note the upper limit (dos_max), to be consistent with the format.
616    ! The use of "E" format is not adequate, for portability of the self-testing procedure.
617    ! header lines depend on the type of DOS (projected etc...) which is output
618 
619    if (.not. iam_master) goto 10
620    call write_extra_headers()
621 
622    if (prtdos==2) then
623      ! E, DOS, IDOS
624      do iene=1,nene
625        write(unitdos, '(f11.5,1x,2(f10.4,1x))') &
626         enemin + (iene-1)*deltaene, min(total_dos(iene,:,1), dos_max), total_dos(iene,:,2)
627      end do
628 
629    else if (prtdos==3) then
630 
631      ! Write E, DOS, IDOS
632      do iene=1,nene
633        write(unt_atsph(0), '(f11.5,1x,2(f10.4,1x))') &
634         enemin + (iene-1)*deltaene, min(eig_dos(iene,1), dos_max), eig_dos(iene,2)
635      end do
636 
637      ! E, DOS(L=1,LMAX), IDOS(L=1,LMAX)
638      ! Here we assume mpsang = 5 in the format.
639      if (paw_dos_flag/=1.or.dtset%pawprtdos==2) then
640        frmt = '(f11.5,1x,5(f9.4,1x),10x,5(f8.2,1x),10x,25(f8.2,1x))'
641        if (bigDOS) frmt = '(f11.5,1x,5(f10.4,1x),10x,5(f8.2,1x),10x,25(f8.2,1x))'
642        ! for extra atoms in vacuum need more precision
643        frmt_extra = '(f11.5,1x,5(f20.16,1x),10x,5(f20.16,1x),10x,25(f20.16,1x))'
644 
645        do iat=1,natsph
646          i1 = (iat-1)*mbesslang+1; i2 = iat*mbesslang
647          if (prtdosm==0) then
648            do iene=1,nene
649              write(unt_atsph(iat), fmt=frmt) enemin + (iene-1)*deltaene, &
650 &              min(total_dos(iene, i1:i2, 1), dos_max), total_dos(iene, i1:i2,2)
651            end do
652          else
653            do iene=1,nene
654              write(unt_atsph(iat), fmt=frmt) enemin + (iene-1)*deltaene, &
655 &              min(total_dos(iene, i1:i2, 1), dos_max),&
656 &              total_dos(iene, i1:i2, 2),&
657 &              min(dos_m(iene,(iat-1)*mbesslang**2+1:iat*mbesslang**2,1), dos_max)
658            end do
659          end if
660        end do
661 
662        ! Extra spheres.
663        do iat=natsph+1,natsph+natsph_extra
664          i1 = (iat-1)*mbesslang+1; i2 = iat*mbesslang
665          if (prtdosm==0) then
666            do iene=1,nene
667              write(unt_atsph(iat), fmt=frmt_extra) enemin + (iene-1)*deltaene, &
668 &             total_dos(iene, i1:i2, 1), &
669 &             total_dos(iene, i1:i2, 2)
670            end do
671          else
672            do iene=1,nene
673              write(unt_atsph(iat), fmt=frmt_extra) enemin + (iene-1)*deltaene, &
674 &             total_dos(iene, i1:i2, 1),&
675 &             total_dos(iene, i1:i2, 2),&
676 &             dos_m(iene,(iat-1)*mbesslang**2+1:iat*mbesslang**2, 1)
677            end do
678          end if
679        end do
680 
681      else
682        frmt = '(f11.5,1x,5(f9.4,1x),3(6x,5f9.4))'
683        if (bigDOS) frmt = '(f11.5,1x,5(f10.4,1x),3(6x,5f10.4))'
684        ! for extra atom spheres in vacuum need more precision
685        frmt_extra = '(f11.5,1x,5(f20.16,1x),3(6x,5f20.16))'
686 
687        do iat=1,natsph
688          i1 = iat*5-4; i2 = iat*5
689          do iene=1,nene
690            write(unt_atsph(iat), fmt=frmt) enemin + (iene-1)*deltaene, &
691 &           min(total_dos(iene,i1:i2,1), dos_max),&
692 &           min(total_dos(iene,i1:i2,1) - dos_paw1(iene,i1:i2,1) + dos_pawt1(iene,i1:i2,1), dos_max),&
693 &           min(dos_paw1(iene,i1:i2,1), dos_max),&
694 &           min(dos_pawt1(iene,i1:i2,1), dos_max)
695          end do
696        end do
697 
698        ! Extra spheres.
699        do iat=natsph+1,natsph+natsph_extra
700          i1 = iat*5-4; i2 = iat*5
701          do iene=1,nene
702            write(unt_atsph(iat), fmt=frmt_extra) enemin + (iene-1)*deltaene, &
703 &            min(total_dos(iene,i1:i2,1), dos_max),&
704 &            min(total_dos(iene,i1:i2,1) - dos_paw1(iene,i1:i2,1) + dos_pawt1(iene,i1:i2,1), dos_max),&
705 &            min(dos_paw1(iene,i1:i2,1), dos_max),&
706 &            min(dos_pawt1(iene,i1:i2,1), dos_max)
707          end do
708        end do
709      end if
710 
711    else if (prtdos==5)then
712      ! E, SPIN-DOS
713      frmt = '(f11.5,1x,7(f9.4,1x),10x,7(f8.2,1x))'
714      if (bigDOS) frmt = '(f11.5,1x,7(f10.4,1x),10x,7(f8.2,1x))'
715      do iene=1,nene
716        write(unitdos, fmt=frmt) enemin + (iene-1)*deltaene, min(total_dos(iene,1:7,1), dos_max), total_dos(iene,1:7,2)
717      end do
718    end if
719 
720 10 continue
721    integral_DOS=sum(total_dos(nene,:,2))
722    write(msg, '(a,es16.8)' ) ' tetrahedron : integrate to',integral_DOS
723    call wrtout(std_out,msg,'COLL')
724  end do ! isppol
725 
726  ! Close files.
727  if (iam_master) then
728    if (any(prtdos == [2, 5])) then
729      close(unitdos)
730    else if (prtdos == 3) then
731      do iat=0,natsph+natsph_extra
732        close(unt_atsph(iat))
733      end do
734      ABI_FREE(unt_atsph)
735    end if
736  end if
737 
738  ABI_FREE(tmp_eigen)
739  ABI_FREE(total_dos)
740  ABI_FREE(wdt)
741  ABI_FREE(eig_dos)
742 
743  if (prtdosm>=1)  then
744    ABI_FREE(dos_m)
745  end if
746 
747  if (paw_dos_flag==1)  then
748    ABI_FREE(dos_paw1)
749    ABI_FREE(dos_pawt1)
750  end if
751 
752  call destroy_tetra(tetra)
753 
754  call cwtime(cpu,wall,gflops,"stop")
755  write(msg,'(2(a,f8.2),a)')" tetrahedron: cpu_time: ",cpu,"[s], walltime: ",wall," [s]"
756  call wrtout(std_out,msg,"PERS")
757 
758 contains
759 
760 subroutine write_extra_headers()
761 
762 
763 !This section has been created automatically by the script Abilint (TD).
764 !Do not modify the following lines by hand.
765 #undef ABI_FUNC
766 #define ABI_FUNC 'write_extra_headers'
767 !End of the abilint section
768 
769  if (nsppol==2) then
770    if(isppol==1) write(msg,'(a,16x,a)')  '#','Spin-up DOS'
771    if(isppol==2) write(msg,'(2a,16x,a)')  ch10,'#','Spin-dn DOS'
772    ! NB: dtset%prtdos == 5 should not happen for nsppol==2
773 
774    if (any(dtset%prtdos == [2, 5])) then
775      write(unitdos, "(a)")trim(msg)
776 
777    else if (dtset%prtdos == 3) then
778      do iat=0,natsph+natsph_extra
779        write(unt_atsph(iat), "(a)")trim(msg)
780      end do
781    end if
782 
783  end if
784 
785  if (prtdos==2) then
786    write(unitdos, '(a)' )'# energy(Ha)     DOS  integrated DOS'
787 
788  else if (prtdos==3) then
789 
790    write(unt_atsph(0), '(a)' )'# energy(Ha)     DOS  integrated DOS'
791 
792    if (paw_dos_flag/=1.or.dtset%pawprtdos==2) then
793      do iat=1,natsph
794        write(unt_atsph(iat), '(3a,i5,a,i5,a,a,es16.6,3a)' ) &
795 &       '# Local DOS (columns 2-6) and integrated local DOS (columns 7-11),',ch10,&
796 &       '# for atom number iat=',iat,'  iatom=',dtset%iatsph(iat),ch10,&
797 &       '# inside sphere of radius ratsph=',dtset%ratsph(dtset%typat(dtset%iatsph(iat))),' Bohr.',ch10,"#"
798 
799        if (dtset%usepaw==1.and.dtset%pawprtdos==2) then
800          write(unt_atsph(iat), '(3a)' ) &
801 &         '# PAW: note that only all-electron on-site part has been used to compute DOS !',ch10,"#"
802        end if
803        if (bigDOS) then
804          write(msg, '(a,a)' ) &
805 &         '# energy(Ha)   l=0       l=1       l=2       l=3       l=4',&
806 &         '    (integral=>)  l=0     l=1     l=2     l=3     l=4'
807        else
808          write(msg, '(a,a)' ) &
809 &         '# energy(Ha)  l=0      l=1      l=2      l=3      l=4',&
810 &         '    (integral=>)  l=0     l=1     l=2     l=3     l=4'
811        end if
812        if (prtdosm>=1) then
813          write(msg, '(7a)' ) trim(msg),'          ',&
814 &         '  lm=0 0',&
815 &         '  lm=1-1  lm=1 0  lm=1 1',&
816 &         '  lm=2-2  lm=2-1  lm=2 0  lm=2 1  lm=2 2',&
817 &         '  lm=3-3  lm=3-2  lm=3-1  lm=3 0  lm=3 1  lm=3 2  lm=3 3',&
818 &         '  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'
819        end if
820        write(unt_atsph(iat), "(a)")trim(msg)
821      end do
822    else
823      do iat=1,natsph
824        write(unt_atsph(iat), '(9a,i5,a,i5,a,a,es16.6,3a)' ) &
825 &       '# Local DOS (columns 2-6),',ch10,&
826 &       '#  plane-waves contrib. to DOS (columns 7-11),',ch10,&
827 &       '#  AE on-site  contrib. to DOS (columns 12-16),',ch10,&
828 &       '# -PS on-site  contrib. to DOS (columns 17-21),',ch10,&
829 &       '# for atom number iat=',iat,'  iatom=',dtset%iatsph(iat),ch10,&
830 &       '# inside sphere of radius ratsph=',dtset%ratsph(dtset%typat(dtset%iatsph(iat))),' Bohr.',ch10,"#"
831        if (bigDOS) then
832          write(msg, '(4a)' ) &
833 &         '#energy(Ha)   l=0       l=1       l=2       l=3       l=4',&
834 &         '       (PW)  l=0       l=1       l=2       l=3       l=4',&
835 &         '      (Phi)  l=0       l=1       l=2       l=3       l=4',&
836 &         '     (tPhi)  l=0       l=1       l=2       l=3       l=4'
837        else
838          write(msg, '(4a)' ) &
839 &         '#energy(Ha)  l=0      l=1      l=2      l=3      l=4',&
840 &         '       (PW) l=0      l=1      l=2      l=3      l=4',&
841 &         '      (Phi) l=0      l=1      l=2      l=3      l=4',&
842 &         '     (tPhi) l=0      l=1      l=2      l=3      l=4'
843        end if
844        write(unt_atsph(iat), "(a)")trim(msg)
845      end do
846    end if
847    do iat=1,natsph_extra
848      write(unt_atsph(natsph+iat), '(3a,i5,2a,es16.6,3a)' ) &
849 &     '# Local DOS (columns 2-6) and integrated local DOS (columns 7-11),',ch10,&
850 &     '# for non-atomic sphere number iat=',iat,ch10,&
851 &     '# of radius ratsph=',dtset%ratsph_extra,' Bohr.',ch10,"#"
852      if (bigDOS) then
853        write(msg, '(a,a)' ) &
854 &       '# energy(Ha)   l=0       l=1       l=2       l=3       l=4',&
855 &       '    (integral=>)  l=0     l=1     l=2     l=3     l=4'
856      else
857        write(msg, '(a,a)' ) &
858 &       '# energy(Ha)  l=0      l=1      l=2      l=3      l=4',&
859 &       '    (integral=>)  l=0     l=1     l=2     l=3     l=4'
860      end if
861      if (prtdosm>=1) then
862        write(msg, '(7a)' ) trim(msg),'          ',&
863 &       '  lm=0 0',&
864 &       '  lm=1-1  lm=1 0  lm=1 1',&
865 &       '  lm=2-2  lm=2-1  lm=2 0  lm=2 1  lm=2 2',&
866 &       '  lm=3-3  lm=3-2  lm=3-1  lm=3 0  lm=3 1  lm=3 2  lm=3 3',&
867 &       '  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'
868      end if
869      write(unt_atsph(natsph+iat), "(a)")trim(msg)
870    end do
871 
872  else if (prtdos==5) then
873    write(unitdos, '(a)' )&
874 &    '# energy(Ha)     DOS up,up  up,dn  dn,up  dn,dn  sigma_x sigma_y sigma_z  and integrated DOS components'
875  end if ! prtdos value
876 
877 end subroutine write_extra_headers
878 
879 end subroutine dos_calcnwrite

m_epjdos/epjdos_free [ Functions ]

[ Top ] [ m_epjdos ] [ Functions ]

NAME

  epjdos_free

FUNCTION

  deallocate memory

PARENTS

      outscfcv

CHILDREN

      cwtime,wrtout

SOURCE

291 subroutine epjdos_free(self)
292 
293 
294 !This section has been created automatically by the script Abilint (TD).
295 !Do not modify the following lines by hand.
296 #undef ABI_FUNC
297 #define ABI_FUNC 'epjdos_free'
298 !End of the abilint section
299 
300  implicit none
301 
302 !Arguments ------------------------------------
303  type(epjdos_t),intent(inout) :: self
304 
305 ! *********************************************************************
306 
307  ! integer
308  if (allocated(self%mlang_type)) then
309    ABI_FREE(self%mlang_type)
310  end if
311 
312  ! real
313  if (allocated(self%fractions)) then
314    ABI_FREE(self%fractions)
315  end if
316  if (allocated(self%fractions_m)) then
317    ABI_FREE(self%fractions_m)
318  end if
319  if (allocated(self%fractions_paw1)) then
320    ABI_FREE(self%fractions_paw1)
321  end if
322  if (allocated(self%fractions_pawt1)) then
323    ABI_FREE(self%fractions_pawt1)
324  end if
325 
326 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

PARENTS

SOURCE

171 type(epjdos_t) function epjdos_new(dtset, psps, pawtab) result(new)
172 
173 
174 !This section has been created automatically by the script Abilint (TD).
175 !Do not modify the following lines by hand.
176 #undef ABI_FUNC
177 #define ABI_FUNC 'epjdos_new'
178 !End of the abilint section
179 
180  implicit none
181 
182 !Arguments ------------------------------------
183  type(dataset_type),intent(in) :: dtset
184  type(pseudopotential_type),intent(in) :: psps
185  type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*psps%usepaw)
186 
187 !Local variables-------------------------------
188 !scalars
189  integer :: ierr,itypat,iat
190 
191 ! *********************************************************************
192 
193  new%nkpt = dtset%nkpt; new%mband = dtset%mband; new%nsppol = dtset%nsppol
194 
195  new%prtdos = dtset%prtdos
196  new%partial_dos_flag = 0
197  if (new%prtdos==2) new%partial_dos_flag = 0 ! Standard DOS with tetra.
198  if (new%prtdos==3) new%partial_dos_flag = 1 ! L-DOS with tetra (prtdosm>0 if LM is wanted in Ylm/Slm basis).
199  if (new%prtdos==4) new%partial_dos_flag = 1 ! L-DOS with gaussian (prtdosm if LM is wanted in Ylm/Slm basis).
200  if (new%prtdos==5) new%partial_dos_flag = 2 ! Spin DOS
201 
202  new%prtdosm=0
203  if (new%partial_dos_flag==1) new%prtdosm=dtset%prtdosm
204  ! paw_dos_flag= 1 if both PAW contributions are evaluated AND stored
205  new%paw_dos_flag=0
206  if (dtset%usepaw==1 .and. new%partial_dos_flag==1 .and. dtset%pawprtdos==1) new%paw_dos_flag=1
207 
208  new%fatbands_flag=0
209  if (dtset%pawfatbnd>0 .and. new%prtdosm==0) new%fatbands_flag=1
210  if (new%prtdosm==1.and.dtset%pawfatbnd>0)then
211 !  because they compute quantities in real and complex harmonics respectively
212    MSG_ERROR('pawfatbnd>0  and prtdosm=1 are not compatible')
213  end if
214 
215  ! mjv : initialization is needed as mbesslang is used for allocation below
216  ! NOTE: 10/5/2010 the whole of this could be looped over ndosfraction,
217  ! to store much less in memory. The DOS is accumulated in an array
218  ! and then printed to file at the end.
219  new%mbesslang = 1
220  if (new%partial_dos_flag==1 .or. new%fatbands_flag==1) then
221 
222    ABI_MALLOC(new%mlang_type, (dtset%ntypat + dtset%natsph_extra))
223    new%mlang_type = 0
224 
225    ! TODO: Could use mbesslang = 4 or compute it from psps/pawtab
226    ! Increment by one (could underestimate if vloc = vlmax)
227    if (dtset%usepaw == 0) then
228      do iat=1,dtset%natsph
229        itypat = dtset%typat(dtset%iatsph(iat))
230        new%mlang_type(itypat) = 1 + maxval(psps%indlmn(1, :, itypat))
231      end do
232    else
233      do iat=1,dtset%natsph
234        itypat= dtset%typat(dtset%iatsph(iat))
235        new%mlang_type(itypat) = 1 + (pawtab(itypat)%l_size - 1) / 2
236      end do
237    end if
238 
239    ! Up to l=g if we have natsph_extra.
240    if (dtset%natsph_extra > 0) new%mlang_type(dtset%ntypat+1:) = 5
241 
242    new%mlang_type = 5  ! This is to preserve the old implementation
243    new%mbesslang = maxval(new%mlang_type)
244    new%ndosfraction = (dtset%natsph + dtset%natsph_extra) * new%mbesslang
245 
246  else if (new%partial_dos_flag == 2) then
247    new%ndosfraction = 7
248 
249  else
250    new%ndosfraction = 1
251    new%mbesslang = 0
252  end if
253 
254  ! Check allocations status as these arrays are not distributed and the wavefunctions are still in memory.
255  ABI_STAT_MALLOC(new%fractions, (dtset%nkpt,dtset%mband,dtset%nsppol,new%ndosfraction), ierr)
256  ABI_CHECK(ierr==0, "out of memory in new%fractions")
257  new%fractions = zero
258 
259  if (new%prtdosm>=1 .or. new%fatbands_flag==1) then
260    ABI_STAT_MALLOC(new%fractions_m,(dtset%nkpt,dtset%mband,dtset%nsppol,new%ndosfraction*new%mbesslang), ierr)
261    ABI_CHECK(ierr==0, "out of memory in new%fractions_m")
262    new%fractions_m = zero
263  end if
264 
265  if (dtset%usepaw==1 .and. new%partial_dos_flag==1) then
266    ABI_STAT_MALLOC(new%fractions_paw1,(dtset%nkpt,dtset%mband,dtset%nsppol,new%ndosfraction), ierr)
267    ABI_CHECK(ierr==0, "out of memory in new%fraction_paw1")
268    ABI_STAT_MALLOC(new%fractions_pawt1,(dtset%nkpt,dtset%mband,dtset%nsppol,new%ndosfraction), ierr)
269    ABI_CHECK(ierr==0, "out of memory in new%fraction_pawt1")
270    new%fractions_paw1 = zero; new%fractions_pawt1 = zero
271  end if
272 
273 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

 94  type,public :: epjdos_t
 95 
 96    integer :: mbesslang
 97    ! Max L+1 used in LM-DOS  (Bessel function expansion)
 98 
 99    integer :: ndosfraction
100    ! Defines the last dimension of the dos arrays.
101    ! Actual value depends on the other variables.
102 
103    integer :: prtdos
104    ! 2 --> Standard DOS with tetra.
105    ! 3 --> L-DOS with tetra (prtdosm>0 if LM is wanted in Ylm/Slm basis).
106    ! 4 --> L-DOS with gaussian (prtdosm if LM is wanted in Ylm/Slm basis).
107    ! 5 --> Spin-DOS
108 
109    integer :: prtdosm
110    ! Option for the m-contributions to the partial DOS
111    ! 1 if LM-projection is done onto complex Ylm
112    ! 2 if LM-projection is done onto real Slm
113 
114    integer :: partial_dos_flag
115 
116    integer :: paw_dos_flag
117    ! 1 if both PAW contributions are evaluated AND stored
118 
119    !integer :: pawfatbnd
120    integer :: fatbands_flag
121 
122    integer :: nkpt, mband, nsppol
123    ! Used to dimension arrays
124 
125    integer,allocatable :: mlang_type(:)
126    ! mlang_type(ntypat + natsph_extra)
127    ! Max L+1 used in LM-DOS for each atom type
128 
129    real(dp),allocatable :: fractions(:,:,:,:)
130    ! fractions(nkpt,mband,nsppol,ndosfraction))
131    ! TODO: replace nsppol with nspden = 1, 2 (nsppol==2) or 4 (nspinor==2)
132 
133    real(dp),allocatable :: fractions_m(:,:,:,:)
134    ! fractions_m(nkpt,mband,nsppol,ndosfraction*mbesslang)
135 
136    real(dp),allocatable :: fractions_paw1(:,:,:,:)
137    ! fractions_paw1(nkpt,mband,nsppol,ndosfraction)
138 
139    real(dp),allocatable :: fractions_pawt1(:,:,:,:)
140    ! fractions_pawt1(nkpt,mband,nsppol,ndosfraction))
141 
142  end type epjdos_t
143 
144  public :: epjdos_new         ! Create new object
145  public :: epjdos_free        ! Free dynamic memory
146 
147  public :: prtfatbands        ! Print PJDOS contributions in xmgrace format.
148  public :: fatbands_ncwrite   ! Write PJDOS contributions to netcdf file.
149 
150 !----------------------------------------------------------------------
151 
152 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

PARENTS

      outscfcv

CHILDREN

      cwtime,wrtout

SOURCE

1815 subroutine fatbands_ncwrite(dos, crystal, ebands, hdr, dtset, psps, pawtab, ncid)
1816 
1817 
1818 !This section has been created automatically by the script Abilint (TD).
1819 !Do not modify the following lines by hand.
1820 #undef ABI_FUNC
1821 #define ABI_FUNC 'fatbands_ncwrite'
1822 !End of the abilint section
1823 
1824  implicit none
1825 
1826 !Arguments ------------------------------------
1827 !scalars
1828  integer,intent(in) :: ncid
1829  type(epjdos_t),intent(in) :: dos
1830  type(crystal_t),intent(in) :: crystal
1831  type(ebands_t),intent(in) :: ebands
1832  type(hdr_type),intent(in) :: hdr
1833  type(dataset_type),intent(in) :: dtset
1834  type(pseudopotential_type),intent(in) :: psps
1835 !arrays
1836  type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw)
1837 
1838 #ifdef HAVE_NETCDF
1839 !Local variables-------------------------------
1840 !scalars
1841  integer :: itype,ncerr,fform
1842  real(dp) :: cpu,wall,gflops
1843  character(len=500) :: msg
1844 !arrays
1845  integer :: lmax_type(crystal%ntypat)
1846 
1847 !*************************************************************************
1848 
1849  ABI_CHECK(dtset%natsph > 0, "natsph <=0")
1850  call cwtime(cpu, wall, gflops, "start")
1851 
1852  fform = fform_from_ext("FATBANDS.nc")
1853  ABI_CHECK(fform /= 0, "Cannot find fform associated to FATBANDS.nc")
1854 
1855  ! Write header, crystal structure and band energies.
1856  NCF_CHECK(hdr_ncwrite(hdr, ncid, fform, nc_define=.True.))
1857  NCF_CHECK(crystal_ncwrite(crystal, ncid))
1858  NCF_CHECK(ebands_ncwrite(ebands, ncid))
1859 
1860  ! Add fatband-specific quantities
1861  ncerr = nctk_def_dims(ncid, [ &
1862    nctkdim_t("natsph", dtset%natsph), &
1863    nctkdim_t("ndosfraction", dos%ndosfraction)], defmode=.True.)
1864  NCF_CHECK(ncerr)
1865 
1866  if (dos%ndosfraction*dos%mbesslang > 0) then
1867    ncerr = nctk_def_dims(ncid, [ &
1868      nctkdim_t("mbesslang", dos%mbesslang), &
1869      nctkdim_t("dos_fractions_m_lastsize", dos%ndosfraction*dos%mbesslang)])
1870    NCF_CHECK(ncerr)
1871  end if
1872  if (dtset%natsph_extra /= 0) then
1873    NCF_CHECK(nctk_def_dims(ncid, [nctkdim_t("natsph_extra", dtset%natsph_extra)]))
1874  end if
1875 
1876  ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "prtdos", "pawprtdos", "prtdosm"])
1877  NCF_CHECK(ncerr)
1878  ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "ratsph_extra"])
1879  NCF_CHECK(ncerr)
1880 
1881  ncerr = nctk_def_arrays(ncid, [&
1882    nctkarr_t("lmax_type", "int", "number_of_atom_species"), &
1883    nctkarr_t("iatsph", "int", "natsph"), &
1884    nctkarr_t("ratsph", "dp", "number_of_atom_species"), &
1885    nctkarr_t("dos_fractions", "dp", "number_of_kpoints, max_number_of_states, number_of_spins, ndosfraction") &
1886  ])
1887  NCF_CHECK(ncerr)
1888 
1889  if (allocated(dos%fractions_m)) then
1890    ncerr = nctk_def_arrays(ncid, &
1891      nctkarr_t("dos_fractions_m", "dp", &
1892                "number_of_kpoints, max_number_of_states, number_of_spins, dos_fractions_m_lastsize"))
1893    NCF_CHECK(ncerr)
1894  end if
1895 
1896  if (allocated(dos%fractions_paw1)) then
1897    ncerr = nctk_def_arrays(ncid, [&
1898      nctkarr_t("dos_fractions_paw1", "dp", "number_of_kpoints, max_number_of_states, number_of_spins, ndosfraction"), &
1899      nctkarr_t("dos_fractions_pawt1", "dp", "number_of_kpoints, max_number_of_states, number_of_spins, ndosfraction") &
1900    ])
1901    NCF_CHECK(ncerr)
1902  end if
1903 
1904  if (dtset%natsph_extra /= 0) then
1905    ncerr = nctk_def_arrays(ncid, [&
1906      nctkarr_t("xredsph_extra", "dp", "number_of_reduced_dimensions, natsph_extra") &
1907    ])
1908    NCF_CHECK(ncerr)
1909  end if
1910 
1911  ! Write variables
1912  NCF_CHECK(nctk_set_datamode(ncid))
1913 
1914  ! scalars
1915  NCF_CHECK(nf90_put_var(ncid, vid("pawprtdos"), dtset%pawprtdos))
1916  NCF_CHECK(nf90_put_var(ncid, vid("prtdos"), dos%prtdos))
1917  NCF_CHECK(nf90_put_var(ncid, vid("prtdosm"), dos%prtdosm))
1918 
1919  ! arrays
1920  if (dtset%usepaw == 1) then
1921    lmax_type = (pawtab(:)%l_size - 1) / 2
1922  else
1923    do itype=1,crystal%ntypat
1924      lmax_type(itype) = maxval(psps%indlmn(1, :, itype))
1925    end do
1926  end if
1927  NCF_CHECK(nf90_put_var(ncid, vid("lmax_type"), lmax_type))
1928  NCF_CHECK(nf90_put_var(ncid, vid("dos_fractions"), dos%fractions))
1929 
1930  if (dos%prtdos == 3) then
1931    NCF_CHECK(nf90_put_var(ncid, vid("iatsph"), dtset%iatsph(1:dtset%natsph)))
1932    NCF_CHECK(nf90_put_var(ncid, vid("ratsph"), dtset%ratsph(1:dtset%ntypat)))
1933    NCF_CHECK(nf90_put_var(ncid, vid("ratsph_extra"), dtset%ratsph_extra))
1934    if (dtset%natsph_extra /= 0) then
1935      NCF_CHECK(nf90_put_var(ncid, vid("xredsph_extra"), dtset%xredsph_extra(:, 1:dtset%natsph_extra)))
1936    end if
1937  end if
1938 
1939  if (allocated(dos%fractions_m)) then
1940    NCF_CHECK(nf90_put_var(ncid, vid("dos_fractions_m"), dos%fractions_m))
1941  end if
1942  if (allocated(dos%fractions_paw1)) then
1943    NCF_CHECK(nf90_put_var(ncid, vid("dos_fractions_paw1"), dos%fractions_paw1))
1944    NCF_CHECK(nf90_put_var(ncid, vid("dos_fractions_pawt1"), dos%fractions_pawt1))
1945  end if
1946 
1947  call cwtime(cpu,wall,gflops,"stop")
1948  write(msg,'(2(a,f8.2),a)')" fatbands_ncwrite: cpu_time: ",cpu,"[s], walltime: ",wall," [s]"
1949  call wrtout(std_out,msg,"PERS")
1950 #endif
1951 
1952 contains
1953  integer function vid(vname)
1954 
1955 !This section has been created automatically by the script Abilint (TD).
1956 !Do not modify the following lines by hand.
1957 #undef ABI_FUNC
1958 #define ABI_FUNC 'vid'
1959 !End of the abilint section
1960 
1961    character(len=*),intent(in) :: vname
1962    vid = nctk_idname(ncid, vname)
1963  end function vid
1964 
1965 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.

PARENTS

      outscfcv

CHILDREN

      cg_getspin,cwtime,destroy_mpi_enreg,getkpgnorm,getph,initmpi_seq
      initylmg,jlspline_free,ph1d3d,recip_ylm,sort_dp,splint,wrtout,xmpi_sum

SOURCE

2023 subroutine partial_dos_fractions(dos,crystal,dtset,eigen,occ,npwarr,kg,cg,mcg,collect,mpi_enreg)
2024 
2025 
2026 !This section has been created automatically by the script Abilint (TD).
2027 !Do not modify the following lines by hand.
2028 #undef ABI_FUNC
2029 #define ABI_FUNC 'partial_dos_fractions'
2030 !End of the abilint section
2031 
2032  implicit none
2033 
2034 !Arguments ------------------------------------
2035 !scalars
2036  integer,intent(in) :: mcg,collect
2037  type(epjdos_t),intent(inout) :: dos
2038  type(MPI_type),intent(inout) :: mpi_enreg
2039  type(dataset_type),intent(in) :: dtset
2040  type(crystal_t),intent(in) :: crystal
2041 !arrays
2042  integer,intent(in) :: kg(3,dtset%mpw*dtset%mkmem),npwarr(dtset%nkpt)
2043  real(dp),intent(in) :: cg(2,mcg)
2044  real(dp),intent(in) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol)
2045  real(dp),intent(in) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol)
2046 
2047 !Local variables-------------------------------
2048 !scalars
2049  logical,parameter :: write_procar = .False.
2050  integer,parameter :: prtsphere0=0 ! do not output all the band by band details for projections.
2051  integer :: shift_b,shift_sk,iat,iatom,iband,ierr,ikpt,ilang,ioffkg,is1, is2, isoff
2052  integer :: ipw,isppol,ixint,mbess,mcg_disk,me_kpt,shift_cg
2053  integer :: mgfft,my_nspinor,n1,n2,n3,natsph_tot,npw_k,nradintmax
2054  integer :: rc_ylm,itypat,nband_k
2055  integer :: abs_shift_b
2056  integer :: unit_procar, ipauli
2057  real(dp),parameter :: bessint_delta = 0.1_dp
2058  real(dp) :: arg,bessarg,bessargmax,kpgmax,rmax
2059  real(dp) :: cpu,wall,gflops
2060  character(len=500) :: msg
2061  character(len=4) :: ikproc_str
2062  character(len=fnlen) :: filename
2063  type(jlspline_t) :: jlspl
2064  type(MPI_type) :: mpi_enreg_seq
2065 !arrays
2066  integer :: iindex(dtset%mpw),nband_tmp(1),npwarr_tmp(1)
2067  integer,allocatable :: iatsph(:),nradint(:),atindx(:),typat_extra(:),kg_k(:,:)
2068  real(dp) :: kpoint(3),spin(3),ylmgr_dum(1)
2069  real(dp) :: xfit(dtset%mpw),yfit(dtset%mpw)
2070  real(dp),allocatable :: ylm_k(:,:)
2071  real(dp),allocatable :: bess_fit(:,:,:)
2072  real(dp),allocatable :: cg_1band(:,:),cg_1kpt(:,:),kpgnorm(:),ph1d(:,:)
2073  real(dp),allocatable :: ph3d(:,:,:),ratsph(:),rint(:),sum_1atom_1ll(:,:,:)
2074  real(dp),allocatable :: sum_1atom_1lm(:,:,:)
2075  real(dp),allocatable :: xred_sph(:,:),znucl_sph(:),phkxred(:,:)
2076  complex(dpc) :: cgcmat(2,2)
2077 
2078 !*************************************************************************
2079 
2080  ! for the moment, only support projection on angular momenta
2081  if (dos%partial_dos_flag /= 1 .and. dos%partial_dos_flag /= 2) then
2082    write(std_out,*) 'Error: partial_dos_fractions only supports angular '
2083    write(std_out,*) ' momentum projection and spinor components for the moment. return to outscfcv'
2084    write(std_out,*) ' partial_dos = ', dos%partial_dos_flag
2085    return
2086  end if
2087 
2088  ! impose all kpoints have same number of bands
2089  do isppol=1,dtset%nsppol
2090    do ikpt=1,dtset%nkpt
2091      if (dtset%nband((isppol-1)*dtset%nkpt + ikpt) /= dtset%mband) then
2092        write(std_out,*) 'Error: partial_dos_fractions wants same number of bands at each kpoint'
2093        write(std_out,*) ' isppol, ikpt = ', isppol,ikpt, dtset%nband((isppol-1)*dtset%nkpt + ikpt), dtset%mband
2094        write(std_out,*) ' all nband = ', dtset%nband
2095        return
2096      end if
2097    end do
2098  end do
2099 
2100  ! Real or complex spherical harmonics?
2101  rc_ylm = 2; if (dos%prtdosm == 2) rc_ylm = 1
2102 
2103  me_kpt = mpi_enreg%me_kpt
2104  my_nspinor = max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
2105 
2106  n1 = dtset%ngfft(1); n2 = dtset%ngfft(2); n3 = dtset%ngfft(3)
2107  mgfft = maxval(dtset%ngfft(1:3))
2108 
2109  call cwtime(cpu, wall, gflops, "start")
2110 
2111  if (write_procar) then
2112    ! open file for each proc, and print header for master node
2113    call int2char4(me_kpt, ikproc_str)
2114    filename = 'PROCAR_'//ikproc_str
2115    if (open_file(filename, msg, newunit=unit_procar, form="formatted", action="write") /= 0) then
2116       MSG_ERROR(msg)
2117    end if
2118    if(mpi_enreg%me==0) then
2119      write (unit_procar,'(a)') 'PROCAR lm decomposed - need to concatenate files in parallel case'
2120      write (unit_procar,'(a,I10,a,I10,a,I10,a)') '# of k-points: ', dtset%nkpt, &
2121        ' # of bands:', dtset%mband, ' # of ions:', dtset%natom, ch10
2122    end if
2123  end if
2124 
2125 !##############################################################
2126 !FIRST CASE: project on angular momenta to get dos parts
2127 !##############################################################
2128 
2129  if (dos%partial_dos_flag == 1) then
2130    natsph_tot = dtset%natsph + dtset%natsph_extra
2131 
2132    ABI_ALLOCATE(iatsph, (natsph_tot))
2133    ABI_ALLOCATE(typat_extra, (natsph_tot))
2134    ABI_ALLOCATE(ratsph, (natsph_tot))
2135    ABI_ALLOCATE(znucl_sph, (natsph_tot))
2136    ABI_ALLOCATE(nradint, (natsph_tot))
2137    ABI_ALLOCATE(atindx, (natsph_tot))
2138    ABI_ALLOCATE(phkxred, (2,natsph_tot))
2139 
2140    ! initialize atindx
2141    do iatom=1,natsph_tot
2142      atindx(iatom) = iatom
2143    end do
2144 
2145    iatsph(1:dtset%natsph) = dtset%iatsph(1:dtset%natsph)
2146    do iatom=1,dtset%natsph
2147      itypat = dtset%typat(iatsph(iatom))
2148      typat_extra(iatom) = itypat
2149      ratsph(iatom) = dtset%ratsph(itypat)
2150      znucl_sph(iatom) = dtset%znucl(itypat)
2151    end do
2152 
2153    ! fictitious atoms are declared with
2154    ! %natsph_extra, %ratsph_extra and %xredsph_extra(3, dtset%natsph_extra)
2155    ! they have atom index (natom + ii) and itype = ntypat + 1
2156    do iatom=1,dtset%natsph_extra
2157      typat_extra(iatom+dtset%natsph) = dtset%ntypat + 1
2158      ratsph(iatom+dtset%natsph) = dtset%ratsph_extra
2159      znucl_sph(iatom+dtset%natsph) = zero
2160      iatsph(iatom+dtset%natsph) = dtset%natom + iatom
2161    end do
2162 
2163    ! init bessel function integral for recip_ylm max ang mom + 1
2164    ABI_ALLOCATE(sum_1atom_1ll,(dtset%nspinor**2,dos%mbesslang,natsph_tot))
2165    ABI_ALLOCATE(sum_1atom_1lm,(dtset%nspinor**2,dos%mbesslang**2,natsph_tot))
2166 
2167    ! Note ecuteff instead of ecut.
2168    kpgmax = sqrt(dtset%ecut * dtset%dilatmx**2)
2169    rmax = zero; bessargmax = zero; nradintmax = 0
2170    do iatom=1,natsph_tot
2171      rmax = max(rmax, ratsph(iatom))
2172      bessarg = ratsph(iatom)*two_pi*kpgmax
2173      bessargmax = max(bessargmax, bessarg)
2174      nradint(iatom) = int (bessarg / bessint_delta) + 1
2175      nradintmax = max(nradintmax,nradint(iatom))
2176    end do
2177    !write(std_out,*)' partial_dos_fractions: rmax=', rmax,' nradintmax: ", nradintmax
2178 !  use same number of grid points to calculate Bessel function and to do the integration later on r
2179 !  and make sure bessargmax is a multiple of bessint_delta
2180    mbess = nradintmax
2181    bessargmax = bessint_delta*mbess
2182 
2183    ABI_ALLOCATE(rint,(nradintmax))
2184    ABI_ALLOCATE(bess_fit,(dtset%mpw,nradintmax,dos%mbesslang))
2185 
2186    ! initialize general Bessel function array on uniform grid xx, from 0 to (2 \pi |k+G|_{max} |r_{max}|)
2187    jlspl = jlspline_new(mbess, bessint_delta, dos%mbesslang)
2188 
2189    ABI_ALLOCATE(xred_sph, (3, natsph_tot))
2190    do iatom=1,dtset%natsph
2191      xred_sph(:,iatom) = crystal%xred(:,iatsph(iatom))
2192    end do
2193    do iatom=1,dtset%natsph_extra
2194      xred_sph(:,dtset%natsph+iatom) = dtset%xredsph_extra(:, iatom)
2195    end do
2196 
2197    ABI_ALLOCATE(ph1d,(2,(2*n1+1 + 2*n2+1 + 2*n3+1)*natsph_tot))
2198    call getph(atindx,natsph_tot,n1,n2,n3,ph1d,xred_sph)
2199 
2200    ! Fake MPI data to be used for sequential call to initylmg.
2201    call initmpi_seq(mpi_enreg_seq)
2202    mpi_enreg_seq%my_natom = dtset%natom
2203 
2204    shift_sk = 0
2205    abs_shift_b =  0 ! offset to allow for automatic update with +1 below
2206    do isppol=1,dtset%nsppol
2207      ioffkg = 0
2208      do ikpt=1,dtset%nkpt
2209        nband_k = dtset%nband((isppol-1)*dtset%nkpt + ikpt)
2210        if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_kpt)) then
2211          abs_shift_b = abs_shift_b + nband_k ! jump the whole kpt in the eig and occ arrays
2212          cycle
2213        end if
2214        npw_k = npwarr(ikpt)
2215        kpoint(:) = dtset%kpt(:,ikpt)
2216 
2217        if (write_procar) then
2218          write (unit_procar,'(a,I7,a,3F12.6,a,F12.6,a)') &
2219            ' k-point ', ikpt, ' : ', kpoint(:), ' weight = ', dtset%wtk(ikpt), ch10
2220        end if
2221 
2222        ! make phkred for all atoms
2223        do iat=1,natsph_tot
2224          arg=two_pi*( kpoint(1)*xred_sph(1,iat) + kpoint(2)*xred_sph(2,iat) + kpoint(3)*xred_sph(3,iat) )
2225          phkxred(1,iat)=cos(arg)
2226          phkxred(2,iat)=sin(arg)
2227        end do
2228 
2229        ABI_MALLOC(kg_k, (3, npw_k))
2230        kg_k = kg(:,ioffkg+1:ioffkg+npw_k)
2231 
2232        ! kpgnorm contains norms only for kpoints used by this processor
2233        ABI_ALLOCATE(kpgnorm, (npw_k))
2234        call getkpgnorm(crystal%gprimd, kpoint, kg_k, kpgnorm, npw_k)
2235 
2236        ! Now get Ylm(k, G) factors: returns "real Ylms", which are real (+m) and
2237        ! imaginary (-m) parts of actual complex Ylm. Yl-m = Ylm*
2238        ! Call initylmg for a single k-point (mind mpi_enreg_seq).
2239        ABI_MALLOC(ylm_k, (npw_k, dos%mbesslang**2))
2240        npwarr_tmp(1) = npw_k; nband_tmp(1) = nband_k
2241        call initylmg(crystal%gprimd,kg_k,kpoint,1,mpi_enreg_seq,dos%mbesslang,&
2242        npw_k,nband_tmp,1,npwarr_tmp,1,0,crystal%rprimd,ylm_k,ylmgr_dum)
2243 
2244        ! get phases exp (2 pi i (k+G).x_tau) in ph3d
2245        ABI_ALLOCATE(ph3d,(2,npw_k,natsph_tot))
2246        call ph1d3d(1,natsph_tot,kg_k,natsph_tot,natsph_tot,npw_k,n1,n2,n3,phkxred,ph1d,ph3d)
2247 
2248        ! get Bessel function factors on array of |k+G|*r distances
2249        ! since we need many r distances and have a large number of different
2250        ! |k+G|, get j_l on uniform grid (above, in array gen_besj),
2251        ! and spline it for each kpt Gvector set.
2252        ! Note that we use the same step based on rmax, this can lead to (hopefully small)
2253        ! inaccuracies when we integrate from 0 up to rmax(iatom)
2254        do ixint=1,nradintmax
2255          rint(ixint) = (ixint-1)*rmax / (nradintmax-1)
2256          do ipw=1,npw_k
2257            xfit(ipw) = two_pi * kpgnorm(ipw) * rint(ixint)
2258            iindex(ipw) = ipw
2259          end do
2260 
2261          call sort_dp(npw_k,xfit,iindex,tol14)
2262          do ilang=1,dos%mbesslang
2263            call splint(mbess, jlspl%xx, jlspl%bess_spl(:,ilang), jlspl%bess_spl_der(:,ilang), npw_k, xfit, yfit)
2264            ! re-order results for different G vectors
2265            do ipw=1,npw_k
2266              bess_fit(iindex(ipw),ixint,ilang) = yfit(ipw)
2267            end do
2268          end do
2269        end do ! ixint
2270 
2271        shift_b = 0
2272        do iband=1,nband_k
2273          abs_shift_b = abs_shift_b + 1
2274          if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband,isppol,me_kpt)) cycle
2275          !write(std_out,*)"in band:",iband
2276          ! TODO: eventually import eig and occ down to here - a pain, but printing outside would imply saving a huge array in memory
2277          if (write_procar) then
2278            write (unit_procar,'(a,I7,a,F12.6,a,F12.6,a)') 'band ', iband, ' # energy ', &
2279              eigen(abs_shift_b), ' # occ. ', occ(abs_shift_b), ch10
2280          end if
2281 
2282          ! Select wavefunction in cg array
2283          shift_cg = shift_sk + shift_b
2284 
2285          call recip_ylm(bess_fit, cg(:,shift_cg+1:shift_cg+my_nspinor*npw_k), dtset%istwfk(ikpt),&
2286 &          mpi_enreg, nradint, nradintmax, dos%mbesslang , dtset%mpw, natsph_tot, typat_extra, dos%mlang_type,&
2287 &          npw_k, dtset%nspinor, ph3d, prtsphere0, rint, ratsph, rc_ylm, sum_1atom_1ll, sum_1atom_1lm,&
2288 &          crystal%ucvol, ylm_k, znucl_sph)
2289          ! on exit the sum_1atom_* have both spinors counted
2290 
2291          ! Accumulate
2292          do iatom=1,natsph_tot
2293            do ilang=1,dos%mbesslang
2294              dos%fractions(ikpt,iband,isppol,dos%mbesslang*(iatom-1) + ilang) &
2295 &             = dos%fractions(ikpt,iband,isppol,dos%mbesslang*(iatom-1) + ilang) &
2296 &             + sum_1atom_1ll(1,ilang,iatom)
2297            end do
2298          end do
2299 
2300          if (dos%prtdosm /= 0) then
2301            do iatom=1,natsph_tot
2302              do ilang=1,dos%mbesslang**2
2303                dos%fractions_m(ikpt,iband,isppol,dos%mbesslang**2*(iatom-1) + ilang) &
2304 &               = dos%fractions_m(ikpt,iband,isppol,dos%mbesslang**2*(iatom-1) + ilang) &
2305 &               + sum_1atom_1lm(1,ilang,iatom)
2306              end do
2307            end do
2308          end if
2309 
2310          ! Increment band, spinor shift
2311          !shift_b = shift_b + npw_k
2312          shift_b = shift_b + my_nspinor*npw_k
2313 
2314          ! now we have both spinor components.
2315          if (write_procar) then
2316            write (unit_procar,'(a)') 'ion      s     py     pz     px    dxy    dyz    dz2    dxz    dx2    tot'
2317            do ipauli= 1,dtset%nspinor**2
2318              ! Contract with Pauli matrices to get projections for this k and band, all atoms and ilang
2319              do iatom = 1, natsph_tot
2320                write (unit_procar, '(I3)', advance='no') iatom
2321                do ilang=1,min(dos%mbesslang**2,9)
2322                  write (unit_procar, '(F7.3)',advance='no') sum_1atom_1lm(ipauli,ilang,iatom)
2323                end do
2324                write (unit_procar, '(F7.3)',advance='yes') sum(sum_1atom_1lm(ipauli,:,iatom))
2325              end do
2326            end do
2327            write (unit_procar,*)
2328         end if
2329 
2330        end do ! band
2331 
2332        ! Increment kpt and (spin, kpt) shifts
2333        ioffkg = ioffkg + npw_k
2334        shift_sk = shift_sk + nband_k*my_nspinor*npw_k
2335 
2336        ABI_FREE(kg_k)
2337        ABI_FREE(kpgnorm)
2338        ABI_FREE(ylm_k)
2339        ABI_DEALLOCATE(ph3d)
2340      end do ! ikpt
2341    end do ! isppol
2342 
2343    ! collect = 1 ==> gather all contributions from different processors
2344    if (collect == 1) then
2345      call xmpi_sum(dos%fractions,mpi_enreg%comm_kpt,ierr)
2346      if (dos%prtdosm /= 0) call xmpi_sum(dos%fractions_m,mpi_enreg%comm_kpt,ierr)
2347 
2348 ! this is now done inside recip_ylm
2349 !     if (mpi_enreg%paral_spinor == 1)then
2350 !       call xmpi_sum(dos%fractions,mpi_enreg%comm_spinor,ierr)
2351 !       if (dos%prtdosm /= 0) call xmpi_sum(dos%fractions_m,mpi_enreg%comm_spinor,ierr)
2352 !     end if
2353    end if
2354 
2355    ABI_DEALLOCATE(atindx)
2356    ABI_DEALLOCATE(bess_fit)
2357    ABI_DEALLOCATE(iatsph)
2358    ABI_DEALLOCATE(typat_extra)
2359    ABI_DEALLOCATE(nradint)
2360    ABI_DEALLOCATE(ph1d)
2361    ABI_DEALLOCATE(phkxred)
2362    ABI_DEALLOCATE(ratsph)
2363    ABI_DEALLOCATE(rint)
2364    ABI_DEALLOCATE(sum_1atom_1ll)
2365    ABI_DEALLOCATE(sum_1atom_1lm)
2366    ABI_DEALLOCATE(xred_sph)
2367    ABI_DEALLOCATE(znucl_sph)
2368 
2369    call jlspline_free(jlspl)
2370    call destroy_mpi_enreg(mpi_enreg_seq)
2371 
2372  !##############################################################
2373  !2ND CASE: project on spinors
2374  !##############################################################
2375 
2376  else if (dos%partial_dos_flag == 2) then
2377 
2378    if (dtset%nsppol /= 1 .or. dtset%nspinor /= 2) then
2379      MSG_WARNING("spinor projection is meaningless if nsppol==2 or nspinor/=2. Not calculating projections.")
2380      return
2381    end if
2382    if (my_nspinor /= 2) then
2383      MSG_WARNING("spinor projection with spinor parallelization is not coded. Not calculating projections.")
2384      return
2385    end if
2386    ABI_CHECK(mpi_enreg%paral_spinor == 0, "prtdos 5 does not support spinor parallelism")
2387 
2388    ! FIXME: We should not allocate such a large chunk of memory!
2389    mcg_disk = dtset%mpw*my_nspinor*dtset%mband
2390    ABI_ALLOCATE(cg_1kpt,(2,mcg_disk))
2391    shift_sk = 0
2392    isppol = 1
2393 
2394    do ikpt=1,dtset%nkpt
2395      nband_k = dtset%nband((isppol-1)*dtset%nkpt + ikpt)
2396      if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_kpt)) cycle
2397      npw_k = npwarr(ikpt)
2398 
2399      cg_1kpt(:,:) = cg(:,shift_sk+1:shift_sk+mcg_disk)
2400      ABI_ALLOCATE(cg_1band,(2,2*npw_k))
2401      shift_b=0
2402      do iband=1,nband_k
2403        if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband,isppol,me_kpt)) cycle
2404 
2405        ! Select wavefunction in cg array
2406        !shift_cg = shift_sk + shift_b
2407        cg_1band(:,:) = cg_1kpt(:,shift_b+1:shift_b+2*npw_k)
2408        call cg_getspin(cg_1band, npw_k, spin, cgcmat=cgcmat)
2409 
2410        ! MG: TODO: imag part of off-diagonal terms is missing.
2411        ! I will add them later on.
2412        do is1=1,2
2413          do is2=1,2
2414            isoff = is2 + (is1-1)*2
2415            dos%fractions(ikpt,iband,isppol,isoff) = dos%fractions(ikpt,iband,isppol,isoff) &
2416 &           + real(cgcmat(is1,is2))
2417          end do
2418        end do
2419 
2420        dos%fractions(ikpt,iband,isppol,5) = dos%fractions(ikpt,iband,isppol,5) + spin(1)
2421        dos%fractions(ikpt,iband,isppol,6) = dos%fractions(ikpt,iband,isppol,6) + spin(2)
2422        dos%fractions(ikpt,iband,isppol,7) = dos%fractions(ikpt,iband,isppol,7) + spin(3)
2423 
2424        shift_b = shift_b + 2*npw_k
2425      end do
2426      ABI_DEALLOCATE(cg_1band)
2427      shift_sk = shift_sk + nband_k*2*npw_k
2428    end do
2429    ABI_DEALLOCATE(cg_1kpt)
2430 
2431    ! Gather all contributions from different processors
2432    if (collect == 1) then
2433      call xmpi_sum(dos%fractions,mpi_enreg%comm_kpt,ierr)
2434      call xmpi_sum(dos%fractions,mpi_enreg%comm_bandfft,ierr)
2435      !below for future use - spinors should not be parallelized for the moment
2436      !if (mpi_enreg%paral_spinor == 1)then
2437      !  call xmpi_sum(dos%fractions,mpi_enreg%comm_spinor,ierr)
2438      !end if
2439    end if
2440 
2441  else
2442    MSG_WARNING('only partial_dos==1 or 2 is coded')
2443  end if
2444 
2445  if (write_procar) close(unit_procar)
2446 
2447  call cwtime(cpu,wall,gflops,"stop")
2448  write(msg,'(2(a,f8.2),a)')" partial_dos_fractions: cpu_time: ",cpu,"[s], walltime: ",wall," [s]"
2449  call wrtout(std_out,msg,"PERS")
2450 
2451 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=informations 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

PARENTS

      outscfcv

CHILDREN

      pawcprj_alloc,pawcprj_free,simp_gen,timab,xmpi_sum

SOURCE

2508 subroutine partial_dos_fractions_paw(dos,cprj,dimcprj,dtset,mcprj,mkmem,mpi_enreg,pawrad,pawtab)
2509 
2510 
2511 !This section has been created automatically by the script Abilint (TD).
2512 !Do not modify the following lines by hand.
2513 #undef ABI_FUNC
2514 #define ABI_FUNC 'partial_dos_fractions_paw'
2515 !End of the abilint section
2516 
2517  implicit none
2518 
2519 !Arguments ------------------------------------
2520 !scalars
2521  integer,intent(in) :: mcprj,mkmem
2522  type(MPI_type),intent(in) :: mpi_enreg
2523  type(dataset_type),intent(in) :: dtset
2524  type(epjdos_t),intent(inout) :: dos
2525 !arrays
2526  integer,intent(in) :: dimcprj(dtset%natom)
2527  type(pawcprj_type),intent(in) :: cprj(dtset%natom,mcprj)
2528  type(pawrad_type),intent(in) :: pawrad(dtset%ntypat)
2529  type(pawtab_type),target,intent(in) :: pawtab(dtset%ntypat)
2530 
2531 !Local variables-------------------------------
2532 !scalars
2533  integer :: bandpp,basis_size,comm_kptband,cplex,fatbands_flag,iat,iatom,iband,ibg,ibsp
2534  integer :: ierr,ikpt,il,ilang,ilmn,iln,im,iorder_cprj,ispinor,isppol,itypat,j0lmn,j0ln
2535  integer :: jl,jlmn,jln,jm,klmn,kln,lmn_size,mbesslang,me_band,me_kpt,my_nspinor
2536  integer :: nband_cprj_k,nband_k,ndosfraction,nprocband,nproc_spkptband,paw_dos_flag,prtdosm
2537  real(dp) :: cpij,one_over_nproc
2538  !character(len=500) :: message
2539 !arrays
2540  integer ,allocatable :: dimcprj_atsph(:)
2541  integer,ABI_CONTIGUOUS pointer :: indlmn(:,:)
2542  real(dp) :: tsec(2)
2543  real(dp),allocatable :: int1(:,:),int2(:,:),int1m2(:,:)
2544  type(pawcprj_type),allocatable :: cprj_k(:,:)
2545 
2546 !******************************************************************************************
2547 
2548  DBG_ENTER("COLL")
2549 
2550  ABI_CHECK(mkmem/=0,"mkmem==0 not supported anymore!")
2551 
2552  fatbands_flag = dos%fatbands_flag
2553  mbesslang = dos%mbesslang
2554  prtdosm = dos%prtdosm
2555  ndosfraction = dos%ndosfraction
2556  paw_dos_flag = dos%paw_dos_flag
2557 
2558 !m-decomposed DOS not compatible with PAW-decomposed DOS
2559  if(prtdosm>=1.and.paw_dos_flag==1) then
2560    MSG_ERROR('m-decomposed DOS not compatible with PAW-decomposed DOS!')
2561  end if
2562 
2563 !Prepare some useful integrals
2564  basis_size=pawtab(1)%basis_size
2565  if (dtset%ntypat>1) then
2566    do itypat=1,dtset%ntypat
2567      basis_size=max(basis_size,pawtab(itypat)%basis_size)
2568    end do
2569  end if
2570  ABI_ALLOCATE(int1  ,(basis_size*(basis_size+1)/2,dtset%natsph))
2571  ABI_ALLOCATE(int2,(basis_size*(basis_size+1)/2,dtset%natsph))
2572  ABI_ALLOCATE(int1m2,(basis_size*(basis_size+1)/2,dtset%natsph))
2573  int1=zero;int2=zero;int1m2=zero
2574  do iat=1,dtset%natsph
2575    iatom=dtset%iatsph(iat)
2576    itypat= dtset%typat(iatom)
2577    do jln=1,pawtab(itypat)%basis_size
2578      j0ln=jln*(jln-1)/2
2579      do iln=1,jln
2580        kln=j0ln+iln
2581        call simp_gen(int1(kln,iat),pawtab(itypat)%phiphj(:,kln),pawrad(itypat))
2582        if (dtset%pawprtdos<2) then
2583          call simp_gen(int2(kln,iat),pawtab(itypat)%tphitphj(:,kln),pawrad(itypat))
2584          int1m2(kln,iat)=int1(kln,iat)-int2(kln,iat)
2585        else
2586          int2(kln,iat)=zero;int1m2(kln,iat)=int1(kln,iat)
2587        end if
2588      end do !iln
2589    end do !jln
2590  end do
2591 
2592 !Antiferro case
2593  if (dtset%nspden==2.and.dtset%nsppol==1.and.dtset%nspinor==1) then
2594    int1m2(:,:)=half*int1m2(:,:)
2595    if (paw_dos_flag==1.or.fatbands_flag==1.or.prtdosm==2) then
2596      int1(:,:)=half*int1(:,:);int2(:,:)=half*int2(:,:)
2597    end if
2598  end if
2599 
2600 !Init parallelism
2601  comm_kptband=mpi_enreg%comm_kptband
2602  nproc_spkptband=xmpi_comm_size(comm_kptband)*mpi_enreg%nproc_spinor
2603  me_kpt=mpi_enreg%me_kpt ; me_band=mpi_enreg%me_band
2604  my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
2605  bandpp=1;if (mpi_enreg%paral_kgb==1) bandpp=mpi_enreg%bandpp
2606 !Check if cprj is distributed over bands
2607  nprocband=my_nspinor*dtset%mband*dtset%mkmem*dtset%nsppol/mcprj
2608  if (nprocband/=mpi_enreg%nproc_band) then
2609    MSG_BUG('wrong mcprj/nproc_band!')
2610  end if
2611 
2612 !Quick hack: in case of parallelism, dos_fractions have already
2613 !  been reduced over MPI processes; they have to be prepared before
2614 !  the next reduction (at the end of the following loop).
2615  if (nproc_spkptband>1) then
2616    one_over_nproc=one/real(nproc_spkptband,kind=dp)
2617 !$OMP  PARALLEL DO COLLAPSE(4) &
2618 !$OMP& DEFAULT(SHARED) PRIVATE(ilang,isppol,iband,ikpt)
2619    do ilang=1,ndosfraction
2620      do isppol=1,dtset%nsppol
2621        do iband=1,dtset%mband
2622          do ikpt=1,dtset%nkpt
2623            dos%fractions(ikpt,iband,isppol,ilang)= &
2624 &           one_over_nproc*dos%fractions(ikpt,iband,isppol,ilang)
2625          end do
2626        end do
2627      end do
2628    end do
2629 !$OMP END PARALLEL DO
2630    if (fatbands_flag==1.or.prtdosm==1.or.prtdosm==2) then
2631 !$OMP  PARALLEL DO COLLAPSE(4) &
2632 !$OMP& DEFAULT(SHARED) PRIVATE(ilang,isppol,iband,ikpt)
2633      do ilang=1,ndosfraction*mbesslang
2634        do isppol=1,dtset%nsppol
2635          do iband=1,dtset%mband
2636            do ikpt=1,dtset%nkpt
2637              dos%fractions_m(ikpt,iband,isppol,ilang)= &
2638 &             one_over_nproc*dos%fractions_m(ikpt,iband,isppol,ilang)
2639            end do
2640          end do
2641        end do
2642      end do
2643 !$OMP END PARALLEL DO
2644    end if
2645  end if
2646 
2647  iorder_cprj=0
2648 
2649 !LOOPS OVER SPINS,KPTS
2650  ibg=0
2651  do isppol =1,dtset%nsppol
2652    do ikpt=1,dtset%nkpt
2653 
2654      nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
2655      if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_kpt)) cycle
2656 
2657      cplex=2;if (dtset%istwfk(ikpt)>1) cplex=1
2658      nband_cprj_k=nband_k/nprocband
2659      ABI_DATATYPE_ALLOCATE(cprj_k,(dtset%natsph,my_nspinor*nband_cprj_k))
2660      ABI_ALLOCATE(dimcprj_atsph,(dtset%natsph))
2661      do iat=1,dtset%natsph
2662        dimcprj_atsph(iat)=dimcprj(dtset%iatsph(iat))
2663      end do
2664      call pawcprj_alloc(cprj_k,0,dimcprj_atsph)
2665      ABI_DEALLOCATE(dimcprj_atsph)
2666 
2667 !    Extract cprj for this k-point.
2668      ibsp=0
2669      do iband=1,nband_cprj_k
2670        do ispinor=1,my_nspinor
2671          ibsp=ibsp+1
2672          do iat=1,dtset%natsph
2673            iatom=dtset%iatsph(iat)
2674            cprj_k(iat,ibsp)%cp(:,:)=cprj(iatom,ibsp+ibg)%cp(:,:)
2675          end do
2676        end do
2677      end do
2678 
2679 !    LOOP OVER ATOMS (natsph_extra is not included on purpose)
2680      do iat=1,dtset%natsph
2681        iatom=dtset%iatsph(iat)
2682        itypat= dtset%typat(iatom)
2683        lmn_size=pawtab(itypat)%lmn_size
2684        indlmn => pawtab(itypat)%indlmn
2685 
2686 !      LOOP OVER BANDS
2687        ibsp=0
2688        do iband=1,nband_k
2689 
2690          if (mod((iband-1)/bandpp,nprocband)/=me_band) cycle
2691          if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband,isppol,me_kpt)) then
2692            ibsp=ibsp+my_nspinor;cycle
2693          end if
2694 
2695          do ispinor=1,my_nspinor
2696            ibsp=ibsp+1
2697 
2698            do ilang=1,mbesslang
2699 
2700              do jlmn=1,lmn_size
2701                jl=indlmn(1,jlmn)
2702                jm=indlmn(2,jlmn)
2703                j0lmn=jlmn*(jlmn-1)/2
2704                do ilmn=1,jlmn
2705                  il=indlmn(1,ilmn)
2706                  im=indlmn(2,ilmn)
2707                  klmn=j0lmn+ilmn
2708                  kln=pawtab(itypat)%indklmn(2,klmn)
2709 
2710                  if (il==ilang-1.and.jl==ilang-1.and.im==jm) then
2711 
2712                    cpij=cprj_k(iat,ibsp)%cp(1,ilmn)*cprj_k(iat,ibsp)%cp(1,jlmn)
2713                    if (cplex==2) cpij=cpij+cprj_k(iat,ibsp)%cp(2,ilmn)*cprj_k(iat,ibsp)%cp(2,jlmn)
2714                    cpij=pawtab(itypat)%dltij(klmn)*cpij
2715 
2716                    dos%fractions(ikpt,iband,isppol,mbesslang*(iat-1)+ilang)=  &
2717 &                   dos%fractions(ikpt,iband,isppol,mbesslang*(iat-1)+ilang) + &
2718 &                   cpij*int1m2(kln,iat)
2719                    if (prtdosm==1) then
2720                      dos%fractions_m(ikpt,iband,isppol,mbesslang**2*(iat-1)+il**2+il+1+im)= &
2721 &                     dos%fractions_m(ikpt,iband,isppol,mbesslang**2*(iat-1)+il**2+il+1+im) + &
2722 &                     cpij*int1m2(kln,iat)
2723                    end if
2724                    if (fatbands_flag==1.or.prtdosm==2) then
2725                      dos%fractions_m(ikpt,iband,isppol,mbesslang**2*(iat-1)+il**2+il+1+im)= &
2726 &                     dos%fractions_m(ikpt,iband,isppol,mbesslang**2*(iat-1)+il**2+il+1+im) + &
2727 &                     cpij*int1(kln,iat)
2728                    end if
2729                    if (paw_dos_flag==1) then
2730                      dos%fractions_paw1(ikpt,iband,isppol,mbesslang*(iat-1)+ilang)=  &
2731 &                     dos%fractions_paw1(ikpt,iband,isppol,mbesslang*(iat-1)+ilang) + &
2732 &                     cpij*int1(kln,iat)
2733                      dos%fractions_pawt1(ikpt,iband,isppol,mbesslang*(iat-1)+ilang)=  &
2734 &                     dos%fractions_pawt1(ikpt,iband,isppol,mbesslang*(iat-1)+ilang) + &
2735 &                     cpij*int2(kln,iat)
2736                    end if
2737 
2738                  end if
2739 
2740                end do !ilmn
2741              end do   !jlmn
2742 
2743            end do ! ilang
2744          end do ! ispinor
2745        end do ! iband
2746 
2747      end do !iatom
2748 
2749      if (mkmem/=0) ibg = ibg + my_nspinor*nband_cprj_k
2750      call pawcprj_free(cprj_k)
2751      ABI_DATATYPE_DEALLOCATE(cprj_k)
2752    end do ! ikpt
2753  end do ! isppol
2754 
2755  ABI_DEALLOCATE(int1)
2756  ABI_DEALLOCATE(int2)
2757  ABI_DEALLOCATE(int1m2)
2758 
2759 !Reduce data in case of parallelism
2760  call timab(48,1,tsec)
2761  call xmpi_sum(dos%fractions,comm_kptband,ierr)
2762  if (prtdosm>=1.or.fatbands_flag==1) then
2763    call xmpi_sum(dos%fractions_m,comm_kptband,ierr)
2764  end if
2765  if (paw_dos_flag==1) then
2766    call xmpi_sum(dos%fractions_paw1,comm_kptband,ierr)
2767    call xmpi_sum(dos%fractions_pawt1,comm_kptband,ierr)
2768  end if
2769  call timab(48,2,tsec)
2770  if (mpi_enreg%paral_spinor==1) then
2771    call xmpi_sum(dos%fractions,mpi_enreg%comm_spinor,ierr)
2772    if (prtdosm>=1.or.fatbands_flag==1) then
2773      call xmpi_sum(dos%fractions_m,mpi_enreg%comm_spinor,ierr)
2774    end if
2775    if (paw_dos_flag==1) then
2776      call xmpi_sum(dos%fractions_paw1, mpi_enreg%comm_spinor,ierr)
2777      call xmpi_sum(dos%fractions_pawt1, mpi_enreg%comm_spinor,ierr)
2778    end if
2779  end if
2780 
2781 !Averaging: A quick hack for m-decomposed LDOS:
2782 !BA: not valid in presence of spin-orbit coupling  !
2783  if (prtdosm==1.and.fatbands_flag==0) then
2784 !  if pawfatbnd is activated, one think in the cubic harmonics basis
2785 !  whereas prtdosm=1 is in the spherical harmonics basis.
2786 !  the following trick is done in order to have everything
2787 !  in the complex spherical basis (not useful for pawfatbnd if we want to
2788 !  have e.g t2g and eg d-orbitals).
2789    do iat=1,dtset%natsph
2790      do il = 0, mbesslang-1
2791        do im = 1, il
2792          dos%fractions_m(:,:,:,mbesslang**2*(iat-1)+il**2+il+1+im) = &
2793 &         (dos%fractions_m(:,:,:,mbesslang**2*(iat-1)+il**2+il+1+im) + &
2794 &         dos%fractions_m(:,:,:,mbesslang**2*(iat-1)+il**2+il+1-im))/2
2795          dos%fractions_m(:,:,:,mbesslang**2*(iat-1)+il**2+il+1-im) = &
2796 &         dos%fractions_m(:,:,:,mbesslang**2*(iat-1)+il**2+il+1+im)
2797        end do
2798      end do
2799    end do !iatom
2800  end if
2801 
2802  DBG_EXIT("COLL")
2803 
2804 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

PARENTS

      outscfcv

CHILDREN

      cwtime,wrtout

SOURCE

1572 subroutine prtfatbands(dos,dtset,ebands,fildata,pawfatbnd,pawtab)
1573 
1574 
1575 !This section has been created automatically by the script Abilint (TD).
1576 !Do not modify the following lines by hand.
1577 #undef ABI_FUNC
1578 #define ABI_FUNC 'prtfatbands'
1579 !End of the abilint section
1580 
1581  implicit none
1582 
1583 !Arguments ------------------------------------
1584 !scalars
1585  integer,intent(in) :: pawfatbnd
1586  type(epjdos_t),intent(in) :: dos
1587  type(ebands_t),intent(in) :: ebands
1588  type(dataset_type),intent(in) :: dtset
1589  character(len=fnlen),intent(in) :: fildata
1590 !arrays
1591  type(pawtab_type),intent(in) :: pawtab(dtset%ntypat)
1592 
1593 !Local variables-------------------------------
1594 !scalars
1595  integer :: iall,il,iat,natsph,inbfatbands,iband,mband,ixfat,isppol,nkpt,lmax,ll,mm
1596  integer :: ikpt,nband_k,ndosfraction,mbesslang
1597  real(dp) :: xfatband,cpu,wall,gflops
1598  character(len=1) :: tag_l,tag_1m,tag_is
1599  character(len=2) :: tag_2m
1600  character(len=10) :: tag_il,tag_at,tag_grace
1601  character(len=1500) :: message
1602  character(len=fnlen) :: tmpfil
1603  type(atomdata_t) :: atom
1604 !arrays
1605  integer,allocatable :: unitfatbands_arr(:,:)
1606  real(dp),allocatable :: eigenvalues(:,:,:)
1607  character(len=2) :: symbol
1608 
1609 !*************************************************************************
1610 
1611  DBG_ENTER("COLL")
1612 
1613  ndosfraction = dos%ndosfraction; mbesslang = dos%mbesslang
1614 
1615  if(dos%prtdosm.ne.0) then
1616    write(message,'(3a)')&
1617 &   'm decomposed dos is activated',ch10, &
1618 &   'Action: deactivate it with prtdosm=0 !'
1619    MSG_ERROR(message)
1620  end if
1621 
1622  if(dtset%nspinor==2) then
1623    MSG_WARNING("Fatbands are not yet available in the case nspinor==2!")
1624  end if
1625 
1626  ABI_CHECK(allocated(dos%fractions_m), "dos%fractions_m is not allocated!")
1627 
1628  natsph=dtset%natsph
1629  nkpt=dtset%nkpt
1630  mband=dtset%mband
1631 
1632  if(natsph>1000) then
1633    write(message,'(3a)')&
1634 &   'Too big number of fat bands!',ch10, &
1635 &   'Action: decrease natsph in input file !'
1636    MSG_ERROR(message)
1637  end if
1638 
1639 !--------------  PRINTING IN LOG
1640  call cwtime(cpu, wall, gflops, "start")
1641  write(message,'(a,a,a,a,i5,a,a,1000i5)') ch10," ***** Print of fatbands activated ****** ",ch10,&
1642 & "  Number of atom: natsph = ",natsph,ch10, &
1643 & "  atoms  are             = ",(dtset%iatsph(iat),iat=1,natsph)
1644  call wrtout(std_out,message,'COLL')
1645  call wrtout(ab_out,message,'COLL')
1646  iall=0;inbfatbands=0
1647 
1648  if(pawfatbnd==1) then
1649    inbfatbands=mbesslang-1
1650    write(message,'(3a)')"  (fatbands are in eV and are given for each value of L)",ch10
1651  else if(pawfatbnd==2) then
1652    write(message,'(3a)')"  (fatbands are in eV and are given for each value of L and M)",ch10
1653    inbfatbands=(mbesslang-1)**2
1654  end if
1655  call wrtout(std_out,message,'COLL')
1656  call wrtout(ab_out,message,'COLL')
1657 
1658  write(message,'(a,e12.5,a,e12.5,a)') "  Fermi energy is ",ebands%fermie*Ha_eV," eV = ",ebands%fermie," Ha"
1659  call wrtout(std_out,message,'COLL')
1660 
1661 !--------------  OPEN AND NAME FILES FOR FATBANDS
1662  ABI_ALLOCATE(unitfatbands_arr,(natsph*inbfatbands,dtset%nsppol))
1663  unitfatbands_arr = -3
1664 
1665  do iat=1,natsph
1666    lmax=(pawtab(dtset%typat(dtset%iatsph(iat)))%l_size-1)/2
1667    call int2char4(dtset%iatsph(iat),tag_at)
1668    ABI_CHECK((tag_at(1:1)/='#'),'Bug: string length too short!')
1669    call atomdata_from_znucl(atom,dtset%znucl(dtset%typat(dtset%iatsph(iat))))
1670    symbol = atom%symbol
1671    do il=1,inbfatbands
1672      iall=iall+1
1673      ll=int(sqrt(float(il-1)))  ! compute l
1674      if(ll.le.lmax) then  ! print only angular momentum included in the PAW data
1675        do isppol=1,dtset%nsppol
1676          write(tag_is,'(i1)')isppol
1677          if(pawfatbnd==1) then
1678            call int2char4(il-1,tag_il)
1679            ABI_CHECK((tag_il(1:1)/='#'),'Bug: string length too short!')
1680            tmpfil = trim(fildata)// &
1681 &           '_at'//trim(tag_at)//'_'//trim(adjustl(symbol))//'_is'//tag_is//'_l'//trim(tag_il)
1682          else if (pawfatbnd==2) then
1683            write(tag_l,'(i1)') ll
1684            mm=il-(ll**2+ll+1)      ! compute m
1685            if(mm<0) write(tag_2m,'(i2)') mm
1686            if(mm>=0) write(tag_1m,'(i1)') mm
1687            if(mm<0) tmpfil = trim(fildata)// &
1688 &           '_at'//trim(tag_at)//'_'//trim(adjustl(symbol))//'_is'//tag_is//'_l'//tag_l//'_m'//tag_2m
1689            if(mm>=0) tmpfil = trim(fildata)// &
1690 &           '_at'//trim(tag_at)//'_'//trim(adjustl(symbol))//'_is'//tag_is//'_l'//tag_l//'_m+'//tag_1m
1691          end if
1692          !unitfatbands_arr(iall,isppol)=tmp_unit+100+iall-1+(natsph*inbfatbands)*(isppol-1)
1693          !open (unit=unitfatbands_arr(iall,isppol),file=trim(tmpfil),status='unknown',form='formatted')
1694          if (open_file(tmpfil, message, newunit=unitfatbands_arr(iall,isppol), status='unknown',form='formatted') /= 0) then
1695            MSG_ERROR(message)
1696          end if
1697 
1698          write(message,'(a,a,a,i4)') 'opened file : ', trim(tmpfil), ' unit', unitfatbands_arr(iall,isppol)
1699          call wrtout(std_out,message,'COLL')
1700          write(message,'(9a)') "# ",ch10,"# ABINIT package : FATBAND file ", ch10,&
1701 &         "# It contains, for each band: the eigenvalues in eV (and the character of the band) as a function of the k-point",&
1702 &         ch10,"# This file can be read with xmgrace (http://plasma-gate.weizmann.ac.il/Grace/)  ",ch10,"#  "
1703          write(unitfatbands_arr(iall,isppol), "(a)")trim(message)
1704          do iband=1,mband
1705            call int2char4(iband-1,tag_grace)
1706            ABI_CHECK((tag_grace(1:1)/='#'),'Bug: string length too short!')
1707            write(message,'(16a)') ch10,"@    s",trim(tag_grace)," line color 1",&
1708 &           ch10,"@    s",trim(tag_grace)," errorbar color 2",&
1709 &           ch10,"@    s",trim(tag_grace)," errorbar riser linewidth 5.0", &
1710 &           ch10,"@    s",trim(tag_grace)," errorbar linestyle 0"
1711            write(unitfatbands_arr(iall,isppol), "(a)")trim(message)
1712          end do  !iband
1713          write(unitfatbands_arr(iall,isppol), '(a,a)') ch10,'@type xydy'
1714        end do   ! isppol
1715      end if ! ll=<lmax
1716    end do   ! il
1717  end do  ! iat
1718 
1719  if(iall.ne.(natsph*inbfatbands)) then
1720    MSG_ERROR("error1 ")
1721  end if
1722 
1723 
1724 
1725 !--------------  WRITE FATBANDS IN FILES
1726  if (pawfatbnd>0) then
1727    ! Store eigenvalues with nkpt as first dimension for efficiency reasons
1728    ABI_ALLOCATE(eigenvalues,(nkpt,mband,dtset%nsppol))
1729    do isppol=1,dtset%nsppol
1730      do ikpt=1,nkpt
1731        nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
1732        do iband=1,mband
1733          eigenvalues(ikpt,iband,isppol)= ebands%eig(iband, ikpt, isppol) - ebands%fermie
1734        end do
1735      end do
1736    end do
1737    iall=0
1738    do iat=1,natsph
1739      lmax=(pawtab(dtset%typat(dtset%iatsph(iat)))%l_size-1)/2
1740      do il=1,inbfatbands
1741        iall=iall+1
1742        ll=int(sqrt(float(il-1)))
1743        if(ll.le.lmax) then
1744          do isppol=1,dtset%nsppol
1745            do iband=1,mband
1746              write(unitfatbands_arr(iall,isppol),'(a,a,i8)') ch10,"# BAND number :",iband
1747              do ikpt=1,nkpt
1748                if(pawfatbnd==1) then
1749                  xfatband=0.d0
1750                  do ixfat=(il-1)**2+1,il**2
1751                    xfatband=xfatband+dos%fractions_m(ikpt,iband,isppol,(iat-1)*mbesslang**2+ixfat)
1752                  end do ! ixfat
1753                else if (pawfatbnd==2) then
1754                  xfatband=dos%fractions_m(ikpt,iband,isppol,(iat-1)*mbesslang**2+il)
1755                end if
1756                write(unitfatbands_arr(iall,isppol),'(i5,e20.5,e20.5)')&
1757                  ikpt-1,eigenvalues(ikpt,iband,isppol)*Ha_eV,xfatband
1758              end do ! ikpt
1759            end do  !iband
1760            write(unitfatbands_arr(iall,isppol),'(a)') '&'
1761            !close(unitfatbands_arr(iall,isppol))
1762          end do  !isppol
1763        end if
1764      end do ! il
1765    end do ! iat
1766    ABI_DEALLOCATE(eigenvalues)
1767  end if
1768 
1769  do isppol=1,size(unitfatbands_arr, dim=2)
1770    do iat=1,size(unitfatbands_arr, dim=1)
1771      if (unitfatbands_arr(iat, isppol) /= -3) close (unitfatbands_arr(iat, isppol))
1772    end do
1773  end do
1774 
1775  ABI_DEALLOCATE(unitfatbands_arr)
1776 
1777  call cwtime(cpu,wall,gflops,"stop")
1778  write(message,'(2(a,f8.2),a)')" prtfatbands: cpu_time: ",cpu,"[s], walltime: ",wall," [s]"
1779  call wrtout(std_out,message,"PERS")
1780 
1781  DBG_EXIT("COLL")
1782 
1783 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.

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.

PARENTS

      m_cut3d,partial_dos_fractions

CHILDREN

      cwtime,wrtout

SOURCE

 934 subroutine recip_ylm (bess_fit, cg_1band, istwfk, mpi_enreg, nradint, nradintmax, mlang,&
 935 &  mpw, natsph, typat_extra, mlang_type, npw_k, nspinor, ph3d, prtsphere, rint, rmax,&
 936 &  rc_ylm, sum_1ll_1atom, sum_1lm_1atom, ucvol, ylm_k, znucl_sph)
 937 
 938 
 939 !This section has been created automatically by the script Abilint (TD).
 940 !Do not modify the following lines by hand.
 941 #undef ABI_FUNC
 942 #define ABI_FUNC 'recip_ylm'
 943 !End of the abilint section
 944 
 945  implicit none
 946 
 947 !Arguments ------------------------------------
 948 !scalars
 949  integer,intent(in) :: istwfk,mlang,mpw,natsph,npw_k,nradintmax
 950  integer,intent(in) :: nspinor
 951  integer,intent(in) :: prtsphere,rc_ylm
 952  real(dp),intent(in) :: ucvol
 953 !arrays
 954  integer,intent(in) :: nradint(natsph),typat_extra(natsph),mlang_type(:)
 955  real(dp),intent(in) :: bess_fit(mpw,nradintmax,mlang),cg_1band(:,:) !(2,my_nspinor*npw_k)
 956  real(dp),intent(in) :: ph3d(2,npw_k,natsph),rint(nradintmax)
 957  real(dp),intent(in) :: rmax(natsph),ylm_k(npw_k,mlang*mlang)
 958  real(dp),intent(in) :: znucl_sph(natsph)
 959  type(MPI_type),intent(in) :: mpi_enreg
 960  real(dp),intent(out) :: sum_1ll_1atom(nspinor**2, mlang, natsph)
 961  real(dp),intent(out) :: sum_1lm_1atom(nspinor**2, mlang*mlang,natsph)
 962 
 963 !Local variables-------------------------------
 964 !scalars
 965  integer :: ilm,iat,ipw,ixint,ll,mm,il,jlm,ierr,lm_size,itypat
 966  integer :: ispinor
 967  integer :: ipauli, is, isp
 968  integer :: my_nspinor
 969  real(dp),parameter :: invsqrt2=one/sqrt2
 970  real(dp) :: sum_all, dr, fact
 971  type(atomdata_t) :: atom
 972  character(len=500) :: msg
 973 !arrays
 974  integer :: ilang(mlang**2)
 975  real(dp) :: c1(2),c2(2)
 976  real(dp) :: sum_1atom(natsph),sum_1ll(mlang),sum_1lm(mlang**2)
 977  real(dp) :: func(nradintmax)
 978  complex(dpc) :: vect(npw_k)
 979  complex(dpc),allocatable :: tmppsia(:,:),tmppsim(:,:),dotc(:)
 980  integer, allocatable :: ispinors(:)
 981  complex(dpc),allocatable :: values(:,:,:,:)
 982 
 983 ! *************************************************************************
 984 
 985  ! Workspace array (used to reduce the number of MPI communications)
 986  ! One could reduce a bit the memory requirement by using non-blocking operations ...
 987  ABI_STAT_MALLOC(values, (nradintmax, nspinor, mlang**2, natsph), ierr)
 988  ABI_CHECK(ierr==0, "oom in values")
 989  values = czero
 990 
 991  my_nspinor = max(1,nspinor/mpi_enreg%nproc_spinor)
 992  ABI_ALLOCATE(tmppsia, (npw_k,my_nspinor))
 993  ABI_ALLOCATE(tmppsim, (npw_k,my_nspinor))
 994  ABI_ALLOCATE(dotc, (my_nspinor))
 995  ABI_ALLOCATE(ispinors, (my_nspinor))
 996  if (my_nspinor == 2) then
 997    ispinors(1) = 1
 998    ispinors(2) = 2
 999  else
1000    ispinors(1) = mpi_enreg%me_spinor+1
1001  end if
1002 
1003  sum_1lm_1atom = zero
1004 
1005  do ll=0,mlang-1
1006    do mm=-ll,ll
1007      ilm = (ll+1)**2-ll+mm
1008      ilang(ilm) = ll+1
1009    end do
1010  end do
1011 
1012  ! Big loop on all atoms
1013  do iat=1,natsph
1014    itypat = typat_extra(iat)
1015    lm_size = mlang_type(itypat) ** 2
1016    dr = rmax(iat) / (nradint(iat)-1)
1017 
1018    ! u(G) e^{i(k+G).Ra}
1019    ! tmppsia = Temporary array for part which depends only on iat
1020    do ispinor=1,my_nspinor
1021      do ipw=1,npw_k
1022        tmppsia(ipw,ispinor) = dcmplx(cg_1band(1,ipw+(ispinor-1)*npw_k),cg_1band(2,ipw+(ispinor-1)*npw_k)) &
1023 &            * dcmplx(ph3d(1,ipw,iat), ph3d(2,ipw,iat))
1024      end do
1025    end do
1026 
1027    ! tmppsim = temporary arrays for part of psi which does not depend on ixint = tmppsia * ylm.
1028    ! could remove this intermediate array to save memory...
1029    ! u(G) Y_LM^*(k+G) e^{i(k+G).Ra}
1030    ! Take into account the fact that ylm_k are REAL spherical harmonics, see initylmg.f
1031    ! For time-reversal states, detailed treatment show that only the real or imaginary
1032    ! part of tmppsia is needed here, depending on l being even or odd: only one of the coef is 1, the other 0
1033    do ilm=1,lm_size
1034      il = ilang(ilm)
1035      ll = ilang(ilm) - 1
1036      mm = ilm - (ll+1)**2 + ll
1037 
1038      select case (rc_ylm)
1039      case (1)
1040        ! to get PDOS for real spherical harmonics, simply multiply here by ylm_k
1041        do ispinor=1,my_nspinor
1042          do ipw=1,npw_k
1043            tmppsim(ipw,ispinor) = tmppsia(ipw,ispinor) * ylm_k(ipw,ilm)
1044          end do
1045        end do
1046 
1047        ! Handle time-reversal
1048        ! TODO: check if time reversal with spinors is special and may need some treatment.
1049        !  normally SOC will simply impose istwfk 1 in appropriate cases (kptopt 4).
1050        if (istwfk /= 1) then
1051          if (mod(ll, 2) == 0) then
1052             tmppsim(:,:) = dcmplx(real(tmppsim(:,:)),zero)
1053          else
1054             tmppsim(:,:) = dcmplx(aimag(tmppsim(:,:)),zero)
1055          end if
1056 ! f2008 version:
1057 !         if (mod(ll, 2) == 0) then
1058 !            tmppsim(:,:)%im = zero
1059 !         else
1060 !            tmppsim(:,:)%re = tmppsim(:,:)%im
1061 !            tmppsim(:,:)%im = zero
1062 !         end if
1063        end if
1064 
1065      case (2)
1066        ! to get PDOS for complex spherical harmonics, build linear combination of real ylm_k
1067        jlm = (ll+1)**2-ll-mm ! index of (l, -m)
1068        if (mm == 0) then
1069          vect(:) = dcmplx(ylm_k(1:npw_k,ilm),zero)
1070        else if (mm > 0) then
1071           !vect(1,:) =  invsqrt2 * ylm_k(1:npw_k,ilm) * (-1)**mm
1072           !vect(2,:) = +invsqrt2 * ylm_k(1:npw_k,jlm) * (-1)**mm
1073           c1 = sy(ll, mm, mm)
1074           c2 = sy(ll,-mm, mm)
1075           vect(:) = dcmplx(c1(1) * ylm_k(1:npw_k,ilm) + c2(1) * ylm_k(1:npw_k,jlm), &
1076 &                         c1(2) * ylm_k(1:npw_k,ilm) + c2(2) * ylm_k(1:npw_k,jlm))
1077 
1078        else if (mm < 0) then
1079           !vect(1,:) =  invsqrt2 * ylm_k(1:npw_k,jlm) !* (-1)**mm
1080           !vect(2,:) = -invsqrt2 * ylm_k(1:npw_k,ilm) !* (-1)**mm
1081           c1 = sy(ll, mm,  mm)
1082           c2 = sy(ll,-mm,  mm)
1083           vect(:) = dcmplx(c1(1) * ylm_k(1:npw_k,ilm) + c2(1) * ylm_k(1:npw_k,jlm),&
1084 &                         c1(2) * ylm_k(1:npw_k,ilm) + c2(2) * ylm_k(1:npw_k,jlm))
1085        end if
1086        vect(:) = dcmplx(real(vect(:)), -aimag(vect(:)))
1087        !vect(:)%im = -vect(:)%im
1088 
1089        if (istwfk == 1) then
1090          do ispinor=1,my_nspinor
1091            do ipw=1,npw_k
1092              tmppsim(ipw,ispinor) = tmppsia(ipw,ispinor) * vect(ipw)
1093            end do
1094          end do
1095        else
1096          ! Handle time-reversal
1097          if (mod(ll, 2) == 0) then
1098            do ispinor=1,my_nspinor
1099              do ipw=1,npw_k
1100                tmppsim(ipw,ispinor) = real(tmppsia(ipw,ispinor)) * vect(ipw)
1101              end do
1102            end do
1103          else
1104            do ispinor=1,my_nspinor
1105              do ipw=1,npw_k
1106                tmppsim(ipw,ispinor) = aimag(tmppsia(ipw,ispinor)) * vect(ipw)
1107              end do
1108            end do
1109          end if
1110        end if
1111 
1112      case default
1113        MSG_ERROR("Wrong value for rc_ylm")
1114      end select
1115 
1116      ! 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 $
1117      ! or more general spinor case integral
1118      !    $ \int_0^{rc} dr r**2 dotc^*_s   \sigma^x_{ss'}   dotc_{s'}
1119      ! where   dotc_s = \sum_G u_s (G) Y_LM^*(k+G) e^{i(k+G).Ra} j_L(|k+G| r)
1120      do ixint=1,nradint(iat)
1121        dotc = czero
1122        do ispinor=1, my_nspinor
1123          do ipw=1,npw_k
1124            dotc(ispinor) = dotc(ispinor) + bess_fit(ipw, ixint, il) * tmppsim(ipw, ispinor)
1125          end do
1126        end do
1127        if (istwfk /= 1) then
1128          dotc = two * dotc
1129          if (istwfk == 2 .and. mpi_enreg%me_g0 == 1) then
1130            dotc(:) = dotc(:) - bess_fit(1, ixint, il) * tmppsim(1, :)
1131          end if
1132        end if
1133 
1134        ! Store results to reduce number of xmpi_sum calls if MPI
1135        do ispinor=1, my_nspinor
1136          values(ixint, ispinors(ispinor), ilm, iat) = dotc(ispinor)
1137        end do
1138      end do ! ixint
1139 
1140    end do ! ilm
1141  end do ! iat
1142 
1143  ABI_DEALLOCATE(tmppsia)
1144  ABI_DEALLOCATE(tmppsim)
1145  ABI_DEALLOCATE(dotc)
1146  ABI_DEALLOCATE(ispinors)
1147 
1148  ! Collect results in comm_pw (data are distributed over plane waves)
1149  call xmpi_sum(values, mpi_enreg%comm_bandfft, ierr)
1150 ! ! Collect results in mpi_enreg%comm_spinor (data are distributed over spinor components)
1151  call xmpi_sum(values, mpi_enreg%comm_spinor, ierr)
1152 
1153  ! Multiply by r**2 and take norm, integrate
1154  do iat=1,natsph
1155    itypat = typat_extra(iat)
1156    lm_size = mlang_type(itypat) ** 2
1157    do ilm=1,lm_size
1158      do ipauli=0,nspinor**2-1
1159        do ixint=1,nradint(iat)
1160          func(ixint) = zero
1161          do is=1,nspinor
1162            do isp=1,nspinor
1163              func(ixint) =  func(ixint) + real(conjg(values(ixint, is, ilm, iat))*pauli_mat(is,isp,ipauli)*&
1164 &                                                    values(ixint, isp, ilm, iat))
1165            end do
1166          end do
1167          func(ixint) = rint(ixint)**2 * func(ixint)
1168        end do
1169        ! Here I should treat the case in which the last point /= rcut
1170        ! NB: indexing is from 1 not 0 for spin matrix components
1171        sum_1lm_1atom(ipauli+1, ilm, iat) = simpson(dr, func(1:nradint(iat)))
1172      end do
1173    end do
1174  end do
1175 
1176  ! Normalize with unit cell volume and include 4pi term coming from Rayleigh expansion.
1177  fact = four_pi**2 / ucvol
1178  sum_1lm_1atom = fact * sum_1lm_1atom
1179  sum_1ll_1atom = zero
1180  do iat=1,natsph
1181    itypat = typat_extra(iat)
1182    lm_size = mlang_type(itypat) ** 2
1183    do ilm=1,lm_size
1184      il = ilang(ilm)
1185      sum_1ll_1atom(:,il, iat) = sum_1ll_1atom(:,il, iat) + sum_1lm_1atom(:,ilm, iat)
1186    end do
1187  end do
1188 
1189  ABI_FREE(values)
1190 
1191  ! Output
1192  if (prtsphere == 1) then
1193    sum_1ll = zero
1194    sum_1lm = zero
1195    sum_1atom = zero
1196    do iat=1,natsph
1197      sum_1atom(iat) = sum(sum_1lm_1atom(1,:,iat))
1198      sum_1ll(:)=sum_1ll(:)+sum_1ll_1atom(1,:,iat)
1199      sum_1lm(:)=sum_1lm(:)+sum_1lm_1atom(1,:,iat)
1200    end do
1201    sum_all = sum(sum_1atom)
1202 
1203    if (rc_ylm == 1) msg = " Angular analysis (real spherical harmonics)"
1204    if (rc_ylm == 2) msg = " Angular analysis (complex spherical harmonics)"
1205    call wrtout(std_out, msg)
1206    do iat=1,natsph
1207      call atomdata_from_znucl(atom, znucl_sph(iat))
1208      call wrtout(std_out, " ")
1209      write(msg,'(a,i3,a,a,a,f10.6)' )' Atom # ',iat, ' is  ',  atom%symbol,', in-sphere charge =',sum_1atom(iat)
1210      call wrtout(std_out, msg)
1211      do ll=0,mlang-1
1212        write(msg,'(a,i1,a,f9.6,a,9f6.3)' )&
1213 &       ' l=',ll,', charge=',sum_1ll_1atom(1,ll+1,iat),&
1214 &       ', m=-l,l splitting:',sum_1lm_1atom(1,1+ll**2:(ll+1)**2,iat)
1215        call wrtout(std_out, msg)
1216      end do ! ll
1217    end do ! iat
1218    write(msg,'(a,a)') ch10,' Sum of angular contributions for all atomic spheres '
1219    call wrtout(std_out, msg)
1220    do ll=0,mlang-1
1221      write(msg,'(a,i1,a,f9.6,a,f9.6)' )&
1222 &     ' l=',ll,', charge =',sum_1ll(ll+1),' proportion =',sum_1ll(ll+1)/sum_all
1223      call wrtout(std_out, msg)
1224    end do
1225    write(msg,'(a,a,f10.6)' ) ch10,' Total over all atoms and l=0 to 4 :',sum_all
1226    call wrtout(std_out, msg)
1227    call wrtout(std_out, " ")
1228  end if
1229 
1230 contains
1231 
1232  function sy(ll, mm, mp)
1233    use  m_paw_sphharm, only : ys
1234    ! Computes the matrix element <Slm|Ylm'>
1235 
1236 !This section has been created automatically by the script Abilint (TD).
1237 !Do not modify the following lines by hand.
1238 #undef ABI_FUNC
1239 #define ABI_FUNC 'sy'
1240 !End of the abilint section
1241 
1242    integer,intent(in) :: ll,mm, mp
1243 
1244    real(dp) :: sy(2)
1245    complex(dpc) :: ys_val
1246 
1247    ! Computes the matrix element <Yl'm'|Slm>
1248    call ys(ll,mp,ll,mm,ys_val)
1249    !call ys(ll,mm,ll,mp,ys_val)
1250    sy(1) = real(ys_val)
1251    sy(2) = -aimag(ys_val)
1252 
1253  end function sy
1254 
1255 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

PARENTS

      m_epjdos

CHILDREN

      cwtime,wrtout

SOURCE

1495 subroutine sphericaldens(fofg,gnorm,nfft,rmax,sphfofg)
1496 
1497 
1498 !This section has been created automatically by the script Abilint (TD).
1499 !Do not modify the following lines by hand.
1500 #undef ABI_FUNC
1501 #define ABI_FUNC 'sphericaldens'
1502 !End of the abilint section
1503 
1504  implicit none
1505 
1506 !Arguments ------------------------------------
1507 !scalars
1508  integer,intent(in) :: nfft
1509  real(dp),intent(in) :: rmax
1510 !arrays
1511  real(dp),intent(in) :: fofg(2,nfft),gnorm(nfft)
1512  real(dp),intent(out) :: sphfofg(2,nfft)
1513 
1514 !Local variables-------------------------------
1515 !scalars
1516  integer :: ifft
1517  real(dp) :: factor,int0yy,rmax_2pi,yy
1518 
1519 ! *************************************************************************
1520 
1521  rmax_2pi=two_pi*rmax
1522  factor=four_pi/(two_pi)**3
1523 
1524  do ifft=1,nfft
1525    if(abs(gnorm(ifft)) < tol12)then
1526      sphfofg(1,ifft)=fofg(1,ifft)*four_pi*third*rmax**3
1527      sphfofg(2,ifft)=fofg(2,ifft)*four_pi*third*rmax**3
1528    else
1529      yy=gnorm(ifft)*rmax_2pi
1530      int0yy=factor*(sin(yy)-yy*cos(yy))/(gnorm(ifft)**3)
1531      sphfofg(1,ifft)=fofg(1,ifft)*int0yy
1532      sphfofg(2,ifft)=fofg(2,ifft)*int0yy
1533    end if
1534  end do
1535 
1536 end subroutine sphericaldens