diff --git a/PW/src/h_psi.f90 b/PW/src/h_psi.f90 index 5a31acf97..abbff7499 100644 --- a/PW/src/h_psi.f90 +++ b/PW/src/h_psi.f90 @@ -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 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) diff --git a/TDDFPT/src/lr_exx_kernel.f90 b/TDDFPT/src/lr_exx_kernel.f90 index e4dbb9fe2..4f6f92850 100644 --- a/TDDFPT/src/lr_exx_kernel.f90 +++ b/TDDFPT/src/lr_exx_kernel.f90 @@ -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