1! 2! Copyright (C) 2001-2013 Quantum ESPRESSO group 3! This file is distributed under the terms of the 4! GNU General Public License. See the file `License' 5! in the root directory of the present distribution, 6! or http://www.gnu.org/copyleft/gpl.txt . 7! 8! 9 10 11 12 subroutine write_wfc_grid_2 13!this subroutine read real wavefunctions from file 14!on the small charge grid, and write on the 15!wavefunction grid in real space 16 17 18 USE kinds, ONLY : DP 19 USE io_files, ONLY : diropn 20 USE io_global, ONLY : stdout 21 USE gvecs, ONLY : doublegrid 22 USE wvfct, ONLY : nbnd 23 USE fft_base, ONLY : dfftp, dffts 24 25 implicit none 26 27 INTEGER, EXTERNAL :: find_free_unit 28 29 INTEGER :: iw,ix,iy,iz,nn 30 REAL(kind=DP), ALLOCATABLE :: tmprealis(:),tmpreal2(:) 31 INTEGER :: iunwfcreal, iunwfcreal2 32 INTEGER :: iqq 33 LOGICAL :: exst 34 INTEGER :: nrxxs2 35 REAL(kind=DP) :: sca 36 37 nrxxs2=(dffts%nr1/2+1)*(dffts%nr2/2+1)*(dffts%nr3/2+1) 38 39 iunwfcreal=find_free_unit() 40 CALL diropn( iunwfcreal, 'real_whole', dffts%nnr, exst ) 41 42 iunwfcreal2=find_free_unit() 43 CALL diropn( iunwfcreal2, 'real_whole-2', nrxxs2, exst ) 44 45 46 allocate(tmprealis(dffts%nnr)) 47 allocate(tmpreal2(nrxxs2)) 48 49 do iw=1,nbnd 50 CALL davcio( tmprealis,dffts%nnr,iunwfcreal,iw,-1) 51 tmpreal2(:)=0.d0 52 iqq=0 53 sca=0.d0 54 do ix=1,dffts%nr1,2 55 do iy=1,dffts%nr2,2 56 do iz=1,dffts%nr3,2 57 iqq=iqq+1 58 nn=(iz-1)*dffts%nr1*dffts%nr2+(iy-1)*dffts%nr1+ix 59 tmpreal2(iqq)=tmprealis(nn) 60 sca=sca+tmprealis(nn)**2.d0!ATTENZIONE 61 enddo 62 enddo 63 enddo 64 !tmpreal2(:)=tmpreal2(:)/(sqrt(sca/dble(iqq))) 65 write(*,*) 'MODULUS', iw,sca/dble(iqq) 66 CALL davcio( tmpreal2,nrxxs2,iunwfcreal2,iw,1) 67 enddo 68 69 deallocate(tmprealis,tmpreal2) 70 close(iunwfcreal) 71 close(iunwfcreal2) 72 return 73 end subroutine 74 75!----------------------------------------------------------------------- 76subroutine matrix_wannier_gamma_big( matsincos, ispin, n_set, itask ) 77 !----------------------------------------------------------------------- 78 ! 79 !this subroutine calculates the terms <Psi_i|exp(iGX)|Psi_j> 80 !in real space for gamma only case 81 82 83 84 USE kinds, ONLY : DP 85 USE cell_base, ONLY: at, alat, tpiba, omega, tpiba2 86 USE constants, ONLY : e2, pi, tpi, fpi 87 USE uspp, ONLY : okvan, nkb 88 USE io_files, ONLY : diropn 89 USE io_global, ONLY : stdout 90 USE gvecs, ONLY : doublegrid 91 ! USE realus, ONLY : qsave, box,maxbox 92 USE wannier_gw, ONLY : becp_gw, expgsave, becp_gw_c, maxiter2,num_nbndv 93 USE ions_base, ONLY : nat, ntyp =>nsp, ityp 94 USE uspp_param, ONLY : lmaxq,upf,nh, nhm 95 USE lsda_mod, ONLY : nspin 96 USE mp_pools, ONLY : me_pool 97 USE mp, ONLY : mp_bcast,mp_barrier,mp_sum 98 USE mp_world, ONLY : world_comm 99 USE fft_base, ONLY : dffts,dfftp 100 USE wvfct, ONLY : nbnd 101 102 implicit none 103 104 105 INTEGER, EXTERNAL :: find_free_unit 106 107 ! 108 INTEGER, INTENT(in) :: ispin!spin polarization considred 109! COMPLEX(dp), INTENT(out) :: matp(nbnd_normal,nbnd_normal,3) 110 REAL(dp), INTENT(out) :: matsincos(nbnd,nbnd,6) 111 INTEGER, INTENT(in) :: n_set !defines the number of states 112 INTEGER, INTENT(in) :: itask !if ==1 consider subspace {C'} 113 114 INTEGER :: iiw,jjw, jw_begin 115 INTEGER :: iw,jw,ir,ix,iy,iz,nn,ii 116 REAL(kind=DP), ALLOCATABLE :: tmprealis(:,:),tmprealjs(:,:), tmpreal(:) 117 COMPLEX(kind=DP), ALLOCATABLE :: tmpexp(:), tmpexp2(:,:) 118 INTEGER :: iunwfcreal2 119 COMPLEX(kind=DP) :: sca,ee,sca1 120 REAL(kind=DP) :: dsgn 121 LOGICAL :: exst 122 INTEGER :: iqq 123 INTEGER :: na, ih, jh, np 124 INTEGER :: ikb, jkb, ijkb0, is 125 INTEGER :: isgn,mdir 126 INTEGER :: nr3s_start, nr3s_end 127 INTEGER :: nr3_start, nr3_end 128 INTEGER :: nbnd_eff 129 130 131 write(stdout,*) 'MATRIX BIG1' 132 FLUSH(stdout) 133 134 iunwfcreal2=find_free_unit() 135 CALL diropn( iunwfcreal2, 'real_whole', dffts%nnr, exst ) 136 137 138 allocate(tmprealis(dffts%nnr,n_set),tmprealjs(dffts%nnr,n_set), tmpreal(dffts%nnr)) 139 allocate(tmpexp2(dffts%nnr,6)) 140 141!set nr3p for not parallel case 142 143#if !defined(__MPI) 144 dfftp%nr3p(1) = dfftp%nr3 145 dffts%nr3p(1) = dffts%nr3 146#endif 147 148 149 150 151 152!set up exponential grid 153 154 tmpexp2(:,:)=(0.d0,0.d0) 155 156#if !defined(__MPI) 157 iqq=0 158 do ix=1,dffts%nr1 159 do iy=1,dffts%nr2 160 do iz=1,dffts%nr3 161 iqq=(iz-1)*(dffts%nr1x*dffts%nr2x)+(iy-1)*dffts%nr1x+ix 162 tmpexp2(iqq,1) = exp(cmplx(0.d0,1.d0)*tpi*real(ix-1)/real(dffts%nr1)) 163 tmpexp2(iqq,2) = exp(cmplx(0.d0,1.d0)*tpi*real(iy-1)/real(dffts%nr2)) 164 tmpexp2(iqq,3) = exp(cmplx(0.d0,1.d0)*tpi*real(iz-1)/real(dffts%nr3)) 165 tmpexp2(iqq,4) = exp(cmplx(0.d0,-1.d0)*tpi*real(ix-1)/real(dffts%nr1)) 166 tmpexp2(iqq,5) = exp(cmplx(0.d0,-1.d0)*tpi*real(iy-1)/real(dffts%nr2)) 167 tmpexp2(iqq,6) = exp(cmplx(0.d0,-1.d0)*tpi*real(iz-1)/real(dffts%nr3)) 168 enddo 169 enddo 170 enddo 171 172 173#else 174 write(stdout,*) 'NRS', dffts%nr1,dffts%nr2,dffts%nr3 175 write(stdout,*) 'NRXS', dffts%nr1x,dffts%nr2x,dffts%nr3x 176 nr3s_start=0 177 nr3s_end =0 178 do ii=1,me_pool + 1 179 nr3s_start=nr3s_end+1 180 nr3s_end=nr3s_end+dffts%nr3p(ii) 181 end do 182 tmpexp2(:,:)=(0.d0,0.d0) 183 do iz=1,dffts%my_nr3p 184 do iy=1,dffts%nr2 185 do ix=1,dffts%nr1 186 iqq=(iz-1)*(dffts%nr1x*dffts%nr2x)+(iy-1)*dffts%nr1+ix 187 tmpexp2(iqq,1) = exp(cmplx(0.d0,1.d0)*tpi*real(ix-1)/real(dffts%nr1)) 188 tmpexp2(iqq,2) = exp(cmplx(0.d0,1.d0)*tpi*real(iy-1)/real(dffts%nr2)) 189 tmpexp2(iqq,3) = exp(cmplx(0.d0,1.d0)*tpi*real(iz+nr3s_start-1-1)/real(dffts%nr3)) 190 tmpexp2(iqq,4) = exp(cmplx(0.d0,-1.d0)*tpi*real(ix-1)/real(dffts%nr1)) 191 tmpexp2(iqq,5) = exp(cmplx(0.d0,-1.d0)*tpi*real(iy-1)/real(dffts%nr2)) 192 tmpexp2(iqq,6) = exp(cmplx(0.d0,-1.d0)*tpi*real(iz+nr3s_start-1-1)/real(dffts%nr3)) 193 enddo 194 enddo 195 enddo 196 197 198#endif 199 200 write(stdout,*) 'Calculate grid' 201 202 203 204 nbnd_eff=num_nbndv(ispin) 205 206 write(stdout,*) 'MATRIX BIG2' 207 FLUSH(stdout) 208 209 do iiw=1,nbnd_eff/n_set+1 210 write(stdout,*) 'MATRIX IIW',iiw 211 FLUSH(stdout) 212 213 do iw=(iiw-1)*n_set+1,min(iiw*n_set,nbnd_eff) 214!read from disk wfc on coarse grid 215 CALL davcio( tmprealis(:,iw-(iiw-1)*n_set),dffts%nnr,iunwfcreal2,iw+(ispin-1)*nbnd,-1) 216 enddo 217!read in iw wfcs 218 do jjw=iiw,nbnd_eff/n_set+1 219 write(stdout,*) 'MATRIX JJW',jjw 220 FLUSH(stdout) 221 222 do jw=(jjw-1)*n_set+1,min(jjw*n_set,nbnd_eff) 223 CALL davcio( tmprealjs(:,jw-(jjw-1)*n_set),dffts%nnr,iunwfcreal2,jw+(ispin-1)*nbnd,-1) 224 enddo 225 !do product 226 227 do iw=(iiw-1)*n_set+1,min(iiw*n_set,nbnd_eff) 228 if(iiw==jjw) then 229 jw_begin=iw 230 else 231 jw_begin=(jjw-1)*n_set+1 232 endif 233 do jw=jw_begin,min(jjw*n_set,nbnd_eff) 234 235 tmpreal(:)=tmprealis(:,iw-(iiw-1)*n_set)*tmprealjs(:,jw-(jjw-1)*n_set) 236 237!put on fine grid 238 239 240!add us part 241 242!calculate matrix element 243 do mdir=1,3 244 sca=0.d0 245 do ir=1,dffts%nnr 246 sca=sca+tmpreal(ir)*tmpexp2(ir,mdir) 247 enddo 248 sca=sca/dble(dffts%nr1*dffts%nr2*dffts%nr3) 249 call mp_sum(sca,world_comm) 250 !call reduce(2,sca) 251 252 matsincos(iw,jw,mdir)=dble(sca) 253 matsincos(jw,iw,mdir)=dble(sca) 254 matsincos(iw,jw,mdir+3)=dimag(sca) 255 matsincos(jw,iw,mdir+3)=dimag(sca) 256 !matp(iw,jw,mdir)=sca 257 !matp(jw,iw,mdir)=sca 258 enddo 259 260 261 enddo 262 enddo 263 enddo 264 enddo 265 266 267 deallocate(tmprealis,tmprealjs) 268 deallocate(tmpexp2) 269 270 write(stdout,*) 'Calculate US' 271 FLUSH(stdout) 272 if(okvan) then 273 allocate(tmpexp(dfftp%nnr)) 274 allocate(expgsave(maxval(nh),maxval(nh),nat,3)) 275 expgsave(:,:,:,:)=0.d0 276 do mdir=1,3 277 278#if !defined(__MPI) 279 if(mdir==1) then 280 do ix=1,dfftp%nr1 281 ee=exp(cmplx(0.d0,1.d0)*tpi*real(ix)/real(dfftp%nr1)) 282 do iy=1,dfftp%nr2 283 do iz=1,dfftp%nr3 284 nn=(iz-1)*dfftp%nr1x*dfftp%nr2+(iy-1)*dfftp%nr1+ix 285 tmpexp(nn)=ee 286 enddo 287 enddo 288 enddo 289 else if(mdir==2) then 290 do iy=1,dfftp%nr2 291 ee=exp(cmplx(0.d0,1.d0)*tpi*real(iy)/real(dfftp%nr2)) 292 do ix=1,dfftp%nr1 293 do iz=1,dfftp%nr3 294 nn=(iz-1)*dfftp%nr1x*dfftp%nr2x+(iy-1)*dfftp%nr1x+ix 295 tmpexp(nn)=ee 296 enddo 297 enddo 298 enddo 299 else if(mdir==3) then 300 do iz=1,dfftp%nr3 301 ee=exp(cmplx(0.d0,1.d0)*tpi*real(iz)/real(dfftp%nr3)) 302 do ix=1,dfftp%nr1 303 do iy=1,dfftp%nr2 304 nn=(iz-1)*dfftp%nr1x*dfftp%nr2x+(iy-1)*dfftp%nr1x+ix 305 tmpexp(nn)=ee 306 enddo 307 enddo 308 enddo 309 endif 310 311#else 312 nr3_start=0 313 nr3_end =0 314 do ii=1,me_pool + 1 315 nr3_start=nr3_end+1 316 nr3_end=nr3_end+dfftp%nr3p(ii) 317 end do 318 319 do iz=1,dfftp%nr3p(me_pool+1) 320 do iy=1,dfftp%nr2 321 do ix=1,dfftp%nr1 322 323 nn=(iz-1)*dfftp%nr1x*dfftp%nr2x+(iy-1)*dfftp%nr1x+ix 324 if(mdir==1) then 325 tmpexp(nn)= exp(cmplx(0.d0,1.d0)*tpi*real(ix-1)/real(dfftp%nr1)) 326 elseif(mdir==2) then 327 tmpexp(nn)= exp(cmplx(0.d0,1.d0)*tpi*real(iy-1)/real(dfftp%nr2)) 328 else 329 tmpexp(nn)= exp(cmplx(0.d0,1.d0)*tpi*real(iz+nr3_start-1-1)/real(dfftp%nr3)) 330 endif 331 enddo 332 enddo 333 enddo 334 335 336#endif 337 338 339 do np = 1, ntyp 340 if ( upf(np)%tvanp ) then 341 do na = 1, nat 342 if ( ityp(na) == np ) then 343 do ih = 1, nh(np) 344 do jh = ih, nh(np) 345 expgsave(ih,jh,na,mdir)=(0.d0,0.d0) 346 !do ir =1,maxbox(na) 347 ! expgsave(ih,jh,na,mdir)=expgsave(ih,jh,na,mdir)+qsave(ih,jh,na)%q(ir)*tmpexp(box(ir,na)) 348 !enddo 349 enddo 350 enddo 351 endif 352 enddo 353 endif 354 enddo 355 356 expgsave(:,:,:,mdir)=expgsave(:,:,:,mdir)*omega/dble(dfftp%nr1*dfftp%nr2*dfftp%nr3) 357 358#if defined(__MPI) 359 ! call reduce (2 *maxval(nh) *maxval(nh)* nat, expgsave(:,:,:,mdir)) 360 call mp_sum( expgsave(:,:,:,mdir),world_comm) 361#endif 362 363 364 do iw=1,nbnd_eff 365 do jw=iw,nbnd_eff 366 367 do is=1, nspin 368 ijkb0 = 0 369 do np = 1, ntyp 370 if ( upf(np)%tvanp ) then 371 do na = 1, nat 372 if ( ityp(na) == np ) then 373 do ih = 1, nh(np) 374 ikb = ijkb0 + ih 375 do jh = 1, nh(np) 376 jkb = ijkb0 + jh 377 if(ih <= jh) then 378 if(itask /= 1) then 379 matsincos(iw,jw,mdir)=matsincos(iw,jw,mdir)+& 380 &dble(expgsave(ih,jh,na,mdir) * becp_gw(ikb,iw,1)*becp_gw(jkb,jw,1)) 381 matsincos(iw,jw,mdir+3)=matsincos(iw,jw,mdir+3)+& 382 &dimag(expgsave(ih,jh,na,mdir) * becp_gw(ikb,iw,1)*becp_gw(jkb,jw,1)) 383 else 384 matsincos(iw,jw,mdir)=matsincos(iw,jw,mdir)+& 385 &dble(expgsave(ih,jh,na,mdir) * becp_gw_c(ikb,iw,1)*becp_gw_c(jkb,jw,1)) 386 matsincos(iw,jw,mdir+3)=matsincos(iw,jw,mdir+3)+& 387 &dimag(expgsave(ih,jh,na,mdir) * becp_gw_c(ikb,iw,1)*becp_gw_c(jkb,jw,1)) 388 endif 389 390 else 391 if(itask /= 1) then 392 !matp(iw,jw,mdir)=matp(iw,jw,mdir)+expgsave(jh,ih,na,mdir) * becp_gw(ikb,iw,1)*becp_gw(jkb,jw,1) 393 matsincos(iw,jw,mdir)=matsincos(iw,jw,mdir)+& 394 &dble(expgsave(jh,ih,na,mdir) * becp_gw(ikb,iw,1)*becp_gw(jkb,jw,1)) 395 matsincos(iw,jw,mdir+3)=matsincos(iw,jw,mdir+3)+& 396 &dimag(expgsave(jh,ih,na,mdir) * becp_gw(ikb,iw,1)*becp_gw(jkb,jw,1)) 397 else 398 !matp(iw,jw,mdir)=matp(iw,jw,mdir)+expgsave(jh,ih,na,mdir) * becp_gw_c(ikb,iw,1)*becp_gw_c(jkb,jw,1) 399 matsincos(iw,jw,mdir)=matsincos(iw,jw,mdir)+& 400 &dble(expgsave(jh,ih,na,mdir) * becp_gw_c(ikb,iw,1)*becp_gw_c(jkb,jw,1)) 401 matsincos(iw,jw,mdir+3)=matsincos(iw,jw,mdir+3)+& 402 &dimag(expgsave(jh,ih,na,mdir) * becp_gw_c(ikb,iw,1)*becp_gw_c(jkb,jw,1)) 403 endif 404 endif 405 enddo 406 enddo 407 ijkb0=ijkb0+nh(np) 408 endif 409 enddo 410 else 411 do na=1,nat 412 if(ityp(na) == np) ijkb0=ijkb0+nh(np) 413 enddo 414 endif 415 enddo 416 enddo 417! matp(jw,iw,mdir)=matp(iw,jw,mdir) 418 matsincos(jw,iw,mdir)=matsincos(iw,jw,mdir) 419 matsincos(jw,iw,mdir+3)=matsincos(iw,jw,mdir+3) 420 enddo 421 enddo 422 enddo 423 424 deallocate(tmpexp) 425 endif 426 427 close(iunwfcreal2) 428 429 return 430 431 end subroutine 432 433 434