1#include "f_defs.h" 2 3!----------------------------------------------------------------------- 4subroutine genwf_kernel(crys,gvec,kg,kgq,syms, & 5 wfnc,wfnv,nspin,ik,ic,iv,indexq,xct,intwfnv,intwfnc) 6! 7! input: crys, gvec, kg, syms types 8! nspin (number of spins) 9! ik (index of k-point in full BZ) 10! ib (index of ck block) 11! output: wfnc (conduction states at point ik, block ib) 12! wfnv (valence states at point ik) 13! 14!----------------------------------------------------------------------- 15 16 use global_m 17 use blas_m 18 use gmap_m 19 use input_utils_m 20 use sort_m 21 use misc_m 22 implicit none 23 24 type (crystal), intent(in) :: crys 25 type (gspace), intent(in) :: gvec 26 type (grid), intent(in) :: kg 27 type (grid), intent(in) :: kgq 28 type (symmetry), intent(in) :: syms 29 type (wavefunction), intent(out) :: wfnc,wfnv 30 integer, intent(in) :: nspin,ik,ic,iv 31 integer, intent(in) :: indexq(kg%nf) 32 type (int_wavefunction), intent(in) :: intwfnc,intwfnv 33 type (xctinfo), intent(in) :: xct 34 35 character :: errmsg*100 36 integer :: cnb,cng,cns,cnsp,vng,vns,vnsp 37 integer :: ig,ii,jj,kk,info 38 integer :: ivown,icown,ikown,invband,incband,ijk 39 real(DP) :: xnorm,qk(3) 40 41 integer, save :: ikold=0 42 integer, save :: ikold2=0 43 integer, save :: vngold=0 44 integer, save :: cngold=0 45 integer, save :: ifirst=0 46 47 integer, allocatable :: isorti(:) 48 integer, allocatable :: ind(:),isort(:) 49 integer, allocatable :: indq(:),isortq(:) 50 integer, allocatable, save :: ind_old(:),isort_old(:) 51 integer, allocatable, save :: indq_old(:),isortq_old(:) 52 integer, allocatable, save :: ind_old2(:),isort_old2(:) 53 integer, allocatable, save :: indq_old2(:),isortq_old2(:) 54 real(DP), allocatable :: ekin(:) 55 SCALAR, allocatable :: ccg(:,:,:),vcg(:,:,:),ph(:),phq(:) 56 SCALAR, allocatable, save :: ph_old(:),phq_old(:) 57 SCALAR, allocatable, save :: ph_old2(:),phq_old2(:) 58 59 PUSH_SUB(genwf_kernel) 60 61!----------------------------------------------------------------------- 62! Deal with the valence wavefunctions 63 64! write(6,*) peinf%inode, 'in genwf',iv,ic,ik 65 66 if (xct%qflag .eq. 0) then 67 else 68 ivown = peinf%ipev(peinf%inode+1,iv,kg%indr(ik)) 69 ikown = peinf%ipek(peinf%inode+1,kg%indr(ik)) 70 endif 71 72 vng = intwfnv%ng(ikown) 73 vns = intwfnv%nspin 74 vnsp = intwfnv%nspinor 75 if (xct%ivpar .eq. 1) then 76 SAFE_ALLOCATE(vcg, (vng,1,vns*vnsp)) 77 else 78 SAFE_ALLOCATE(vcg, (vng,xct%nvb_co,vns*vnsp)) 79 endif 80 81 SAFE_ALLOCATE(indq, (vng)) 82 SAFE_ALLOCATE(phq, (vng)) 83 SAFE_ALLOCATE(isortq, (gvec%ng)) 84 85 if (ik .ne. ikold .and. ik .ne. ikold2 .and. ifirst .ne. 0) then 86 SAFE_DEALLOCATE(indq_old2) 87 SAFE_DEALLOCATE(phq_old2) 88 SAFE_DEALLOCATE(isortq_old2) 89 SAFE_ALLOCATE(indq_old2, (vngold)) 90 SAFE_ALLOCATE(phq_old2, (vngold)) 91 SAFE_ALLOCATE(isortq_old2, (gvec%ng)) 92 93 indq_old2 = indq_old 94 phq_old2 = phq_old 95 isortq_old2 = isortq_old 96 97 SAFE_DEALLOCATE(indq_old) 98 SAFE_DEALLOCATE(phq_old) 99 SAFE_DEALLOCATE(isortq_old) 100 SAFE_ALLOCATE(indq_old, (vng)) 101 SAFE_ALLOCATE(phq_old, (vng)) 102 SAFE_ALLOCATE(isortq_old, (gvec%ng)) 103 endif 104 if (ifirst .eq. 0) then 105 SAFE_ALLOCATE(indq_old, (vng)) 106 SAFE_ALLOCATE(phq_old, (vng)) 107 SAFE_ALLOCATE(isortq_old, (gvec%ng)) 108 SAFE_ALLOCATE(indq_old2, (vng)) 109 SAFE_ALLOCATE(phq_old2, (vng)) 110 SAFE_ALLOCATE(isortq_old2, (gvec%ng)) 111 endif 112 113! Initalize parameters in variable wfnv 114 115 wfnv%ng=vng 116 wfnv%nband=1 117 wfnv%nspin=vns 118 wfnv%nspinor=vnsp 119 if (intwfnv%nspin.ne.nspin) then 120 write(errmsg,'(2a,2i2)') 'spin number mismatch: ', nspin, vns 121 call die(errmsg, only_root_writes = .true.) 122 endif 123 124 if (xct%ivpar .eq. 1) then 125 SAFE_ALLOCATE(wfnv%cg, (wfnv%ng,1,wfnv%nspin*wfnv%nspinor)) 126 vcg(1:vng,1,:) = intwfnv%cg(1:vng,ivown,:) 127 else 128 SAFE_ALLOCATE(wfnv%cg, (wfnv%ng,xct%nvb_co,wfnv%nspin*wfnv%nspinor)) 129 vcg(1:vng,:,:) = intwfnv%cg(1:vng,ivown:ivown+xct%nvb_co-1,:) 130 endif 131 SAFE_ALLOCATE(wfnv%isort, (gvec%ng)) 132 133 isortq(:) = intwfnv%isort(:,ikown) 134 135! Compute inverse index array of Fourier components around rk-kpoint 136 137 if (ik .ne. ikold .and. ik .ne. ikold2) then 138 SAFE_ALLOCATE(isorti, (gvec%ng)) 139 isorti(:)=0 140 do ii=1,wfnv%ng 141 isorti(isortq(ii))=ii 142 enddo 143 endif 144 145! Compute index array of Fourier components around fk-kpoint 146 147 if (ik .ne. ikold .and. ik .ne. ikold2) then 148 if (peinf%verb_debug) then 149 write(6,*) 'ikneikold',peinf%inode,ik,ikold 150 endif 151 SAFE_ALLOCATE(ekin, (gvec%ng)) 152 if (xct%qflag .eq. 0) then 153 qk(1:3)=kgq%f(1:3,indexq(ik)) 154 elseif (xct%qflag .eq. 2) then 155 qk(1:3)=kg%f(1:3,indexq(ik)) 156 else 157 qk(1:3)=kg%f(1:3,ik) 158 endif 159 call kinetic_energies(gvec, crys%bdot, ekin, qvec = qk) 160 call sortrx(gvec%ng,ekin,isortq) 161 SAFE_DEALLOCATE(ekin) 162 isortq_old(:)=isortq(:) 163 elseif (ik .eq. ikold) then 164 isortq(:)=isortq_old(:) 165 else 166 isortq(:)=isortq_old2(:) 167 endif 168 169! Find ind and ph relating wavefunctions in fk to rk-kpoint 170 171 if (ik .ne. ikold .and. ik .ne. ikold2) then 172 indq=0 173 phq=ZERO 174 if (xct%qflag .eq. 0) then 175 else 176 call gmap(gvec,syms,wfnv%ng,kg%itran(ik), & 177 kg%kg0(:,ik),isortq,isorti,indq,phq,.true.) 178 endif 179 SAFE_DEALLOCATE(isorti) 180 indq_old(:)=indq(:) 181 phq_old(:)=phq(:) 182 else if (ik .eq. ikold) then 183 indq(:)=indq_old(:) 184 phq(:)=phq_old(:) 185 else 186 indq(:)=indq_old2(:) 187 phq(:)=phq_old2(:) 188 endif 189 190! Compute and renormalize valence wavefunctions in fk-kpoint 191 192 if (xct%ivpar .eq. 1) invband = 1 193 if (xct%ivpar .eq. 0) invband = xct%nvb_co 194 195 do ijk = 1, invband 196 do kk=1,wfnv%nspin 197 do jj=1,wfnv%nspinor 198 do ii=1,wfnv%ng 199 if (indq(ii) .gt. 0) then 200 wfnv%cg(ii,ijk,kk*jj)=phq(ii)*vcg(indq(ii),ijk,kk*jj) 201 else 202 wfnv%cg(ii,ijk,kk*jj) = ZERO 203 endif 204 enddo !ii 205 enddo ! jj 206 call checknorm('wfnv%cg',ijk,ik,wfnv%ng,kk,wfnv%nspinor,wfnv%cg(1:wfnv%ng,ijk,:)) 207 enddo ! kk 208 enddo ! ijk 209 vcg=wfnv%cg 210 211 212 wfnv%cg=vcg 213 wfnv%isort=isortq 214 215 216 SAFE_DEALLOCATE(vcg) 217 SAFE_DEALLOCATE(indq) 218 SAFE_DEALLOCATE(phq) 219 SAFE_DEALLOCATE(isortq) 220 221!----------------------------------------------------------------------- 222! Deal with the conduction wavefunctions 223 224 icown = peinf%ipec(peinf%inode+1,ic,kg%indr(ik)) 225 ikown = peinf%ipek(peinf%inode+1,kg%indr(ik)) 226 227 cng = intwfnc%ng(ikown) 228 cnb = 1 229 cns = intwfnc%nspin 230 cnsp = intwfnc%nspinor 231 232 if (xct%icpar .eq. 1) then 233 SAFE_ALLOCATE(ccg, (cng,1,cns*cnsp)) 234 else 235 SAFE_ALLOCATE(ccg, (cng,xct%ncb_co,cns*cnsp)) 236 endif 237 238 SAFE_ALLOCATE(ind, (cng)) 239 SAFE_ALLOCATE(ph, (cng)) 240 SAFE_ALLOCATE(isort, (gvec%ng)) 241 242 if (ik .ne. ikold .and. ik .ne. ikold2 .and. ifirst .ne. 0) then 243 SAFE_DEALLOCATE(ind_old2) 244 SAFE_DEALLOCATE(ph_old2) 245 SAFE_DEALLOCATE(isort_old2) 246 SAFE_ALLOCATE(ind_old2, (cngold)) 247 SAFE_ALLOCATE(ph_old2, (cngold)) 248 SAFE_ALLOCATE(isort_old2, (gvec%ng)) 249 250 ind_old2 = ind_old 251 ph_old2 = ph_old 252 isort_old2 = isort_old 253 254 SAFE_DEALLOCATE(ind_old) 255 SAFE_DEALLOCATE(ph_old) 256 SAFE_DEALLOCATE(isort_old) 257 SAFE_ALLOCATE(ind_old, (cng)) 258 SAFE_ALLOCATE(ph_old, (cng)) 259 SAFE_ALLOCATE(isort_old, (gvec%ng)) 260 endif 261 if (ifirst .eq. 0) then 262 SAFE_ALLOCATE(ind_old, (cng)) 263 SAFE_ALLOCATE(ph_old, (cng)) 264 SAFE_ALLOCATE(isort_old, (gvec%ng)) 265 SAFE_ALLOCATE(ind_old2, (cng)) 266 SAFE_ALLOCATE(ph_old2, (cng)) 267 SAFE_ALLOCATE(isort_old2, (gvec%ng)) 268 endif 269 270 wfnc%ng=cng 271 wfnc%nband=1 272 wfnc%nspin=cns 273 wfnc%nspinor=cnsp 274 275 if (cns.ne.nspin) then 276 write(errmsg,'(2a,2i2)') 'spin number mismatch: ', nspin, cns 277 call die(errmsg, only_root_writes = .true.) 278 endif 279 280 SAFE_ALLOCATE(wfnc%isort, (gvec%ng)) 281 282 isort(:)=intwfnc%isort(:,ikown) 283 if (xct%icpar .eq. 1) then 284 SAFE_ALLOCATE(wfnc%cg, (wfnc%ng,1,wfnc%nspin*wfnc%nspinor)) 285 ccg(1:cng,1,:) = intwfnc%cg(1:cng,icown,:) 286 else 287 SAFE_ALLOCATE(wfnc%cg, (wfnc%ng,xct%ncb_co,wfnc%nspin*wfnc%nspinor)) 288 ccg(1:cng,:,:) = intwfnc%cg(1:cng,icown:icown+xct%ncb_co-1,:) 289 endif 290 291! JRD: Below is now necessary because kg might be different from kgq 292! Compute inverse index array of Fourier components around rk-kpoint 293 294 if (ik .ne. ikold .and. ik .ne. ikold2) then 295 SAFE_ALLOCATE(isorti, (gvec%ng)) 296 isorti(:)=0 297 do ii=1,wfnc%ng 298 isorti(isort(ii))=ii 299 enddo 300 endif 301 302! Compute index array of Fourier components around fk-kpoint 303 304 if (ik .ne. ikold .and. ik .ne. ikold2) then 305 SAFE_ALLOCATE(ekin, (gvec%ng)) 306 call kinetic_energies(gvec, crys%bdot, ekin, qvec = kg%f(:, ik)) 307 call sortrx(gvec%ng,ekin,isort) 308 SAFE_DEALLOCATE(ekin) 309 isort_old(:) = isort(:) 310 else if (ik .eq. ikold) then 311 isort(:)=isort_old(:) 312 else 313 isort(:)=isort_old2(:) 314 endif 315 316! Find ind and ph relating wavefunctions in fk to rk-kpoint 317 318 if (ik .ne. ikold .and. ik .ne. ikold2) then 319 ind=0 320 ph=ZERO 321 call gmap(gvec,syms,wfnc%ng,kg%itran(ik), & 322 kg%kg0(:,ik),isort,isorti,ind,ph,.true.) 323 SAFE_DEALLOCATE(isorti) 324 ind_old(:)=ind(:) 325 ph_old(:)=ph(:) 326 else if (ik .eq. ikold) then 327 ind(:)=ind_old(:) 328 ph(:)=ph_old(:) 329 else 330 ind(:)=ind_old2(:) 331 ph(:)=ph_old2(:) 332 endif 333 334! Compute and renormalize conduction wavefunctions 335 336 if (xct%icpar .eq. 1) incband =1 337 if (xct%icpar .eq. 0) incband = xct%ncb_co 338 339 do ijk =1, incband 340 do kk=1,wfnc%nspin 341 do jj=1,wfnc%nspinor 342 do ii=1,wfnc%ng 343 if (ind(ii) .gt. 0) then 344 wfnc%cg(ii,ijk,kk*jj)=ph(ii)*ccg(ind(ii),ijk,kk*jj) 345 else 346 wfnc%cg(ii,ijk,kk*jj)=ZERO 347 endif 348 enddo 349 enddo ! jj 350 call checknorm('wfnc%cg',ijk,ik,wfnc%ng,kk,wfnc%nspinor,wfnc%cg(1:wfnc%ng,ijk,:)) 351 enddo ! kk 352 enddo ! ijk 353 ccg=wfnc%cg 354 355 356 wfnc%cg=ccg 357 wfnc%isort=isort 358 359 SAFE_DEALLOCATE(ccg) 360 SAFE_DEALLOCATE(ind) 361 SAFE_DEALLOCATE(ph) 362 SAFE_DEALLOCATE(isort) 363 364 if (ik .ne. ikold .and. ik .ne. ikold2) then 365 ikold2=ikold 366 ikold=ik 367 cngold=cng 368 vngold=vng 369 endif 370 371 ifirst=-1 372 373 POP_SUB(genwf_kernel) 374 375 return 376end subroutine genwf_kernel 377