1!========================================================================== 2! 3! Routines: 4! 5! (1) genwf_mpi() Originally By (JRD) Last Modified 11/2009 (JRD) 6! 7! On entry: 8! qq = current q-vector 9! rk = current k point in irr. zone 10! irk = its index 11! 12! On exit: 13! vwfn%ev and vwfn%zv hold eigenvalues and wavefunctions (valence) 14! cwfn%ec and cwfn%zc hold eigenvalues and wavefunctions (conduction) 15! 16! with proper phases and ordering for the k-point rk (given the 17! data for irr. zone) 18! 19! subroutine generates valence-band wavefunctions for rk(irk) 20! and conduction-band wavefunctions for rk(irk) from the 21! wavefunctions available for rk in the irreducible wedge of the 22! bz 23! 24! i rk rk(irk) 25! o c ngv,...,ev valence-band wavefunctions for rk+q 26! and associated data 27! o c ngc,...,ec conduction-band wavefunctions for rk 28! and associated data 29! 30!========================================================================== 31 32#include "f_defs.h" 33 34module genwf_mpi_m 35 36 use global_m 37 use find_kpt_match_m 38 use gmap_m 39 use input_utils_m 40 use misc_m 41 use sort_m 42 use timing_m, only: timing => epsilon_timing 43 44 implicit none 45 46 private 47 integer, allocatable :: GK(:) 48 logical :: GK_allocated=.false. 49 public :: genwf_mpi 50 51contains 52 53subroutine genwf_mpi(syms,gvec,crys,kp,kpq,irk,rk,qq,vwfn,pol,cwfn,use_wfnq,intwfnv,intwfnvq,intwfnc,ivin) 54 type (symmetry), intent(in) :: syms 55 type (gspace), intent(in) :: gvec 56 type (crystal), intent(in) :: crys 57 type (kpoints), target, intent(in) :: kp 58 type (kpoints), target, intent(in) :: kpq 59 integer, intent(in) :: irk 60 real(DP), intent(in) :: rk(3) 61 real(DP), intent(in) :: qq(3) 62 type (valence_wfns), intent(inout) :: vwfn 63 type (polarizability), intent(in) :: pol 64 type (conduction_wfns), intent(inout) :: cwfn 65 logical, intent(in) :: use_wfnq 66 type (int_wavefunction), intent(in) :: intwfnv 67 type (int_wavefunction), intent(in) :: intwfnvq 68 type (int_wavefunction), intent(in) :: intwfnc 69 integer, intent(in) :: ivin 70 71 integer :: itval,jband 72 integer :: ng 73 integer, allocatable :: isortc(:) 74 real(DP), allocatable :: eig(:,:) 75 SCALAR, allocatable :: zin(:,:),zinc(:,:) 76 SCALAR, allocatable :: zintemp(:,:) 77 type(kpoints), pointer :: kp_point 78 79 character :: tmpstr*120 80 integer :: itqq,i 81 integer :: n,ig,ispin,iband 82 integer :: naddc 83 integer :: ikrkq,kgq(3),kgqq(3) 84 integer, allocatable :: ind(:),isorti(:) 85 real(DP), allocatable :: xnorm(:) 86 real(DP) :: qk(3),rkq(3) 87 real(DP), allocatable :: ekin(:) 88 real(DP) :: rkmatch(3) 89 90 SCALAR, allocatable :: ph(:) 91 92 SCALAR, allocatable, save :: ph_old(:) 93 94 integer, allocatable, save :: ind_old(:) 95 integer, save :: irk_old=0 96 97 PUSH_SUB(genwf_mpi) 98 99 if (.not.GK_allocated) then 100 SAFE_ALLOCATE(GK, (gvec%ng)) 101 call get_GK_array_from_gvecs(gvec%ng, gvec%components, GK) 102 GK_allocated = .true. 103 endif 104 SAFE_ALLOCATE(isortc, (gvec%ng)) 105 SAFE_ALLOCATE(eig, (cwfn%nband,kp%nspin)) 106 SAFE_ALLOCATE(xnorm, (kp%nspin)) 107 108 if(peinf%inode.eq.0) call timing%start(timing%genwf_val) 109 110! rkq = rk + qq (i.e. rkq = current kpoint in irr. zone + qq) 111 112 rkq(1:3) = rk(1:3) + qq(1:3) 113 114! Should we look for the valence WFNs in WFNq or WFN? 115 if (use_wfnq) then 116 kp_point => kpq 117 else 118 kp_point => kp 119 endif 120 121! We used to assume WFNq contained the needed q->0 point, 122! but it is better to unfold the shifted grid for more flexibility. 123 124 call find_kpt_match(kp_point, syms, rkq, ikrkq, itqq, kgqq) 125 126 if(ikrkq == 0) then 127 if(peinf%inode == 0) write(0,'(a,3f12.6,/,a,3f12.6)') 'rk = ', rk(:), 'qq = ', qq(:) 128 write(tmpstr,'(a,3f8.3,a)') 'genwf_mpi: No match for rkq point:',rkq,' in file ' 129 if(use_wfnq) then 130 write(tmpstr,'(a,a)') TRUNC(tmpstr), ' WFNq' 131 else 132 write(tmpstr,'(a,a)') TRUNC(tmpstr), ' WFN' 133 endif 134 call die(tmpstr, only_root_writes = .true.) 135 endif 136 137 rkmatch(1:3) = kp_point%rk(1:3, ikrkq) 138 vwfn%idx_kp = ikrkq 139 140! if(peinf%inode.eq.0) then 141 itval=vwfn%nband+pol%ncrit 142 if (.not.use_wfnq) then 143 ng = intwfnv%ng(ikrkq) 144 if (irk .ne. irk_old) then 145 eig(1:itval,1:kp%nspin)=kp_point%el(1+vwfn%ncore_excl:itval+vwfn%ncore_excl,ikrkq,1:kp%nspin) 146 isortc(1:ng)=intwfnv%isort(1:ng,ikrkq) 147 endif 148 qk(:)=intwfnv%qk(:,ikrkq) 149 else 150 ng = intwfnvq%ng(ikrkq) 151 if (irk .ne. irk_old) then 152 eig(1:itval,1:kp%nspin)=kp_point%el(1+vwfn%ncore_excl:itval+vwfn%ncore_excl,ikrkq,1:kp%nspin) 153 isortc(1:ng)=intwfnvq%isort(1:ng,ikrkq) 154 endif 155 qk(:)=intwfnvq%qk(:,ikrkq) 156 endif 157! endif 158 159 SAFE_ALLOCATE(zintemp, (ng,kp%nspin*kp%nspinor)) 160 SAFE_ALLOCATE(zin, (ng,kp%nspin*kp%nspinor)) 161 162 jband = (ikrkq-1)*peinf%nvownactual + ivin 163 zintemp=0D0 164 if (use_wfnq) then 165 zintemp(1:ng,1:kp%nspin*kp%nspinor)=intwfnvq%cg(1:ng,jband,1:kp%nspin*kp%nspinor) 166 else 167 zintemp(1:ng,1:kp%nspin*kp%nspinor)=intwfnv%cg(1:ng,jband,1:kp%nspin*kp%nspinor) 168 endif 169 170 zin(:,:)=zintemp(:,:) 171 SAFE_DEALLOCATE(zintemp) 172 173! Check kpoint 174 175 if(any(abs(rkmatch(1:3) - qk(1:3)) .gt. TOL_Small)) call die('genwf_mpi: rkmatch') 176 177 if (irk .ne. irk_old) then 178 179! Compute kinetic energies for rkq+g 180 181 SAFE_ALLOCATE(ekin, (gvec%ng)) 182 if (peinf%inode == 0) call timing%start(timing%genwf_ekin) 183 call kinetic_energies(gvec, crys%bdot, ekin, qvec = rkq) 184 if (peinf%inode == 0) call timing%stop(timing%genwf_ekin) 185 186! sort array ekin to ascending order, store indices in array vwfn%isort 187! WARNING: one should not assume that in the case of 188! q-->0 the orders as provided below and as read in from 189! WFNq file is the same (sorting may change....) 190! ===> need to call gmap also for q-->0 (see below) 191! EKC: We initialize vwfn%isort to the appropriate array 192! before reading in. this way we do not get zeroes in the array 193! these are the valence wave-functions that do not need 194! to be changed 195 196 if(peinf%inode.eq.0) call timing%start(timing%genwf_sort) 197 call sort_threaded(gvec%ng, ekin, vwfn%isort, GK) 198 !call sortrx(gvec%ng,ekin,vwfn%isort,gvec=gvec%components) 199 if(peinf%inode.eq.0) call timing%stop(timing%genwf_sort) 200 201 SAFE_ALLOCATE(isorti, (gvec%ng)) 202 do i=1,gvec%ng 203 isorti(vwfn%isort(i))=i 204 enddo 205 do i=1,ng 206 isorti(isortc(i))=i 207 enddo 208 209 vwfn%ngv=ng 210 211! SIB: put read eigenvalues into vwfn%ev(band,spin). 212 213 SAFE_ALLOCATE(vwfn%ev, ((vwfn%nband+pol%ncrit),kp%nspin)) 214 vwfn%ev(1:(vwfn%nband+pol%ncrit),:) = eig(1:(vwfn%nband+pol%ncrit),:) 215 216! JRD: Map planewave components for rk+q, to those of rk 217! (even for q--> 0) 218! 219! SIB: get phases (ph) and indices (ind) for g-vectors 220! gvec%components(:,vwfn%isort(1:vwfn%ngv))+kgqq 221 222 SAFE_ALLOCATE(ind, (vwfn%ngv)) 223 SAFE_ALLOCATE(ph, (vwfn%ngv)) 224 call gmap(gvec,syms,vwfn%ngv,itqq,kgqq,vwfn%isort,isorti,ind,ph,.true.) 225 if (irk_old .ne. 0) then 226 SAFE_DEALLOCATE(ph_old) 227 SAFE_DEALLOCATE(ind_old) 228 endif 229 SAFE_ALLOCATE(ph_old, (vwfn%ngv)) 230 SAFE_ALLOCATE(ind_old, (vwfn%ngv)) 231 ph_old(:) = ph(:) 232 ind_old(:) = ind(:) 233 else 234 235 SAFE_ALLOCATE(ph, (vwfn%ngv)) 236 SAFE_ALLOCATE(ind, (vwfn%ngv)) 237 ph(:)=ph_old(:) 238 ind(:)=ind_old(:) 239 240 endif ! irk = irk_old 241 SAFE_DEALLOCATE(eig) 242 243 SAFE_ALLOCATE(vwfn%zv, (vwfn%ngv,kp%nspin*kp%nspinor)) 244 245! XAV: vwfn%zv(ig) corresponds really to the 246! vwfn%isort(ig) G-vector (between 1 and ng) 247! The subroutine gmap assumes that, as read from WFNq or WFN, 248! zin(ig) corresponds really to isortc(ig) G-vector !!! 249! 250! SIB: checks that zin vectors have norm greater than 1e-6, and then 251! normalizes them to have unit square modulus. 252 253 do ispin=1,kp%nspin*kp%nspinor 254 vwfn%zv(:,ispin)=zin(ind(:),ispin)*ph(:) 255 enddo 256 257 258! BAB: we check if the norm differs appreciably from unity. 259! there is no longer a need to further normalize the vector 260 if (peinf%check_norms) then 261 do ispin=1,kp%nspin 262 call compute_norm(xnorm(ispin),ispin,vwfn%ngv,kp%nspinor,vwfn%zv(:,:)) 263 if(abs(xnorm(ispin) - 1d0) > TOL_Small) then 264 write(0,*) 'Bad norm: sym op=',itqq,'ik=',ikrkq,'ispin=',ispin,'norm=',xnorm(ispin) 265 call die('genwf_mpi: Bad norm') 266 endif 267 enddo 268 endif 269 270! End calculation of valence-band wavefunctions 271 272 SAFE_DEALLOCATE(zin) 273 SAFE_DEALLOCATE(ind) 274 SAFE_DEALLOCATE(ph) 275 SAFE_DEALLOCATE(xnorm) 276 if(peinf%inode.eq.0) call timing%stop(timing%genwf_val) 277 278 !FHJ: only generate conduction WFNs if ivin=1 279 if(irk .ne. irk_old) then 280 281!------------------------------------------------------------------------ 282! Generate conduction-band wavefunctions for rk 283! find rk, r, g0 such that rk=r(rk)+g0 284 285! SIB: This seems redundant, but find a k-point and symmetry so that 286! rk = sym%mtrx(:,:,itqq)*kp%rk(irk_,:) + kgq, where kgq is integer 3-vec 287 288 if(peinf%inode.eq.0) call timing%start(timing%genwf_cond) 289 290 call find_kpt_match(kp, syms, rk, ikrkq, itqq, kgq) 291 cwfn%idx_kp = ikrkq 292 293 if(ikrkq == 0) call die('genwf_mpi: kgq mismatch') 294 295! Write out rk, it and kgq 296 297 if (peinf%verb_debug .and. peinf%inode==0) then 298 write(6,7000) (rk(i),i=1,3),ikrkq,(kp%rk(i,ikrkq),i=1,3),itqq,(kgq(i),i=1,3) 2997000 format(1x,'rk=',3f7.3,1x,'irk_=',i5,1x,' rk=',3f7.3,1x,'it=',i5,1x,'kg0=',3i3) 300 endif 301 302! SIB: if we already found this k-point last time, get its qk, ng, 303! and isortc(:). Otherwise, skip ikrkq-1 records, and read in information (qk,cwfn%ec,ng,isortc), 304 305 SAFE_ALLOCATE(cwfn%ec, (cwfn%nband,kp%nspin)) 306 307 ng = intwfnc%ng(ikrkq) 308 cwfn%ec(1:cwfn%nband,1:kp%nspin)=kp%el(1+vwfn%ncore_excl:cwfn%nband+vwfn%ncore_excl,ikrkq,1:kp%nspin) 309 qk(:)=intwfnc%qk(:,ikrkq) 310 isortc(1:ng)=intwfnc%isort(1:ng,ikrkq) 311 312! Check kpoint (again ... boring...) 313! Check that kp%rk(:,ikrkq) = qk (why it wouldn`t is a mystery!) 314 315 if(any(abs(kp%rk(1:3, ikrkq) - qk(1:3)) .gt. TOL_Small)) & 316 call die('genwf_mpi: qk mismatch') 317 318 cwfn%ngc=ng 319 320! Compute inverse to isort 321! NOTE: isortc orders |kp%rk+G|^2 322! It is not necessarily the same order as |rk+G|^2 323! (particularly if umklapp, ie kgq non zero) 324 325 if(peinf%inode.eq.0) call timing%start(timing%genwf_ekin) 326 call kinetic_energies(gvec, crys%bdot, ekin, qvec = rk) 327 if(peinf%inode.eq.0) call timing%stop(timing%genwf_ekin) 328 329! Sort array ekin to ascending order 330! store indices in array isort 331 332 if(peinf%inode.eq.0) call timing%start(timing%genwf_sort) 333 call sort_threaded(gvec%ng, ekin, cwfn%isort, GK) 334 !call sortrx(gvec%ng,ekin,cwfn%isort,gvec=gvec%components) 335 if(peinf%inode.eq.0) call timing%stop(timing%genwf_sort) 336 337 do i=1,ng 338 isorti(isortc(i))=i 339 enddo 340 SAFE_DEALLOCATE(ekin) 341 342! map planewave components for rk to those of rk 343! compute phases 344! We do not the isorti related to kp%rk BUT the cwfn%isort related to rk 345 346 SAFE_ALLOCATE(ind, (cwfn%ngc)) 347 SAFE_ALLOCATE(ph, (cwfn%ngc)) 348 call gmap(gvec,syms,cwfn%ngc,itqq,kgq,cwfn%isort,isorti,ind,ph,.true.) 349 SAFE_DEALLOCATE(isorti) 350 351! generate conduction-band wavefunctions 352! loop over wavefunctions 353! read conduction-band from file one by one 354 355 SAFE_ALLOCATE(cwfn%zc, (peinf%ncownactual*cwfn%ngc,kp%nspin*kp%nspinor)) 356 SAFE_ALLOCATE(zinc, (cwfn%ngc,kp%nspin*kp%nspinor)) 357 358 do n=1,peinf%ncownactual 359 360 jband= (ikrkq-1)*peinf%ncownactual + n 361 362 iband = intwfnc%cbi(jband) 363 zinc(1:ng,1:kp%nspin*kp%nspinor)=intwfnc%cg(1:ng,jband,1:kp%nspin*kp%nspinor) 364 365 if(iband .ne. peinf%invindexc(n)) call die('genwf_mpi: invindexc mismatch') 366 367 naddc=(n-1)*cwfn%ngc 368 369! 370! Loop over components of wfns 371! note that only conduction-band wfns are stored- 372! they start in the 1st position with the state nvband+1 373! 374 375 do ig=1,cwfn%ngc 376 cwfn%zc(naddc+ig,:)=ph(ig)*zinc(ind(ig),:) 377 enddo 378 379 380 enddo ! n (cond-bands per node [ncownactual] loop) 381 382 SAFE_DEALLOCATE(zinc) 383 SAFE_DEALLOCATE(ind) 384 SAFE_DEALLOCATE(ph) 385 386! end generation of conduction-band wavefunctions for rk 387 388 if(peinf%inode.eq.0) call timing%stop(timing%genwf_cond) 389 390 endif ! (irk .ne. irk_old) 391 392 irk_old = irk 393 SAFE_DEALLOCATE(isortc) 394 395 POP_SUB(genwf_mpi) 396 397 return 398end subroutine genwf_mpi 399 400end module genwf_mpi_m 401