1!----------------------------------------------------------------------- 2subroutine read_export (pp_file,kunit,uspp_spsi, ascii, single_file, raw) 3 !----------------------------------------------------------------------- 4 ! 5 use kinds, ONLY : DP 6 use pwcom 7 use control_flags, ONLY : gamma_only 8 use becmod, ONLY : bec_type, becp, calbec, & 9 allocate_bec_type, deallocate_bec_type 10! use symme, ONLY : nsym, s, invsym, sname, irt, ftau 11! use symme, ONLY : nsym, s, invsym, irt, ftau 12! use char, ONLY : sname 13! sname is in symme which is now outside pwcom 14 use uspp, ONLY : nkb, vkb 15 use wavefunctions, ONLY : evc 16 use io_files, ONLY : prefix, iunwfc, nwordwfc, iunsat, nwordatwfc 17 use io_files, ONLY : pseudo_dir, psfile 18 use io_global, ONLY : ionode, stdout 19 USE ions_base, ONLY : atm, nat, ityp, tau, nsp 20 use mp_pools, ONLY : nproc_pool, my_pool_id, intra_pool_comm, inter_pool_comm 21 USE mp_world, ONLY : world_comm, nproc 22 use mp, ONLY : mp_sum, mp_max 23 use ldaU, ONLY : lda_plus_u 24 USE gvecw, ONLY : ecutwfc 25 USE cell_base, ONLY : bg 26 USE gvect, ONLY : ngm, ngm_g,g, gstart, ig_l2g, mill 27! USE pwcom, ONLY : igk_k 28 USE cell_base, ONLY : tpiba2 29 30 implicit none 31 32 CHARACTER(5), PARAMETER :: fmt_name="QEXPT" 33 CHARACTER(5), PARAMETER :: fmt_version="1.1.0" 34 35 integer, intent(in) :: kunit 36 character(80), intent(in) :: pp_file 37 logical, intent(in) :: uspp_spsi, ascii, single_file, raw 38 39 integer :: i, j, k, ig, ik, ibnd, na, ngg,ig_, ierr 40 integer, allocatable :: kisort(:) 41 real(DP) :: xyz(3), tmp(3) 42 integer :: npool, nkbl, nkl, nkr, npwx_g 43 integer :: ike, iks, npw_g, ispin, local_pw 44 integer, allocatable :: ngk_g( : ) 45 integer, allocatable :: itmp_g( :, : ) 46 real(DP),allocatable :: rtmp_g( :, : ) 47 real(DP),allocatable :: rtmp_gg( : ) 48 integer, allocatable :: itmp1( : ) 49 integer, allocatable :: igwk( :, : ) 50 integer, allocatable :: l2g_new( : ) 51 integer, allocatable :: igk_l2g( :, : ) 52 53 54 real(DP) :: wfc_scal 55 logical :: twf0, twfm 56 complex(DP), allocatable :: sevc (:,:) 57 58 CHARACTER(LEN=256) :: outdir 59 60 write(stdout,*) "nkstot=", nkstot 61 62 IF( nkstot > 0 ) THEN 63 64 IF( ( kunit < 1 ) .OR. ( MOD( nkstot, kunit ) /= 0 ) ) & 65 CALL errore( ' write_export ',' wrong kunit ', 1 ) 66 67 IF( ( nproc_pool > nproc ) .OR. ( MOD( nproc, nproc_pool ) /= 0 ) ) & 68 CALL errore( ' write_export ',' nproc_pool ', 1 ) 69 70 ! find out the number of pools 71 npool = nproc / nproc_pool 72 73 ! find out number of k points blocks 74 nkbl = nkstot / kunit 75 76 ! k points per pool 77 nkl = kunit * ( nkbl / npool ) 78 79 ! find out the reminder 80 nkr = ( nkstot - nkl * npool ) / kunit 81 82 ! Assign the reminder to the first nkr pools 83 IF( my_pool_id < nkr ) nkl = nkl + kunit 84 85 ! find out the index of the first k point in this pool 86 iks = nkl * my_pool_id + 1 87 IF( my_pool_id >= nkr ) iks = iks + nkr * kunit 88 89 ! find out the index of the last k point in this pool 90 ike = iks + nkl - 1 91 92 END IF 93 94 write(stdout,*) "after first init" 95 96 ! find out the global number of G vectors: ngm_g 97 ngm_g = ngm 98 call mp_sum( ngm_g , intra_pool_comm) 99 100 ! collect all G vectors across processors within the pools 101 ! and compute their modules 102 ! 103 allocate( itmp_g( 3, ngm_g ) ) 104 allocate( rtmp_g( 3, ngm_g ) ) 105 allocate( rtmp_gg( ngm_g ) ) 106 107 itmp_g = 0 108 do ig = 1, ngm 109 itmp_g( 1, ig_l2g( ig ) ) = mill(1, ig ) 110 itmp_g( 2, ig_l2g( ig ) ) = mill(2, ig ) 111 itmp_g( 3, ig_l2g( ig ) ) = mill(3, ig ) 112 end do 113 call mp_sum( itmp_g , intra_pool_comm) 114 ! 115 ! here we are in crystal units 116 rtmp_g(1:3,1:ngm_g) = REAL( itmp_g(1:3,1:ngm_g) ) 117 ! 118 ! go to cartesian units (tpiba) 119 call cryst_to_cart( ngm_g, rtmp_g, bg , 1 ) 120 ! 121 ! compute squared moduli 122 do ig = 1, ngm_g 123 rtmp_gg(ig) = rtmp_g(1,ig)**2 + rtmp_g(2,ig)**2 + rtmp_g(3,ig)**2 124 enddo 125 deallocate( rtmp_g ) 126 127 ! build the G+k array indexes 128 allocate ( igk_l2g ( npwx, nks ) ) 129 allocate ( kisort( npwx ) ) 130 do ik = 1, nks 131 kisort = 0 132 npw = npwx 133 call gk_sort (xk (1, ik+iks-1), ngm, g, ecutwfc / tpiba2, npw, kisort(1), g2kin) 134 ! 135 ! mapping between local and global G vector index, for this kpoint 136 ! 137 DO ig = 1, npw 138 ! 139 igk_l2g(ig,ik) = ig_l2g( kisort(ig) ) 140 ! 141 END DO 142 ! 143 igk_l2g( npw+1 : npwx, ik ) = 0 144 ! 145 ngk (ik) = npw 146 end do 147 deallocate (kisort) 148 149 ! compute the global number of G+k vectors for each k point 150 allocate( ngk_g( nkstot ) ) 151 ngk_g = 0 152 ngk_g( iks:ike ) = ngk( 1:nks ) 153 CALL mp_sum( ngk_g,world_comm ) 154 155 ! compute the Maximum G vector index among all G+k and processors 156 npw_g = MAXVAL( igk_l2g(:,:) ) 157 CALL mp_max( npw_g,world_comm ) 158 159 ! compute the Maximum number of G vector among all k points 160 npwx_g = MAXVAL( ngk_g( 1:nkstot ) ) 161 162 deallocate(rtmp_gg) 163 164 allocate( igwk( npwx_g,nkstot ) ) 165 166 write(stdout,*) "after g stuff" 167 168! wfc grids 169 170 DO ik = 1, nkstot 171 igwk(:,ik) = 0 172 ! 173 ALLOCATE( itmp1( npw_g ), STAT= ierr ) 174 IF ( ierr/=0 ) CALL errore('pw_export','allocating itmp1', ABS(ierr) ) 175 itmp1 = 0 176 ! 177 IF( ik >= iks .AND. ik <= ike ) THEN 178 DO ig = 1, ngk( ik-iks+1 ) 179 itmp1( igk_l2g( ig, ik-iks+1 ) ) = igk_l2g( ig, ik-iks+1 ) 180 END DO 181 END IF 182 ! 183 CALL mp_sum( itmp1 ,world_comm) 184 ! 185 ngg = 0 186 DO ig = 1, npw_g 187 IF( itmp1( ig ) == ig ) THEN 188 ngg = ngg + 1 189 igwk( ngg , ik) = ig 190 END IF 191 END DO 192 IF( ngg /= ngk_g( ik ) ) THEN 193 WRITE( stdout,*) ' ik, ngg, ngk_g = ', ik, ngg, ngk_g( ik ) 194 END IF 195 ! 196 DEALLOCATE( itmp1 ) 197 ! 198 ENDDO 199 ! 200 deallocate( itmp_g ) 201 202 write(stdout,*)"after wfc waves" 203 204#ifdef __MPI 205 call poolrecover (et, nbnd, nkstot, nks) 206#endif 207 208 wfc_scal = 1.0d0 209 twf0 = .true. 210 twfm = .false. 211 212 do ik = 1, nkstot 213 local_pw = 0 214 IF( (ik >= iks) .AND. (ik <= ike) ) THEN 215 216 call davcio (evc, 2*nwordwfc, iunwfc, (ik-iks+1), - 1) 217 !IF ( lda_plus_u ) CALL davcio( swfcatom, nwordatwfc, iunsat, (ik-iks+1), -1 ) 218 local_pw = ngk(ik-iks+1) 219 220 ENDIF 221 222 223 allocate(l2g_new(local_pw)) 224 225 l2g_new = 0 226 do ig = 1, local_pw 227 ngg = igk_l2g(ig,ik-iks+1) 228 do ig_ = 1, ngk_g(ik) 229 if(ngg == igwk(ig_,ik)) then 230 l2g_new(ig) = ig_ 231 exit 232 end if 233 end do 234 end do 235 236 237 ispin = isk( ik ) 238 ! WRITE(0,*) ' ### ', ik,nkstot,iks,ike,kunit,nproc,nproc_pool 239 deallocate(l2g_new) 240 end do 241 ! 242 243 write(stdout,*) "after davcio" 244 245 ! If specified and if USPP are used the wfcs S_psi are written 246 ! | spsi_nk > = \hat S | psi_nk > 247 ! where S is the overlap operator of US PP 248 ! 249 IF ( uspp_spsi .AND. nkb > 0 ) THEN 250 251 ALLOCATE( sevc(npwx,nbnd), STAT=ierr ) 252 IF (ierr/=0) CALL errore( ' read_export ',' Unable to allocate SEVC ', ABS(ierr) ) 253 254 CALL allocate_bec_type (nkb,nbnd,becp) 255 256 do ik = 1, nkstot 257 258 local_pw = 0 259 IF( (ik >= iks) .AND. (ik <= ike) ) THEN 260 261 CALL gk_sort (xk (1, ik+iks-1), ngm, g, ecutwfc / tpiba2, npw, igk_k(1,ik), g2kin) 262 CALL davcio (evc, 2*nwordwfc, iunwfc, (ik-iks+1), - 1) 263 264 CALL init_us_2(npw, igk_k(1,ik), xk(1, ik), vkb) 265 local_pw = ngk(ik-iks+1) 266 267 IF ( gamma_only ) THEN 268 if(nkb>0) CALL calbec ( ngk_g(ik), vkb, evc, becp ) 269 ELSE 270 CALL calbec ( npw, vkb, evc, becp ) 271 ENDIF 272 CALL s_psi(npwx, npw, nbnd, evc, sevc) 273 ENDIF 274 275 ALLOCATE(l2g_new(local_pw)) 276 277 l2g_new = 0 278 DO ig = 1, local_pw 279 ngg = igk_l2g(ig,ik-iks+1) 280 DO ig_ = 1, ngk_g(ik) 281 IF(ngg == igwk(ig_,ik)) THEN 282 l2g_new(ig) = ig_ 283 EXIT 284 ENDIF 285 ENDDO 286 ENDDO 287 288 ispin = isk( ik ) 289 DEALLOCATE(l2g_new) 290 ENDDO 291 292 DEALLOCATE( sevc, STAT=ierr ) 293 IF ( ierr/= 0 ) CALL errore('read_export','Unable to deallocate SEVC',ABS(ierr)) 294 CALL deallocate_bec_type ( becp ) 295 ENDIF 296 297 DEALLOCATE( igk_l2g ) 298 DEALLOCATE( igwk ) 299 DEALLOCATE ( ngk_g ) 300 301end subroutine read_export 302