1!=============================================================================
2!
3! Utilities:
4!
5! (1) wfn_rho_vxc_info    Originally By DAS      Last Modified 12/12/2011 (DAS)
6!
7!     Prints the contents of the header of a WFN, RHO, or VXC file in
8!     a human-readable format.
9!
10!==============================================================================
11
12#include "f_defs.h"
13
14program wfn_rho_vxc_info
15
16  use global_m
17  use wfn_rho_vxc_io_m
18  use wfn_io_hdf5_m
19#ifdef HDF5
20  use hdf5
21#endif
22  implicit none
23
24  type(mf_header_t) :: mf
25  character*256 :: infile, usage, unit_string, length_label, energy_label
26  integer :: nargs, ii, iat, ik, isym, is, ierr
27  logical :: is_atomic, use_hdf5
28  real(DP) :: length_factor, energy_factor
29
30  usage = 'Usage: wfn_rho_vxc_info.x wfn[.h5] [optional units: eVA/au]'
31
32! Get file names from command-line arguments
33
34  nargs = iargc()
35
36  if (nargs < 1 .or. nargs > 2) then
37    call die(usage)
38  endif
39
40  call getarg(1, infile)
41
42  if(nargs == 2) then
43    call getarg(2, unit_string)
44    if(trim(unit_string) == 'eVA') then
45      is_atomic = .false.
46    else if(trim(unit_string) == 'au') then
47      is_atomic = .true.
48    else
49      call die("Unknown unit specifier " // trim(unit_string) // &
50        ". Choose 'eVA' for eV/Angstrom, or 'au' for atomic units.")
51    endif
52  else
53    is_atomic = .true.
54  endif
55
56  if(is_atomic) then
57    length_label = "bohr"
58    length_factor = 1
59    energy_label = "Ry"
60    energy_factor = 1
61  else
62    length_label = "Ang"
63    length_factor = BOHR
64    energy_label = "eV"
65    energy_factor = RYD
66  endif
67
68  use_hdf5 = .false.
69#ifdef HDF5
70  use_hdf5 = index(infile, '.h5') == len(TRUNC(infile)) - 2
71#endif
72  if (.not.use_hdf5) then
73    call open_file(unit=11, file=TRUNC(infile), form='unformatted', status='old')
74    call read_mf_header(11, mf, warn=.false., dont_warn_kgrid=.true.)
75    call close_file(11)
76  else
77#ifdef HDF5
78    call h5open_f(ierr)
79    call read_hdf5_mf_header(TRUNC(infile), mf)
80    call h5close_f(ierr)
81#endif
82  endif
83
84  write(6,'(a)') '====== GENERAL ====='
85  write(6,'(a,a)')  'Type: ', mf%sheader
86  if(mf%iflavor == 1) then
87    write(6,'(a)') 'Flavor: real'
88  else
89    write(6,'(a)') 'Flavor: complex'
90  endif
91  write(6,'(a,a)')  'Date created: ', mf%sdate
92  write(6,'(a,a)')  'Time created: ', mf%stime
93  write(6,'(a,i1)') 'Number of spins: ', mf%kp%nspin
94
95  write(6,'(a)') '====== G-VECTORS ====='
96  write(6,'(a,i8)') 'Number of G-vectors: ', mf%gvec%ng
97  write(6,'(a,f12.6,a)') 'Charge density cutoff: ', mf%gvec%ecutrho, ' Ry'
98  write(6,'(a,3i8)') 'FFT grid: ', mf%gvec%FFTgrid(1:3)
99  if(mf%sheader == 'WFN') then
100    write(6,'(a,i8)') 'Max number of wfn G-vectors: ', mf%kp%ngkmax
101    write(6,'(a,f12.6,a)') 'Wavefunction cutoff: ', mf%kp%ecutwfc, ' Ry'
102  endif
103
104  write(6,'(a)') '====== ATOMS ====='
105  write(6,'(a,i6)') 'Number of atoms: ', mf%crys%nat
106  write(6,'(2a10,a30)') 'Index', 'Species', 'Coordinates (' // trim(length_label) // ')'
107  do iat = 1, mf%crys%nat
108    write(6,'(2i10,3f12.6)') iat, mf%crys%atyp(iat), mf%crys%apos(1:3, iat) * length_factor
109  enddo
110
111  write(6,'(a)') '====== LATTICE ====='
112  write(6,'(a,f12.6)') 'Cell volume (real space, ' // trim(length_label) // '^3): ', &
113    mf%crys%celvol * length_factor**3
114  write(6,'(a,f12.6)') 'Lattice constant (real space, ' // trim(length_label) // '): ', &
115    mf%crys%alat  * length_factor
116  write(6,'(a)') 'Lattice vectors (real space):'
117  do ii = 1, 3
118    write(6,'(10x,3f12.6)') mf%crys%avec(1:3, ii)
119  enddo
120  write(6,'(a)') 'Metric (real space, ' // trim(length_label) // '^2):'
121  do ii = 1, 3
122    write(6,'(10x,3f12.6)') mf%crys%adot(1:3, ii) * length_factor**2
123  enddo
124
125  write(6,'(a,f12.6)') 'Cell volume (reciprocal space, ' // trim(length_label) // '^-3): ', &
126    mf%crys%recvol / length_factor**3
127  write(6,'(a,f12.6)') 'Lattice constant (reciprocal space, ' // trim(length_label) // '^-1): ', &
128    mf%crys%blat / length_factor
129  write(6,'(a)') 'Lattice vectors (reciprocal space):'
130  do ii = 1, 3
131    write(6,'(10x,3f12.6)') mf%crys%bvec(1:3, ii)
132  enddo
133  write(6,'(a)') 'Metric (reciprocal space, ' // trim(length_label) // '^-2):'
134  do ii = 1, 3
135    write(6,'(10x,3f12.6)') mf%crys%bdot(1:3, ii) / length_factor**2
136  enddo
137
138  write(6,'(a)') '====== SYMMETRIES ====='
139  write(6,'(a,i2)') 'Number of symmetries: ', mf%syms%ntran
140  if(mf%syms%cell_symmetry == 0) then
141    write(6,'(a)') 'Symmetry type: cubic'
142  else
143    write(6,'(a)') 'Symmetry type: hexagonal'
144  endif
145  write(6,'(a7,a31,12x,a33)') 'Index', 'Rotation matrix', 'Fractional translations'
146  do isym = 1, mf%syms%ntran
147    write(6,'(i5,1x,a,2x,3(3i4,2x),3f12.6)') isym, ':', &
148      mf%syms%mtrx(1:3, 1:3, isym), mf%syms%tnp(1:3, isym)
149  enddo
150
151  if(mf%sheader == 'WFN') then
152    write(6,'(a)') '====== K-POINTS ====='
153    write(6,'(a,i8)') 'Number of k-points: ', mf%kp%nrk
154    write(6,'(a,i8)') 'Number of bands: ', mf%kp%mnband
155    write(6,'(a,3i4)') 'k-grid: ', mf%kp%kgrid(1:3)
156    write(6,'(a,3f12.6)') 'k-shifts: ', mf%kp%shift(1:3)
157    write(6,'(a)') '[ifmin = lowest occupied band, ifmax = highest occupied band, for each spin]'
158    write(6,'(a8,a30,6x,a12,a22)',advance='no') 'Index', 'Coordinates (crystal)', 'Weight', 'Number of G-vectors'
159    if(mf%kp%nspin == 1) then
160      write(6,'(2a10)') 'ifmin', 'ifmax'
161    else
162      write(6,'(4a10)') 'ifmin1', 'ifmax1', 'ifmin2', 'ifmax2'
163    endif
164    do ik = 1, mf%kp%nrk
165      write(6,'(i8,4f12.6,i22,4i10)') ik, mf%kp%rk(1:3, ik), mf%kp%w(ik), mf%kp%ngk(ik), &
166        (mf%kp%ifmin(ik, is), mf%kp%ifmax(ik, is), is = 1, mf%kp%nspin)
167    enddo
168
169    write(6,'(a)') '====== ENERGIES (' // trim(energy_label) // ')/OCCUPATIONS ====='
170    do is = 1, mf%kp%nspin
171      if(mf%kp%nspin > 1) write(6,'(a,i2)') 'Spin ', is
172      do ik = 1, mf%kp%nrk
173        write(6,'(a,i6)') 'k-point ', ik
174        write(6,'(9999999f14.6)') mf%kp%el(:, ik, is) * energy_factor
175        write(6,'(9999999f14.6)') mf%kp%occ(:, ik, is)
176      enddo
177    enddo
178  endif
179
180  call dealloc_header_type(mf%sheader, mf%crys, mf%kp)
181
182end program wfn_rho_vxc_info
183