1!
2! SIB: This routine looks the same as input().
3! Except it reads from WFNq and
4! writes to iunit_v='INT_VWFQ' and only valence bands.
5! And k-point information is read into kpq.
6!
7!     SUBROUTINE READS CRYSTAL DATA AND WAVEFUNCTIONS FROM TAPE26
8!     AND PARAMETERS FOR POLARIZABILITY CALCULATION FROM TAPE5
9!     TAPE10 (OUTPUT TAPE) IS INITIALIZED
10
11#include "f_defs.h"
12
13module input_q_m
14
15  use global_m
16  use eqpcor_m
17  use input_utils_m
18  use misc_m
19  use scissors_m
20  use wfn_rho_vxc_io_m
21  use io_utils_m
22#ifdef HDF5
23  use hdf5
24#endif
25  use wfn_io_hdf5_m
26
27  implicit none
28
29  private
30  public :: input_q
31
32contains
33
34subroutine input_q(gvec,kpq,cwfn,vwfn,pol,intwfnvq)
35  type (gspace), intent(in) :: gvec
36  type (kpoints), intent(out) :: kpq
37  type (conduction_wfns), intent(in) :: cwfn
38  type (valence_wfns), intent(in) :: vwfn
39  type (polarizability), intent(inout) :: pol
40  type (int_wavefunction), intent(out) :: intwfnvq
41
42  type (crystal) :: crys
43  type (symmetry) :: syms
44
45  real(DP) :: vcell
46  integer :: dummygvec(1, 1)
47  character :: fncor*32
48
49  character(len=3) :: sheader
50  integer :: iflavor
51  type(gspace) :: gvecq
52
53  PUSH_SUB(input_q)
54
55  sheader = 'WFN'
56  iflavor = 0
57  if(pol%os_hdf5) then
58  else
59    if(peinf%inode == 0) then
60      write(6,'(a)') ' Reading header of WFNq'
61      call open_file(26,file='WFNq',form='unformatted',status='old')
62    endif
63    call read_binary_header_type(26, sheader, iflavor, kpq, gvecq, syms, crys, &
64      dont_warn_kgrid=pol%subsample, warn=.false.)
65    call read_binary_gvectors(26, gvecq%ng, gvecq%ng, dummygvec, dont_read = .true.)
66  endif
67  call check_trunc_kpts(pol%icutv, kpq)
68
69  call scissors_shift(kpq, pol%scis)
70  call get_volume(vcell,crys%bdot)
71  if (abs(crys%celvol-vcell).gt.TOL_Small) then
72    call die('volume mismatch')
73  endif
74
75!-----------------------------------------------------------------
76! If a quasi-particle correction file exists, read the corrected
77! quasiparticle energies from file (in eV)
78
79  if (pol%eqp_corrections) then
80    fncor='eqp_q.dat'
81    call eqpcor(fncor,peinf%inode,peinf%npes,kpq,1+vwfn%ncore_excl,vwfn%nband+vwfn%ncore_excl+pol%ncrit,0,0,kpq%el,kpq%el,&
82                 kpq%el,1,0)
83  endif
84
85  call find_efermi(pol%rfermi, pol%efermi, pol%efermi_input, kpq, kpq%mnband, 1+vwfn%ncore_excl, &
86    "shifted grid", should_search = .false., should_update = .false., write7 = .true.)
87
88  if (peinf%inode==0) then
89    if(any (kpq%ifmax(:,:) < vwfn%nband+vwfn%ncore_excl .or. kpq%ifmax(:,:) > vwfn%nband +vwfn%ncore_excl+ pol%ncrit)) then
90      write(0,'(a,i6,a,i6,a)') 'epsilon.inp says there are ', vwfn%nband, ' fully occupied bands and ', &
91        pol%ncrit, ' partially occupied.'
92      write(0,'(a,2i6)') 'This is inconsistent with highest bands in WFNq file; min, max = ', minval(kpq%ifmax), maxval(kpq%ifmax)
93      call die("band_occupation, number_partial_occup, and WFNq inconsistent.")
94    endif
95
96    if(maxval(kpq%ifmax) - minval(kpq%ifmax) > pol%ncrit) then
97      write(0,'(a,i6,a)') 'epsilon.inp says there are ', pol%ncrit, ' partially occupied bands.'
98      write(0,'(a,i6)') 'This is less than the number partially occupied in WFNq file: ', maxval(kpq%ifmax) - minval(kpq%ifmax)
99      call die("number_partial_occup and WFNq inconsistent.")
100    endif
101  endif
102
103  if(pol%os_hdf5) then
104  else
105    call read_wavefunctions(kpq, gvec, pol, cwfn, vwfn, intwfnvq)
106    if(peinf%inode.eq.0) then
107      call close_file(26)
108    endif
109  endif
110  POP_SUB(input_q)
111
112  return
113end subroutine input_q
114
115
116subroutine read_wavefunctions(kpq, gvec, pol, cwfn, vwfn, intwfnvq)
117  type (kpoints), intent(in) :: kpq
118  type (gspace), intent(in) :: gvec
119  type (polarizability), intent(in) :: pol
120  type (conduction_wfns), intent(in) :: cwfn
121  type (valence_wfns), intent(in) :: vwfn
122  type (int_wavefunction), intent(out) :: intwfnvq
123
124  integer, allocatable :: isort(:)
125  SCALAR, allocatable :: zc(:,:)
126
127  character :: filenamevq*20
128  integer :: i,i2,j,k,ik,iiii,is
129  integer :: iunit_v
130  real(DP) :: qk(3)
131  logical :: dont_read
132
133  type(gspace) :: gvec_kpt
134  type(progress_info) :: prog_info !< a user-friendly progress report
135
136  PUSH_SUB(read_wavefunctions)
137
138  SAFE_ALLOCATE(intwfnvq%ng, (kpq%nrk))
139  SAFE_ALLOCATE(intwfnvq%isort, (kpq%ngkmax,kpq%nrk))
140  SAFE_ALLOCATE(intwfnvq%cg, (kpq%ngkmax,kpq%nrk*peinf%nvownactual,kpq%nspin*kpq%nspinor))
141  SAFE_ALLOCATE(intwfnvq%qk, (3,kpq%nrk))
142
143  call progress_init(prog_info, 'reading wavefunctions (WFNq)', 'state', kpq%nrk*(vwfn%nband+pol%ncrit))
144  do ik=1,kpq%nrk
145    qk(1:3) = kpq%rk(1:3, ik)
146    SAFE_ALLOCATE(gvec_kpt%components, (3, kpq%ngk(ik)))
147
148    call read_binary_gvectors(26, kpq%ngk(ik), kpq%ngk(ik), gvec_kpt%components)
149
150    SAFE_ALLOCATE(isort, (kpq%ngk(ik)))
151    do i = 1, kpq%ngk(ik)
152      call findvector(isort(i), gvec_kpt%components(:, i), gvec)
153      if (isort(i) == 0) then
154        if(peinf%inode == 0) write(0,*) 'ik = ', ik, 'ig = ', i, 'gvec = ', gvec_kpt%components(:, i)
155        call die('input_q: could not find gvec')
156      endif
157    enddo
158    SAFE_DEALLOCATE_P(gvec_kpt%components)
159
160    intwfnvq%ng(ik)=kpq%ngk(ik)
161    intwfnvq%isort(1:kpq%ngk(ik),ik)=isort(1:kpq%ngk(ik))
162    intwfnvq%qk(:,ik)=qk(:)
163!
164! SIB:  loop on max number of bands, and proc 0 reads the wave function
165! from unit 26, checks normalization, and if the band is less than
166! cwfn%nband (# of bands in total) **AND** is a valence band (so
167! its index is <= vwfn%nband), then it is written to iunit_v.
168!
169    SAFE_ALLOCATE(zc, (kpq%ngk(ik), kpq%nspinor*kpq%nspin))
170
171    do i=1,kpq%mnband
172
173      dont_read = (i > cwfn%nband .or. i <= vwfn%ncore_excl)
174      if(.not. dont_read) dont_read = i > vwfn%nband+vwfn%ncore_excl+pol%ncrit
175      call read_binary_data(26, kpq%ngk(ik), kpq%ngk(ik), kpq%nspin*kpq%nspinor, zc, dont_read = dont_read)
176
177      ! FHJ: the following lines were introduced in r6294 and are supposed to
178      ! be a shortcut if we are past the last band of the last k-point. However,
179      ! in light of a previous bug (#223), this feature is commented out for now.
180      !! FHJ: shortcut if this is past the last band of the last k-point
181      !if (dont_read .and. ik==kpq%nrk) exit
182      if (.not.dont_read) then
183        ! DVF: recall that we redefined the number of valence bands to exclude the
184        ! core states. So, we have to subtract ncore_excl right here because ib is
185        ! referenced to the full wavefunction file including all the core states.
186        i2=i-vwfn%ncore_excl
187        call progress_step(prog_info, (ik-1)*(vwfn%nband+pol%ncrit) + i)
188        if (peinf%inode == 0) then
189          do is = 1, kpq%nspin
190            call checknorm('WFNq',i,ik,kpq%ngk(ik),is,kpq%nspinor,zc(:,:))
191          enddo
192        endif
193
194        if (peinf%doiownv(i2)) then
195          iiii=peinf%indexv(i2)+(ik-1)*peinf%nvownactual
196          intwfnvq%cg(1:kpq%ngk(ik),iiii,1:kpq%nspin*kpq%nspinor)=zc(1:kpq%ngk(ik),1:kpq%nspinor*kpq%nspin)
197        endif
198      endif
199    enddo
200    SAFE_DEALLOCATE(isort)
201    SAFE_DEALLOCATE(zc)
202  enddo                     ! end loop over k+q points
203  call progress_free(prog_info)
204
205  POP_SUB(read_wavefunctions)
206  return
207
208end subroutine read_wavefunctions
209
210end module input_q_m
211