diff --git a/src/GT/GG0T0pp.f90 b/src/GT/GG0T0pp.f90 index c0012d4..b5453a8 100644 --- a/src/GT/GG0T0pp.f90 +++ b/src/GT/GG0T0pp.f90 @@ -1,4 +1,4 @@ -subroutine GG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE, & +subroutine GG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, & linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF) ! Perform one-shot calculation with a T-matrix self-energy (G0T0) @@ -13,7 +13,7 @@ subroutine GG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d logical,intent(in) :: doACFDT logical,intent(in) :: exchange_kernel logical,intent(in) :: doXBS - logical,intent(in) :: dophBSE + logical,intent(in) :: dophBSE,dophBSE2 logical,intent(in) :: doppBSE logical,intent(in) :: TDA_T logical,intent(in) :: TDA @@ -37,7 +37,7 @@ subroutine GG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d ! Local variables - logical :: print_T = .false. + logical :: print_T = .true. integer :: nOO integer :: nVV @@ -71,8 +71,7 @@ subroutine GG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d write(*,*)'* Generalized G0T0pp Calculation *' write(*,*)'**********************************' write(*,*) - - + ! TDA for T if(TDA_T) then @@ -104,22 +103,23 @@ subroutine GG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d ! Compute linear response allocate(Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO)) - + Bpp(:,:) = 0d0 + Cpp(:,:) = 0d0 + Dpp(:,:) = 0d0 call ppGLR_C(nOrb,nC,nO,nV,nR,nVV,1d0,eHF,ERI,Cpp) call ppGLR_D(nOrb,nC,nO,nV,nR,nOO,1d0,eHF,ERI,Dpp) if(.not.TDA_T) call ppGLR_B(nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) - + call ppGLR(TDA_T,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA) deallocate(Bpp,Cpp,Dpp) if(print_T) call print_excitation_energies('ppRPA@GHF','2p',nVV,Om1) - if(print_T) call print_excitation_energies('ppRPA@FHF','2h',nOO,Om2) + if(print_T) call print_excitation_energies('ppRPA@GHF','2h',nOO,Om2) !---------------------------------------------- ! Compute excitation densities !---------------------------------------------- - call GGTpp_excitation_density(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2) !---------------------------------------------- @@ -170,6 +170,24 @@ subroutine GG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d call print_GG0T0pp(nOrb,nO,eHF,ENuc,EGHF,Sig,Z,eGT,EcGM,EcRPA) +!---------------------------------------------- +! ppBSE calculation +!---------------------------------------------- + + if(doppBSE) then + + call GGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nOO,nVV, & + ERI,dipole_int,eHF,eGT,EcBSE) + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp@RHF correlation energy = ',EcBSE,' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp@RHF total energy = ',ENuc + EGHF + EcBSE,' au' + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + + end if + ! Testing zone if(dotest) then diff --git a/src/GT/GGT.f90 b/src/GT/GGT.f90 new file mode 100644 index 0000000..02e2c1a --- /dev/null +++ b/src/GT/GGT.f90 @@ -0,0 +1,117 @@ +subroutine GGT(dotest,doG0T0,doevGT,doqsGT,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE, & + TDA_T,TDA,dBSE,dTDA,linearize,eta,doSRG,nNuc,ZNuc,rNuc,ENuc,nBas,nBas2,nC,nO,nV,nR,nS,EGHF,S,X,T,V,Hc, & + ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,eHF) + +! GT module + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: dotest + + logical,intent(in) :: doG0T0 + logical,intent(in) :: doevGT + logical,intent(in) :: doqsGT + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + logical,intent(in) :: doACFDT + logical,intent(in) :: exchange_kernel + logical,intent(in) :: doXBS + logical,intent(in) :: dophBSE + logical,intent(in) :: dophBSE2 + logical,intent(in) :: TDA_T + logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA + logical,intent(in) :: doppBSE + logical,intent(in) :: linearize + double precision,intent(in) :: eta + logical,intent(in) :: doSRG + + integer,intent(in) :: nNuc + double precision,intent(in) :: ZNuc(nNuc) + double precision,intent(in) :: rNuc(nNuc,ncart) + double precision,intent(in) :: ENuc + + integer,intent(in) :: nBas + integer,intent(in) :: nBas2 + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + + double precision,intent(in) :: EGHF + double precision,intent(in) :: eHF(nBas2) + double precision,intent(in) :: cHF(nBas2,nBas2) + double precision,intent(in) :: PHF(nBas2,nBas2) + double precision,intent(in) :: S(nBas2,nBas2) + double precision,intent(in) :: T(nBas2,nBas2) + double precision,intent(in) :: V(nBas2,nBas2) + double precision,intent(in) :: Hc(nBas2,nBas2) + double precision,intent(in) :: X(nBas2,nBas2) + double precision,intent(in) :: ERI_AO(nBas2,nBas2,nBas2,nBas2) + double precision,intent(in) :: ERI(nBas2,nBas2,nBas2,nBas2) + double precision,intent(in) :: dipole_int_AO(nBas2,nBas2,ncart) + double precision,intent(in) :: dipole_int(nBas2,nBas2,ncart) + +! Local variables + + double precision :: start_GT ,end_GT ,t_GT + +!------------------------------------------------------------------------ +! Perform G0T0 calculatiom +!------------------------------------------------------------------------ + + if(doG0T0) then + call wall_time(start_GT) + call GG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, & + linearize,eta,doSRG,nBas2,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF) + call wall_time(end_GT) + + t_GT = end_GT - start_GT + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0 = ',t_GT,' seconds' + write(*,*) + + end if + +!------------------------------------------------------------------------ +! Perform evGT calculation +!------------------------------------------------------------------------ + + ! if(doevGT) then + + ! call wall_time(start_GT) + ! call evGGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, & + ! linearize,eta,doSRG,nBas2,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF) + ! call wall_time(end_GT) + + ! t_GT = end_GT - start_GT + ! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGT = ',t_GT,' seconds' + ! write(*,*) + + ! end if + +!------------------------------------------------------------------------ +! Perform qsGT calculation +!------------------------------------------------------------------------ + + ! if(doqsGT) then + + ! call wall_time(start_GT) + ! call qsGGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, & + ! eta,doSRG,nNuc,ZNuc,rNuc,ENuc,nBas,nBas2,nC,nO,nV,nR,nS,EGHF,S,X,T,V,Hc,ERI_AO,ERI, & + ! dipole_int_AO,dipole_int,PHF,cHF,eHF) + ! call wall_time(end_GT) + + ! t_GT = end_GT - start_GT + ! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGT = ',t_GT,' seconds' + ! write(*,*) + + ! end if + +end subroutine diff --git a/src/GT/GGT_Tmatrix.f90 b/src/GT/GGT_Tmatrix.f90 new file mode 100644 index 0000000..909f2ab --- /dev/null +++ b/src/GT/GGT_Tmatrix.f90 @@ -0,0 +1,74 @@ +subroutine GGT_Tmatrix(nOrb,nC,nO,nV,nR,nOO,nVV,lambda,ERI,eGF,Om1,rho1,Om2,rho2,T) + +! Compute the T-matrix tensor elements + + implicit none + include 'parameters.h' + + ! Input variables + + integer,intent(in) :: nOrb + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nOO + integer,intent(in) :: nVV + double precision,intent(in) :: lambda + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: eGF(nOrb) + double precision,intent(in) :: Om1(nVV) + double precision,intent(in) :: rho1(nOrb,nOrb,nVV) + double precision,intent(in) :: Om2(nOO) + double precision,intent(in) :: rho2(nOrb,nOrb,nOO) + + ! Local variables + + integer :: p,q,r,s + integer :: c,d,k,l + integer :: kl,cd + + ! Output variables + + double precision,intent(out) :: T(nOrb,nOrb,nOrb,nOrb) + + ! Initialization + T(:,:,:,:) = 0d0 + +! Start by building the tensor elements of T +! This is probabbly not a good idea because this tensor is really large + !$OMP PARALLEL & + !$OMP SHARED(nC,nO,nOrb,nR,T,ERI,rho1,rho2,Om1,Om2) & + !$OMP PRIVATE(p,q,r,s,c,d,cd,k,l,kl) & + !$OMP DEFAULT(NONE) + !$OMP DO + do s=nC+1,nOrb-nR + do r=nC+1,nOrb-nR + do q=nC+1,nOrb-nR + do p=nC+1,nOrb-nR + T(p,q,r,s) = ERI(p,q,r,s) - ERI(p,q,s,r) + + cd = 0 + do c = nO+1, nOrb-nR + do d = c+1, nOrb-nR + cd = cd + 1 + T(p,q,r,s) = T(p,q,r,s) - rho1(p,q,cd) * rho1(r,s,cd) / Om1(cd) + end do ! d + end do ! c + + kl = 0 + do k = nC+1, nO + do l = k+1, nO + kl = kl + 1 + T(p,q,r,s) = T(p,q,r,s) + rho2(p,q,kl) * rho2(r,s,kl) / Om2(kl) + enddo + enddo + + enddo + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + +end subroutine diff --git a/src/GT/GGTpp_excitation_density.f90 b/src/GT/GGTpp_excitation_density.f90 index 0550f25..b388077 100644 --- a/src/GT/GGTpp_excitation_density.f90 +++ b/src/GT/GGTpp_excitation_density.f90 @@ -1,4 +1,4 @@ -subroutine GGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2) +subroutine GGTpp_excitation_density(nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2) ! Compute excitation densities for T-matrix self-energy @@ -12,7 +12,6 @@ subroutine GGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho ! Input variables - integer,intent(in) :: ispin integer,intent(in) :: nBas integer,intent(in) :: nC integer,intent(in) :: nO diff --git a/src/GT/GGTpp_ppBSE.f90 b/src/GT/GGTpp_ppBSE.f90 new file mode 100644 index 0000000..5f8a861 --- /dev/null +++ b/src/GT/GGTpp_ppBSE.f90 @@ -0,0 +1,127 @@ +subroutine GGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nOO,nVV, & + ERI,dipole_int,eT,eGT,EcBSE) + +! Compute the Bethe-Salpeter excitation energies with the T-matrix kernel + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: TDA_T + logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA + + double precision,intent(in) :: eta + integer,intent(in) :: nOrb + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + + integer,intent(in) :: nOO + integer,intent(in) :: nVV + + double precision,intent(in) :: eT(nOrb) + double precision,intent(in) :: eGT(nOrb) + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart) + +! Local variables + + double precision :: EcRPA(nspin) + double precision,allocatable :: Bpp(:,:),Cpp(:,:),Dpp(:,:) + double precision,allocatable :: Om1(:), Om2(:) + double precision,allocatable :: X1(:,:), X2(:,:) + double precision,allocatable :: Y1(:,:), Y2(:,:) + double precision,allocatable :: rho1(:,:,:), rho2(:,:,:) + double precision,allocatable :: KB_sta(:,:),KC_sta(:,:),KD_sta(:,:) + double precision,allocatable :: T(:,:,:,:) + +! Output variables + + double precision,intent(out) :: EcBSE + +!---------------------------------------------- +! Compute linear response +!---------------------------------------------- + allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO)) + allocate(Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO)) + Bpp(:,:) = 0d0 + Cpp(:,:) = 0d0 + Dpp(:,:) = 0d0 + call ppGLR_C(nOrb,nC,nO,nV,nR,nVV,1d0,eT,ERI,Cpp) + call ppGLR_D(nOrb,nC,nO,nV,nR,nOO,1d0,eT,ERI,Dpp) + if(.not.TDA_T) call ppGLR_B(nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) + + call ppGLR(TDA_T,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA) + + deallocate(Bpp,Cpp,Dpp) + +!---------------------------------------------- +! Compute excitation densities +!---------------------------------------------- + allocate(rho1(nOrb,nOrb,nVV),rho2(nOrb,nOrb,nOO)) + call GGTpp_excitation_density(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2) + + deallocate(X1,Y1,X2,Y2) + +!---------------------------------------------- +! Compute new ppRPA block +!---------------------------------------------- + + allocate(Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO)) + + call ppGLR_C(nOrb,nC,nO,nV,nR,nVV,1d0,eGT,ERI,Cpp) + call ppGLR_D(nOrb,nC,nO,nV,nR,nOO,1d0,eGT,ERI,Dpp) + if(.not.TDA_T) call ppGLR_B(nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) + +!---------------------------------------------- +! Compute T matrix tensor +!---------------------------------------------- + + allocate(T(nOrb,nOrb,nOrb,nOrb)) + + call GGT_Tmatrix(nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGT,Om1,rho1,Om2,rho2,T) + +!---------------------------------------------- +! Compute kernels +!---------------------------------------------- + + allocate(KB_sta(nVV,nOO),KC_sta(nVV,nVV),KD_sta(nOO,nOO)) + + call GGTpp_ppBSE_static_kernel_C(eta,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGT, & + Om1,rho1,Om2,rho2,T,KC_sta) + call GGTpp_ppBSE_static_kernel_D(eta,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGT, & + Om1,rho1,Om2,rho2,T,KD_sta) + if(.not.TDA_T) call GGTpp_ppBSE_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGT, & + Om1,rho1,Om2,rho2,T,KB_sta) + + deallocate(Om1,Om2,rho1,rho2) +! Deallocate the 4-tensor T + deallocate(T) + +!---------------------------------------------- +! Diagonalize ppBSE +!---------------------------------------------- + + Bpp(:,:) = Bpp(:,:) + KB_sta(:,:) + Cpp(:,:) = Cpp(:,:) + KC_sta(:,:) + Dpp(:,:) = Dpp(:,:) + KD_sta(:,:) + + allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO)) + + call ppGLR(TDA_T,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE) + + call ppLR_transition_vectors(.true.,nOrb,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2) + + !----------------------------------------------------! + ! Compute the dynamical screening at the ppBSE level ! + !----------------------------------------------------! + + ! TODO + + deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta) + +end subroutine diff --git a/src/GT/GGTpp_ppBSE_static_kernel_B.f90 b/src/GT/GGTpp_ppBSE_static_kernel_B.f90 index 2afe373..e70bcee 100644 --- a/src/GT/GGTpp_ppBSE_static_kernel_B.f90 +++ b/src/GT/GGTpp_ppBSE_static_kernel_B.f90 @@ -1,42 +1,48 @@ -subroutine GGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,nVVx,lambda,Om1,rho1,Om2,rho2,TB) +subroutine GGTpp_ppBSE_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nOO,nVV,lambda,ERI,eGF,Om1,rho1,Om2,rho2,T,KB_sta) -! Compute the VVOO block of the static T-matrix +! Compute the VVOO of the T-matrix static pp kernel implicit none include 'parameters.h' ! Input variables - integer,intent(in) :: ispin double precision,intent(in) :: eta - integer,intent(in) :: nBas + integer,intent(in) :: nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR integer,intent(in) :: nOO integer,intent(in) :: nVV - integer,intent(in) :: nOOx - integer,intent(in) :: nVVx double precision,intent(in) :: lambda + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: eGF(nOrb) double precision,intent(in) :: Om1(nVV) - double precision,intent(in) :: rho1(nBas,nBas,nVV) + double precision,intent(in) :: rho1(nOrb,nOrb,nVV) double precision,intent(in) :: Om2(nOO) - double precision,intent(in) :: rho2(nBas,nBas,nOO) + double precision,intent(in) :: rho2(nOrb,nOrb,nOO) + double precision,intent(in) :: T(nOrb,nOrb,nOrb,nOrb) ! Local variables - double precision :: chi - double precision :: eps - integer :: i,j,a,b,ij,ab,cd,kl + double precision :: dem,num + integer :: p,q,r,s,e,m + integer :: i,j,a,b + integer :: ij,ab,kl,cd ! Output variables - double precision,intent(out) :: TB(nVVx,nOOx) + double precision,intent(out) :: KB_sta(nVV,nOO) + +! Initialization + KB_sta(:,:) = 0d0 +! Computing the kernel +! This is the same code as for the GF(2) kernel with elements T instead of ERI ab = 0 - do a=nO+1,nBas-nR - do b=a+1,nBas-nR + do a=nO+1,nOrb-nR + do b=a+1,nOrb-nR ab = ab + 1 ij = 0 @@ -44,25 +50,23 @@ subroutine GGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n do j=i+1,nO ij = ij + 1 - chi = 0d0 - - do cd=1,nVV - eps = + Om1(cd) - chi = chi + rho1(a,b,cd)*rho1(i,j,cd)*eps/(eps**2 + eta**2) + do m=nC+1,nO + do e=nO+1,nOrb-nR + + dem = eGF(m) - eGF(e) + num = (T(a,m,i,e) - T(a,m,e,i)) * (T(e,b,m,j) - T(e,b,j,m)) + num = num + (T(a,e,i,m) - T(a,e,m,i)) * (T(m,b,e,j) - T(m,b,j,e)) + num = num - (T(b,m,i,e) - T(b,m,e,i)) * (T(e,a,m,j) - T(e,a,j,m)) + num = num - (T(b,e,i,m) - T(b,e,m,i)) * (T(m,a,e,j) - T(m,a,j,e)) + KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2) + + end do end do - do kl=1,nOO - eps = - Om2(kl) - chi = chi + rho2(a,b,kl)*rho2(i,j,kl)*eps/(eps**2 + eta**2) - end do - - TB(ab,ij) = lambda*chi - end do end do end do end do - - + end subroutine diff --git a/src/GT/GGTpp_ppBSE_static_kernel_C.f90 b/src/GT/GGTpp_ppBSE_static_kernel_C.f90 index 0ce97d4..3808d90 100644 --- a/src/GT/GGTpp_ppBSE_static_kernel_C.f90 +++ b/src/GT/GGTpp_ppBSE_static_kernel_C.f90 @@ -1,67 +1,73 @@ -subroutine GGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,nVVx,lambda,Om1,rho1,Om2,rho2,TC) +subroutine GGTpp_ppBSE_static_kernel_C(eta,nOrb,nC,nO,nV,nR,nOO,nVV,lambda,ERI,eGF,Om1,rho1,Om2,rho2,T,KC_sta) -! Compute the VVVV block of the static T-matrix +! Compute the VVVV block of the T-matrix static pp kernels implicit none include 'parameters.h' ! Input variables - integer,intent(in) :: ispin double precision,intent(in) :: eta - integer,intent(in) :: nBas + integer,intent(in) :: nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR integer,intent(in) :: nOO integer,intent(in) :: nVV - integer,intent(in) :: nOOx - integer,intent(in) :: nVVx double precision,intent(in) :: lambda + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: eGF(nOrb) double precision,intent(in) :: Om1(nVV) - double precision,intent(in) :: rho1(nBas,nBas,nVV) + double precision,intent(in) :: rho1(nOrb,nOrb,nVV) double precision,intent(in) :: Om2(nOO) - double precision,intent(in) :: rho2(nBas,nBas,nOO) + double precision,intent(in) :: rho2(nOrb,nOrb,nOO) + double precision,intent(in) :: T(nOrb,nOrb,nOrb,nOrb) + ! Local variables - double precision,external :: Kronecker_delta - double precision :: chi - double precision :: eps - integer :: a,b,c,d,ab,cd,ef,mn + double precision :: dem,num + integer :: p,q,r,s,e,m + integer :: a,b,c,d,k,l + integer :: ab,kl,cd ! Output variables - double precision,intent(out) :: TC(nVVx,nVVx) + double precision,intent(out) :: KC_sta(nVV,nVV) +! Initialization + KC_sta(:,:) = 0d0 + +! Computing the kernel +! This is the same code as for the GF(2) kernel with elements T instead of ERI ab = 0 - do a=nO+1,nBas-nR - do b=a+1,nBas-nR + do a=nO+1,nOrb-nR + do b=a+1,nOrb-nR ab = ab + 1 - + cd = 0 - do c=nO+1,nBas-nR - do d=c+1,nBas-nR + do c=nO+1,nOrb-nR + do d=c+1,nOrb-nR cd = cd + 1 - - chi = 0d0 - - do ef=1,nVV - eps = + Om1(ef) - chi = chi + rho1(a,b,ef)*rho1(c,d,ef)*eps/(eps**2 + eta**2) + + do m=nC+1,nO + do e=nO+1,nOrb-nR + + dem = eGF(m) - eGF(e) + num = (T(a,m,c,e) - T(a,m,e,c)) * (T(e,b,m,d) - T(e,b,d,m)) + num = num + (T(a,e,c,m) - T(a,e,m,c)) * (T(m,b,e,d) - T(m,b,d,e)) + num = num - (T(b,m,c,e) - T(b,m,e,c)) * (T(e,a,m,d) - T(e,a,d,m)) + num = num - (T(b,e,c,m) - T(b,e,m,c)) * (T(m,a,e,d) - T(m,a,d,e)) + + KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2) + + end do end do - - do mn=1,nOO - eps = - Om2(mn) - chi = chi + rho2(a,b,mn)*rho2(c,d,mn)*eps/(eps**2 + eta**2) - end do - - TC(ab,cd) = lambda*chi - + end do end do - + end do end do diff --git a/src/GT/GGTpp_ppBSE_static_kernel_D.f90 b/src/GT/GGTpp_ppBSE_static_kernel_D.f90 index b920a76..15e2e4f 100644 --- a/src/GT/GGTpp_ppBSE_static_kernel_D.f90 +++ b/src/GT/GGTpp_ppBSE_static_kernel_D.f90 @@ -1,39 +1,45 @@ -subroutine GGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,nVVx,lambda,Om1,rho1,Om2,rho2,TD) +subroutine GGTpp_ppBSE_static_kernel_D(eta,nOrb,nC,nO,nV,nR,nOO,nVV,lambda,ERI,eGF,Om1,rho1,Om2,rho2,T,KD_sta) -! Compute the OOOO block of the static T-matrix +! Compute the OOOO block of the T-matrix static pp kernel implicit none include 'parameters.h' ! Input variables - integer,intent(in) :: ispin double precision,intent(in) :: eta - integer,intent(in) :: nBas + integer,intent(in) :: nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR integer,intent(in) :: nOO integer,intent(in) :: nVV - integer,intent(in) :: nOOx - integer,intent(in) :: nVVx double precision,intent(in) :: lambda + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: eGF(nOrb) double precision,intent(in) :: Om1(nVV) - double precision,intent(in) :: rho1(nBas,nBas,nVV) + double precision,intent(in) :: rho1(nOrb,nOrb,nVV) double precision,intent(in) :: Om2(nOO) - double precision,intent(in) :: rho2(nBas,nBas,nOO) + double precision,intent(in) :: rho2(nOrb,nOrb,nOO) + double precision,intent(in) :: T(nOrb,nOrb,nOrb,nOrb) ! Local variables - double precision :: chi - double precision :: eps - integer :: i,j,k,l,ij,kl,ef,mn + double precision :: dem,num + integer :: p,q,r,s,e,m + integer :: i,j,k,l + integer :: ij,kl,cd ! Output variables - double precision,intent(out) :: TD(nOOx,nOOx) + double precision,intent(out) :: KD_sta(nOO,nOO) + +! Initialization + KD_sta(:,:) = 0d0 +! Computing the kernel +! This is the same code as for the GF(2) kernel with elements T instead of ERI ij = 0 do i=nC+1,nO do j=i+1,nO @@ -44,24 +50,20 @@ subroutine GGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n do l=k+1,nO kl = kl + 1 - chi = 0d0 - - do ef=1,nVV - eps = + Om1(ef) - chi = chi + rho1(i,j,ef)*rho1(k,l,ef)*eps/(eps**2 + eta**2) + do m=nC+1,nO + do e=nO+1,nOrb-nR + + dem = eGF(m) - eGF(e) + num = T(i,m,k,e) * T(e,j,m,l) + T(i,e,k,m) * T(m,j,e,l) + num = num - T(j,m,k,e) * T(e,i,m,l) - T(j,e,k,m) * T(m,i,e,l) + + KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2) + end do end do - - do mn=1,nOO - eps = - Om2(mn) - chi = chi + rho2(i,j,mn)*rho2(k,l,mn)*eps/(eps**2 + eta**2) - end do - - TD(ij,kl) = lambda*chi - + end do end do - end do end do - + end subroutine diff --git a/src/GT/RG0T0pp.f90 b/src/GT/RG0T0pp.f90 index ac10827..6794419 100644 --- a/src/GT/RG0T0pp.f90 +++ b/src/GT/RG0T0pp.f90 @@ -212,6 +212,7 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d call ppLR_D(isp_T,nOrb,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp) call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T)) +! call ppLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nOOs,nVVs,dipole_int,Om1s,X1s,Y1s,Om2s,X2s,Y2s) deallocate(Bpp,Cpp,Dpp) @@ -224,6 +225,7 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d call ppLR_D(isp_T,nOrb,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp) call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T)) +! call ppLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,dipole_int,Om1t,X1t,Y1t,Om2t,X2t,Y2t) deallocate(Bpp,Cpp,Dpp) diff --git a/src/GT/RGT.f90 b/src/GT/RGT.f90 index c18d794..9490fa2 100644 --- a/src/GT/RGT.f90 +++ b/src/GT/RGT.f90 @@ -1,6 +1,3 @@ - -! --- - subroutine RGT(dotest, doG0T0pp, doevGTpp, doqsGTpp, doufG0T0pp, doG0T0eh, doevGTeh, doqsGTeh, & maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, dophBSE2, & doppBSE, TDA_T, TDA, dBSE, dTDA, singlet, triplet, linearize, eta, regularize, & diff --git a/src/GT/RGTpp_SigC.f90 b/src/GT/RGTpp_SigC.f90 index 5f69f7f..1edd7bb 100644 --- a/src/GT/RGTpp_SigC.f90 +++ b/src/GT/RGTpp_SigC.f90 @@ -42,12 +42,12 @@ double precision function RGTpp_SigC(p,w,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVV do cd=1,nVVs eps = w + e(i) - Om1s(cd) - RGTpp_SigC = RGTpp_SigC + rho1s(p,i,cd)**2*eps/(eps**2 + eta**2) + RGTpp_SigC = RGTpp_SigC + (1d0/2d0)*rho1s(p,i,cd)**2*eps/(eps**2 + eta**2) end do do cd=1,nVVt eps = w + e(i) - Om1t(cd) - RGTpp_SigC = RGTpp_SigC + rho1t(p,i,cd)**2*eps/(eps**2 + eta**2) + RGTpp_SigC = RGTpp_SigC + (3d0/2d0)*rho1t(p,i,cd)**2*eps/(eps**2 + eta**2) end do end do @@ -60,12 +60,12 @@ double precision function RGTpp_SigC(p,w,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVV do kl=1,nOOs eps = w + e(a) - Om2s(kl) - RGTpp_SigC = RGTpp_SigC + rho2s(p,a,kl)**2*eps/(eps**2 + eta**2) + RGTpp_SigC = RGTpp_SigC + (1d0/2d0)*rho2s(p,a,kl)**2*eps/(eps**2 + eta**2) end do do kl=1,nOOt eps = w + e(a) - Om2t(kl) - RGTpp_SigC = RGTpp_SigC + rho2t(p,a,kl)**2*eps/(eps**2 + eta**2) + RGTpp_SigC = RGTpp_SigC + (3d0/2d0)*rho2t(p,a,kl)**2*eps/(eps**2 + eta**2) end do end do diff --git a/src/GT/RGTpp_dSigC.f90 b/src/GT/RGTpp_dSigC.f90 index 7e0404f..7e5cac4 100644 --- a/src/GT/RGTpp_dSigC.f90 +++ b/src/GT/RGTpp_dSigC.f90 @@ -42,31 +42,30 @@ double precision function RGTpp_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV do cd=1,nVVs eps = w + e(i) - Om1s(cd) - RGTpp_dSigC = RGTpp_dSigC - rho1s(p,i,cd)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + RGTpp_dSigC = RGTpp_dSigC - (1d0/2d0)*rho1s(p,i,cd)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do do cd=1,nVVt eps = w + e(i) - Om1t(cd) - RGTpp_dSigC = RGTpp_dSigC - rho1t(p,i,cd)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + RGTpp_dSigC = RGTpp_dSigC - (3d0/2d0)*rho1t(p,i,cd)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do end do !---------------------------------------------- ! Virtual part of the T-matrix self-energy -! !---------------------------------------------- do a=nO+1,nBas-nR do kl=1,nOOs eps = w + e(a) - Om2s(kl) - RGTpp_dSigC = RGTpp_dSigC - rho2s(p,a,kl)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + RGTpp_dSigC = RGTpp_dSigC - (1d0/2d0)*rho2s(p,a,kl)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do do kl=1,nOOt eps = w + e(a) - Om2t(kl) - RGTpp_dSigC = RGTpp_dSigC - rho2t(p,a,kl)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + RGTpp_dSigC = RGTpp_dSigC - (3d0/2d0)*rho2t(p,a,kl)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do end do diff --git a/src/GT/RGTpp_excitation_density.f90 b/src/GT/RGTpp_excitation_density.f90 index 8a853f1..bddd33d 100644 --- a/src/GT/RGTpp_excitation_density.f90 +++ b/src/GT/RGTpp_excitation_density.f90 @@ -66,7 +66,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho cd = cd + 1 rho1(p,q,ab) = rho1(p,q,ab) & + (ERI(p,q,c,d) + ERI(p,q,d,c))*X1(cd,ab)/ & - sqrt(1d0 + Kronecker_delta(c,d))/sqrt(2d0) + sqrt(1d0 + Kronecker_delta(c,d)) end do end do @@ -76,7 +76,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho kl = kl + 1 rho1(p,q,ab) = rho1(p,q,ab) & + (ERI(p,q,k,l) + ERI(p,q,l,k))*Y1(kl,ab)/ & - sqrt(1d0 + Kronecker_delta(k,l))/sqrt(2d0) + sqrt(1d0 + Kronecker_delta(k,l)) end do end do @@ -93,7 +93,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho cd = cd + 1 rho2(p,q,ij) = rho2(p,q,ij) & + (ERI(p,q,c,d) + ERI(p,q,d,c))*X2(cd,ij)/ & - sqrt(1d0 + Kronecker_delta(c,d))/sqrt(2d0) + sqrt(1d0 + Kronecker_delta(c,d)) end do end do @@ -103,7 +103,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho kl = kl + 1 rho2(p,q,ij) = rho2(p,q,ij) & + (ERI(p,q,k,l) + ERI(p,q,l,k))*Y2(kl,ij)/ & - sqrt(1d0 + Kronecker_delta(k,l))/sqrt(2d0) + sqrt(1d0 + Kronecker_delta(k,l)) end do end do @@ -148,7 +148,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho cd = cd + 1 rho1(p,q,ab) = rho1(p,q,ab) & - + sqrt(3d0/2d0)*(ERI(p,q,c,d) - ERI(p,q,d,c))*X1(cd,ab) + + (ERI(p,q,c,d) - ERI(p,q,d,c))*X1(cd,ab) end do ! d end do ! c @@ -159,7 +159,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho kl = kl + 1 rho1(p,q,ab) = rho1(p,q,ab) & - + sqrt(3d0/2d0)*(ERI(p,q,k,l) - ERI(p,q,l,k))*Y1(kl,ab) + + (ERI(p,q,k,l) - ERI(p,q,l,k))*Y1(kl,ab) end do ! l end do ! k end do ! b @@ -179,7 +179,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho cd = cd + 1 rho2(p,q,ij) = rho2(p,q,ij) & - + sqrt(3d0/2d0)*(ERI(p,q,c,d) - ERI(p,q,d,c))*X2(cd,ij) + + (ERI(p,q,c,d) - ERI(p,q,d,c))*X2(cd,ij) end do ! d end do ! c @@ -190,7 +190,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho kl = kl + 1 rho2(p,q,ij) = rho2(p,q,ij) & - + sqrt(3d0/2d0)*(ERI(p,q,k,l) - ERI(p,q,l,k))*Y2(kl,ij) + + (ERI(p,q,k,l) - ERI(p,q,l,k))*Y2(kl,ij) end do ! l end do ! k end do ! j @@ -249,8 +249,8 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho deallocate(ERI_1, ERI_2) - rho1 = rho1*sqrt(3d0/2d0) - rho2 = rho2*sqrt(3d0/2d0) + rho1 = rho1 + rho2 = rho2 endif endif diff --git a/src/GT/RGTpp_ppBSE.f90 b/src/GT/RGTpp_ppBSE.f90 index 48a2788..1751075 100644 --- a/src/GT/RGTpp_ppBSE.f90 +++ b/src/GT/RGTpp_ppBSE.f90 @@ -34,7 +34,7 @@ subroutine RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR, ! Local variables - integer :: ispin + integer :: ispin, isp_T double precision :: EcRPA(nspin) double precision,allocatable :: Bpp(:,:),Cpp(:,:),Dpp(:,:) @@ -46,67 +46,75 @@ subroutine RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR, double precision,allocatable :: X2s(:,:),X2t(:,:) double precision,allocatable :: Y2s(:,:),Y2t(:,:) double precision,allocatable :: rho2s(:,:,:),rho2t(:,:,:) - double precision,allocatable :: TBs(:,:),TCs(:,:),TDs(:,:) - double precision,allocatable :: TBt(:,:),TCt(:,:),TDt(:,:) + double precision,allocatable :: KB_sta(:,:) + double precision,allocatable :: KC_sta(:,:) + double precision,allocatable :: KD_sta(:,:) + double precision,allocatable :: Taaaa(:,:,:,:),Tabab(:,:,:,:),Tbaab(:,:,:,:) ! Output variables double precision,intent(out) :: EcBSE(nspin) -!------------------------------------! -! Compute T-matrix for singlet block ! -!------------------------------------! - - ispin = 1 +!--------------------------------- +! Compute ppRPA excitation density +!--------------------------------- - allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs),Om2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), & - Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) - - if(.not.TDA_T) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp) - call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVVs,1d0,eT,ERI,Cpp) - call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOOs,1d0,eT,ERI,Dpp) - - call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin)) - - deallocate(Bpp,Cpp,Dpp) - allocate(TBs(nVVs,nOOs),TCs(nVVs,nVVs),TDs(nOOs,nOOs)) - - if(.not.TDA_T) call RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOs,nVVs,1d0, & - Om1s,rho1s,Om2s,rho2s,TBs) - call RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOs,nVVs,1d0, & - Om1s,rho1s,Om2s,rho2s,TCs) - call RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOs,nVVs,1d0, & - Om1s,rho1s,Om2s,rho2s,TDs) + isp_T = 1 + allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) - deallocate(Om1s,X1s,Y1s,Om2s,X2s,Y2s) + if(.not.TDA_T) call ppLR_B(isp_T,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp) + call ppLR_C(isp_T,nBas,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp) + call ppLR_D(isp_T,nBas,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp) -!------------------------------------! -! Compute T-matrix for triplet block ! -!------------------------------------! + allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs)) + allocate(Om2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs)) + call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T)) +! call ppLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nOOs,nVVs,dipole_int,Om1s,X1s,Y1s,Om2s,X2s,Y2s) - ispin = 2 + + allocate(rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs)) + call RGTpp_excitation_density(isp_T,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s) - allocate(Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt),Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), & - Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) + deallocate(X1s,Y1s,X2s,Y2s,Bpp,Cpp,Dpp) - if(.not.TDA_T) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp) - call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVVt,1d0,eT,ERI,Cpp) - call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOOt,1d0,eT,ERI,Dpp) + isp_T = 2 + allocate(Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt)) + allocate(Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt)) + allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) - call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin)) + if(.not.TDA_T) call ppLR_B(isp_T,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp) + call ppLR_C(isp_T,nBas,nC,nO,nV,nR,nVVt,1d0,eGT,ERI,Cpp) + call ppLR_D(isp_T,nBas,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp) - deallocate(Bpp,Cpp,Dpp) - allocate(TBt(nVVt,nOOt),TCt(nVVt,nVVt),TDt(nOOt,nOOt)) + call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T)) +! call ppLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,dipole_int,Om1t,X1t,Y1t,Om2t,X2t,Y2t) - if(.not.TDA_T) call RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,nOOs,nVVs,1d0, & - Om1t,rho1t,Om2t,rho2t,TBt) - call RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,nOOs,nVVs,1d0, & - Om1t,rho1t,Om2t,rho2t,TCt) - call RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,nOOs,nVVs,1d0, & - Om1t,rho1t,Om2t,rho2t,TDt) + allocate(rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt)) + call RGTpp_excitation_density(isp_T,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t) - deallocate(Om1t,X1t,Y1t,Om2t,X2t,Y2t) + deallocate(X1t,Y1t,X2t,Y2t,Bpp,Cpp,Dpp) + +!--------------------------------- +! Compute T matrix elements +!--------------------------------- + + ! This correspond to the alpha-alpha-alpha-alpha elements + isp_T = 1 + allocate(Taaaa(nBas,nBas,nBas,nBas)) + call RGT_Tmatrix(isp_T,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,1d0,ERI,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t,Taaaa) + ! This correspond to the alpha-alpha-alpha-alpha elements + isp_T = 2 + allocate(Tabab(nBas,nBas,nBas,nBas)) + call RGT_Tmatrix(isp_T,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,1d0,ERI,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t,Tabab) + + ! This correspond to the alpha-alpha-alpha-alpha elements + isp_T = 3 + allocate(Tbaab(nBas,nBas,nBas,nBas)) + call RGT_Tmatrix(isp_T,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,1d0,ERI,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t,Tbaab) + + deallocate(Om1s,Om2s,Om1t,Om2t,rho1s,rho2s,rho1t,rho2t) + !------------------! ! Singlet manifold ! !------------------! @@ -114,23 +122,35 @@ subroutine RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR, if(singlet) then ispin = 1 + + allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs)) + allocate(Om2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs)) - allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs),Om2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), & - Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) + ! Compute BSE excitation energies + + allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) + allocate(KB_sta(nVVs,nOOs),KC_sta(nVVs,nVVs),KD_sta(nOOs,nOOs)) if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp) call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp) call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp) - Bpp(:,:) = Bpp(:,:) - TBs(:,:) - TBt(:,:) - Cpp(:,:) = Cpp(:,:) - TCs(:,:) - TCt(:,:) - Dpp(:,:) = Dpp(:,:) - TDs(:,:) - TDt(:,:) + if(.not.TDA_T) call RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT, & + Taaaa,Tabab,Tbaab,KB_sta) + call RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT, & + Taaaa,Tabab,Tbaab,KC_sta) + call RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT, & + Taaaa,Tabab,Tbaab,KD_sta) + + Bpp(:,:) = Bpp(:,:) + KB_sta(:,:) + Cpp(:,:) = Cpp(:,:) + KC_sta(:,:) + Dpp(:,:) = Dpp(:,:) + KD_sta(:,:) call ppLR(TDA,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcBSE(ispin)) call ppLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nOOs,nVVs,dipole_int,Om1s,X1s,Y1s,Om2s,X2s,Y2s) - deallocate(Om1s,X1s,Y1s,Om2s,X2s,Y2s,Bpp,Cpp,Dpp) + deallocate(Om1s,X1s,Y1s,Om2s,X2s,Y2s,Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta) end if @@ -141,21 +161,33 @@ subroutine RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR, if(triplet) then ispin = 2 - + EcBSE(ispin) = 0d0 - allocate(Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt),Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), & - Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) + allocate(Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt)) + allocate(Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt)) + + ! Compute BSE excitation energies + + allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) + allocate(KB_sta(nVVt,nOOt),KC_sta(nVVt,nVVt),KD_sta(nOOt,nOOt)) - if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp) - call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp) - call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp) + if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp) + call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVVt,1d0,eGT,ERI,Cpp) + call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp) - Bpp(:,:) = Bpp(:,:) + TBs(:,:) - TBt(:,:) - Cpp(:,:) = Cpp(:,:) + TCs(:,:) - TCt(:,:) - Dpp(:,:) = Dpp(:,:) + TDs(:,:) - TDt(:,:) + if(.not.TDA_T) call RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT, & + Taaaa,Tabab,Tbaab,KB_sta) + call RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT, & + Taaaa,Tabab,Tbaab,KC_sta) + call RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT, & + Taaaa,Tabab,Tbaab,KD_sta) + + Bpp(:,:) = Bpp(:,:) + KB_sta(:,:) + Cpp(:,:) = Cpp(:,:) + KC_sta(:,:) + Dpp(:,:) = Dpp(:,:) + KD_sta(:,:) - call ppLR(TDA,nOOs,nVVs,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcBSE(ispin)) + call ppLR(TDA,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcBSE(ispin)) call ppLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,dipole_int,Om1t,X1t,Y1t,Om2t,X2t,Y2t) diff --git a/src/GT/RGTpp_ppBSE_static_kernel_B.f90 b/src/GT/RGTpp_ppBSE_static_kernel_B.f90 index 9b84f3f..de14f2c 100644 --- a/src/GT/RGTpp_ppBSE_static_kernel_B.f90 +++ b/src/GT/RGTpp_ppBSE_static_kernel_B.f90 @@ -1,4 +1,4 @@ -subroutine RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,nVVx,lambda,Om1,rho1,Om2,rho2,TB) +subroutine RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,lambda,eGF,Taaaa,Tabab,Tbaab,KB_sta) ! Compute the VVOO block of the static T-matrix @@ -16,13 +16,11 @@ subroutine RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n integer,intent(in) :: nR integer,intent(in) :: nOO integer,intent(in) :: nVV - integer,intent(in) :: nOOx - integer,intent(in) :: nVVx double precision,intent(in) :: lambda - double precision,intent(in) :: Om1(nVV) - double precision,intent(in) :: rho1(nBas,nBas,nVV) - double precision,intent(in) :: Om2(nOO) - double precision,intent(in) :: rho2(nBas,nBas,nOO) + double precision,intent(in) :: eGF(nBas) + double precision,intent(in) :: Taaaa(nBas,nBas,nBas,nBas) + double precision,intent(in) :: Tabab(nBas,nBas,nBas,nBas) + double precision,intent(in) :: Tbaab(nBas,nBas,nBas,nBas) ! Local variables @@ -32,8 +30,12 @@ subroutine RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n ! Output variables - double precision,intent(out) :: TB(nVVx,nOOx) + double precision,intent(out) :: KB_sta(nVV,nOO) +! Initialization + + KB_sta(:,:) = 0d0 + !===============! ! singlet block ! !===============! @@ -53,16 +55,16 @@ subroutine RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n chi = 0d0 do cd=1,nVV - eps = + Om1(cd) - chi = chi + rho1(a,b,cd)*rho1(i,j,cd)*eps/(eps**2 + eta**2) + eps = 0d0 + chi = chi + 0d0 end do do kl=1,nOO - eps = - Om2(kl) - chi = chi + rho2(a,b,kl)*rho2(i,j,kl)*eps/(eps**2 + eta**2) + eps = 0d0 + chi = chi + 0d0 end do - TB(ab,ij) = lambda*chi + KB_sta(ab,ij) = lambda*chi end do end do @@ -91,54 +93,16 @@ subroutine RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n chi = 0d0 do cd=1,nVV - eps = + Om1(cd) - chi = chi + rho1(a,b,cd)*rho1(i,j,cd)*eps/(eps**2 + eta**2) + eps = 0d0 + chi = chi + 0d0 end do do kl=1,nOO - eps = - Om2(kl) - chi = chi + rho2(a,b,kl)*rho2(i,j,kl)*eps/(eps**2 + eta**2) + eps = 0d0 + chi = chi + 0d0 end do - TB(ab,ij) = lambda*chi - - end do - end do - - end do - end do - - end if - -!==================! -! alpha-beta block ! -!==================! - - if(ispin == 3) then - - ab = 0 - do a=nO+1,nBas-nR - do b=nO+1,nBas-nR - ab = ab + 1 - - ij = 0 - do i=nC+1,nO - do j=nC+1,nO - ij = ij + 1 - - chi = 0d0 - - do cd=1,nVV - eps = + Om1(cd) - chi = chi + rho1(a,b,cd)*rho1(i,j,cd)*eps/(eps**2 + eta**2) - end do - - do kl=1,nOO - eps = - Om2(kl) - chi = chi + rho2(a,b,kl)*rho2(i,j,kl)*eps/(eps**2 + eta**2) - end do - - TB(ab,ij) = lambda*chi + KB_sta(ab,ij) = lambda*chi end do end do diff --git a/src/GT/RGTpp_ppBSE_static_kernel_C.f90 b/src/GT/RGTpp_ppBSE_static_kernel_C.f90 index d6089ae..c8b36e4 100644 --- a/src/GT/RGTpp_ppBSE_static_kernel_C.f90 +++ b/src/GT/RGTpp_ppBSE_static_kernel_C.f90 @@ -1,4 +1,4 @@ -subroutine RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,nVVx,lambda,Om1,rho1,Om2,rho2,TC) +subroutine RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,lambda,eGF,Taaaa,Tabab,Tbaab,KC_sta) ! Compute the VVVV block of the static T-matrix @@ -16,25 +16,28 @@ subroutine RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n integer,intent(in) :: nR integer,intent(in) :: nOO integer,intent(in) :: nVV - integer,intent(in) :: nOOx - integer,intent(in) :: nVVx double precision,intent(in) :: lambda - double precision,intent(in) :: Om1(nVV) - double precision,intent(in) :: rho1(nBas,nBas,nVV) - double precision,intent(in) :: Om2(nOO) - double precision,intent(in) :: rho2(nBas,nBas,nOO) + double precision,intent(in) :: eGF(nBas) + double precision,intent(in) :: Taaaa(nBas,nBas,nBas,nBas) + double precision,intent(in) :: Tabab(nBas,nBas,nBas,nBas) + double precision,intent(in) :: Tbaab(nBas,nBas,nBas,nBas) ! Local variables double precision,external :: Kronecker_delta + double precision :: dem,num double precision :: chi double precision :: eps - integer :: a,b,c,d,ab,cd,ef,mn + integer :: a,b,c,d,ab,cd,ef,mn,m,e ! Output variables - double precision,intent(out) :: TC(nVVx,nVVx) + double precision,intent(out) :: KC_sta(nVV,nVV) +! Initialization + + KC_sta(:,:) = 0d0 + !===============! ! singlet block ! !===============! @@ -54,18 +57,16 @@ subroutine RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n chi = 0d0 do ef=1,nVV - eps = + Om1(ef) - chi = chi + rho1(a,b,ef)*rho1(c,d,ef)*eps/(eps**2 + eta**2) & - + rho1(a,b,ef)*rho1(d,c,ef)*eps/(eps**2 + eta**2) + eps = 0d0 + chi = chi + 0d0 end do do mn=1,nOO - eps = - Om2(mn) - chi = chi + rho2(a,b,mn)*rho2(c,d,mn)*eps/(eps**2 + eta**2) & - + rho2(a,b,mn)*rho2(d,c,mn)*eps/(eps**2 + eta**2) + eps = 0d0 + chi = chi + 0d0 end do - TC(ab,cd) = 0.5d0*lambda*chi/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) + KC_sta(ab,cd) = 0.5d0*lambda*chi/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) end do end do @@ -91,58 +92,15 @@ subroutine RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n do d=c+1,nBas-nR cd = cd + 1 - chi = 0d0 - - do ef=1,nVV - eps = + Om1(ef) - chi = chi + rho1(a,b,ef)*rho1(c,d,ef)*eps/(eps**2 + eta**2) + do m=nC+1,nO + do e=nO+1,nBas-nR + dem = eGF(m) - eGF(e) + num = 2d0*(Taaaa(a,m,c,e)*Taaaa(e,b,m,d) + Tabab(a,m,c,e)*Tabab(e,b,m,d)) + + KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2) + end do end do - do mn=1,nOO - eps = - Om2(mn) - chi = chi + rho2(a,b,mn)*rho2(c,d,mn)*eps/(eps**2 + eta**2) - end do - - TC(ab,cd) = lambda*chi - - end do - end do - - end do - end do - - end if - -!==================! -! alpha-beta block ! -!==================! - - if(ispin == 3) then - - ab = 0 - do a=nO+1,nBas-nR - do b=nO+1,nBas-nR - ab = ab + 1 - - cd = 0 - do c=nO+1,nBas-nR - do d=nO+1,nBas-nR - cd = cd + 1 - - chi = 0d0 - - do ef=1,nVV - eps = + Om1(ef) - chi = chi + rho1(a,b,ef)*rho1(c,d,ef)*eps/(eps**2 + eta**2) - end do - - do mn=1,nOO - eps = - Om2(mn) - chi = chi + rho2(a,b,mn)*rho2(c,d,mn)*eps/(eps**2 + eta**2) - end do - - TC(ab,cd) = lambda*chi - end do end do diff --git a/src/GT/RGTpp_ppBSE_static_kernel_D.f90 b/src/GT/RGTpp_ppBSE_static_kernel_D.f90 index 3a53998..d3c91f4 100644 --- a/src/GT/RGTpp_ppBSE_static_kernel_D.f90 +++ b/src/GT/RGTpp_ppBSE_static_kernel_D.f90 @@ -1,4 +1,4 @@ -subroutine RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,nVVx,lambda,Om1,rho1,Om2,rho2,TD) +subroutine RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,lambda,eGF,Taaaa,Tabab,Tbaab,KD_sta) ! Compute the OOOO block of the static T-matrix @@ -16,24 +16,26 @@ subroutine RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n integer,intent(in) :: nR integer,intent(in) :: nOO integer,intent(in) :: nVV - integer,intent(in) :: nOOx - integer,intent(in) :: nVVx double precision,intent(in) :: lambda - double precision,intent(in) :: Om1(nVV) - double precision,intent(in) :: rho1(nBas,nBas,nVV) - double precision,intent(in) :: Om2(nOO) - double precision,intent(in) :: rho2(nBas,nBas,nOO) - -! Local variables + double precision,intent(in) :: eGF(nBas) + double precision,intent(in) :: Taaaa(nBas,nBas,nBas,nBas) + double precision,intent(in) :: Tabab(nBas,nBas,nBas,nBas) + double precision,intent(in) :: Tbaab(nBas,nBas,nBas,nBas) + double precision,external :: Kronecker_delta + double precision :: dem,num double precision :: chi double precision :: eps - integer :: i,j,k,l,ij,kl,ef,mn + integer :: i,j,k,l,ij,kl,ef,mn,m,e ! Output variables - double precision,intent(out) :: TD(nOOx,nOOx) + double precision,intent(out) :: KD_sta(nOO,nOO) +! Initialization + + KD_sta(:,:) = 0d0 + !===============! ! singlet block ! !===============! @@ -50,20 +52,31 @@ subroutine RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n do l=k,nO kl = kl + 1 - chi = 0d0 - - do ef=1,nVV - eps = + Om1(ef) - chi = chi + rho1(i,j,ef)*rho1(k,l,ef)*eps/(eps**2 + eta**2) + do m=nC+1,nO + do e=nO+1,nBas-nR + dem = eGF(m) - eGF(e) + ! Wabab_{ijkl} + num = Taaaa(i,m,k,e)*Tabab(e,j,m,l) + Tabab(i,m,k,e)*Taaaa(e,j,m,l) + KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2) + + num = Taaaa(i,e,k,m)*Tabab(m,j,e,l) + Tabab(i,e,k,m)*Taaaa(m,j,e,l) + KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2) + + num = Taaaa(j,m,k,e)*Tabab(e,i,m,l) + Tabab(j,m,k,e)*Taaaa(e,i,m,l) + KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2) + + num = Taaaa(j,e,k,m)*Tabab(m,i,e,l) + Tabab(j,e,k,m)*Taaaa(m,i,e,l) + KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2) + + num = Tabab(i,m,k,e)*Tbaab(e,j,m,l) + Tbaab(i,e,k,m)*Tabab(m,j,e,l) + KD_sta(ij,kl) = KD_sta(ij,kl) - num*dem/(dem**2 + eta**2) + + num = Tabab(j,m,k,e)*Tbaab(e,i,m,l) + Tbaab(j,e,k,m)*Tabab(m,i,e,l) + KD_sta(ij,kl) = KD_sta(ij,kl) - num*dem/(dem**2 + eta**2) + end do end do - - do mn=1,nOO - eps = - Om2(mn) - chi = chi + rho2(i,j,mn)*rho2(k,l,mn)*eps/(eps**2 + eta**2) - end do - - TD(ij,kl) = lambda*chi - + + KD_sta(ij,kl) = KD_sta(ij,kl) / sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) end do end do @@ -88,20 +101,28 @@ subroutine RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n do l=k+1,nO kl = kl + 1 - chi = 0d0 - - do ef=1,nVV - eps = + Om1(ef) - chi = chi + rho1(i,j,ef)*rho1(k,l,ef)*eps/(eps**2 + eta**2) + do m=nC+1,nO + do e=nO+1,nBas-nR + dem = eGF(m) - eGF(e) + num = Taaaa(i,m,k,e)*Taaaa(e,j,m,l) + Tabab(i,m,k,e)*Tabab(e,j,m,l) + + KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2) + + num = Taaaa(i,e,k,m)*Taaaa(m,j,e,l) + Tabab(i,e,k,m)*Tabab(m,j,e,l) + + KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2) + + num = Taaaa(j,m,k,e)*Taaaa(e,i,m,l) + Tabab(j,m,k,e)*Tabab(e,i,m,l) + + KD_sta(ij,kl) = KD_sta(ij,kl) - num*dem/(dem**2 + eta**2) + + num = Taaaa(j,e,k,m)*Taaaa(m,i,e,l) + Tabab(j,e,k,m)*Tabab(m,i,e,l) + + KD_sta(ij,kl) = KD_sta(ij,kl) - num*dem/(dem**2 + eta**2) + + end do end do - - do mn=1,nOO - eps = - Om2(mn) - chi = chi + rho2(i,j,mn)*rho2(k,l,mn)*eps/(eps**2 + eta**2) - end do - - TD(ij,kl) = lambda*chi - + end do end do @@ -110,43 +131,4 @@ subroutine RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n end if -!==================! -! alpha-beta block ! -!==================! - - if(ispin == 3) then - - ij = 0 - do i=nC+1,nO - do j=nC+1,nO - ij = ij + 1 - - kl = 0 - do k=nC+1,nO - do l=nC+1,nO - kl = kl + 1 - - chi = 0d0 - - do ef=1,nVV - eps = + Om1(ef) - chi = chi + rho1(i,j,ef)*rho1(k,l,ef)*eps/(eps**2 + eta**2) - end do - - do mn=1,nOO - eps = - Om2(mn) - chi = chi + rho2(i,j,mn)*rho2(k,l,mn)*eps/(eps**2 + eta**2) - end do - - TD(ij,kl) = lambda*chi - - end do - end do - - end do - end do - - - end if - end subroutine diff --git a/src/GT/RGTpp_self_energy.f90 b/src/GT/RGTpp_self_energy.f90 index 12f5d1e..86dbc34 100644 --- a/src/GT/RGTpp_self_energy.f90 +++ b/src/GT/RGTpp_self_energy.f90 @@ -48,14 +48,14 @@ subroutine RGTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1s,rho do cd=1,nVVs eps = e(p) + e(i) - Om1s(cd) - num = rho1s(p,i,cd)*rho1s(q,i,cd) + num = (1d0/2d0)*rho1s(p,i,cd)*rho1s(q,i,cd) Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2) if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do do cd=1,nVVt eps = e(p) + e(i) - Om1t(cd) - num = rho1t(p,i,cd)*rho1t(q,i,cd) + num = (3d0/2d0)*rho1t(p,i,cd)*rho1t(q,i,cd) Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2) if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do @@ -74,14 +74,14 @@ subroutine RGTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1s,rho do kl=1,nOOs eps = e(p) + e(a) - Om2s(kl) - num = rho2s(p,a,kl)*rho2s(q,a,kl) + num = (1d0/2d0)*rho2s(p,a,kl)*rho2s(q,a,kl) Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2) if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do do kl=1,nOOt eps = e(p) + e(a) - Om2t(kl) - num = rho2t(p,a,kl)*rho2t(q,a,kl) + num = (3d0/2d0)*rho2t(p,a,kl)*rho2t(q,a,kl) Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2) if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do @@ -99,13 +99,13 @@ subroutine RGTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1s,rho do cd=1,nVVs eps = e(i) + e(j) - Om1s(cd) - num = rho1s(i,j,cd)*rho1s(i,j,cd) + num = (1d0/2d0)*rho1s(i,j,cd)*rho1s(i,j,cd) EcGM = EcGM + num*eps/(eps**2 + eta**2) end do do cd=1,nVVt eps = e(i) + e(j) - Om1t(cd) - num = rho1t(i,j,cd)*rho1t(i,j,cd) + num = (3d0/2d0)*rho1t(i,j,cd)*rho1t(i,j,cd) EcGM = EcGM + num*eps/(eps**2 + eta**2) end do @@ -117,13 +117,13 @@ subroutine RGTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1s,rho do kl=1,nOOs eps = e(a) + e(b) - Om2s(kl) - num = rho2s(a,b,kl)*rho2s(a,b,kl) + num = (1d0/2d0)*rho2s(a,b,kl)*rho2s(a,b,kl) EcGM = EcGM - num*eps/(eps**2 + eta**2) end do do kl=1,nOOt eps = e(a) + e(b) - Om2t(kl) - num = rho2t(a,b,kl)*rho2t(a,b,kl) + num = (3d0/2d0)*rho2t(a,b,kl)*rho2t(a,b,kl) EcGM = EcGM - num*eps/(eps**2 + eta**2) end do diff --git a/src/GT/RGTpp_self_energy_diag.f90 b/src/GT/RGTpp_self_energy_diag.f90 index d2c0c1c..5a1fcc8 100644 --- a/src/GT/RGTpp_self_energy_diag.f90 +++ b/src/GT/RGTpp_self_energy_diag.f90 @@ -48,14 +48,14 @@ subroutine RGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1 do cd=1,nVVs eps = e(p) + e(i) - Om1s(cd) - num = rho1s(p,i,cd)**2 + num = (1d0/2d0)*rho1s(p,i,cd)**2 Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do do cd=1,nVVt eps = e(p) + e(i) - Om1t(cd) - num = rho1t(p,i,cd)**2 + num = (3d0/2d0)*rho1t(p,i,cd)**2 Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do @@ -72,14 +72,14 @@ subroutine RGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1 do kl=1,nOOs eps = e(p) + e(a) - Om2s(kl) - num = rho2s(p,a,kl)**2 + num = (1d0/2d0)*rho2s(p,a,kl)**2 Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do do kl=1,nOOt eps = e(p) + e(a) - Om2t(kl) - num = rho2t(p,a,kl)**2 + num = (3d0/2d0)*rho2t(p,a,kl)**2 Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do @@ -96,13 +96,13 @@ subroutine RGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1 do cd=1,nVVs eps = e(i) + e(j) - Om1s(cd) - num = rho1s(i,j,cd)**2 + num = (1d0/2d0)*rho1s(i,j,cd)**2 EcGM = EcGM + num*eps/(eps**2 + eta**2) end do do cd=1,nVVt eps = e(i) + e(j) - Om1t(cd) - num = rho1t(i,j,cd)**2 + num = (3d0/2d0)*rho1t(i,j,cd)**2 EcGM = EcGM + num*eps/(eps**2 + eta**2) end do @@ -114,13 +114,13 @@ subroutine RGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1 do kl=1,nOOs eps = e(a) + e(b) - Om2s(kl) - num = rho2s(a,b,kl)**2 + num = (1d0/2d0)*rho2s(a,b,kl)**2 EcGM = EcGM - num*eps/(eps**2 + eta**2) end do do kl=1,nOOt eps = e(a) + e(b) - Om2t(kl) - num = rho2t(a,b,kl)**2 + num = (3d0/2d0)*rho2t(a,b,kl)**2 EcGM = EcGM - num*eps/(eps**2 + eta**2) end do diff --git a/src/GW/evRGW.f90 b/src/GW/evRGW.f90 index a2d07c2..3d89073 100644 --- a/src/GW/evRGW.f90 +++ b/src/GW/evRGW.f90 @@ -85,7 +85,7 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop ! SRG regularization - flow = 500d0 + flow = 100d0 if(doSRG) then diff --git a/src/GW/qsRGW.f90 b/src/GW/qsRGW.f90 index a8e4ce2..43f5d9f 100644 --- a/src/GW/qsRGW.f90 +++ b/src/GW/qsRGW.f90 @@ -314,7 +314,7 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop ! Deallocate memory deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,Z,Om,XpY,XmY,rho,err,err_diis,F_diis) - + ! Perform BSE calculation if(dophBSE) then diff --git a/src/LR/ppGLR.f90 b/src/LR/ppGLR.f90 index 5ec5dac..c47ff9c 100644 --- a/src/LR/ppGLR.f90 +++ b/src/LR/ppGLR.f90 @@ -49,7 +49,6 @@ subroutine ppGLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA) allocate(M(nPP,nPP),Z(nPP,nPP),Om(nPP)) ! Hermitian case for TDA - if(TDA) then X1(:,:) = +Cpp(:,:) @@ -72,7 +71,7 @@ subroutine ppGLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA) M( 1:nVV,nVV+1:nPP) = + Bpp(1:nVV,1:nOO) M(nVV+1:nPP, 1:nVV) = - transpose(Bpp(1:nVV,1:nOO)) -! if((nOO .eq. 0) .or. (nVV .eq. 0)) then + ! if((nOO .eq. 0) .or. (nVV .eq. 0)) then ! Diagonalize the pp matrix @@ -82,36 +81,36 @@ subroutine ppGLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA) call sort_ppRPA(nOO,nVV,nPP,Om,Z,Om1,X1,Y1,Om2,X2,Y2) -! else + ! else -! thr_d = 1d-6 ! to determine if diagonal elements of L.T x R are close enouph to 1 -! thr_nd = 1d-6 ! to determine if non-diagonal elements of L.T x R are close enouph to 1 -! thr_deg = 1d-8 ! to determine if two eigenvectors are degenerate or not -! imp_bio = .True. ! impose bi-orthogonality -! verbose = .False. -! call diagonalize_nonsym_matrix(N, M, Z, Om, thr_d, thr_nd, thr_deg, imp_bio, verbose) -! -! do i = 1, nOO -! Om2(i) = Om(i) -! do j = 1, nVV -! X2(j,i) = Z(j,i) -! enddo -! do j = 1, nOO -! Y2(j,i) = Z(nVV+j,i) -! enddo -! enddo -! -! do i = 1, nVV -! Om1(i) = Om(nOO+i) -! do j = 1, nVV -! X1(j,i) = M(j,nOO+i) -! enddo -! do j = 1, nOO -! Y1(j,i) = M(nVV+j,nOO+i) -! enddo -! enddo + ! thr_d = 1d-6 ! to determine if diagonal elements of L.T x R are close enouph to 1 + ! thr_nd = 1d-6 ! to determine if non-diagonal elements of L.T x R are close enouph to 1 + ! thr_deg = 1d-8 ! to determine if two eigenvectors are degenerate or not + ! imp_bio = .True. ! impose bi-orthogonality + ! verbose = .False. + ! call diagonalize_nonsym_matrix(Npp, M, Z, Om, thr_d, thr_nd, thr_deg, imp_bio, verbose) + + ! do i = 1, nOO + ! Om2(i) = Om(i) + ! do j = 1, nVV + ! X2(j,i) = Z(j,i) + ! enddo + ! do j = 1, nOO + ! Y2(j,i) = Z(nVV+j,i) + ! enddo + ! enddo + + ! do i = 1, nVV + ! Om1(i) = Om(nOO+i) + ! do j = 1, nVV + ! X1(j,i) = M(j,nOO+i) + ! enddo + ! do j = 1, nOO + ! Y1(j,i) = M(nVV+j,nOO+i) + ! enddo + ! enddo -! endif + ! endif end if diff --git a/src/LR/ppLR.f90 b/src/LR/ppLR.f90 index 362fcd7..4d417b0 100644 --- a/src/LR/ppLR.f90 +++ b/src/LR/ppLR.f90 @@ -82,36 +82,36 @@ subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA) call sort_ppRPA(nOO,nVV,nPP,Om,Z,Om1,X1,Y1,Om2,X2,Y2) -! else + ! else -! thr_d = 1d-6 ! to determine if diagonal elements of L.T x R are close enouph to 1 -! thr_nd = 1d-6 ! to determine if non-diagonal elements of L.T x R are close enouph to 1 -! thr_deg = 1d-8 ! to determine if two eigenvectors are degenerate or not -! imp_bio = .True. ! impose bi-orthogonality -! verbose = .False. -! call diagonalize_nonsym_matrix(N, M, Z, Om, thr_d, thr_nd, thr_deg, imp_bio, verbose) -! -! do i = 1, nOO -! Om2(i) = Om(i) -! do j = 1, nVV -! X2(j,i) = Z(j,i) -! enddo -! do j = 1, nOO -! Y2(j,i) = Z(nVV+j,i) -! enddo -! enddo -! -! do i = 1, nVV -! Om1(i) = Om(nOO+i) -! do j = 1, nVV -! X1(j,i) = M(j,nOO+i) -! enddo -! do j = 1, nOO -! Y1(j,i) = M(nVV+j,nOO+i) -! enddo -! enddo + ! thr_d = 1d-6 ! to determine if diagonal elements of L.T x R are close enouph to 1 + ! thr_nd = 1d-6 ! to determine if non-diagonal elements of L.T x R are close enouph to 1 + ! thr_deg = 1d-8 ! to determine if two eigenvectors are degenerate or not + ! imp_bio = .True. ! impose bi-orthogonality + ! verbose = .False. + ! call diagonalize_nonsym_matrix(Npp, M, Z, Om, thr_d, thr_nd, thr_deg, imp_bio, verbose) + + ! do i = 1, nOO + ! Om2(i) = Om(i) + ! do j = 1, nVV + ! X2(j,i) = Z(j,i) + ! enddo + ! do j = 1, nOO + ! Y2(j,i) = Z(nVV+j,i) + ! enddo + ! enddo + + ! do i = 1, nVV + ! Om1(i) = Om(nOO+i) + ! do j = 1, nVV + ! X1(j,i) = M(j,nOO+i) + ! enddo + ! do j = 1, nOO + ! Y1(j,i) = M(nVV+j,nOO+i) + ! enddo + ! enddo -! endif + ! endif end if diff --git a/src/QuAcK/GQuAcK.f90 b/src/QuAcK/GQuAcK.f90 index 52ca828..22441fd 100644 --- a/src/QuAcK/GQuAcK.f90 +++ b/src/QuAcK/GQuAcK.f90 @@ -1,11 +1,12 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & dodrCCD,dorCCD,docrCCD,dolCCD,dophRPA,dophRPAx,docrRPA,doppRPA, & - doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, & + doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2,doG0T0pp,doevGTpp,doqsGTpp, & nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO, & maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,reg_MP, & maxSCF_CC,max_diis_CC,thresh_CC, & TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, & - maxSCF_GW,max_diis_GW,thresh_GW,TDA_W,lin_GW,reg_GW,eta_GW, & + maxSCF_GW,max_diis_GW,thresh_GW,TDA_W,lin_GW,reg_GW,eta_GW, & + maxSCF_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, & dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS) implicit none @@ -23,6 +24,7 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do logical,intent(in) :: dophRPA,dophRPAx,docrRPA,doppRPA logical,intent(in) :: doG0F2,doevGF2,doqsGF2 logical,intent(in) :: doG0W0,doevGW,doqsGW + logical,intent(in) :: doG0T0pp,doevGTpp,doqsGTpp integer,intent(in) :: nNuc,nBas integer,intent(in) :: nC @@ -62,12 +64,17 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do logical,intent(in) :: TDA_W,lin_GW,reg_GW double precision,intent(in) :: eta_GW + integer,intent(in) :: maxSCF_GT,max_diis_GT + double precision,intent(in) :: thresh_GT + logical,intent(in) :: TDA_T,lin_GT,reg_GT + double precision,intent(in) :: eta_GT + logical,intent(in) :: dophBSE,dophBSE2,doppBSE,dBSE,dTDA logical,intent(in) :: doACFDT,exchange_kernel,doXBS ! Local variables - logical :: doMP,doCC,doRPA,doGF,doGW + logical :: doMP,doCC,doRPA,doGF,doGW,doGT double precision :: start_HF ,end_HF ,t_HF double precision :: start_stab ,end_stab ,t_stab @@ -77,6 +84,7 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do double precision :: start_RPA ,end_RPA ,t_RPA double precision :: start_GF ,end_GF ,t_GF double precision :: start_GW ,end_GW ,t_GW + double precision :: start_GT ,end_GT ,t_GT double precision,allocatable :: cHF(:,:),eHF(:),PHF(:,:),FHF(:,:) double precision :: EGHF @@ -288,6 +296,27 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GW = ',t_GW,' seconds' write(*,*) + end if + +!-----------------! +! T-matrix module ! +!-----------------! + + doGT = doG0T0pp .or. doevGTpp .or. doqsGTpp + + if(doGT) then + call wall_time(start_GT) + call GGT(dotest,doG0T0pp,doevGTpp,doqsGTpp, & + maxSCF_GT,thresh_GT,max_diis_GT,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE, & + TDA_T,TDA,dBSE,dTDA,lin_GT,eta_GT,reg_GT,nNuc,ZNuc,rNuc,ENuc, & + nBas,nBas2,nC,nO,nV,nR,nS,EGHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO, & + dipole_int_MO,PHF,cHF,eHF) + call wall_time(end_GT) + + t_GT = end_GT - start_GT + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GT = ',t_GT,' seconds' + write(*,*) + end if end subroutine diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 160ee7e..55a8b79 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -237,15 +237,15 @@ program QuAcK !--------------------------! ! Generalized QuAcK branch ! !--------------------------! - if(doGQuAcK) & call GQuAcK(doGtest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & dodrCCD,dorCCD,docrCCD,dolCCD,dophRPA,dophRPAx,docrRPA,doppRPA, & - doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, & + doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2,doG0T0pp,doevGTpp,doqsGTpp, & nNuc,nBas,sum(nC),sum(nO),sum(nV),sum(nR),ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO, & maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,reg_MP, & maxSCF_CC,max_diis_CC,thresh_CC,TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, & maxSCF_GW,max_diis_GW,thresh_GW,TDA_W,lin_GW,reg_GW,eta_GW, & + maxSCF_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, & dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS) !-----------!