Skip to content

GitLab

Commit 3693bc1c authored by Matteo Gerosa's avatar Matteo Gerosa
Browse files

gamma_only and k cases treated in same subroutine calc_exx2

parent 9d3e7f84
Show whitespace changes
Inline Side-by-side
Wfreq/calc_exx2.f90
View file @ 3693bc1c
......@@ -11,197 +11,7 @@
! Marco Govoni
!
!-----------------------------------------------------------------------
SUBROUTINE calc_exx2_gamma( sigma_exx, nb1, nb2 )
!-----------------------------------------------------------------------
!
! store in sigma_exx(n,iks) = < n,iks | V_exx | n,iks > n = nb1, nb2
!
USE kinds, ONLY : DP
USE mp_global, ONLY : inter_image_comm,intra_bgrp_comm,my_image_id
USE mp, ONLY : mp_sum,mp_bcast
USE io_global, ONLY : stdout, ionode
USE io_files, ONLY : nwordwfc, iunwfc
USE scf, ONLY : rho, rho_core, rhog_core
USE gvect, ONLY : g,nl,gstart,ngm_g,ngm
USE gvecw, ONLY : gcutw
USE cell_base, ONLY : tpiba2,omega,tpiba,at,alat
USE fft_base, ONLY : dfftp,dffts
USE fft_interfaces, ONLY : fwfft, invfft
USE constants, ONLY : tpi,fpi,rytoev,e2
USE pwcom, ONLY : npw,npwx,et,nks,current_spin,isk,xk,nbnd,lsda,igk_k,g2kin,nspin,current_k,ngk
USE fft_at_gamma, ONLY : single_invfft_gamma,single_fwfft_gamma
USE fft_at_k, ONLY : single_invfft_k,single_fwfft_k
USE wavefunctions_module, ONLY : evc,psic,psic_nc
USE westcom, ONLY : iuwfc,lrwfc,npwq,nbnd_occ
USE control_flags, ONLY : gamma_only
USE noncollin_module, ONLY : noncolin,npol
USE buffers, ONLY : get_buffer
USE funct, ONLY : dft_is_hybrid,init_dft_exxrpa,stop_exx
USE exx, ONLY : x_gamma_extrapolation,exxdiv_treatment,exx_grid_init,exx_div_check,&
&deallocate_exx,exxinit,vexx,exx_grid_initialized
USE uspp, ONLY : vkb,nkb
USE bar, ONLY : bar_type,start_bar_type,update_bar_type,stop_bar_type
USE io_push, ONLY : io_push_bar
USE class_idistribute, ONLY : idistribute
USE coulomb_vcut_module, ONLY : vcut_init, vcut_type, vcut_info, &
vcut_get, vcut_spheric_get, vcut_destroy
USE types_coulomb, ONLY : pot3D
!
IMPLICIT NONE
!
! I/O
!
INTEGER, INTENT(IN) :: nb1, nb2
REAL(DP),INTENT(OUT) :: sigma_exx( nb1:nb2, nks)
!
! Workspace
!
COMPLEX(DP),ALLOCATABLE :: pertg(:),pertr(:),pertr_nc(:,:)
REAL(DP), EXTERNAL :: DDOT
INTEGER :: ib,iv,i1,ir,iks,ig,iv_glob
INTEGER :: nbndval
TYPE(idistribute) :: vband
TYPE(bar_type) :: barra
INTEGER :: barra_load
!
WRITE(stdout,'(5x,a)') ' '
CALL io_push_bar()
WRITE(stdout,'(5x,a)') '(X)-Sigma'
CALL io_push_bar()
!
ALLOCATE( pertg( ngm ) )
! IF(noncolin) THEN
! ALLOCATE( pertr_nc( dffts%nnr, npol ) )
! ELSE
ALLOCATE( pertr( dffts%nnr ) )
! ENDIF
!
CALL pot3D%init('Dense',.FALSE.,'gb')
!
! Set to zero
!
sigma_exx = 0._DP
!
barra_load = nks * ( nb2-nb1 + 1 )
CALL start_bar_type( barra, 'sigmax', barra_load )
!
! LOOP
!
DO iks = 1, nks ! KPOINT-SPIN
!
! ... Set k-point, spin, kinetic energy, needed by Hpsi
!
current_k = iks
IF ( lsda ) current_spin = isk(iks)
call g2_kin( iks )
!
! ... More stuff needed by the hamiltonian: nonlocal projectors
!
IF ( nkb > 0 ) CALL init_us_2( ngk(iks), igk_k(1,iks), xk(1,iks), vkb )
!npw = ngk(iks)
!
! ... read in wavefunctions from the previous iteration
!
IF(nks>1) THEN
!iuwfc = 20
!lrwfc = nbnd * npwx * npol
!!CALL get_buffer( evc, nwordwfc, iunwfc, iks )
IF(my_image_id==0) CALL get_buffer( evc, lrwfc, iuwfc, iks )
!CALL mp_bcast(evc,0,inter_image_comm)
!CALL davcio(evc,lrwfc,iuwfc,iks,-1)
CALL mp_bcast(evc,0,inter_image_comm)
ENDIF
! !
! ! ... Needed for LDA+U
! !
! IF ( nks > 1 .AND. lda_plus_u .AND. (U_projection .NE. 'pseudo') ) &
! CALL get_buffer ( wfcU, nwordwfcU, iunhub, iks )
! !
! current_k = iks
! current_spin = isk(iks)
! !
! CALL gk_sort(xk(1,iks),ngm,g,gcutw,npw,igk,g2kin)
! g2kin=g2kin*tpiba2
! !
! ! reads unperturbed wavefuctions psi_k in G_space, for all bands
! !
! !
! CALL init_us_2 (npw, igk, xk (1, iks), vkb)
!
nbndval = nbnd_occ(iks)
!
vband = idistribute()
CALL vband%init(nbndval,'i','nbndval',.FALSE.)
!
DO ib = nb1, nb2
!
!IF(gamma_only) THEN
CALL single_invfft_gamma(dffts,npw,npwx,evc(1,ib),psic,'Wave')
!ELSEIF(noncolin) THEN
! CALL single_invfft_k(dffts,npw,npwx,evc(1 ,ib),psic_nc(1,1),'Wave',igk_k(1,current_k))
! CALL single_invfft_k(dffts,npw,npwx,evc(1+npwx,ib),psic_nc(1,2),'Wave',igk_k(1,current_k))
! ELSE
! CALL single_invfft_k(dffts,npw,npwx,evc(1,ib),psic,'Wave',igk_k(1,current_k))
!ENDIF
!
DO iv = 1, vband%nloc
!
iv_glob = vband%l2g(iv)
!
! Bring it to R-space
!IF(gamma_only) THEN
CALL single_invfft_gamma(dffts,npw,npwx,evc(1,iv_glob),pertr,'Wave')
DO ir=1,dffts%nnr
pertr(ir)=psic(ir)*pertr(ir)
ENDDO
CALL single_fwfft_gamma(dffts,ngm,ngm,pertr,pertg,'Dense')
!ELSEIF(noncolin) THEN
! CALL single_invfft_k(dffts,npw,npwx,evc(1 ,iv_glob),pertr_nc(1,1),'Wave',igk_k(1,current_k))
! CALL single_invfft_k(dffts,npw,npwx,evc(1+npwx,iv_glob),pertr_nc(1,2),'Wave',igk_k(1,current_k))
! DO ir=1,dffts%nnr
! pertr_nc(ir,1)=DCONJG(psic_nc(ir,1))*pertr_nc(ir,1)+DCONJG(psic_nc(ir,2))*pertr_nc(ir,2)
! ENDDO
! CALL single_fwfft_k(dffts,ngm,ngm,pertr_nc(1,1),pertg,'Dense') ! no igk
!ELSE
! CALL single_invfft_k(dffts,npw,npwx,evc(1,iv_glob),pertr,'Wave',igk_k(1,current_k))
! DO ir=1,dffts%nnr
! pertr(ir)=DCONJG(psic(ir))*pertr(ir)
! ENDDO
! CALL single_fwfft_k(dffts,ngms,ngms,pertr,pertg,'Smooth') ! no igk
!ENDIF
!
DO ig = 1,ngm
pertg(ig) = pertg(ig) * pot3D%sqvc(ig)
ENDDO
sigma_exx( ib, iks ) = sigma_exx( ib, iks ) - 2._DP * DDOT( 2*ngm, pertg(1), 1, pertg(1), 1) / omega
!IF(gstart==2) sigma_exx( ib, iks ) = sigma_exx( ib, iks ) + REAL( pertg(1), KIND = DP )**2 / omega
IF( ib == iv_glob .AND. gstart == 2 ) sigma_exx( ib, iks ) = sigma_exx( ib, iks ) - pot3D%div
!
ENDDO
!
CALL update_bar_type( barra, 'sigmax', 1 )
!
ENDDO ! ib
!
ENDDO ! iks
!
CALL stop_bar_type( barra, 'sigmax' )
!
CALL mp_sum( sigma_exx, intra_bgrp_comm )
CALL mp_sum( sigma_exx, inter_image_comm )
!
DEALLOCATE( pertg )
! DEALLOCATE( mysqvc )
! IF( noncolin ) THEN
! DEALLOCATE( pertr_nc )
! ELSE
DEALLOCATE( pertr )
! ENDIF
!
END SUBROUTINE
!
!-----------------------------------------------------------------------
SUBROUTINE calc_exx2_k( sigma_exx, nb1, nb2 )
SUBROUTINE calc_exx2( sigma_exx, nb1, nb2 )
!-----------------------------------------------------------------------
!
! store in sigma_exx(n,iks) = < n,iks | V_exx | n,iks > n = nb1, nb2
......@@ -258,7 +68,7 @@ SUBROUTINE calc_exx2_k( sigma_exx, nb1, nb2 )
TYPE(bar_type) :: barra
INTEGER :: barra_load
LOGICAL :: l_gammaq
REAL(DP) :: g0(3)
REAL(DP) :: g0(3), peso
!
WRITE(stdout,'(5x,a)') ' '
CALL io_push_bar()
......@@ -266,8 +76,13 @@ SUBROUTINE calc_exx2_k( sigma_exx, nb1, nb2 )
CALL io_push_bar()
!
ALLOCATE( pertg( ngm ) )
IF (gamma_only) THEN
peso = 2._DP
ELSE
peso = 1._DP
ALLOCATE( phase(dffts%nnr) )
ALLOCATE( evckmq(npwx*npol,nbnd) )
ENDIF
IF(noncolin) THEN
ALLOCATE( pertr_nc( dffts%nnr, npol ) )
ELSE
......@@ -328,16 +143,14 @@ SUBROUTINE calc_exx2_k( sigma_exx, nb1, nb2 )
! !
! CALL init_us_2 (npw, igk, xk (1, iks), vkb)
!
nbndval = nbnd_occ(iks)
!
vband = idistribute()
CALL vband%init(nbndval,'i','nbndval',.FALSE.)
!
DO ib = nb1, nb2
!
sigma_exx(ib,iks) = 0._DP
!
IF (noncolin) THEN
IF (gamma_only) THEN
CALL single_invfft_gamma(dffts,npw,npwx,evc(1,ib),psic,'Wave')
ELSEIF (noncolin) THEN
CALL single_invfft_k(dffts,npw,npwx,evc(1 ,ib),psic_nc(1,1),'Wave',igk_k(1,current_k))
CALL single_invfft_k(dffts,npw,npwx,evc(npwx+1,ib),psic_nc(1,2),'Wave',igk_k(1,current_k))
ELSE
......@@ -346,25 +159,38 @@ SUBROUTINE calc_exx2_k( sigma_exx, nb1, nb2 )
!
DO iq = 1, q_grid%np
!
IF (gamma_only) THEN
l_gammaq = .TRUE.
CALL pot3D%init('Dense',.FALSE.,'gb')
nbndval = nbnd_occ(iks)
ELSE
l_gammaq = q_grid%l_pIsGamma(iq)
!
CALL k_grid%find( k_grid%p_cart(:,ik) - q_grid%p_cart(:,iq), is, 'cart', ikqs, g0 )
!
CALL pot3D%init('Dense',.FALSE.,'gb',iq)
!
CALL k_grid%find( k_grid%p_cart(:,ik) - q_grid%p_cart(:,iq), is, 'cart', ikqs, g0 )
CALL compute_phase( g0, 'cart', phase )
!
nbndval = nbnd_occ(ikqs)
npwkq = ngk(ikqs)
!
IF ( my_image_id == 0 ) CALL get_buffer( evckmq, lrwfc, iuwfc, ikqs )
CALL mp_bcast( evckmq, 0, inter_image_comm )
ENDIF
!
vband = idistribute()
CALL vband%init(nbndval,'i','nbndval',.FALSE.)
!
DO iv = 1, vband%nloc
!
iv_glob = vband%l2g(iv)
!
! Bring it to R-space
IF(noncolin) THEN
IF (gamma_only) THEN
CALL single_invfft_gamma(dffts,npw,npwx,evc(1,iv_glob),pertr,'Wave')
DO ir=1,dffts%nnr
pertr(ir)=psic(ir)*pertr(ir)
ENDDO
CALL single_fwfft_gamma(dffts,ngm,ngm,pertr,pertg,'Dense')
ELSEIF(noncolin) THEN
CALL single_invfft_k(dffts,npwkq,npwx,evckmq(1 ,iv_glob),pertr_nc(1,1),'Wave',igk_k(1,ikqs))
CALL single_invfft_k(dffts,npwkq,npwx,evckmq(1+npwx,iv_glob),pertr_nc(1,2),'Wave',igk_k(1,ikqs))
DO ir=1,dffts%nnr
......@@ -382,7 +208,7 @@ SUBROUTINE calc_exx2_k( sigma_exx, nb1, nb2 )
DO ig = 1,ngm
pertg(ig) = pertg(ig) * pot3D%sqvc(ig)
ENDDO
sigma_exx( ib, iks ) = sigma_exx( ib, iks ) - DDOT( 2*ngm, pertg(1), 1, pertg(1), 1)/omega*q_grid%weight(iq)
sigma_exx( ib, iks ) = sigma_exx( ib, iks ) - peso*DDOT( 2*ngm, pertg(1), 1, pertg(1), 1)/omega*q_grid%weight(iq)
!IF(gstart==2) sigma_exx( ib, iks ) = sigma_exx( ib, iks ) + REAL( pertg(1), KIND = DP )**2 / omega
IF( ib == iv_glob .AND. gstart == 2 .AND. l_gammaq ) sigma_exx( ib, iks ) = sigma_exx( ib, iks ) - pot3D%div
!
......@@ -407,5 +233,9 @@ SUBROUTINE calc_exx2_k( sigma_exx, nb1, nb2 )
ELSE
DEALLOCATE( pertr )
ENDIF
IF (.NOT.gamma_only) THEN
DEALLOCATE( phase )
DEALLOCATE( evckmq )
ENDIF
!
END SUBROUTINE
Wfreq/solve_hf.f90
View file @ 3693bc1c
......@@ -74,7 +74,7 @@ SUBROUTINE solve_hf_gamma( )
!
! Get SIGMA EXX
!
CALL calc_exx2_gamma(sigma_exx, qp_bandrange(1), qp_bandrange(2))
CALL calc_exx2(sigma_exx, qp_bandrange(1), qp_bandrange(2))
!
! Get SIGMA X
!
......@@ -109,7 +109,7 @@ SUBROUTINE solve_hf_gamma( )
nbndval = MIN( MAXVAL( nbnd_occ(:) ), nbnd )
ALLOCATE(sigma_exx_all_occupied(nbndval,nks))
!
CALL calc_exx2_gamma( sigma_exx_all_occupied, 1, nbndval )
CALL calc_exx2( sigma_exx_all_occupied, 1, nbndval )
!
exx_etot = 0._DP
DO iks = 1, nks
......@@ -180,7 +180,7 @@ SUBROUTINE solve_hf_k( )
!
! Get SIGMA EXX
!
CALL calc_exx2_k(sigma_exx, qp_bandrange(1), qp_bandrange(2))
CALL calc_exx2(sigma_exx, qp_bandrange(1), qp_bandrange(2))
!
! Get SIGMA X
!
......@@ -215,7 +215,7 @@ SUBROUTINE solve_hf_k( )
nbndval = MIN( MAXVAL( nbnd_occ(:) ), nbnd )
ALLOCATE(sigma_exx_all_occupied(nbndval,k_grid%nps))
!
CALL calc_exx2_k( sigma_exx_all_occupied, 1, nbndval )
CALL calc_exx2( sigma_exx_all_occupied, 1, nbndval )
!
exx_etot = 0._DP
DO iks = 1, k_grid%nps
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment