[Q-e-developers] extension of PP/src/sumpdos.f90
Guido Fratesi
guido.fratesi at unimi.it
Thu Apr 7 15:48:32 CEST 2016
Dear QE developers,
the program PP/src/sumpdos.f90 currently allows one to sum projected
density of states, integrated over the full Brillouin zone, as computed
by PP/src/projwfc.x.
I have modified it slightly so as to sum also the k-resolved DOS files
(always from projwfc, with kresolveddos=.true.).
[notice that the pyton script in PP/tools/sum_states.py is useless in
this respect because it also sums over k]
I tested it and it works both for standard and k-resolved files.
The output format has been made more similar to that of the files
written by projwfc.x and read by the program.
Is it sufficient to attach it here, so as to include it in the QE
distribution?
Thank you for your kind attention.
Best regards,
Guido
--
Guido Fratesi
Dipartimento di Fisica
Universita` degli Studi di Milano
Via Celoria 16, 20133 Milano, Italy
Phone: +39 02 503 17348
email: guido.fratesi at unimi.it
web: https://sites.google.com/site/guidofratesi/
-------------- next part --------------
!
! Copyright (C) 2005 Andrea Ferretti
! Modified 2016 Guido Fratesi
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
PROGRAM sumpdos
!
IMPLICIT NONE
!
! AUTHOR: Andrea Ferretti
! edited: Guido Fratesi (to sum K-resolved DOS files)
!
! this program reads and sum pdos from different
! files (which are related to different atoms)
!
! file names are read from stdin
! USAGE: sumpdos <file1> ... <fileN>
!
INTEGER :: ngrid ! dimension of the energy grid
INTEGER :: nfile ! number of files to sum
INTEGER :: nspin ! number of spin_component
CHARACTER(256), ALLOCATABLE :: file(:) ! names of the files to sum
CHARACTER(256) :: filein
CHARACTER(10) :: cdum, str1, str2, str3
LOGICAL :: exist, kresolveddos
REAL :: efermi = 0.0d0 ! translate the input grid
REAL, ALLOCATABLE :: pdos(:,:,:)
REAL, ALLOCATABLE :: egrid(:)
REAL, ALLOCATABLE :: mysum(:,:)
INTEGER :: ios, ierr, iarg, ie, isp, ifile, i
INTEGER :: ik, nk, iun
!**************************************************************
! User should supply input values here
!
efermi = 0.0d0
!**************************************************************
!
! get the number of arguments (i.e. the number of files)
!
nfile = command_argument_count ()
IF ( nfile == 0 ) THEN
WRITE(0,"( 'No file to sum' )")
STOP
ENDIF
CALL get_command_argument ( 1, str1 )
!
SELECT CASE ( trim(str1) )
CASE ( "-h" )
!
! write the manual
!
WRITE(0,"(/,'USAGE: sumpdos [-h] [-f <filein>] [<file1> ... <fileN>]', /, &
&' Sum the pdos from the file specified in input and write the sum ', /, &
&' to stdout', /, &
&' -h : write this manual',/, &
&' -f <filein> : takes the list of pdos files from <filein> ', /, &
&' (one per line) instead of command line',/, &
&' <fileM> : the M-th pdos file', &
& / )")
STOP
!
CASE ( "-f" )
!
! read file names from file
!
CALL get_command_argument ( 2, filein )
IF ( len_trim(filein) == 0 ) CALL errore('sumpdos','provide filein name',2)
INQUIRE( FILE=trim(filein), EXIST=exist )
IF (.not. exist) CALL errore('sumpdos','file '//trim(filein)//' does not exist',3)
OPEN( 10, FILE=trim(filein), IOSTAT=ios )
IF (ios/=0) CALL errore('sumpdos','opening '//trim(filein),abs(ios))
!
! get the number of non-empty lines in the file
! (which is assumed to be the number of files to sum)
!
ios = 0
nfile = 0
!
DO WHILE ( ios == 0 )
nfile = nfile + 1
READ(10, *, IOSTAT=ios ) cdum
IF ( ios ==0 .and. len_trim(cdum)==0 ) nfile = nfile -1
ENDDO
nfile = nfile -1
!
IF (nfile ==0 ) CALL errore('sumpdos','no file to sum in '//trim(filein),4)
!
ALLOCATE( file(nfile), STAT=ierr )
IF (ierr/=0) CALL errore('sumpdos','allocating FILE',abs(ierr))
!
REWIND(10)
DO i = 1, nfile
file(i) = ' '
DO WHILE( len_trim(file(i)) == 0 )
READ(10,*, IOSTAT=ios) file(i)
IF (ios /=0 ) CALL errore('sumpdos','reading from '//trim(filein),i)
ENDDO
ENDDO
CASE DEFAULT
!
! get the names of the files
!
ALLOCATE( file(nfile), STAT=ierr )
IF (ierr/=0) CALL errore('sumpdos','allocating FILE',abs(ierr))
DO iarg = 1, nfile
CALL get_command_argument ( iarg, file(iarg) )
ENDDO
END SELECT
!
! open the first file and get data about spin
! and grid dimensions
!
INQUIRE( FILE=trim(file(1)), EXIST=exist )
IF (.not. exist) CALL errore('sumpdos','file '//trim(file(1))//' does not exist',3)
!
WRITE(0,"('Reading dimensions from file: ',a)") trim(file(1))
!
OPEN(10, FILE=trim(file(1)), IOSTAT=ios)
IF (ios/=0) CALL errore("sumpdos", "error opening "//trim(file(1)), 1)
!
! try to understand if it is k-resolved
!
kresolveddos=.false.
READ(10,*, IOSTAT=ios) cdum, str1
IF (ios/=0) CALL errore("sumpdos", "reading first line of "//trim(file(1)), 1)
IF ( trim(str1) == 'ik' ) THEN
kresolveddos=.true.
WRITE(0,"('Summing K-resolved DOS files')")
END IF
REWIND(10)
!
! try to understand if we have 1 or 2 spin
!
IF (kresolveddos) THEN
READ(10,*, IOSTAT=ios) cdum, cdum, cdum, cdum, str1, str2
ELSE
READ(10,*, IOSTAT=ios) cdum, cdum, cdum, str1, str2
END IF
IF (ios/=0) CALL errore("sumpdos", "reading first line of "//trim(file(1)), 1)
!
IF ( trim(str1) == 'ldos(E)' ) THEN
nspin = 1
ELSEIF ( trim(str1) == 'ldosup(E)' .and. trim(str2) == 'ldosdw(E)' ) THEN
nspin = 2
ELSE
CALL errore("sumpdos", "wrong fmf in the first line of "//trim(file(1)), 1)
ENDIF
!
! determine the dimension fo the energy mesh
! no further control will be done on the consistency of the energy
! grid of each file
!
ie = 0
DO WHILE ( .true. )
IF (kresolveddos) THEN
READ( 10, *, IOSTAT=ios ) nk
IF ( nk>1 ) exit
IF ( ios /= 0 ) exit
ELSE
READ( 10, *, IOSTAT=ios )
IF ( ios /= 0 ) exit
END IF
ie = ie + 1
ENDDO
ngrid = ie
!
! K-resolved DOS: continue reading to get the number of K-points
IF (kresolveddos) THEN
DO WHILE ( .true. )
READ( 10, *, IOSTAT=ios ) nk
IF ( ios /= 0 ) exit
END DO
ELSE
nk=1
END IF
!
CLOSE(10)
!
! allocations
!
ALLOCATE( pdos( ngrid, nspin, nfile), STAT=ierr )
IF (ierr/=0) CALL errore("sumpdos", "allocating pdos", ierr)
ALLOCATE( mysum( ngrid, nspin), STAT=ierr )
IF (ierr/=0) CALL errore("sumpdos", "allocating mysum", ierr)
ALLOCATE( egrid( ngrid) )
IF (ierr/=0) CALL errore("sumpdos", "allocating egrid", ierr)
!
! get data
!
WRITE(0,"('Reading the following ',i5,' files: ')") nfile
!
DO ifile = 1, nfile
iun=10+ifile
!
INQUIRE( FILE=trim(file(ifile)), EXIST=exist )
IF (.not. exist) &
CALL errore('sumpdos','file '//trim(file(ifile))//' does not exist',ifile)
!
WRITE(0,"(2x,'Reading file: ',a)") trim(file(ifile))
OPEN(iun, FILE=trim(file(ifile)), IOSTAT=ios)
IF (ios/=0) CALL errore("sumpdos", "error opening "//trim(file(ifile)), ios )
!
END DO
!
DO ik = 1, nk
DO ifile = 1, nfile
iun=10+ifile
!
READ(iun,*, IOSTAT=ios)
IF (ios/=0) &
CALL errore("sumpdos", "reading first line in "//trim(file(ifile)), ios )
!
! egrid is overwritten every time
!
DO ie = 1, ngrid
IF (kresolveddos) THEN
READ(iun, *, IOSTAT=ios ) cdum, egrid(ie), pdos(ie, 1:nspin, ifile)
ELSE
READ(iun, *, IOSTAT=ios ) egrid(ie), pdos(ie, 1:nspin, ifile)
END IF
IF (ios/=0) &
CALL errore("sumpdos", "reading first line in "//trim(file(ifile)), ie )
ENDDO
ENDDO
!
! perform the sum and write
!
IF (kresolveddos) THEN
IF ( ik == 1 ) THEN
IF ( nspin == 1 ) THEN
WRITE(6,"('# ik E (eV) pdos(E) ')")
ELSEIF ( nspin == 2) THEN
WRITE(6,"('# ik E (eV) pdosup(E) pdosdw(E) ')")
ELSE
CALL errore("sumpdos", "really sure NSPIN /= 1 or 2 ???", 3 )
ENDIF
ELSE
WRITE(6,*)
END IF
ELSE
IF ( nspin == 1 ) THEN
WRITE(6,"('# E (eV) pdos(E) ')")
ELSEIF ( nspin == 2) THEN
WRITE(6,"('# E (eV) pdosup(E) pdosdw(E) ')")
ELSE
CALL errore("sumpdos", "really sure NSPIN /= 1 or 2 ???", 3 )
ENDIF
END IF
mysum = 0.0d0
DO ie=1,ngrid
DO isp=1,nspin
mysum(ie,isp) = sum( pdos(ie,isp,:) )
ENDDO
IF (kresolveddos) THEN
WRITE(6,"(i5,' ',f7.3,2e11.3)") ik, egrid(ie) - efermi, mysum(ie,1:nspin)
ELSE
WRITE(6,"(f7.3,2e11.3)") egrid(ie) - efermi, mysum(ie,1:nspin)
END IF
ENDDO
END DO
!
! clean
!
DEALLOCATE( file, STAT=ierr )
IF (ierr/=0) CALL errore("sumpdos", "deallocating file", ierr)
DEALLOCATE( pdos, STAT=ierr )
IF (ierr/=0) CALL errore("sumpdos", "deallocating pdos", ierr)
DEALLOCATE( mysum, STAT=ierr )
IF (ierr/=0) CALL errore("sumpdos", "deallocating mysum", ierr)
DEALLOCATE( egrid, STAT=ierr )
IF (ierr/=0) CALL errore("sumpdos", "deallocating egrid", ierr)
DO ifile = 1, nfile
iun=10+ifile
CLOSE (iun)
END DO
CONTAINS
!*************************************************
SUBROUTINE errore(routine, msg, ierr)
!*************************************************
IMPLICIT NONE
CHARACTER(*), INTENT(in) :: routine, msg
INTEGER, INTENT(in) :: ierr
!
WRITE( UNIT = 0, FMT = '(/,1X,78("*"))')
WRITE( UNIT = 0, &
FMT = '(5X,"from ",A," : error #",I10)' ) routine, ierr
WRITE( UNIT = 0, FMT = '(5X,A)' ) msg
WRITE( UNIT = 0, FMT = '(1X,78("*"),/)' )
!
STOP
RETURN
END SUBROUTINE errore
END PROGRAM sumpdos
More information about the developers
mailing list