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