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