mirror of https://gitlab.com/QEF/q-e.git
optimized usage of Hybrid functionals in TDDFPT codes
reorganized do loops in k1d_term_gamma functions in h_psi vexx is skipped for TDDFPT using Gamma k2d_term_gamma function became k2d_vexx_term_gamma calculates k2d and vexx terms on the same cycle of the loop
This commit is contained in:
Oscar Baseggio
parent
3f8097e34e
commit
1b19c2af05
|
|
@ -28,21 +28,17 @@ SUBROUTINE h_psi( lda, n, m, psi, hpsi )
|
|||
USE noncollin_module, ONLY : npol
|
||||
USE funct, ONLY : exx_is_active
|
||||
USE mp_bands, ONLY : use_bgrp_in_hpsi, inter_bgrp_comm
|
||||
USE mp, ONLY : mp_allgather, mp_size, &
|
||||
mp_type_create_column_section, mp_type_free
|
||||
USE mp, ONLY : mp_sum
|
||||
!
|
||||
IMPLICIT NONE
|
||||
!
|
||||
INTEGER, INTENT(IN) :: lda, n, m
|
||||
COMPLEX(DP), INTENT(IN) :: psi(lda*npol,m)
|
||||
COMPLEX(DP), INTENT(OUT) :: hpsi(lda*npol,m)
|
||||
COMPLEX(DP), INTENT(OUT) :: hpsi(lda*npol,m)
|
||||
!
|
||||
INTEGER :: m_start, m_end
|
||||
INTEGER :: column_type
|
||||
INTEGER, ALLOCATABLE :: recv_counts(:), displs(:)
|
||||
INTEGER :: m_start, m_end
|
||||
!
|
||||
CALL start_clock( 'h_psi_bgrp' ); !write (*,*) 'start h_psi_bgrp'; FLUSH(6)
|
||||
|
||||
! band parallelization with non-distributed bands is performed if
|
||||
! 1. enabled (variable use_bgrp_in_hpsi must be set to .T.)
|
||||
! 2. exact exchange is not active (if it is, band parallelization is already
|
||||
|
|
@ -51,18 +47,12 @@ SUBROUTINE h_psi( lda, n, m, psi, hpsi )
|
|||
!
|
||||
IF (use_bgrp_in_hpsi .AND. .NOT. exx_is_active() .AND. m > 1) THEN
|
||||
! use band parallelization here
|
||||
ALLOCATE( recv_counts(mp_size(inter_bgrp_comm)), displs(mp_size(inter_bgrp_comm)) )
|
||||
CALL divide_all(inter_bgrp_comm,m,m_start,m_end,recv_counts,displs)
|
||||
CALL mp_type_create_column_section(hpsi(1,1), 0, lda*npol, lda*npol, column_type)
|
||||
!
|
||||
hpsi(:,:) = (0.d0,0.d0)
|
||||
CALL divide(inter_bgrp_comm,m,m_start,m_end)
|
||||
! Check if there at least one band in this band group
|
||||
IF (m_end >= m_start) &
|
||||
CALL h_psi_( lda, n, m_end-m_start+1, psi(1,m_start), hpsi(1,m_start) )
|
||||
CALL mp_allgather(hpsi, column_type, recv_counts, displs, inter_bgrp_comm)
|
||||
!
|
||||
CALL mp_type_free( column_type )
|
||||
DEALLOCATE( recv_counts )
|
||||
DEALLOCATE( displs )
|
||||
CALL mp_sum(hpsi,inter_bgrp_comm)
|
||||
ELSE
|
||||
! don't use band parallelization here
|
||||
CALL h_psi_( lda, n, m, psi, hpsi )
|
||||
|
|
@ -99,7 +89,7 @@ SUBROUTINE h_psi_( lda, n, m, psi, hpsi )
|
|||
USE ldaU, ONLY : lda_plus_u, U_projection
|
||||
USE gvect, ONLY : gstart
|
||||
USE funct, ONLY : dft_is_meta
|
||||
USE control_flags, ONLY : gamma_only
|
||||
USE control_flags, ONLY : gamma_only, tddfpt
|
||||
USE noncollin_module, ONLY: npol, noncolin
|
||||
USE realus, ONLY : real_space, &
|
||||
invfft_orbital_gamma, fwfft_orbital_gamma, calbec_rs_gamma, add_vuspsir_gamma, &
|
||||
|
|
@ -116,28 +106,17 @@ SUBROUTINE h_psi_( lda, n, m, psi, hpsi )
|
|||
COMPLEX(DP), INTENT(IN) :: psi(lda*npol,m)
|
||||
COMPLEX(DP), INTENT(OUT) :: hpsi(lda*npol,m)
|
||||
!
|
||||
INTEGER :: ipol, ibnd
|
||||
INTEGER :: ipol, ibnd, incr
|
||||
REAL(dp) :: ee
|
||||
!
|
||||
CALL start_clock( 'h_psi' ); !write (*,*) 'start h_psi';FLUSH(6)
|
||||
!
|
||||
! ... Here we add the kinetic energy (k+G)^2 psi and clean up garbage
|
||||
!
|
||||
!$omp parallel do
|
||||
DO ibnd = 1, m
|
||||
hpsi (1:n, ibnd) = g2kin (1:n) * psi (1:n, ibnd)
|
||||
IF (n<lda) hpsi (n+1:lda, ibnd) = (0.0_dp, 0.0_dp)
|
||||
IF ( noncolin ) THEN
|
||||
hpsi (lda+1:lda+n, ibnd) = g2kin (1:n) * psi (lda+1:lda+n, ibnd)
|
||||
IF (n<lda) hpsi (lda+n+1:lda+lda, ibnd) = (0.0_dp, 0.0_dp)
|
||||
END IF
|
||||
END DO
|
||||
!$omp end parallel do
|
||||
|
||||
hpsi (:, 1:m) = (0.0_dp, 0.0_dp)
|
||||
|
||||
CALL start_clock( 'h_psi:pot' ); !write (*,*) 'start h_pot';FLUSH(6)
|
||||
!
|
||||
! ... Here the product with the local potential V_loc psi
|
||||
!
|
||||
CALL start_clock( 'h_psi:pot' ); !write (*,*) 'start h_pot';FLUSH(6)
|
||||
!
|
||||
IF ( gamma_only ) THEN
|
||||
!
|
||||
IF ( real_space .and. nkb > 0 ) then
|
||||
|
|
@ -145,10 +124,12 @@ SUBROUTINE h_psi_( lda, n, m, psi, hpsi )
|
|||
! ... real-space algorithm
|
||||
! ... fixme: real_space without beta functions does not make sense
|
||||
!
|
||||
IF ( dffts%has_task_groups ) &
|
||||
CALL errore( 'h_psi', 'task_groups not implemented with real_space', 1 )
|
||||
|
||||
DO ibnd = 1, m, 2
|
||||
IF ( dffts%has_task_groups ) then
|
||||
incr = 2 * fftx_ntgrp(dffts)
|
||||
ELSE
|
||||
incr = 2
|
||||
ENDIF
|
||||
DO ibnd = 1, m, incr
|
||||
! ... transform psi to real space -> psic
|
||||
CALL invfft_orbital_gamma(psi,ibnd,m)
|
||||
! ... compute becp%r = < beta|psi> from psic in real space
|
||||
|
|
@ -180,9 +161,11 @@ SUBROUTINE h_psi_( lda, n, m, psi, hpsi )
|
|||
! ... real-space algorithm
|
||||
! ... fixme: real_space without beta functions does not make sense
|
||||
!
|
||||
IF ( dffts%has_task_groups ) &
|
||||
CALL errore( 'h_psi', 'task_groups not implemented with real_space', 1 )
|
||||
|
||||
IF ( dffts%has_task_groups ) then
|
||||
incr = fftx_ntgrp(dffts)
|
||||
ELSE
|
||||
incr = 1
|
||||
ENDIF
|
||||
DO ibnd = 1, m
|
||||
! ... transform psi to real space -> psic
|
||||
CALL invfft_orbital_k(psi,ibnd,m)
|
||||
|
|
@ -220,6 +203,15 @@ SUBROUTINE h_psi_( lda, n, m, psi, hpsi )
|
|||
!
|
||||
CALL stop_clock( 'h_psi:pot' )
|
||||
!
|
||||
! ... Here we add the kinetic energy (k+G)^2 psi
|
||||
!
|
||||
DO ibnd = 1, m
|
||||
hpsi (1:n, ibnd) = hpsi(1:n, ibnd) + g2kin (1:n) * psi (1:n, ibnd)
|
||||
IF ( noncolin ) THEN
|
||||
hpsi (lda+1:lda+n, ibnd) = hpsi(lda+1:lda+n,ibnd) + g2kin (1:n) * psi (lda+1:lda+n, ibnd)
|
||||
END IF
|
||||
END DO
|
||||
!
|
||||
if (dft_is_meta()) call h_psi_meta (lda, n, m, psi, hpsi)
|
||||
!
|
||||
! ... Here we add the Hubbard potential times psi
|
||||
|
|
@ -236,7 +228,7 @@ SUBROUTINE h_psi_( lda, n, m, psi, hpsi )
|
|||
!
|
||||
! ... Here the exact-exchange term Vxx psi
|
||||
!
|
||||
IF ( exx_is_active() ) THEN
|
||||
IF ( exx_is_active() .and..not. (tddfpt .and. gamma_only)) THEN
|
||||
IF ( use_ace) THEN
|
||||
IF (gamma_only) THEN
|
||||
CALL vexxace_gamma(lda,m,psi,ee,hpsi)
|
||||
|
|
|
|||
|
|
@ -36,7 +36,7 @@ MODULE lr_exx_kernel
|
|||
USE lr_variables, ONLY : gamma_only, lr_verbosity
|
||||
USE realus, ONLY : invfft_orbital_gamma, fwfft_orbital_gamma,&
|
||||
& invfft_orbital_k, fwfft_orbital_k
|
||||
USE wavefunctions, ONLY : psic
|
||||
USE wavefunctions_module, ONLY : psic
|
||||
USE cell_base, ONLY : omega
|
||||
USE exx_base, ONLY : g2_convolution
|
||||
USE exx, ONLY : exxalfa, npwt, gt, dfftt
|
||||
|
|
@ -391,6 +391,10 @@ SUBROUTINE lr_exx_kernel_noint ( evc, int_vect )
|
|||
COMPLEX(DP), ALLOCATABLE :: psic_all(:), temppsic_all(:), temppsic(:)
|
||||
!
|
||||
INTEGER :: ibnd_start_gamma, ibnd_end_gamma
|
||||
|
||||
LOGICAL :: exst
|
||||
|
||||
|
||||
! Offset ibnd_start and ibnd_ed to avoid conflicts with gamma_tricks for
|
||||
! even values
|
||||
!
|
||||
|
|
@ -451,9 +455,9 @@ SUBROUTINE lr_exx_kernel_noint ( evc, int_vect )
|
|||
w2=0.0d0
|
||||
ENDIF
|
||||
! Update the container with the actual interaction for this band(s).
|
||||
revc_int(1:dfftt%nnr,:)= revc_int(1:dfftt%nnr,:) &
|
||||
& -1.d0 * scale * k1d_term_gamma(w1,w2,psic,fac,ibnd) &
|
||||
& +1.d0 * scale * k2d_term_gamma(w1,w2,psic,fac,ibnd)
|
||||
revc_int(1:dfftt%nnr,:)= revc_int(1:dfftt%nnr,:) &
|
||||
& -1.d0 * scale * k1d_term_gamma(w1,w2,psic,fac,ibnd,evc(:,:,1)) &
|
||||
& +1.d0 * scale * k2d_vexx_term_gamma(w1,w2,psic,fac,ibnd,.false.)
|
||||
!
|
||||
ENDDO
|
||||
!
|
||||
|
|
@ -684,8 +688,8 @@ SUBROUTINE lr_exx_kernel_int ( orbital, ibnd, nbnd, ikk )
|
|||
!
|
||||
IF (.NOT.ltammd) THEN
|
||||
revc_int(1:dfftt%nnr,:)= revc_int(1:dfftt%nnr,:)&
|
||||
& -1.d0 * scale * k1d_term_gamma(w1,w2,psic,fac,ibnd) &
|
||||
& -1.d0 * scale * k2d_term_gamma(w1,w2,psic,fac,ibnd)
|
||||
& -1.d0 * scale * k1d_term_gamma(w1,w2,psic,fac,ibnd,orbital) &
|
||||
& -1.d0 * scale * k2d_vexx_term_gamma(w1,w2,psic,fac,ibnd,.true.)
|
||||
ELSE
|
||||
!
|
||||
! Slightly different interaction in the Tamm--Dancoff case.
|
||||
|
|
@ -694,7 +698,7 @@ SUBROUTINE lr_exx_kernel_int ( orbital, ibnd, nbnd, ikk )
|
|||
! case.
|
||||
!
|
||||
revc_int(1:dfftt%nnr,:)= revc_int(1:dfftt%nnr,:)&
|
||||
& -2.d0 * scale * k1d_term_gamma(w1,w2,psic,fac,ibnd)
|
||||
& -2.d0 * scale * k1d_term_gamma(w1,w2,psic,fac,ibnd,orbital)
|
||||
ENDIF
|
||||
!
|
||||
ELSE
|
||||
|
|
@ -786,94 +790,6 @@ SUBROUTINE lr_exx_kernel_int ( orbital, ibnd, nbnd, ikk )
|
|||
|
||||
END SUBROUTINE lr_exx_kernel_int
|
||||
!
|
||||
FUNCTION k1d_term_gamma(w1, w2, psi, fac_in, ibnd) RESULT (psi_int)
|
||||
!------------------------------------------------------------------
|
||||
!
|
||||
! This routine computes the K^1d term, Eq (21) from Eq Ref (1).
|
||||
! As the integral in this term remains the same throughout the
|
||||
! chain it can also be calculated once for each combination of
|
||||
! bands and stored in k1d_pot().
|
||||
!
|
||||
! psi contains two bands of |a> with w1, w2 the associated weights
|
||||
! fac_in contains the interaction W(G) and ibnd the band index n'.
|
||||
!
|
||||
! (1) Rocca, Lu and Galli, J. Chem. Phys., 133, 164109 (2010)
|
||||
!------------------------------------------------------------------
|
||||
!
|
||||
|
||||
IMPLICIT NONE
|
||||
|
||||
COMPLEX(KIND=DP),DIMENSION(:), INTENT(IN) :: psi
|
||||
REAL(KIND=DP),DIMENSION(:), INTENT(IN) :: fac_in
|
||||
REAL(kind=DP), INTENT(IN) :: w1, w2
|
||||
REAL(KIND=DP), ALLOCATABLE :: psi_int(:,:)
|
||||
INTEGER, INTENT(IN) :: ibnd
|
||||
!
|
||||
! Workspaces
|
||||
!
|
||||
INTEGER :: ibnd2, is, npw_, ngm_, nnr_
|
||||
!
|
||||
npw_=npwt
|
||||
ngm_=dfftt%ngm
|
||||
nnr_=dfftt%nnr
|
||||
!
|
||||
ALLOCATE(psi_int(nnr_, nbnd))
|
||||
psi_int = 0.d0
|
||||
!
|
||||
!
|
||||
DO ibnd2=1,nbnd,1
|
||||
!
|
||||
!
|
||||
! Set up the pseudo density for this Hartree like interaction.
|
||||
!
|
||||
vhart(:,:) = (0.d0,0.d0)
|
||||
pseudo_dens_c(:) = (0.d0,0.d0)
|
||||
!
|
||||
! The following code is to check if the value of ibnd2 is even or odd
|
||||
! and therefore whether the REAL or IMAGINARY part of red_revc0 is to be
|
||||
! used. This is because red_revc0 is stored using gamma_tricks.
|
||||
!
|
||||
IF (MOD(ibnd2,2)==1) THEN
|
||||
pseudo_dens_c(1:nnr_) = CMPLX( w1*DBLE(red_revc0(1:nnr_,ibnd,1)) *&
|
||||
& DBLE(red_revc0(1:nnr_,ibnd2,1)), &
|
||||
& w2*AIMAG(red_revc0(1:nnr_,ibnd, 1)) *&
|
||||
& DBLE(red_revc0(1:nnr_,ibnd2,1)), kind=DP )
|
||||
ELSE
|
||||
pseudo_dens_c(1:nnr_) = CMPLX( w1*DBLE(red_revc0(1:nnr_,ibnd,1)) *&
|
||||
&AIMAG(red_revc0(1:nnr_,ibnd2-1,1)),&
|
||||
& w2*AIMAG(red_revc0(1:nnr_,ibnd,1)) *&
|
||||
&AIMAG(red_revc0(1:nnr_,ibnd2-1,1)), kind=DP )
|
||||
ENDIF
|
||||
!
|
||||
CALL fwfft ('Rho', pseudo_dens_c, dfftt)
|
||||
!
|
||||
! hartree contribution is computed in reciprocal space
|
||||
!
|
||||
DO is = 1, nspin
|
||||
!
|
||||
vhart(dfftt%nl(1:ngm_),is) =&
|
||||
& pseudo_dens_c(dfftt%nl(1:ngm_)) *&
|
||||
& fac_in(1:ngm_)
|
||||
IF (gamma_only) vhart(dfftt%nlm(1:ngm_),is) = &
|
||||
& pseudo_dens_c(dfftt%nlm(1:ngm_)) *&
|
||||
& fac_in(1:ngm_)
|
||||
!
|
||||
! and transformed back to real space
|
||||
!
|
||||
CALL invfft ('Rho', vhart (:, is), dfftt)
|
||||
!
|
||||
ENDDO
|
||||
!
|
||||
! Finally return the interaction
|
||||
!
|
||||
psi_int(1:nnr_,ibnd2) = psi_int(1:nnr_,ibnd2) &
|
||||
& +DBLE(vhart(1:nnr_, 1)) * DBLE(psi(1:nnr_)) &
|
||||
& +AIMAG(vhart(1:nnr_,1)) * AIMAG(psi(1:nnr_))
|
||||
ENDDO
|
||||
!
|
||||
!
|
||||
END FUNCTION k1d_term_gamma
|
||||
!
|
||||
FUNCTION k1d_term_k(w1, psi, fac_in, ibnd, ik,ikq) RESULT (psi_int)
|
||||
!------------------------------------------------------------------
|
||||
!
|
||||
|
|
@ -938,87 +854,6 @@ FUNCTION k1d_term_k(w1, psi, fac_in, ibnd, ik,ikq) RESULT (psi_int)
|
|||
!
|
||||
END FUNCTION k1d_term_k
|
||||
!
|
||||
FUNCTION k2d_term_gamma(w1, w2, psi, fac_in, ibnd) RESULT (psi_int)
|
||||
!------------------------------------------------------------------
|
||||
!
|
||||
! This routine computes the K^2d term, Eq (22) from Eq Ref (1).
|
||||
! psi contains two bands of |b> with w1, w2 the associated weights
|
||||
! fac_in contains the interaction W(G) and ibnd the band index n'.
|
||||
!
|
||||
! (1) Rocca, Lu and Galli, J. Chem. Phys., 133, 164109 (2010)
|
||||
!------------------------------------------------------------------
|
||||
!
|
||||
|
||||
IMPLICIT NONE
|
||||
|
||||
COMPLEX(KIND=DP),DIMENSION(:), INTENT(IN) :: psi
|
||||
REAL(KIND=DP),DIMENSION(:), INTENT(IN) :: fac_in
|
||||
REAL(kind=DP), INTENT(IN) :: w1, w2
|
||||
INTEGER, INTENT(IN) :: ibnd
|
||||
REAL(KIND=DP), ALLOCATABLE :: psi_int(:,:)
|
||||
!
|
||||
! Workspaces
|
||||
!
|
||||
INTEGER :: ibnd2, is, npw_,ngm_, nnr_
|
||||
!
|
||||
nnr_ = dfftt%nnr
|
||||
ngm_ = dfftt%ngm
|
||||
npw_ = npwt
|
||||
!
|
||||
ALLOCATE(psi_int(nnr_, nbnd))
|
||||
psi_int = 0.d0
|
||||
!
|
||||
!
|
||||
DO ibnd2=1,nbnd,1
|
||||
!
|
||||
! Set up the pseudo density for this Hartree like interaction.
|
||||
!
|
||||
vhart(:,:) = (0.d0,0.d0)
|
||||
pseudo_dens_c(:) = (0.d0,0.d0)
|
||||
!
|
||||
! The following code is to check if the value of ibnd2 is even or odd
|
||||
! and therefore whether the REAL or IMAGINARY part of red_revc0 is to be
|
||||
! used. This is because red_revc0 is stored using gamma_tricks.
|
||||
!
|
||||
IF (MOD(ibnd2,2)==1) THEN
|
||||
pseudo_dens_c(1:nnr_) = CMPLX( &
|
||||
& w1*DBLE(psi(1:nnr_))*DBLE(red_revc0(1:nnr_,ibnd2,1)),&
|
||||
& w2*AIMAG(psi(1:nnr_))*DBLE(red_revc0(1:nnr_,ibnd2,1)), kind=DP )
|
||||
ELSE
|
||||
pseudo_dens_c(1:nnr_) = CMPLX( &
|
||||
& w1*DBLE(psi(1:nnr_))*AIMAG(red_revc0(1:nnr_,ibnd2-1,1)),&
|
||||
& w2*AIMAG(psi(1:nnr_))*AIMAG(red_revc0(1:nnr_,ibnd2-1,1)), kind=DP )
|
||||
ENDIF
|
||||
!
|
||||
CALL fwfft ('Rho', pseudo_dens_c, dfftt)
|
||||
!
|
||||
! hartree contribution is computed in reciprocal space
|
||||
!
|
||||
DO is = 1, nspin
|
||||
!
|
||||
vhart(dfftt%nl(1:ngm_),is) = pseudo_dens_c(dfftt%nl(1:ngm_)) * fac_in(1:ngm_)
|
||||
IF (gamma_only) vhart(dfftt%nlm(1:ngm_),is) = &
|
||||
& pseudo_dens_c(dfftt%nlm(1:ngm_)) *&
|
||||
& fac_in(1:ngm_)
|
||||
!
|
||||
! and transformed back to real space
|
||||
!
|
||||
CALL invfft ('Rho', vhart (:, is), dfftt)
|
||||
!
|
||||
ENDDO
|
||||
!
|
||||
!
|
||||
! Finally return the interaction
|
||||
!
|
||||
psi_int(1:nnr_,ibnd2) = psi_int(1:nnr_,ibnd2) &
|
||||
& +DBLE(vhart(1:nnr_,1)) * DBLE(red_revc0(1:nnr_,ibnd,1)) &
|
||||
& +AIMAG(vhart(1:nnr_,1)) * AIMAG(red_revc0(1:nnr_,ibnd,1))
|
||||
!
|
||||
ENDDO
|
||||
!
|
||||
!
|
||||
END FUNCTION k2d_term_gamma
|
||||
!
|
||||
FUNCTION k2d_term_k(w1, psi, fac_in, ibnd, ik, ikq) RESULT (psi_int)
|
||||
!------------------------------------------------------------------
|
||||
!
|
||||
|
|
@ -1158,4 +993,465 @@ SUBROUTINE fwfft_orbital_custom_gamma(orbital, ibnd, nbnd, npwt, dfftt)
|
|||
!
|
||||
END SUBROUTINE fwfft_orbital_custom_gamma
|
||||
|
||||
SUBROUTINE invfft_orbital_ibnd2_gamma(orbital, psitemp, ibnd2, npwt, dfftt)
|
||||
|
||||
USE kinds, ONLY : DP
|
||||
USE fft_types, ONLY : fft_type_descriptor
|
||||
|
||||
IMPLICIT NONE
|
||||
|
||||
COMPLEX(DP), INTENT(IN) :: orbital(:,:)
|
||||
INTEGER, INTENT(IN) :: ibnd2, npwt
|
||||
TYPE(fft_type_descriptor), INTENT(IN) :: dfftt
|
||||
COMPLEX(DP), INTENT(OUT) :: psitemp(:)
|
||||
!
|
||||
psitemp=(0.0_dp, 0.0_dp)
|
||||
!
|
||||
psitemp(dfftt%nl(1:npwt)) = orbital(1:npwt,ibnd2)
|
||||
psitemp(dfftt%nlm(1:npwt)) = CONJG(orbital(1:npwt,ibnd2))
|
||||
!
|
||||
CALL invfft ('Wave', psitemp, dfftt)
|
||||
!
|
||||
RETURN
|
||||
!
|
||||
END SUBROUTINE invfft_orbital_ibnd2_gamma
|
||||
|
||||
FUNCTION k1d_term_gamma(w1, w2, psi, fac_in, ibnd, orbital) RESULT (psi_int)
|
||||
!------------------------------------------------------------------
|
||||
!
|
||||
! This routine computes the K^1d term, Eq (21) from Eq Ref (1).
|
||||
! As the integral in this term remains the same throughout the
|
||||
! chain it can also be calculated once for each combination of
|
||||
! bands and stored in k1d_pot().
|
||||
!
|
||||
! psi contains two bands of |a> with w1, w2 the associated weights
|
||||
! fac_in contains the interaction W(G) and ibnd the band index n'.
|
||||
!
|
||||
! vhart is symmetric integral so calculate only nbnd(nbnd+1)/2
|
||||
! terms (Baseggio 2018)
|
||||
!
|
||||
! (1) Rocca, Lu and Galli, J. Chem. Phys., 133, 164109 (2010)
|
||||
!------------------------------------------------------------------
|
||||
!
|
||||
USE io_global, ONLY : stdout
|
||||
|
||||
|
||||
IMPLICIT NONE
|
||||
|
||||
COMPLEX(KIND=DP),DIMENSION(:), INTENT(IN) :: psi
|
||||
REAL(KIND=DP),DIMENSION(:), INTENT(IN) :: fac_in
|
||||
REAL(kind=DP), INTENT(IN) :: w1, w2
|
||||
REAL(KIND=DP), ALLOCATABLE :: psi_int(:,:)
|
||||
INTEGER, INTENT(IN) :: ibnd
|
||||
COMPLEX(DP), INTENT(IN) :: orbital(:,:)
|
||||
COMPLEX(DP) :: psitemp(dfftt%nnr)
|
||||
REAL(KIND=DP) :: w3
|
||||
!
|
||||
! Workspaces
|
||||
!
|
||||
INTEGER :: ibnd2, is, npw_, ngm_, nnr_
|
||||
INTEGER :: nrec
|
||||
INTEGER :: idxvhartr, idxvharti
|
||||
!
|
||||
npw_=npwt
|
||||
ngm_=dfftt%ngm
|
||||
nnr_=dfftt%nnr
|
||||
!
|
||||
ALLOCATE(psi_int(nnr_, nbnd))
|
||||
psi_int = 0.d0
|
||||
!
|
||||
!
|
||||
IF (ibnd < nbnd) then
|
||||
!
|
||||
vhart(:,:) = (0.d0,0.d0)
|
||||
pseudo_dens_c(:) = (0.d0,0.d0)
|
||||
!
|
||||
! calculate vhart for couples ibnd,ibnd and ibnd+1,ibnd
|
||||
!
|
||||
pseudo_dens_c(1:nnr_) = CMPLX( w1*DBLE(red_revc0(1:nnr_,ibnd,1)) *&
|
||||
& DBLE(red_revc0(1:nnr_,ibnd,1)), &
|
||||
& w2*AIMAG(red_revc0(1:nnr_,ibnd, 1)) *&
|
||||
& DBLE(red_revc0(1:nnr_,ibnd,1)), kind=DP )
|
||||
!
|
||||
CALL fwfft ('Rho', pseudo_dens_c, dfftt)
|
||||
!
|
||||
DO is = 1, nspin
|
||||
!
|
||||
vhart(dfftt%nl(1:ngm_),is) =&
|
||||
& pseudo_dens_c(dfftt%nl(1:ngm_)) *&
|
||||
& fac_in(1:ngm_)
|
||||
IF (gamma_only) vhart(dfftt%nlm(1:ngm_),is) = &
|
||||
& pseudo_dens_c(dfftt%nlm(1:ngm_)) *&
|
||||
& fac_in(1:ngm_)
|
||||
!
|
||||
! and transformed back to real space
|
||||
!
|
||||
CALL invfft ('Rho', vhart (:, is), dfftt)
|
||||
ENDDO
|
||||
!
|
||||
! Finally return the interaction psi_int for terms:
|
||||
! ibnd, ibnd, ibnd
|
||||
! ibnd+1, ibnd+1, ibnd
|
||||
! ibnd, ibnd, ibnd+1
|
||||
!
|
||||
psi_int(1:nnr_,ibnd) = psi_int(1:nnr_,ibnd) &
|
||||
& + DBLE(vhart(1:nnr_, 1)) * DBLE(psi(1:nnr_)) &
|
||||
& + AIMAG(vhart(1:nnr_,1)) * AIMAG(psi(1:nnr_))
|
||||
! & + w1*DBLE(vhart(1:nnr_, 1)) * DBLE(psi(1:nnr_)) &
|
||||
! & + w2*AIMAG(vhart(1:nnr_,1)) * AIMAG(psi(1:nnr_))
|
||||
!
|
||||
psi_int(1:nnr_,ibnd+1) = psi_int(1:nnr_,ibnd+1) &
|
||||
& + AIMAG(vhart(1:nnr_,1)) * DBLE(psi(1:nnr_))
|
||||
! & + w1* AIMAG(vhart(1:nnr_,1)) * DBLE(psi(1:nnr_))
|
||||
!
|
||||
! calculate vhart for couple ibnd+1,ibnd+1
|
||||
!
|
||||
vhart(:,:) = (0.d0,0.d0)
|
||||
pseudo_dens_c(:) = (0.d0,0.d0)
|
||||
!
|
||||
pseudo_dens_c(1:nnr_) = CMPLX( w2*AIMAG(red_revc0(1:nnr_,ibnd,1)) *&
|
||||
& AIMAG(red_revc0(1:nnr_,ibnd,1)), 0.0d0, kind=DP )
|
||||
!
|
||||
CALL fwfft ('Rho', pseudo_dens_c, dfftt)
|
||||
!
|
||||
DO is = 1, nspin
|
||||
!
|
||||
vhart(dfftt%nl(1:ngm_),is) =&
|
||||
& pseudo_dens_c(dfftt%nl(1:ngm_)) *&
|
||||
& fac_in(1:ngm_)
|
||||
IF (gamma_only) vhart(dfftt%nlm(1:ngm_),is) = &
|
||||
& pseudo_dens_c(dfftt%nlm(1:ngm_)) *&
|
||||
& fac_in(1:ngm_)
|
||||
!
|
||||
! and transformed back to real space
|
||||
!
|
||||
CALL invfft ('Rho', vhart (:, is), dfftt)
|
||||
ENDDO
|
||||
!
|
||||
! Finally return the interaction psi_int for term:
|
||||
! ibnd+1, ibnd+1, ibnd+1
|
||||
!
|
||||
psi_int(1:nnr_,ibnd+1) = psi_int(1:nnr_,ibnd+1) &
|
||||
& + DBLE(vhart(1:nnr_,1)) * AIMAG(psi(1:nnr_))
|
||||
! & + w2* DBLE(vhart(1:nnr_,1)) * AIMAG(psi(1:nnr_))
|
||||
!
|
||||
!
|
||||
! start second loop over bands
|
||||
!
|
||||
DO ibnd2=1,ibnd-1,1
|
||||
!
|
||||
! calculate vhart for couples ibnd,ibnd2 and ibnd+1,ibnd2
|
||||
!
|
||||
vhart(:,:) = (0.d0,0.d0)
|
||||
pseudo_dens_c(:) = (0.d0,0.d0)
|
||||
!
|
||||
! The following code is to check if the value of ibnd2 is even or odd
|
||||
! and therefore whether the REAL or IMAGINARY part of red_revc0 is to be
|
||||
! used. This is because red_revc0 is stored using gamma_tricks.
|
||||
!
|
||||
IF (MOD(ibnd2,2)==1) THEN
|
||||
pseudo_dens_c(1:nnr_) = CMPLX( w1*DBLE(red_revc0(1:nnr_,ibnd,1)) *&
|
||||
& DBLE(red_revc0(1:nnr_,ibnd2,1)), &
|
||||
& w2*AIMAG(red_revc0(1:nnr_,ibnd, 1)) *&
|
||||
& DBLE(red_revc0(1:nnr_,ibnd2,1)), kind=DP )
|
||||
ELSE
|
||||
pseudo_dens_c(1:nnr_) = CMPLX( w1*DBLE(red_revc0(1:nnr_,ibnd,1)) *&
|
||||
&AIMAG(red_revc0(1:nnr_,ibnd2-1,1)),&
|
||||
&w2*AIMAG(red_revc0(1:nnr_,ibnd,1)) *&
|
||||
&AIMAG(red_revc0(1:nnr_,ibnd2-1,1)), kind=DP )
|
||||
ENDIF
|
||||
!
|
||||
CALL fwfft ('Rho', pseudo_dens_c, dfftt)
|
||||
!
|
||||
! hartree contribution is computed in reciprocal space
|
||||
!
|
||||
DO is = 1, nspin
|
||||
!
|
||||
vhart(dfftt%nl(1:ngm_),is) =&
|
||||
& pseudo_dens_c(dfftt%nl(1:ngm_)) *&
|
||||
& fac_in(1:ngm_)
|
||||
IF (gamma_only) vhart(dfftt%nlm(1:ngm_),is) = &
|
||||
& pseudo_dens_c(dfftt%nlm(1:ngm_)) *&
|
||||
& fac_in(1:ngm_)
|
||||
!
|
||||
! and transformed back to real space
|
||||
!
|
||||
CALL invfft ('Rho', vhart (:, is), dfftt)
|
||||
!
|
||||
ENDDO
|
||||
!
|
||||
! Finally return the interaction psi_int for terms:
|
||||
! ibnd, ibnd, ibnd2
|
||||
! ibnd+1, ibnd+1, ibnd2
|
||||
! ibnd2, ibnd2, ibnd
|
||||
! ibnd2, ibnd2, ibnd+1
|
||||
!
|
||||
psi_int(1:nnr_,ibnd2) = psi_int(1:nnr_,ibnd2) &
|
||||
& + DBLE(vhart(1:nnr_, 1)) * DBLE(psi(1:nnr_)) &
|
||||
& + AIMAG(vhart(1:nnr_,1)) * AIMAG(psi(1:nnr_))
|
||||
! & + w1*DBLE(vhart(1:nnr_, 1)) * DBLE(psi(1:nnr_)) &
|
||||
! & + w2*AIMAG(vhart(1:nnr_,1)) * AIMAG(psi(1:nnr_))
|
||||
!
|
||||
CALL invfft_orbital_ibnd2_gamma(orbital(:,:), psitemp, ibnd2, npw_, dfftt)
|
||||
!
|
||||
! w3 = wg(ibnd2,1)/omega
|
||||
!
|
||||
psi_int(1:nnr_,ibnd) = psi_int(1:nnr_,ibnd) &
|
||||
& + DBLE(vhart(1:nnr_, 1)) * DBLE(psitemp(1:nnr_))
|
||||
! & + w3*DBLE(vhart(1:nnr_, 1)) * DBLE(psitemp(1:nnr_))
|
||||
!
|
||||
psi_int(1:nnr_,ibnd+1) = psi_int(1:nnr_,ibnd+1) &
|
||||
& + AIMAG(vhart(1:nnr_, 1)) * DBLE(psitemp(1:nnr_))
|
||||
! & + w3*AIMAG(vhart(1:nnr_, 1)) * DBLE(psitemp(1:nnr_))
|
||||
!
|
||||
ENDDO
|
||||
!
|
||||
ELSE
|
||||
!
|
||||
! calculate vhart for couple ibnd,ibnd
|
||||
!
|
||||
vhart(:,:) = (0.d0,0.d0)
|
||||
pseudo_dens_c(:) = (0.d0,0.d0)
|
||||
!
|
||||
pseudo_dens_c(1:nnr_) = CMPLX( w1*DBLE(red_revc0(1:nnr_,ibnd,1)) *&
|
||||
& DBLE(red_revc0(1:nnr_,ibnd,1)), 0.0d0, kind=DP )
|
||||
!
|
||||
CALL fwfft ('Rho', pseudo_dens_c, dfftt)
|
||||
!
|
||||
DO is = 1, nspin
|
||||
!
|
||||
vhart(dfftt%nl(1:ngm_),is) =&
|
||||
& pseudo_dens_c(dfftt%nl(1:ngm_)) *&
|
||||
& fac_in(1:ngm_)
|
||||
IF (gamma_only) vhart(dfftt%nlm(1:ngm_),is) = &
|
||||
& pseudo_dens_c(dfftt%nlm(1:ngm_)) *&
|
||||
& fac_in(1:ngm_)
|
||||
!
|
||||
! and transformed back to real space
|
||||
!
|
||||
CALL invfft ('Rho', vhart (:, is), dfftt)
|
||||
ENDDO
|
||||
!
|
||||
! Finally return the interaction psi_int for term:
|
||||
! ibnd, ibnd, ibnd
|
||||
!
|
||||
psi_int(1:nnr_,ibnd) = psi_int(1:nnr_,ibnd) &
|
||||
& + DBLE(vhart(1:nnr_,1)) * DBLE(psi(1:nnr_))
|
||||
! & + w1* DBLE(vhart(1:nnr_,1)) * DBLE(psi(1:nnr_))
|
||||
!
|
||||
! start second loop over bands
|
||||
!
|
||||
DO ibnd2=1,ibnd-1,1
|
||||
!
|
||||
!
|
||||
! Set up the pseudo density for this Hartree like interaction.
|
||||
! calculate vhart for couple ibnd,ibnd2
|
||||
!
|
||||
vhart(:,:) = (0.d0,0.d0)
|
||||
pseudo_dens_c(:) = (0.d0,0.d0)
|
||||
!
|
||||
! The following code is to check if the value of ibnd2 is even or odd
|
||||
! and therefore whether the REAL or IMAGINARY part of red_revc0 is to be
|
||||
! used. This is because red_revc0 is stored using gamma_tricks.
|
||||
!
|
||||
IF (MOD(ibnd2,2)==1) THEN
|
||||
pseudo_dens_c(1:nnr_) = CMPLX( w1*DBLE(red_revc0(1:nnr_,ibnd,1)) *&
|
||||
& DBLE(red_revc0(1:nnr_,ibnd2,1)), 0.0d0, kind=DP )
|
||||
ELSE
|
||||
pseudo_dens_c(1:nnr_) = CMPLX( w1*DBLE(red_revc0(1:nnr_,ibnd,1)) *&
|
||||
& AIMAG(red_revc0(1:nnr_,ibnd2-1,1)), 0.0d0, kind=DP )
|
||||
ENDIF
|
||||
!
|
||||
CALL fwfft ('Rho', pseudo_dens_c, dfftt)
|
||||
!
|
||||
! hartree contribution is computed in reciprocal space
|
||||
!
|
||||
DO is = 1, nspin
|
||||
!
|
||||
vhart(dfftt%nl(1:ngm_),is) =&
|
||||
& pseudo_dens_c(dfftt%nl(1:ngm_)) *&
|
||||
& fac_in(1:ngm_)
|
||||
IF (gamma_only) vhart(dfftt%nlm(1:ngm_),is) = &
|
||||
& pseudo_dens_c(dfftt%nlm(1:ngm_)) *&
|
||||
& fac_in(1:ngm_)
|
||||
!
|
||||
! and transformed back to real space
|
||||
!
|
||||
CALL invfft ('Rho', vhart (:, is), dfftt)
|
||||
!
|
||||
ENDDO
|
||||
!
|
||||
! Finally return the interaction psi_int for terms:
|
||||
! ibnd, ibnd, ibnd2
|
||||
! ibn2, ibnd2, ibnd
|
||||
!
|
||||
psi_int(1:nnr_,ibnd2) = psi_int(1:nnr_,ibnd2) &
|
||||
& + DBLE(vhart(1:nnr_, 1)) * DBLE(psi(1:nnr_))
|
||||
! & + w1*DBLE(vhart(1:nnr_, 1)) * DBLE(psi(1:nnr_))
|
||||
!
|
||||
CALL invfft_orbital_ibnd2_gamma(orbital(:,:), psitemp, ibnd2, npw_, dfftt)
|
||||
!
|
||||
! w3 = wg(ibnd2,1)/omega
|
||||
!
|
||||
psi_int(1:nnr_,ibnd) = psi_int(1:nnr_,ibnd) &
|
||||
& + DBLE(vhart(1:nnr_, 1)) * DBLE(psitemp(1:nnr_))
|
||||
! & + w3*DBLE(vhart(1:nnr_, 1)) * DBLE(psitemp(1:nnr_))
|
||||
!
|
||||
ENDDO
|
||||
|
||||
ENDIF
|
||||
!
|
||||
!
|
||||
END FUNCTION k1d_term_gamma
|
||||
|
||||
FUNCTION k2d_vexx_term_gamma(w1, w2, psi, fac_in, ibnd, interaction) RESULT (psi_int)
|
||||
!------------------------------------------------------------------
|
||||
!
|
||||
! This routine computes the K^2d term, Eq (22) from Eq Ref (1).
|
||||
! psi contains two bands of |b> with w1, w2 the associated weights
|
||||
! fac_in contains the interaction W(G) and ibnd the band index n'.
|
||||
!
|
||||
! Also computes the Vexx term, Eq (15) from Eq Ref (2).
|
||||
! So in gamma calculations, h_psi soubrite doesn't call anymore
|
||||
! vexx routine
|
||||
!
|
||||
! (1) Rocca, Lu and Galli, J. Chem. Phys., 133, 164109 (2010)
|
||||
! (2) Ge, Binnie, Rocca, Gebauer, Baroni, Comp. Phys. Comm.,
|
||||
! 185, 2080 (2014)
|
||||
!------------------------------------------------------------------
|
||||
!
|
||||
|
||||
IMPLICIT NONE
|
||||
|
||||
COMPLEX(KIND=DP),DIMENSION(:), INTENT(IN) :: psi
|
||||
REAL(KIND=DP),DIMENSION(:), INTENT(IN) :: fac_in
|
||||
REAL(kind=DP), INTENT(IN) :: w1, w2
|
||||
INTEGER, INTENT(IN) :: ibnd
|
||||
REAL(KIND=DP), ALLOCATABLE :: psi_int(:,:)
|
||||
LOGICAL, INTENT(IN) :: interaction
|
||||
!
|
||||
! Workspaces
|
||||
!
|
||||
INTEGER :: ibnd2, is, npw_,ngm_, nnr_
|
||||
!
|
||||
nnr_ = dfftt%nnr
|
||||
ngm_ = dfftt%ngm
|
||||
npw_ = npwt
|
||||
!
|
||||
ALLOCATE(psi_int(nnr_, nbnd))
|
||||
psi_int = 0.d0
|
||||
!
|
||||
!
|
||||
DO ibnd2=1,nbnd,1
|
||||
!
|
||||
! Set up the pseudo density for this Hartree like interaction.
|
||||
!
|
||||
vhart(:,:) = (0.d0,0.d0)
|
||||
pseudo_dens_c(:) = (0.d0,0.d0)
|
||||
!
|
||||
! The following code is to check if the value of ibnd2 is even or odd
|
||||
! and therefore whether the REAL or IMAGINARY part of red_revc0 is to be
|
||||
! used. This is because red_revc0 is stored using gamma_tricks.
|
||||
!
|
||||
IF (MOD(ibnd2,2)==1) THEN
|
||||
pseudo_dens_c(1:nnr_) = CMPLX( &
|
||||
& w1*DBLE(psi(1:nnr_))*DBLE(red_revc0(1:nnr_,ibnd2,1)),&
|
||||
& w2*AIMAG(psi(1:nnr_))*DBLE(red_revc0(1:nnr_,ibnd2,1)), kind=DP )
|
||||
!
|
||||
CALL fwfft ('Rho', pseudo_dens_c, dfftt)
|
||||
!
|
||||
! hartree contribution is computed in reciprocal space
|
||||
!
|
||||
DO is = 1, nspin
|
||||
!
|
||||
vhart(dfftt%nl(1:ngm_),is) = pseudo_dens_c(dfftt%nl(1:ngm_)) * fac_in(1:ngm_)
|
||||
IF (gamma_only) vhart(dfftt%nlm(1:ngm_),is) = &
|
||||
& pseudo_dens_c(dfftt%nlm(1:ngm_)) *&
|
||||
& fac_in(1:ngm_)
|
||||
!
|
||||
! and transformed back to real space
|
||||
!
|
||||
CALL invfft ('Rho', vhart (:, is), dfftt)
|
||||
!
|
||||
ENDDO
|
||||
!
|
||||
!
|
||||
! Finally return the interaction
|
||||
!
|
||||
psi_int(1:nnr_,ibnd2) = psi_int(1:nnr_,ibnd2) &
|
||||
& +DBLE(vhart(1:nnr_,1)) * DBLE(red_revc0(1:nnr_,ibnd,1)) &
|
||||
& +AIMAG(vhart(1:nnr_,1)) * AIMAG(red_revc0(1:nnr_,ibnd,1))
|
||||
!
|
||||
IF (interaction) then
|
||||
psi_int(1:nnr_,ibnd) = psi_int(1:nnr_,ibnd) &
|
||||
& +DBLE(vhart(1:nnr_,1)) * DBLE(red_revc0(1:nnr_,ibnd2,1))
|
||||
IF (ibnd < nbnd) then
|
||||
psi_int(1:nnr_,ibnd+1) = psi_int(1:nnr_,ibnd+1) &
|
||||
& +AIMAG(vhart(1:nnr_,1)) * DBLE(red_revc0(1:nnr_,ibnd2,1))
|
||||
ENDIF
|
||||
ELSE
|
||||
psi_int(1:nnr_,ibnd) = psi_int(1:nnr_,ibnd) &
|
||||
& -DBLE(vhart(1:nnr_,1)) * DBLE(red_revc0(1:nnr_,ibnd2,1))
|
||||
IF (ibnd < nbnd) then
|
||||
psi_int(1:nnr_,ibnd+1) = psi_int(1:nnr_,ibnd+1) &
|
||||
& -AIMAG(vhart(1:nnr_,1)) * DBLE(red_revc0(1:nnr_,ibnd2,1))
|
||||
ENDIF
|
||||
ENDIF
|
||||
!
|
||||
ELSE
|
||||
pseudo_dens_c(1:nnr_) = CMPLX( &
|
||||
& w1*DBLE(psi(1:nnr_))*AIMAG(red_revc0(1:nnr_,ibnd2-1,1)),&
|
||||
& w2*AIMAG(psi(1:nnr_))*AIMAG(red_revc0(1:nnr_,ibnd2-1,1)), kind=DP )
|
||||
!
|
||||
CALL fwfft ('Rho', pseudo_dens_c, dfftt)
|
||||
!
|
||||
! hartree contribution is computed in reciprocal space
|
||||
!
|
||||
DO is = 1, nspin
|
||||
!
|
||||
vhart(dfftt%nl(1:ngm_),is) = pseudo_dens_c(dfftt%nl(1:ngm_)) * fac_in(1:ngm_)
|
||||
IF (gamma_only) vhart(dfftt%nlm(1:ngm_),is) = &
|
||||
& pseudo_dens_c(dfftt%nlm(1:ngm_)) *&
|
||||
& fac_in(1:ngm_)
|
||||
!
|
||||
! and transformed back to real space
|
||||
!
|
||||
CALL invfft ('Rho', vhart (:, is), dfftt)
|
||||
!
|
||||
ENDDO
|
||||
!
|
||||
!
|
||||
! Finally return the interaction
|
||||
!
|
||||
psi_int(1:nnr_,ibnd2) = psi_int(1:nnr_,ibnd2) &
|
||||
& +DBLE(vhart(1:nnr_,1)) * DBLE(red_revc0(1:nnr_,ibnd,1)) &
|
||||
& +AIMAG(vhart(1:nnr_,1)) * AIMAG(red_revc0(1:nnr_,ibnd,1))
|
||||
!
|
||||
!
|
||||
! and interaction for Vexx term
|
||||
!
|
||||
IF (interaction) then
|
||||
psi_int(1:nnr_,ibnd) = psi_int(1:nnr_,ibnd) &
|
||||
& +DBLE(vhart(1:nnr_,1)) * AIMAG(red_revc0(1:nnr_,ibnd2-1,1))
|
||||
IF (ibnd < nbnd) then
|
||||
psi_int(1:nnr_,ibnd+1) = psi_int(1:nnr_,ibnd+1) &
|
||||
& +AIMAG(vhart(1:nnr_,1)) * AIMAG(red_revc0(1:nnr_,ibnd2-1,1))
|
||||
ENDIF
|
||||
ELSE
|
||||
psi_int(1:nnr_,ibnd) = psi_int(1:nnr_,ibnd) &
|
||||
& -DBLE(vhart(1:nnr_,1)) * AIMAG(red_revc0(1:nnr_,ibnd2-1,1))
|
||||
IF (ibnd < nbnd) then
|
||||
psi_int(1:nnr_,ibnd+1) = psi_int(1:nnr_,ibnd+1) &
|
||||
& -AIMAG(vhart(1:nnr_,1)) * AIMAG(red_revc0(1:nnr_,ibnd2-1,1))
|
||||
ENDIF
|
||||
ENDIF
|
||||
!
|
||||
ENDIF
|
||||
!
|
||||
!
|
||||
ENDDO
|
||||
!
|
||||
!
|
||||
END FUNCTION k2d_vexx_term_gamma
|
||||
|
||||
END MODULE lr_exx_kernel
|
||||
|
|
|
|||
Loading…
Reference in New Issue