1!
2! Copyright (C) 2005 Andrea Ferretti
3!      Modified 2016 Guido Fratesi
4! This file is distributed under the terms of the
5! GNU General Public License. See the file `License'
6! in the root directory of the present distribution,
7! or http://www.gnu.org/copyleft/gpl.txt .
8!
9PROGRAM sumpdos
10  !
11  IMPLICIT NONE
12  !
13  ! AUTHOR: Andrea Ferretti
14  ! edited: Guido Fratesi (to sum K-resolved DOS files)
15  !
16  ! this program reads and sum pdos from different
17  ! files (which are related to different atoms)
18  !
19  ! file names are read from stdin
20  ! USAGE: sumpdos <file1> ... <fileN>
21  !
22  INTEGER             :: ngrid              ! dimension of the energy grid
23  INTEGER             :: nfile              ! number of files to sum
24  INTEGER             :: nspin              ! number of spin_component
25
26
27  CHARACTER(256), ALLOCATABLE    :: file(:) ! names of the files to sum
28  CHARACTER(256)      :: filein
29  CHARACTER(10)       :: cdum, str1, str2, str3
30
31  LOGICAL             :: exist, kresolveddos
32  REAL                :: efermi = 0.0d0       ! translate the input grid
33  REAL, ALLOCATABLE   :: pdos(:,:,:)
34  REAL, ALLOCATABLE   :: egrid(:)
35  REAL, ALLOCATABLE   :: mysum(:,:)
36
37  INTEGER :: ios, ierr, iarg, ie, isp, ifile, i
38  INTEGER :: ik, nk, iun
39
40  !**************************************************************
41  ! User should supply input values here
42  !
43  efermi = 0.0d0
44
45  !**************************************************************
46
47  !
48  ! get the number of arguments (i.e. the number of files)
49  !
50  nfile = command_argument_count ()
51  IF ( nfile == 0 ) THEN
52     WRITE(0,"( 'No file to sum' )")
53     STOP
54  ENDIF
55
56  CALL get_command_argument ( 1, str1 )
57  !
58  SELECT CASE ( trim(str1) )
59  CASE ( "-h" )
60     !
61     ! write the manual
62     !
63     WRITE(0,"(/,'USAGE: sumpdos [-h] [-f <filein>] [<file1> ... <fileN>]', /, &
64          &'  Sum the pdos from the file specified in input and write the sum ', /, &
65          &'  to stdout', /, &
66          &'     -h           : write this manual',/, &
67          &'     -f <filein>  : takes the list of pdos files from <filein> ', /, &
68          &'                    (one per line) instead of command line',/, &
69          &'     <fileM>      : the M-th pdos file', &
70          & / )")
71     STOP
72     !
73  CASE ( "-f" )
74     !
75     ! read file names from file
76     !
77     CALL get_command_argument ( 2, filein )
78     IF ( len_trim(filein) == 0 ) CALL errore('sumpdos','provide filein name',2)
79
80     INQUIRE( FILE=trim(filein), EXIST=exist )
81     IF (.not. exist) CALL errore('sumpdos','file '//trim(filein)//' does not exist',3)
82     OPEN( 10, FILE=trim(filein), IOSTAT=ios )
83     IF (ios/=0) CALL errore('sumpdos','opening '//trim(filein),abs(ios))
84
85     !
86     ! get the number of non-empty lines in the file
87     ! (which is assumed to be the number of files to sum)
88     !
89     ios = 0
90     nfile = 0
91     !
92     DO WHILE ( ios == 0 )
93        nfile = nfile + 1
94        READ(10, *, IOSTAT=ios ) cdum
95        IF ( ios ==0 .and. len_trim(cdum)==0 ) nfile = nfile -1
96     ENDDO
97     nfile = nfile -1
98
99     !
100     IF (nfile ==0 ) CALL errore('sumpdos','no file to sum in '//trim(filein),4)
101     !
102     ALLOCATE( file(nfile), STAT=ierr )
103     IF (ierr/=0) CALL errore('sumpdos','allocating FILE',abs(ierr))
104     !
105     REWIND(10)
106
107     DO i = 1, nfile
108        file(i) = ' '
109        DO WHILE( len_trim(file(i)) == 0 )
110           READ(10,*, IOSTAT=ios) file(i)
111           IF (ios /=0 ) CALL errore('sumpdos','reading from '//trim(filein),i)
112        ENDDO
113     ENDDO
114
115  CASE DEFAULT
116
117     !
118     ! get the names of the files
119     !
120     ALLOCATE( file(nfile), STAT=ierr )
121     IF (ierr/=0) CALL errore('sumpdos','allocating FILE',abs(ierr))
122     DO iarg = 1, nfile
123        CALL get_command_argument ( iarg, file(iarg) )
124     ENDDO
125
126  END SELECT
127
128  !
129  ! open the first file and get data about spin
130  ! and grid dimensions
131  !
132  INQUIRE( FILE=trim(file(1)), EXIST=exist )
133  IF (.not. exist) CALL errore('sumpdos','file '//trim(file(1))//' does not exist',3)
134  !
135  WRITE(0,"('Reading dimensions from file: ',a)") trim(file(1))
136  !
137  OPEN(10, FILE=trim(file(1)), IOSTAT=ios)
138  IF (ios/=0) CALL errore("sumpdos", "error opening "//trim(file(1)), 1)
139  !
140  ! try to understand if it is k-resolved
141  !
142  kresolveddos=.false.
143  READ(10,*, IOSTAT=ios) cdum, str1
144  IF (ios/=0) CALL errore("sumpdos", "reading first line of "//trim(file(1)), 1)
145  IF ( trim(str1) == 'ik' ) THEN
146     kresolveddos=.true.
147     WRITE(0,"('Summing K-resolved DOS files')")
148  END IF
149  REWIND(10)
150  !
151  ! try to understand if we have 1 or 2 spin
152  !
153  IF (kresolveddos) THEN
154     READ(10,*, IOSTAT=ios) cdum, cdum, cdum, cdum, str1, str2
155  ELSE
156     READ(10,*, IOSTAT=ios) cdum, cdum, cdum, str1, str2
157  END IF
158  IF (ios/=0) CALL errore("sumpdos", "reading first line of "//trim(file(1)), 1)
159  !
160  IF ( trim(str1) == 'ldos(E)' ) THEN
161     nspin = 1
162  ELSEIF ( trim(str1) == 'ldosup(E)' .and.  trim(str2) == 'ldosdw(E)' ) THEN
163     nspin = 2
164  ELSE
165     CALL errore("sumpdos", "wrong fmf in the first line of "//trim(file(1)), 1)
166  ENDIF
167  !
168  ! determine the dimension fo the energy mesh
169  ! no further control will be done on the consistency of the energy
170  ! grid of each file
171  !
172  ie = 0
173  DO WHILE ( .true.  )
174     IF (kresolveddos) THEN
175        READ( 10, *, IOSTAT=ios ) nk
176        IF ( nk>1 ) exit
177        IF ( ios /= 0 ) exit
178     ELSE
179        READ( 10, *, IOSTAT=ios )
180        IF ( ios /= 0 ) exit
181     END IF
182     ie = ie + 1
183  ENDDO
184  ngrid = ie
185  !
186  ! K-resolved DOS: continue reading to get the number of K-points
187  IF (kresolveddos) THEN
188     DO WHILE ( .true. )
189        READ( 10, *, IOSTAT=ios ) nk
190        IF ( ios /= 0 ) exit
191     END DO
192  ELSE
193     nk=1
194  END IF
195  !
196  CLOSE(10)
197
198  !
199  ! allocations
200  !
201  ALLOCATE( pdos( ngrid, nspin, nfile), STAT=ierr )
202  IF (ierr/=0) CALL errore("sumpdos", "allocating pdos", ierr)
203  ALLOCATE( mysum( ngrid, nspin), STAT=ierr )
204  IF (ierr/=0) CALL errore("sumpdos", "allocating mysum", ierr)
205  ALLOCATE( egrid( ngrid) )
206  IF (ierr/=0) CALL errore("sumpdos", "allocating egrid", ierr)
207
208
209  !
210  ! get data
211  !
212  WRITE(0,"('Reading the following ',i5,' files: ')") nfile
213  !
214  DO ifile = 1, nfile
215     iun=10+ifile
216     !
217     INQUIRE( FILE=trim(file(ifile)), EXIST=exist )
218     IF (.not. exist) &
219          CALL errore('sumpdos','file '//trim(file(ifile))//' does not exist',ifile)
220     !
221     WRITE(0,"(2x,'Reading file: ',a)") trim(file(ifile))
222     OPEN(iun, FILE=trim(file(ifile)), IOSTAT=ios)
223     IF (ios/=0) CALL errore("sumpdos", "error opening "//trim(file(ifile)), ios )
224     !
225  END DO
226  !
227  DO ik = 1, nk
228     DO ifile = 1, nfile
229        iun=10+ifile
230        !
231        READ(iun,*, IOSTAT=ios)
232        IF (ios/=0) &
233             CALL errore("sumpdos", "reading first line in "//trim(file(ifile)), ios )
234        !
235        ! egrid is overwritten every time
236        !
237        DO ie = 1, ngrid
238           IF (kresolveddos) THEN
239              READ(iun, *, IOSTAT=ios ) cdum, egrid(ie), pdos(ie, 1:nspin, ifile)
240           ELSE
241              READ(iun, *, IOSTAT=ios ) egrid(ie), pdos(ie, 1:nspin, ifile)
242           END IF
243           IF (ios/=0) &
244                CALL errore("sumpdos", "reading first line in "//trim(file(ifile)), ie )
245        ENDDO
246     ENDDO
247
248     !
249     ! perform the sum and write
250     !
251     IF (kresolveddos) THEN
252        IF ( ik == 1 ) THEN
253           IF ( nspin == 1 ) THEN
254              WRITE(6,"('# ik   E (eV)  pdos(E) ')")
255           ELSEIF ( nspin == 2) THEN
256              WRITE(6,"('# ik   E (eV)  pdosup(E)  pdosdw(E) ')")
257           ELSE
258              CALL errore("sumpdos", "really sure NSPIN /= 1 or 2 ???", 3 )
259           ENDIF
260        ELSE
261           WRITE(6,*)
262        END IF
263     ELSE
264        IF ( nspin == 1 ) THEN
265           WRITE(6,"('# E (eV) pdos(E) ')")
266        ELSEIF ( nspin == 2) THEN
267           WRITE(6,"('# E (eV) pdosup(E)  pdosdw(E) ')")
268        ELSE
269           CALL errore("sumpdos", "really sure NSPIN /= 1 or 2 ???", 3 )
270        ENDIF
271     END IF
272
273     mysum = 0.0d0
274     DO ie=1,ngrid
275        DO isp=1,nspin
276           mysum(ie,isp) = sum( pdos(ie,isp,:) )
277        ENDDO
278        IF (kresolveddos) THEN
279           WRITE(6,"(i5,' ',f7.3,2e11.3)") ik, egrid(ie) - efermi, mysum(ie,1:nspin)
280        ELSE
281           WRITE(6,"(f7.3,2e11.3)") egrid(ie) - efermi, mysum(ie,1:nspin)
282        END IF
283     ENDDO
284
285  END DO
286
287  !
288  ! clean
289  !
290  DEALLOCATE( file, STAT=ierr )
291  IF (ierr/=0) CALL errore("sumpdos", "deallocating file", ierr)
292  DEALLOCATE( pdos, STAT=ierr )
293  IF (ierr/=0) CALL errore("sumpdos", "deallocating pdos", ierr)
294  DEALLOCATE( mysum, STAT=ierr )
295  IF (ierr/=0) CALL errore("sumpdos", "deallocating mysum", ierr)
296  DEALLOCATE( egrid, STAT=ierr )
297  IF (ierr/=0) CALL errore("sumpdos", "deallocating egrid", ierr)
298  DO ifile = 1, nfile
299     iun=10+ifile
300     CLOSE (iun)
301  END DO
302
303
304CONTAINS
305
306  !*************************************************
307  SUBROUTINE errore(routine, msg, ierr)
308    !*************************************************
309    IMPLICIT NONE
310    CHARACTER(*),    INTENT(in) :: routine, msg
311    INTEGER,         INTENT(in) :: ierr
312
313    !
314    WRITE( UNIT = 0, FMT = '(/,1X,78("*"))')
315    WRITE( UNIT = 0, &
316         FMT = '(5X,"from ",A," : error #",I10)' ) routine, ierr
317    WRITE( UNIT = 0, FMT = '(5X,A)' ) msg
318    WRITE( UNIT = 0, FMT = '(1X,78("*"),/)' )
319    !
320    STOP
321    RETURN
322  END SUBROUTINE errore
323
324END PROGRAM sumpdos
325
326