1MODULE contract_w 2 3USE kinds, ONLY : DP,sgl 4use bse_basic_structures, ONLY : ii_mat 5 6SAVE 7 8type(ii_mat) :: iimat_contract 9REAL(kind=sgl), POINTER :: vphipizeta_save(:,:) 10REAL(kind=sgl), POINTER :: vww_save(:,:) 11COMPLEX(kind=sgl), POINTER :: vphipizeta_save_g(:,:) 12COMPLEX(kind=sgl), POINTER :: vww_save_g(:,:) 13INTEGER, POINTER :: vpmax_ii(:),vpmax_ii_start(:),vpmax_ii_end(:) 14INTEGER :: vpmax_tot 15 16CONTAINS 17 18subroutine free_memory_contrac_w 19 use bse_wannier, ONLY : l_gtrick 20 implicit none 21 if(.not.l_gtrick) then 22 deallocate(vphipizeta_save) 23 else 24 deallocate(vphipizeta_save_g) 25 endif 26 deallocate(vpmax_ii) 27 deallocate(vpmax_ii_start,vpmax_ii_end) 28 if(.not.l_gtrick) then 29 deallocate(vww_save) 30 else 31 deallocate(vww_save_g) 32 endif 33end subroutine free_memory_contrac_w 34 35subroutine contract_w_build(fc) 36! this subroutine computes the w part of the direct term of the exc Hamiltonian 37 38USE wvfct, ONLY : npw 39USE fft_custom_gwl 40use bse_basic_structures 41use exciton 42USE gvect 43use bse_wannier, ONLY: l_truncated_coulomb, & 44 truncation_radius, l_gtrick 45USE constants, ONLY : e2, fpi 46USE cell_base, ONLY : tpiba,omega,tpiba2 47!USE io_files, ONLY : find_free_unit, prefix, diropn 48USE io_files, ONLY : prefix, diropn 49USE wavefunctions, ONLY : psic 50USE io_global, ONLY : stdout, ionode, ionode_id 51USE mp_world, ONLY : mpime, nproc 52USE mp_pools, ONLY: intra_pool_comm 53USE mp_wave, ONLY : mergewf,splitwf 54USE polarization 55USE lsda_mod, ONLY :nspin 56USE io_global, ONLY : stdout,ionode 57USE mp, ONLY :mp_barrier 58USE mp_world, ONLY : world_comm 59 60 61 62 63 64implicit none 65INTEGER, EXTERNAL :: find_free_unit 66logical :: debug=.false. 67 68type(bse_z) :: z 69type(polaw) :: pw 70 71type(fft_cus) :: fc 72type(ii_mat) :: iimat 73 74 75 76REAL(kind=DP), ALLOCATABLE :: fac(:) 77COMPLEX(kind=DP), ALLOCATABLE :: p_basis(:,:) 78COMPLEX(kind=DP), ALLOCATABLE :: p_basis_t(:,:) 79REAL(kind=DP), ALLOCATABLE :: p_basis_r(:,:) 80REAL(kind=DP), ALLOCATABLE :: zvphi(:) 81REAL(kind=DP), ALLOCATABLE :: zvv(:) 82COMPLEX(kind=DP), ALLOCATABLE :: evc_g(:) 83 84INTEGER ::iungprod 85INTEGER :: ig,ii,iv,ispin 86REAL(kind=DP) :: qq 87LOGICAL :: exst 88 89 90INTEGER :: vpmax,k 91REAL(kind=DP), allocatable :: zp(:,:) 92REAL(kind=DP), allocatable :: pizeta(:,:) 93REAL(kind=DP), allocatable :: vphipizeta(:,:) 94 95INTEGER :: kilobytes 96 97 98call start_clock('direct_w_exc') 99!CALL memstat( kilobytes ) 100!write(stdout,*) 'memory0', kilobytes 101!FLUSH(stdout) 102 103 104 105 106! read iimat 107call initialize_imat(iimat) 108 109do ispin=1,nspin 110 call read_iimat(iimat,ispin) 111enddo 112 113! read z terms 114call initialize_bse_z(z) 115call read_z(1,iimat,z) 116 117FLUSH( stdout ) 118 119! get Coulomb potential 120allocate(fac(npw)) 121if(l_truncated_coulomb) then 122 do ig=1,npw 123 qq = g(1,ig)**2.d0 + g(2,ig)**2.d0 + g(3,ig)**2.d0 124 if (qq > 1.d-8) then 125 fac(ig)=(e2*fpi/(tpiba2*qq))*(1.d0-dcos(dsqrt(qq)*truncation_radius*tpiba)) 126 else 127 fac(ig)=e2*fpi*(truncation_radius**2.d0/2.d0) 128 endif 129 enddo 130 fac(:)=fac(:)/omega 131else 132 133 fac(:)=0.d0 134 fac(1:npw)=vg_q(1:npw) 135endif 136 137 138! read polarization basis and multiply per V 139 140iungprod = find_free_unit() 141allocate(p_basis(npw,z%numw_prod)) 142CALL diropn( iungprod, 'wiwjwfc_red', npw*2, exst ) 143 144do ii=1,z%numw_prod 145 call davcio(p_basis(:,ii),npw*2,iungprod,ii,-1) 146 p_basis(1:npw,ii)=p_basis(1:npw,ii)*dcmplx(fac(1:npw)) 147enddo 148 149call mp_barrier(world_comm) 150 151close(iungprod) 152!CALL memstat( kilobytes ) 153!write(stdout,*) 'memory1', kilobytes 154!FLUSH(stdout) 155 156! FFT to real space (dual grid) 157allocate(p_basis_t(fc%npwt,z%numw_prod)) 158allocate(p_basis_r(fc%nrxxt,z%numw_prod)) 159allocate(evc_g(fc%ngmt_g )) 160 161 162if(fc%dual_t==4.d0) then 163 p_basis_t(1:fc%npwt,1:z%numw_prod)=p_basis(1:npw,1:z%numw_prod) 164else 165 call reorderwfp_col(z%numw_prod,npw,fc%npwt,p_basis(1,1),p_basis_t(1,1),npw,fc%npwt, & 166 & ig_l2g,fc%ig_l2gt,fc%ngmt_g,mpime, nproc,intra_pool_comm ) 167 168! do ii=1,z%numw_prod 169! call mergewf(p_basis(:,ii),evc_g,a_in%npw,ig_l2g,mpime,nproc,ionode_id,intra_pool_comm) 170! call splitwf(p_basis_t(:,ii),evc_g,fc%npwt,fc%ig_l2gt,mpime,nproc,ionode_id,intra_pool_comm) 171! enddo 172endif 173 174deallocate(evc_g) 175deallocate(p_basis) 176 177call start_clock('direct_w_cft3t') 178do ii=1,z%numw_prod,2 179 psic(1:fc%nrxxt)=(0.d0,0.d0) 180 if (ii==z%numw_prod) then 181 psic(fc%nlt(1:fc%npwt)) = p_basis_t(1:fc%npwt,ii) 182 psic(fc%nltm(1:fc%npwt)) = CONJG( p_basis_t(1:fc%npwt,ii) ) 183 else 184 psic(fc%nlt(1:fc%npwt))=p_basis_t(1:fc%npwt,ii)+(0.d0,1.d0)*p_basis_t(1:fc%npwt,ii+1) 185 psic(fc%nltm(1:fc%npwt))=CONJG(p_basis_t(1:fc%npwt,ii))+(0.d0,1.d0)*CONJG(p_basis_t(1:fc%npwt,ii+1)) 186 endif 187 CALL cft3t( fc, psic, fc%nr1t, fc%nr2t, fc%nr3t, fc%nrx1t, fc%nrx2t, fc%nrx3t, 2 ) 188 p_basis_r(1:fc%nrxxt,ii)= DBLE(psic(1:fc%nrxxt)) 189 if(ii/=z%numw_prod) p_basis_r(1:fc%nrxxt,ii+1)= DIMAG(psic(1:fc%nrxxt)) 190enddo 191call stop_clock('direct_w_cft3t') 192 193deallocate(p_basis_t) 194 195!read P 196call initialize_polaw(pw) 197call read_polaw_global(0, pw) 198 199 200call mp_barrier(world_comm) 201 202!CALL memstat( kilobytes ) 203!write(stdout,*) 'memory2', kilobytes 204!FLUSH(stdout) 205 206 207! allocate tmp matrix 208 209 210!compute line by line the output excitonic vector 211 212!!!!!!!!!!!!!!!!!dgemm subroutine!!!!!!!!!!!!!!!!!!!!! 213call start_clock('direct_w_dgemv') 214!write(stdout,*) 'memory2',z%numw_prod,iimat%np_max 215allocate(zp(z%numw_prod,iimat%np_max)) 216!!allocate(pizeta(z%numw_prod,iimat%np_max)) 217!!allocate(vphipizeta(fc%nrxxt,iimat%np_max)) 218!calculate vpmax_tot 219vpmax_tot=0 220allocate(vpmax_ii( iimat%numb_v)) 221allocate(vpmax_ii_start( iimat%numb_v)) 222allocate(vpmax_ii_end( iimat%numb_v)) 223do iv=1, iimat%numb_v 224 225 vpmax=0 226 vpmax_ii_start(iv)=vpmax_tot+1 227 do ii=1, iimat%np_max 228 if (iimat%iimat(ii,iv)==0) cycle 229 vpmax=vpmax+1 230 enddo 231 vpmax_ii(iv)=vpmax 232 vpmax_tot=vpmax_tot+vpmax 233 vpmax_ii_end(iv)=vpmax_tot 234enddo 235write(stdout,*) 'VPHIZETA_SAVE :', fc%nrxxt,vpmax_tot 236FLUSH(stdout) 237if(.not.l_gtrick) then 238 allocate(vphipizeta_save(fc%nrxxt,vpmax_tot)) 239else 240 allocate(vphipizeta_save_g(fc%npwt,vpmax_tot)) 241endif 242write(stdout,*) 'VPHIZETA_SAVE :', fc%nrxxt,vpmax_tot 243! 244!CALL memstat( kilobytes ) 245!write(stdout,*) 'memory3', kilobytes 246!FLUSH(stdout) 247do iv=1, iimat%numb_v 248 zp(1:z%numw_prod,1:iimat%np_max)=0.d0 249 vpmax=0 250 251 call start_clock('dgemv1') 252 do ii=1, iimat%np_max 253 if (iimat%iimat(ii,iv)==0) cycle 254 vpmax=vpmax+1 255 do k=1, z%numw_prod 256 zp(k,vpmax)=z%z(k,ii,iv)!ATTENZIONE era zp(ii) 257 enddo 258 enddo 259 260 if(vpmax>0) then 261 allocate(pizeta(z%numw_prod,vpmax)) 262 allocate(vphipizeta(fc%nrxxt,vpmax)) 263 else 264 allocate(pizeta(z%numw_prod,1)) 265 allocate(vphipizeta(fc%nrxxt,1)) 266 endif 267 call stop_clock('dgemv1') 268 269 270 call start_clock('dgemv2') 271 272 273 call dgemm('N','N', z%numw_prod,vpmax, z%numw_prod,1.d0,pw%pw,z%numw_prod,zp,& 274 z%numw_prod,0.d0,pizeta,z%numw_prod) 275 call stop_clock('dgemv2') 276 277 278 call start_clock('dgemv3') 279 call dgemm('N','N', fc%nrxxt, vpmax, z%numw_prod,1.d0,p_basis_r(1,1),fc%nrxxt,& 280 pizeta(1,1), z%numw_prod, 0.d0, vphipizeta(1,1),fc%nrxxt) 281 call stop_clock('dgemv3') 282 if(.not.l_gtrick) then 283 vphipizeta_save(1:fc%nrxxt,vpmax_ii_start(iv):vpmax_ii_end(iv))= real(vphipizeta(1:fc%nrxxt,1:vpmax)) 284 else 285 do ii=1,vpmax,2 286 if(ii==vpmax) then 287 psic(1:fc%nrxxt)=dcmplx(vphipizeta(1:fc%nrxxt,ii),0.d0) 288 else 289 psic(1:fc%nrxxt)=dcmplx(vphipizeta(1:fc%nrxxt,ii),vphipizeta(1:fc%nrxxt,ii+1)) 290 endif 291 CALL cft3t(fc, psic, fc%nr1t, fc%nr2t, fc%nr3t, fc%nrx1t, fc%nrx2t, fc%nrx3t, -2 ) 292 if(ii==vpmax) then 293 vphipizeta_save_g(1:fc%npwt,vpmax_ii_start(iv)+ii-1)=cmplx(psic(fc%nlt(1:fc%npwt))) 294 else 295 vphipizeta_save_g(1:fc%npwt,vpmax_ii_start(iv)+ii-1)=& 296 &cmplx(0.5d0*(psic(fc%nlt(1:fc%npwt))+conjg( psic(fc%nltm(1:fc%npwt))))) 297 vphipizeta_save_g(1:fc%npwt,vpmax_ii_start(iv)+ii-1+1)=& 298 &cmplx((0.d0,-0.5d0)*(psic(fc%nlt(1:fc%npwt)) - conjg(psic(fc%nltm(1:fc%npwt))))) 299 endif 300 enddo 301 endif 302 303 ! 304! 305 deallocate(pizeta) 306 deallocate(vphipizeta) 307call mp_barrier(world_comm) 308enddo 309call stop_clock('direct_w_dgemv') 310 311deallocate(zp) 312deallocate(p_basis_r) 313 314call free_bse_z(z) 315call free_memory_polaw(pw) 316call free_imat(iimat) 317 318 319!CALL memstat( kilobytes ) 320!write(stdout,*) 'memory4', kilobytes 321!FLUSH(stdout) 322call stop_clock('direct_w_exc') 323 324return 325end subroutine contract_w_build 326 327 328 329subroutine contract_w_apply(a_in,fc,a_out) 330! this subroutine computes the w part of the direct term of the exc Hamiltonian 331 332USE fft_custom_gwl 333use bse_basic_structures 334use exciton 335USE gvect 336use bse_wannier, ONLY: l_truncated_coulomb, & 337 truncation_radius,l_gtrick 338USE constants, ONLY : e2, fpi 339USE cell_base, ONLY : tpiba,omega,tpiba2 340!USE io_files, ONLY : find_free_unit, prefix, diropn 341USE io_files, ONLY : prefix, diropn 342USE wavefunctions, ONLY : psic 343USE io_global, ONLY : stdout, ionode, ionode_id 344USE mp_world, ONLY : mpime, nproc 345USE mp_pools, ONLY: intra_pool_comm 346USE mp_wave, ONLY : mergewf,splitwf 347USE polarization 348USE lsda_mod, ONLY :nspin 349USE io_global, ONLY : stdout,ionode 350USE mp, ONLY :mp_barrier 351USE mp_world, ONLY : world_comm 352 353 354 355 356 357implicit none 358INTEGER, EXTERNAL :: find_free_unit 359 360type(polaw) :: pw 361type(exc):: a_in 362type(exc):: a_out 363type(exc_r):: a_in_rt 364type(exc_r):: a_tmp_rt 365type(fft_cus) :: fc 366 367 368 369 370 371INTEGER ::iungprod 372INTEGER :: ig,ii,iv,ispin 373REAL(kind=DP) :: qq 374LOGICAL :: exst 375 376 377INTEGER ::k 378 379 380logical debug 381 382call start_clock('direct_w_contract') 383 384debug=.false. 385! read iimat 386 387!FFT the input excitonic vector to real space (dual grid) 388call initialize_exc_r(a_in_rt) 389call fft_a_exc(a_in,fc,a_in_rt) 390 391call mp_barrier(world_comm) 392 393! allocate tmp matrix 394call initialize_exc_r(a_tmp_rt) 395a_tmp_rt%nrxxt=fc%nrxxt 396a_tmp_rt%numb_v=a_in%numb_v 397a_tmp_rt%label=12 398allocate(a_tmp_rt%ar(a_tmp_rt%nrxxt,a_tmp_rt%numb_v)) 399 400!compute line by line the output excitonic vector 401 402!!!!!!!!!!!!!!!!!dgemm subroutine!!!!!!!!!!!!!!!!!!!!! 403call start_clock('contract_w_dgemv') 404 405 406! 407a_tmp_rt%ar(1:a_tmp_rt%nrxxt,1:a_tmp_rt%numb_v) =0.d0 408! 409do iv=1, a_in%numb_v 410 411 if(.not.l_gtrick) then 412 do ii=1,vpmax_ii(iv) 413 a_tmp_rt%ar(1:fc%nrxxt,iv)= & 414 &a_tmp_rt%ar(1:fc%nrxxt,iv)+a_in_rt%ar(1:fc%nrxxt,iimat_contract%iimat(ii,iv))*& 415 &dble(vphipizeta_save(1:fc%nrxxt,vpmax_ii_start(iv)-1+ii)) 416 enddo 417 else 418 do ii=1,vpmax_ii(iv),2 419 psic(1:fc%nrxxt)=(0.d0,0.d0) 420 if (ii==vpmax_ii(iv)) then 421 psic(fc%nlt(1:fc%npwt)) = dcmplx(vphipizeta_save_g(1:fc%npwt,vpmax_ii_start(iv)-1+ii)) 422 psic(fc%nltm(1:fc%npwt)) = dcmplx(CONJG( vphipizeta_save_g(1:fc%npwt,vpmax_ii_start(iv)-1+ii) )) 423 else 424 psic(fc%nlt(1:fc%npwt))= dcmplx(vphipizeta_save_g(1:fc%npwt,vpmax_ii_start(iv)-1+ii))+& 425 &(0.0,1.0)* dcmplx(vphipizeta_save_g(1:fc%npwt,vpmax_ii_start(iv)-1+ii+1)) 426 psic(fc%nltm(1:fc%npwt))=DCMPLX(CONJG( vphipizeta_save_g(1:fc%npwt,vpmax_ii_start(iv)-1+ii))+& 427 &(0.0,1.0)*CONJG( vphipizeta_save_g(1:fc%npwt,vpmax_ii_start(iv)-1+ii+1))) 428 endif 429 CALL cft3t( fc, psic, fc%nr1t, fc%nr2t, fc%nr3t, fc%nrx1t, fc%nrx2t, fc%nrx3t, 2 ) 430 a_tmp_rt%ar(1:fc%nrxxt,iv)= & 431 &a_tmp_rt%ar(1:fc%nrxxt,iv)+a_in_rt%ar(1:fc%nrxxt,iimat_contract%iimat(ii,iv))*& 432 &dble(psic(1:fc%nrxxt)) 433 if (ii/=vpmax_ii(iv)) then 434 a_tmp_rt%ar(1:fc%nrxxt,iv)= & 435 &a_tmp_rt%ar(1:fc%nrxxt,iv)+a_in_rt%ar(1:fc%nrxxt,iimat_contract%iimat(ii+1,iv))*& 436 &dimag(psic(1:fc%nrxxt)) 437 438 endif 439 enddo 440 endif 441 442 443 444enddo 445call stop_clock('contract_w_dgemv') 446 447 448 449 450call free_memory_exc_a_r(a_in_rt) 451 452 453 454call start_clock('wdirect_fftback') 455!FFT back to provide the output excitonic wave vector in G-space 456call fftback_a_exc(a_tmp_rt,fc,a_out) 457call stop_clock('wdirect_fftback') 458 459call free_memory_exc_a_r(a_tmp_rt) 460 461FLUSH( stdout ) 462call stop_clock('direct_w_contract') 463 464return 465end subroutine contract_w_apply 466 467 468subroutine contract_v_build(fc) 469!TO BE CALLED AFTER CONTRACT_W_BUILD 470!IMPORTANT: REQUIRES IIMAT_CONTRACT 471! computes the v part of the direct term of the exc Hamiltonian 472USE fft_custom_gwl 473use bse_basic_structures 474use exciton 475USE wavefunctions, ONLY : psic 476USE gvect, ONLY : ig_l2g 477USE io_global, ONLY : stdout, ionode, ionode_id 478USE mp_world, ONLY : mpime, nproc 479USE mp_pools, ONLY: intra_pool_comm 480USE mp_wave, ONLY : mergewf,splitwf 481USE io_global, ONLY : stdout,ionode 482USE lsda_mod, ONLY :nspin 483USE gvect, ONLY : gstart 484USE mp, ONLY : mp_sum 485USE mp_world, ONLY : world_comm 486USE wvfct, ONLY : npw 487USE bse_wannier, ONLY : l_gtrick 488 489implicit none 490 491 492type(vww_prod) :: vww 493type(fft_cus) :: fc 494 495 496COMPLEX(kind=DP), allocatable :: vwwg_t(:,:) 497COMPLEX(kind=DP), ALLOCATABLE :: evc_g(:) 498 499 500!real(kind=dp), allocatable :: phivwwr(:,:) 501COMPLEX(kind=DP) :: csca 502 503integer ii, iv,ispin, iimax 504 505logical debug 506 507write(stdout,*) 'Routine contract_v_build' 508 509debug=.false. 510 511!if(debug) then 512! if(ionode) write(stdout,*) 'Direct_v_exc #1' 513!endif 514write(stdout,*) 'VWW_SAVE :', fc%nrxxt,vpmax_tot 515if(.not.l_gtrick) then 516 allocate(vww_save(fc%nrxxt,vpmax_tot)) 517else 518 allocate(vww_save_g(fc%npwt,vpmax_tot)) 519endif 520 521 522call initialize_vww_prod(vww) 523call read_vww_prod(1,iimat_contract%numb_v,npw,iimat_contract%np_max,iimat_contract,vww) 524 525!if(debug) then 526! if(ionode) write(stdout,*) 'Direct_v_exc #6' 527!endif 528 529! for every element iv of the excitonic wavefunction vector, here we FFT all 530! the available v*w_iv*w_ivp(G) products, multiply by a_in_rt%ar(:,ivp) 531! sum over ivp, and FFT back 532 533allocate(vwwg_t(fc%npwt,iimat_contract%np_max)) 534!allocate(phivwwr(fc%nrxxt,a_in%numb_v)) 535allocate(evc_g(fc%ngmt_g )) 536 537 538write(stdout,*) 'ATT1' 539do iv=1,iimat_contract%numb_v 540 541 vwwg_t(1:fc%npwt,1:iimat_contract%np_max)=dcmplx(0.d0,0.d0) 542 iimax=0 543 do ii=1,iimat_contract%np_max 544 if (iimat_contract%iimat(ii,iv)>0) then 545 iimax=iimax+1 546 else 547 exit 548 endif 549 enddo 550 if(iimax>0) then 551 call reorderwfp_col(iimax,vww%npw,fc%npwt,vww%vww(1,1,iv),vwwg_t, vww%npw,fc%npwt, & 552 & ig_l2g,fc%ig_l2gt,fc%ngmt_g,mpime, nproc,intra_pool_comm ) 553 endif 554 555 if(.not.l_gtrick) then 556 do ii=1,vpmax_ii(iv),2 557 558 559 psic(1:fc%nrxxt)=(0.d0,0.d0) 560 if (ii==vpmax_ii(iv)) then 561 psic(fc%nlt(1:fc%npwt)) = vwwg_t(1:fc%npwt,ii) 562 psic(fc%nltm(1:fc%npwt)) = CONJG( vwwg_t(1:fc%npwt,ii) ) 563 else 564 psic(fc%nlt(1:fc%npwt))=vwwg_t(1:fc%npwt,ii)+(0.d0,1.d0)*vwwg_t(1:fc%npwt,ii+1) 565 psic(fc%nltm(1:fc%npwt))=CONJG(vwwg_t(1:fc%npwt,ii))+(0.d0,1.d0)*CONJG(vwwg_t(1:fc%npwt,ii+1)) 566 endif 567 CALL cft3t( fc, psic, fc%nr1t, fc%nr2t, fc%nr3t, fc%nrx1t, fc%nrx2t, fc%nrx3t, 2 ) 568 vww_save(1:fc%nrxxt,vpmax_ii_start(iv)+ii-1)=real(psic(1:fc%nrxxt)) 569 if (ii/=vpmax_ii(iv)) then 570 vww_save(1:fc%nrxxt,vpmax_ii_start(iv)+ii)=aimag(psic(1:fc%nrxxt)) 571 endif 572 573 enddo 574 else 575 do ii=1,vpmax_ii(iv) 576 vww_save_g(1:fc%npwt,vpmax_ii_start(iv)+ii-1)=cmplx(vwwg_t(1:fc%npwt,ii)) 577 end do 578 endif 579enddo 580 581!if(debug) then 582! if(ionode) write(stdout,*) 'Direct_v_exc #7' 583!endif 584 585FLUSH( stdout ) 586 587 588call free_vww_prod(vww) 589deallocate(vwwg_t) 590!deallocate(vwwr_t) 591deallocate(evc_g) 592 593 594end subroutine 595 596 597subroutine contract_v_apply(a_in,fc,a_out) 598! computes the v part of the direct term of the exc Hamiltonian 599USE fft_custom_gwl 600use bse_basic_structures 601use exciton 602USE wavefunctions, ONLY : psic 603USE gvect, ONLY : ig_l2g 604USE io_global, ONLY : stdout, ionode, ionode_id 605USE mp_world, ONLY : mpime, nproc 606USE mp_pools, ONLY: intra_pool_comm 607USE mp_wave, ONLY : mergewf,splitwf 608USE io_global, ONLY : stdout,ionode 609USE lsda_mod, ONLY :nspin 610USE gvect, ONLY : gstart 611USE mp, ONLY : mp_sum 612USE mp_world, ONLY : world_comm 613USE bse_wannier, ONLY : l_gtrick 614 615 616implicit none 617type(exc), intent(in) :: a_in 618type(exc):: a_out 619type(exc_r):: a_in_rt 620type(exc_r):: a_tmp_rt 621 622 623type(fft_cus) :: fc 624 625 626COMPLEX(kind=DP), allocatable :: vwwg_t(:,:) 627COMPLEX(kind=DP), ALLOCATABLE :: evc_g(:) 628 629COMPLEX(kind=DP) :: csca 630 631integer ii, iv,ispin, iimax 632 633logical debug 634 635call start_clock('direct_v_contract') 636debug=.false. 637 638call initialize_exc_r(a_tmp_rt) 639a_tmp_rt%nrxxt=fc%nrxxt 640a_tmp_rt%numb_v=a_in%numb_v 641a_tmp_rt%label=12 642allocate(a_tmp_rt%ar(a_tmp_rt%nrxxt,a_tmp_rt%numb_v)) 643 644 645! FFT a_in to real space (dual grid) 646call initialize_exc_r(a_in_rt) 647call fft_a_exc(a_in,fc,a_in_rt) 648 649 650 651 652 653a_tmp_rt%ar(1:a_tmp_rt%nrxxt,1:a_tmp_rt%numb_v) =0.d0 654do iv=1,a_in%numb_v 655 656 657 if(.not.l_gtrick) then 658 do ii=1,vpmax_ii(iv) 659 a_tmp_rt%ar(1:fc%nrxxt,iv)= & 660 &a_tmp_rt%ar(1:fc%nrxxt,iv)+DBLE(vww_save(1:fc%nrxxt,vpmax_ii_start(iv)+ii-1))*& 661 &a_in_rt%ar(1:a_in_rt%nrxxt,iimat_contract%iimat(ii,iv)) 662 enddo 663 else 664 do ii=1,vpmax_ii(iv),2 665 psic(1:fc%nrxxt)=(0.d0,0.d0) 666 if (ii==vpmax_ii(iv)) then 667 psic(fc%nlt(1:fc%npwt)) = dcmplx(vww_save_g(1:fc%npwt,vpmax_ii_start(iv)+ii-1)) 668 psic(fc%nltm(1:fc%npwt)) = dcmplx(CONJG( vww_save_g(1:fc%npwt,vpmax_ii_start(iv)+ii-1) )) 669 else 670 psic(fc%nlt(1:fc%npwt))=dcmplx(vww_save_g(1:fc%npwt,vpmax_ii_start(iv)+ii-1)+& 671 &(0.0,1.0)*vww_save_g(1:fc%npwt,vpmax_ii_start(iv)+ii-1+1)) 672 psic(fc%nltm(1:fc%npwt))=dcmplx(CONJG(vww_save_g(1:fc%npwt,vpmax_ii_start(iv)+ii-1))& 673 &+(0.0,1.0)*CONJG(vww_save_g(1:fc%npwt,vpmax_ii_start(iv)+ii-1+1))) 674 endif 675 call start_clock('d_v_fft') 676 CALL cft3t( fc, psic, fc%nr1t, fc%nr2t, fc%nr3t, fc%nrx1t, fc%nrx2t, fc%nrx3t, 2 ) 677 call stop_clock('d_v_fft') 678 a_tmp_rt%ar(1:fc%nrxxt,iv)= & 679 &a_tmp_rt%ar(1:fc%nrxxt,iv)+DBLE(psic(1:fc%nrxxt))*& 680 &a_in_rt%ar(1:a_in_rt%nrxxt,iimat_contract%iimat(ii,iv)) 681 if (ii/=vpmax_ii(iv)) then 682 a_tmp_rt%ar(1:fc%nrxxt,iv)= & 683 &a_tmp_rt%ar(1:fc%nrxxt,iv)+DIMAG(psic(1:fc%nrxxt))*& 684 &a_in_rt%ar(1:a_in_rt%nrxxt,iimat_contract%iimat(ii+1,iv)) 685 endif 686 enddo 687 endif 688enddo 689 690 691 692call free_memory_exc_a_r(a_in_rt) 693 694call fftback_a_exc(a_tmp_rt,fc,a_out) 695 696! free memory 697call free_memory_exc_a_r(a_tmp_rt) 698 699 700 701call stop_clock('direct_v_contract') 702end subroutine 703 704 705 706 707 708 709 710END MODULE contract_w 711