1! 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt. 6! See Docs/Contributors.txt for a list of contributors. 7! 8 9program readwfx 10 11 use m_getopts 12 use f2kcli 13 14! This program READWF reads a the coefficients of wavefunctions 15! on the expansion of the atomic orbitals basis set, as written 16! by SIESTA (unformatted) and writes them in ascii, user-friendly 17! form. 18 19! Written by P. Ordejon, June 2003 20! Updated to WFSX format and given more features by A. Garcia. 21 22! threshold to plot the coefficients of the wavefunctions. If the 23! norm of the weight of a wavefunction on a given orbital is smaller 24! than this threshold then it is not printed. This is useful for 25! large systems, with a very large number of basis orbitals. 26 27implicit none 28 29integer, parameter :: dp = selected_real_kind(10,100) 30integer, parameter :: sp = selected_real_kind(5,10) 31 32 character(len=200) :: opt_arg 33 character(len=10) :: opt_name 34 integer :: nargs, iostat, n_opts, nlabels 35 36 integer :: io, wfs_u, nk, nspin_blocks, ik, ispin, is, idummy, & 37 number_of_wfns, iw, indwf, j, nuotot, jj, nspin_flag, is0 38 39 character(len=256) :: wfsx_file, outfile = "" 40 41 integer, allocatable, dimension(:) :: iaorb,iphorb,cnfigfio 42 character(len=20), allocatable, dimension(:) :: symfio,labelfis 43 44 real(sp), allocatable, dimension(:,:) :: psi 45 logical :: gamma, non_coll 46 real(dp) :: k(3), eigval, norm, wk 47 48 real(dp) :: threshold = 0.0_dp 49 50 real(dp) :: emin = -huge(1.0_dp) 51 real(dp) :: emax = huge(1.0_dp) 52 logical :: emin_given = .false. 53 logical :: emax_given = .false. 54 integer :: min_band = -huge(1) 55 integer :: max_band = huge(1) 56 logical :: min_band_set = .false. 57 logical :: max_band_set = .false. 58 logical :: energies_only = .false. 59 60 real(dp) :: min_eigval, max_eigval 61 real(dp) :: min_eigval_in_band_set, max_eigval_in_band_set 62 integer :: nwfmin, nwfmax, wfs_spin_flag 63 64! Process options 65! 66 n_opts = 0 67 do 68 call getopts('hlt:e:m:E:M:b:B:o:',opt_name,opt_arg,n_opts,iostat) 69 if (iostat /= 0) exit 70 select case(opt_name) 71 case ('o') 72 read(opt_arg,*) outfile 73 case ('m', 'e') 74 emin_given = .true. 75 read(opt_arg,*) emin 76 case ('M', 'E') 77 emax_given = .true. 78 read(opt_arg,*) emax 79 case ('b') 80 read(opt_arg,*) min_band 81 case ('B') 82 read(opt_arg,*) max_band 83 case ('t') 84 read(opt_arg,*) threshold 85 case ('l') 86 energies_only = .true. 87 case ('h') 88 call manual() 89 STOP 90 case ('?',':') 91 write(0,*) "Invalid option: ", opt_arg(1:1) 92 write(0,*) "Use -h option for manual" 93 write(0,*) "" 94 call manual() 95 STOP 96 end select 97 enddo 98 99 nargs = command_argument_count() 100 nlabels = nargs - n_opts + 1 101 if (nlabels /= 1) then 102 write(0,*) "Use -h option for manual" 103 write(0,*) "" 104 call manual() 105 STOP 106 endif 107 108 call get_command_argument(n_opts,value=wfsx_file,status=iostat) 109 if ( iostat /= 0 ) then 110 stop "Cannot get WFSX file" 111 end if 112 113 threshold = threshold**2 ! We use the square later on 114 115 wfs_u = 10 116 io = 6 117 118 open(wfs_u, file=wfsx_file, form='unformatted', status='old' ) 119 120 rewind(wfs_u) 121 read(wfs_u) nk, gamma 122 123 read(wfs_u) wfs_spin_flag ! 1, 2, or 4 124 non_coll = (wfs_spin_flag >= 4) 125 read(wfs_u) nuotot 126 read(wfs_u) !! Symbols, etc 127 128 129 if (non_coll) then 130 nspin_blocks = 1 131 else 132 nspin_blocks = wfs_spin_flag 133 endif 134 135 !------------------------------------- 136 137 nwfmax = -huge(1) 138 nwfmin = huge(1) 139 min_eigval = huge(1.0_dp) 140 max_eigval = -huge(1.0_dp) 141 min_eigval_in_band_set = huge(1.0_dp) 142 max_eigval_in_band_set = -huge(1.0_dp) 143 144 do ik=1,nk 145 do is=1,nspin_blocks 146 147 read(wfs_u) idummy, k(1:3), wk 148 if (idummy /= ik) stop "ik index mismatch in WFS file" 149 read(wfs_u) is0 150 read(wfs_u) number_of_wfns 151 nwfmax = max(nwfmax,number_of_wfns) 152 nwfmin = min(nwfmin,number_of_wfns) 153 154 do iw=1,number_of_wfns 155 read(wfs_u) indwf 156 read(wfs_u) eigval 157 min_eigval = min(min_eigval,eigval) 158 max_eigval = max(max_eigval,eigval) 159 ! 160 ! 161 if ((iw>=min_band).and.(iw<=max_band)) then 162 min_eigval_in_band_set = min(min_eigval_in_band_set,eigval) 163 max_eigval_in_band_set = max(max_eigval_in_band_set,eigval) 164 endif 165 read(wfs_u) 166 enddo 167 enddo 168 enddo 169 170 print "(a,2i5)", "Minimum/Maximum number of wfs per k-point: ", nwfmin, nwfmax 171 print "(a,2f12.4)", "Min_eigval, max_eigval on WFS file: ", & 172 min_eigval, max_eigval 173 174 print "(a,2f12.4)", "Min_eigval, max_eigval in band set : ", & 175 min_eigval_in_band_set, max_eigval_in_band_set 176 177 if (energies_only) STOP ! energies only 178 179 if (outfile /= "") then 180 open(io, file=outfile, form='formatted', status='unknown' ) 181 endif 182 183 rewind (wfs_u) 184 185 read(wfs_u) nk, gamma 186 read(wfs_u) nspin_flag 187 read(wfs_u) nuotot 188 189 if (nspin_flag == 4) then 190 non_coll = .true. 191 nspin_blocks = 1 192 else 193 non_coll = .false. 194 nspin_blocks = nspin_flag 195 endif 196 197 if (gamma) then 198 if (non_coll) then 199 allocate(psi(4,nuotot)) 200 else 201 allocate(psi(1,nuotot)) 202 endif 203 else 204 if (non_coll) then 205 allocate(psi(4,nuotot)) 206 else 207 allocate(psi(2,nuotot)) 208 endif 209 endif 210 211 allocate(iaorb(nuotot), labelfis(nuotot), & 212 iphorb(nuotot), cnfigfio(nuotot), & 213 symfio(nuotot)) 214 215 read(wfs_u) (iaorb(j),labelfis(j), & 216 iphorb(j), cnfigfio(j), symfio(j), j=1,nuotot) 217 218 write(io,*) 219 write(io,'(a22,2x,i6)') 'Nr of k-points = ',nk 220 if (non_coll) then 221 write(io,'(a,2x,i6)') 'Nr of Spins blocks ' // & 222 '(non-collinear) = ', nspin_blocks 223 else 224 write(io,'(a,2x,i6)') 'Nr of Spins blocks = ',nspin_blocks 225 endif 226 write(io,'(a,2x,i6)') 'Nr of basis orbs = ',nuotot 227 write(io,*) 228 229 kpoints: do ik = 1,nk 230 spin_blocks: do ispin = 1,nspin_blocks 231 232 read(wfs_u) idummy, k(1),k(2),k(3) 233 if (idummy .ne. ik) stop 'error in index of k-point' 234 ! if ik not in k-range CYCLE kpoints 235 read(wfs_u) idummy 236 if (.not. non_coll) then 237 if (idummy .ne. ispin) stop 'error in index of spin' 238 ! if ispin not in requested spin-channel CYCLE spin_blocks 239 endif 240 read(wfs_u) number_of_wfns 241 242 write(io,*) 243 write(io,fmt='(a,2x,i6,2x,3f10.6)',advance="no") 'k-point = ',ik, k(1),k(2),k(3) 244 if (.not. non_coll) then 245 write(io,fmt='(2x,a,2x,i1)',advance="no") 'Spin component = ',ispin 246 endif 247 write(io,fmt='(2x,a,1x,i6)') 'Num. wavefunctions = ',number_of_wfns 248 249! Loop over wavefunctions 250 251 do iw = 1,number_of_wfns 252 253 read(wfs_u) indwf 254 read(wfs_u) eigval 255 if ( (eigval < emin) .or. (eigval > emax) .or. & 256 (iw < min_band) .or. (iw > max_band)) then 257 read(wfs_u) ! skip over wfn data 258 CYCLE 259 endif 260 261 write(io,'(a22,2x,i6,2x,a,f10.6)') & 262 'Wavefunction = ', indwf, 'Eigval (eV) = ', eigval 263 264 if (threshold >= 100*100) then 265 read(wfs_u) ! Skip wf data 266 CYCLE ! Early exit to print only energies 267 endif 268 269 write(io,"(60('-'))") 270 if (non_coll) then 271 write(io,'(a72)') & 272 ' Atom Species Orb-global Orb-in-atom '// & 273 ' .Orb-type [Re(psi) Im(psi)]Up ' // & 274 ' [Re(psi) Im(psi)]Down ' 275 else 276 write(io,'(a72)') & 277 ' Atom Species Orb-global Orb-in-atom ' // & 278 ' .Orb-type Re(psi) Im(psi)' 279 endif 280 281 read(wfs_u) (psi(1:,j), j=1,nuotot) 282 do jj = 1,nuotot 283 norm = dot_product(psi(1:,jj),psi(1:,jj)) ! valid for all 284 if (norm .lt. threshold) CYCLE 285 write(io, & 286 '(i6,5x,a10,1x,i10,2x,i3,2x,i1,a20,2(2(f10.6),2x))') & 287 iaorb(jj),labelfis(jj),jj, iphorb(jj), cnfigfio(jj), & 288 symfio(jj), psi(1:,jj) 289 290 enddo 291 292 write(io,"(60('-'))") 293 294 enddo 295 enddo spin_blocks 296 enddo kpoints 297 298 close (wfs_u) 299 300 CONTAINS 301 subroutine manual 302 303 write(6,*) "Usage: readwfx [ options ] WFSX_FILE" 304 write(6,*) "Options:" 305 write(6,*) " -h: print manual " 306 write(6,*) " -l: print summary of energy information " 307 write(6,*) " -o File : set output file (default: use stdout) " 308 write(6,*) " " 309 write(6,*) " Selection of states to be printed: by eigenvalue range or band index: " 310 write(6,*) " " 311 write(6,*) " -m Min_e : set lower bound of eigenvalue range " 312 write(6,*) " -M Max_e : set upper bound of eigenvalue range " 313 write(6,*) " -b Min_band : set minimum band index to be used " 314 write(6,*) " -B Max_band : set maximum band index to be used " 315 write(6,*) " " 316 write(6,*) " -t threshold : set threshold for orbital contributions to be printed " 317 write(6,*) " (use large value to get only energy information) " 318 end subroutine manual 319 320 end program readwfx 321