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