diff --git a/src/GF/RG0F2.f90 b/src/GF/RG0F2.f90 index 8ed0a38..73e47f6 100644 --- a/src/GF/RG0F2.f90 +++ b/src/GF/RG0F2.f90 @@ -110,7 +110,7 @@ subroutine RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize, if(doppBSE) then - call RGF2_ppBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE) + call RGF2_ppBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE) EcBSE(2) = 3d0*EcBSE(2) diff --git a/src/GT/GG0T0pp.f90 b/src/GT/GG0T0pp.f90 index c0012d4..7be6092 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,23 @@ 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@GHF correlation energy = ',EcBSE,' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp@GHF 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..d38f584 --- /dev/null +++ b/src/GT/GGT_Tmatrix.f90 @@ -0,0 +1,66 @@ +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 :: 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,nOO,nVV,T,ERI,rho1,rho2,Om1,Om2) & + !$OMP PRIVATE(p,q,r,s,cd,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) + + do cd=1,nVV + T(p,q,r,s) = T(p,q,r,s) - rho1(p,q,cd)*rho1(r,s,cd)/Om1(cd) + end do + + do kl=1,nOO + T(p,q,r,s) = T(p,q,r,s) + rho2(p,q,kl)*rho2(r,s,kl)/Om2(kl) + end do + + 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..30f1d2b --- /dev/null +++ b/src/GT/GGTpp_ppBSE.f90 @@ -0,0 +1,124 @@ +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 + 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)) + + if(.not.TDA_T) call ppGLR_B(nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) + 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) + + 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)) + + if(.not.TDA) call ppGLR_B(nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) + 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) + +!---------------------------------------------- +! 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)) + + if(.not.TDA) call GGTpp_ppBSE_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGT,Om1,rho1,Om2,rho2,T,KB_sta) + 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) + + 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,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..848cbfd 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,21 @@ 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(e,b,m,j) + T(a,e,i,m) * T(m,b,e,j) + num = num - T(b,m,i,e) * T(e,a,m,j) - T(b,e,i,m) * T(m,a,e,j) + 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..c0c3866 100644 --- a/src/GT/GGTpp_ppBSE_static_kernel_C.f90 +++ b/src/GT/GGTpp_ppBSE_static_kernel_C.f90 @@ -1,67 +1,70 @@ -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(e,b,m,d) + T(a,e,c,m) * T(m,b,e,d) + num = num - T(b,m,c,e) * T(e,a,m,d) - T(b,e,c,m) * T(m,a,e,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 diff --git a/src/GT/GGTpp_ppBSE_static_kernel_D.f90 b/src/GT/GGTpp_ppBSE_static_kernel_D.f90 index b920a76..5052c6b 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,21 @@ 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 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) - do ef=1,nVV - eps = + Om1(ef) - chi = chi + rho1(i,j,ef)*rho1(k,l,ef)*eps/(eps**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/GGTpp_self_energy_diag.f90 b/src/GT/GGTpp_self_energy_diag.f90 index c34bb73..1f34638 100644 --- a/src/GT/GGTpp_self_energy_diag.f90 +++ b/src/GT/GGTpp_self_energy_diag.f90 @@ -81,7 +81,7 @@ subroutine GGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rh do cd=1,nVV eps = e(i) + e(j) - Om1(cd) - num = rho1(i,j,cd)**2 + num = 0.5d0*rho1(i,j,cd)**2 EcGM = EcGM + num*eps/(eps**2 + eta**2) end do @@ -93,7 +93,7 @@ subroutine GGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rh do kl=1,nOO eps = e(a) + e(b) - Om2(kl) - num = rho2(a,b,kl)**2 + num = 0.5d0*rho2(a,b,kl)**2 EcGM = EcGM - num*eps/(eps**2 + eta**2) end do 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/RGT_Tmatrix.f90 b/src/GT/RGT_Tmatrix.f90 new file mode 100644 index 0000000..6c34ce9 --- /dev/null +++ b/src/GT/RGT_Tmatrix.f90 @@ -0,0 +1,157 @@ +subroutine RGT_Tmatrix(isp_T,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,lambda,ERI,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t,T) + +! Compute the T-matrix tensor elements + + implicit none + include 'parameters.h' + + ! Input variables + + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nOOs, nOOt + integer,intent(in) :: nVVs, nVVt + integer,intent(in) :: isp_T + double precision,intent(in) :: lambda + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: Om1s(nVVs) + double precision,intent(in) :: rho1s(nBas,nBas,nVVs) + double precision,intent(in) :: Om2s(nOOs) + double precision,intent(in) :: rho2s(nBas,nBas,nOOs) + double precision,intent(in) :: Om1t(nVVt) + double precision,intent(in) :: rho1t(nBas,nBas,nVVt) + double precision,intent(in) :: Om2t(nOOt) + double precision,intent(in) :: rho2t(nBas,nBas,nOOt) + + ! Local variables + + double precision,external :: Kronecker_delta + integer :: p,q,r,s + integer :: c,d,k,l + integer :: kl,cd + + ! Output variables + + double precision,intent(out) :: T(nBas,nBas,nBas,nBas) + + ! Initialization + + T(:,:,:,:) = 0d0 + + ! Elements aaaa + + if(isp_T == 1) then + + !$OMP PARALLEL & + !$OMP SHARED(nC,nO,nBas,nR,T,ERI,nOOt,nVVt,rho1t,rho2t,Om1t,Om2t) & + !$OMP PRIVATE(p,q,r,s,cd,kl) & + !$OMP DEFAULT(NONE) + !$OMP DO + do s=nC+1,nBas-nR + do r=nC+1,nBas-nR + do q=nC+1,nBas-nR + do p=nC+1,nBas-nR + + T(p,q,r,s) = ERI(p,q,r,s) - ERI(p,q,s,r) + + do cd=1,nVVt + T(p,q,r,s) = T(p,q,r,s) - rho1t(p,q,cd) * rho1t(r,s,cd) / Om1t(cd) + end do ! cd + + do kl=1,nOOt + T(p,q,r,s) = T(p,q,r,s) + rho2t(p,q,kl) * rho2t(r,s,kl) / Om2t(kl) + enddo ! kl + + enddo ! p + enddo ! q + enddo ! r + enddo ! s + !$OMP END DO + !$OMP END PARALLEL + + endif + + ! Elements abab + + if(isp_T == 2) then + !$OMP PARALLEL & + !$OMP SHARED(nC,nO,nBas,nR,T,ERI,nOOs,nOOt,nVVs,nVVt,rho1s,rho2s,Om1s,Om2s,rho1t,rho2t,Om1t,Om2t) & + !$OMP PRIVATE(p,q,r,s,cd,kl) & + !$OMP DEFAULT(NONE) + !$OMP DO + do s=nC+1,nBas-nR + do r=nC+1,nBas-nR + do q=nC+1,nBas-nR + do p=nC+1,nBas-nR + + T(p,q,r,s) = ERI(p,q,r,s) + + do cd=1,nVVs + T(p,q,r,s) = T(p,q,r,s) - 0.5d0 * rho1s(p,q,cd) * rho1s(r,s,cd) / Om1s(cd) + end do ! cd + + do cd=1,nVVt + T(p,q,r,s) = T(p,q,r,s) - 0.5d0 * rho1t(p,q,cd) * rho1t(r,s,cd) / Om1t(cd) + end do ! cd + + do kl=1,nOOs + T(p,q,r,s) = T(p,q,r,s) + 0.5d0 * rho2s(p,q,kl) * rho2s(r,s,kl) / Om2s(kl) + enddo ! kl + + do kl=1,nOOt + T(p,q,r,s) = T(p,q,r,s) + 0.5d0 * rho2t(p,q,kl) * rho2t(r,s,kl) / Om2t(kl) + enddo ! kl + + enddo ! p + enddo ! q + enddo ! r + enddo ! s + !$OMP END DO + !$OMP END PARALLEL + + endif + + ! Elements baab + + if(isp_T == 3) then + !$OMP PARALLEL & + !$OMP SHARED(nC,nO,nBas,nR,T,ERI,nOOs,nOOt,nVVs,nVVt,rho1s,rho2s,Om1s,Om2s,rho1t,rho2t,Om1t,Om2t) & + !$OMP PRIVATE(p,q,r,s,cd,kl) & + !$OMP DEFAULT(NONE) + !$OMP DO + do s=nC+1,nBas-nR + do r=nC+1,nBas-nR + do q=nC+1,nBas-nR + do p=nC+1,nBas-nR + + T(p,q,r,s) = - ERI(p,q,s,r) + + do cd=1,nVVs + T(p,q,r,s) = T(p,q,r,s) + 0.5d0 * rho1s(p,q,cd) * rho1s(r,s,cd) / Om1s(cd) + end do ! cd + + do cd=1,nVVt + T(p,q,r,s) = T(p,q,r,s) - 0.5d0 * rho1t(p,q,cd) * rho1t(r,s,cd) / Om1t(cd) + end do ! cd + + do kl=1,nOOs + T(p,q,r,s) = T(p,q,r,s) - 0.5d0 * rho2s(p,q,kl) * rho2s(r,s,kl) / Om2s(kl) + enddo ! kl + + do kl=1,nOOt + T(p,q,r,s) = T(p,q,r,s) + 0.5d0 * rho2t(p,q,kl) * rho2t(r,s,kl) / Om2t(kl) + enddo ! kl + + enddo ! p + enddo ! q + enddo ! r + enddo ! s + !$OMP END DO + !$OMP END PARALLEL + + endif + +end subroutine 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..787b650 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,90 @@ 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) + ! Singlet contribution - deallocate(Om1s,X1s,Y1s,Om2s,X2s,Y2s) + isp_T = 1 -!------------------------------------! -! Compute T-matrix for triplet block ! -!------------------------------------! + allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) - ispin = 2 + 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,eT,ERI,Cpp) + call ppLR_D(isp_T,nBas,nC,nO,nV,nR,nOOs,1d0,eT,ERI,Dpp) - 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(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs)) + allocate(Om2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs)) - 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) + 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) + + allocate(rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs)) - call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin)) + call RGTpp_excitation_density(isp_T,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s) - deallocate(Bpp,Cpp,Dpp) - allocate(TBt(nVVt,nOOt),TCt(nVVt,nVVt),TDt(nOOt,nOOt)) + deallocate(X1s,Y1s,X2s,Y2s,Bpp,Cpp,Dpp) - 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) + ! Triplet contribution - deallocate(Om1t,X1t,Y1t,Om2t,X2t,Y2t) + isp_T = 2 + allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) + + 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,eT,ERI,Cpp) + call ppLR_D(isp_T,nBas,nC,nO,nV,nR,nOOt,1d0,eT,ERI,Dpp) + + allocate(Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt)) + allocate(Om2t(nOOt),X2t(nVVt,nOOt),Y2t(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) + + 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(X1t,Y1t,X2t,Y2t,Bpp,Cpp,Dpp) + +!--------------------------------- +! Compute T matrix elements +!--------------------------------- + + ! Elements aaaa + + 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) + + ! Elements abab + + 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) + + ! Elements baab + + 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 +137,32 @@ 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) 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 +173,30 @@ 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) 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..d96a236 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,24 +16,26 @@ 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 - double precision :: chi - double precision :: eps - integer :: i,j,a,b,ij,ab,cd,kl + double precision,external :: Kronecker_delta + double precision :: dem,num + integer :: i,j,a,b,ij,ab,cd,kl,m,e ! Output variables - double precision,intent(out) :: TB(nVVx,nOOx) + double precision,intent(out) :: KB_sta(nVV,nOO) +! Initialization + + KB_sta(:,:) = 0d0 + !===============! ! singlet block ! !===============! @@ -50,19 +52,31 @@ subroutine RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n do j=i,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,nBas-nR + dem = eGF(m) - eGF(e) + ! Wabab_{ijkl} + num = Taaaa(a,m,i,e)*Tabab(e,b,m,j) + Tabab(a,m,i,e)*Taaaa(e,b,m,j) + KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2) + + num = Taaaa(a,e,i,m)*Tabab(m,b,e,j) + Tabab(a,e,i,m)*Taaaa(m,b,e,j) + KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2) + + num = Taaaa(b,m,i,e)*Tabab(e,a,m,j) + Tabab(b,m,i,e)*Taaaa(e,a,m,j) + KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2) + + num = Taaaa(b,e,i,m)*Tabab(m,a,e,j) + Tabab(b,e,i,m)*Taaaa(m,a,e,j) + KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2) + + num = Tbaab(a,m,i,e)*Tbaab(e,b,m,j) + Tbaab(a,e,i,m)*Tbaab(m,b,e,j) + KB_sta(ab,ij) = KB_sta(ab,ij) - num*dem/(dem**2 + eta**2) + + num = Tbaab(b,m,i,e)*Tbaab(e,a,m,j) + Tbaab(b,e,i,m)*Tbaab(m,a,e,j) + 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 + KB_sta(ab,ij) = KB_sta(ab,ij) / sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j))) end do end do @@ -88,58 +102,25 @@ subroutine RGTpp_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,nBas-nR + dem = eGF(m) - eGF(e) + + num = Taaaa(a,m,i,e)*Taaaa(e,b,m,j) + Tabab(a,m,i,e)*Tabab(e,b,m,j) + KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2) + + num = Taaaa(a,e,i,m)*Taaaa(m,b,e,j) + Tabab(a,e,i,m)*Tabab(m,b,e,j) + KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2) + + num = Taaaa(b,m,i,e)*Taaaa(e,a,m,j) + Tabab(b,m,i,e)*Tabab(e,a,m,j) + KB_sta(ab,ij) = KB_sta(ab,ij) - num*dem/(dem**2 + eta**2) + + num = Taaaa(b,e,i,m)*Taaaa(m,a,e,j) + Tabab(b,e,i,m)*Tabab(m,a,e,j) + 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 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 - 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..a9ae81b 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,26 @@ 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 :: chi - double precision :: eps - integer :: a,b,c,d,ab,cd,ef,mn + double precision :: dem,num + 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 ! !===============! @@ -51,21 +52,31 @@ subroutine RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,n do d=c,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) & - + rho1(a,b,ef)*rho1(d,c,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(a,m,c,e)*Tabab(e,b,m,d) + Tabab(a,m,c,e)*Taaaa(e,b,m,d) + KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2) + + num = Taaaa(a,e,c,m)*Tabab(m,b,e,d) + Tabab(a,e,c,m)*Taaaa(m,b,e,d) + KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2) + + num = Taaaa(b,m,c,e)*Tabab(e,a,m,d) + Tabab(b,m,c,e)*Taaaa(e,a,m,d) + KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2) + + num = Taaaa(b,e,c,m)*Tabab(m,a,e,d) + Tabab(b,e,c,m)*Taaaa(m,a,e,d) + KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2) + + num = Tbaab(a,m,c,e)*Tbaab(e,b,m,d) + Tbaab(a,e,c,m)*Tbaab(m,b,e,d) + KC_sta(ab,cd) = KC_sta(ab,cd) - num*dem/(dem**2 + eta**2) + + num = Tbaab(b,m,c,e)*Tbaab(e,a,m,d) + Tbaab(b,e,c,m)*Tbaab(m,a,e,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) & - + rho2(a,b,mn)*rho2(d,c,mn)*eps/(eps**2 + eta**2) - end do - - TC(ab,cd) = 0.5d0*lambda*chi/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) + KC_sta(ab,cd) = KC_sta(ab,cd) / sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) end do end do @@ -91,58 +102,25 @@ 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 = 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) + + num = Taaaa(a,e,c,m)*Taaaa(m,b,e,d) + Tabab(a,e,c,m)*Tabab(m,b,e,d) + KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2) + + num = Taaaa(b,m,c,e)*Taaaa(e,a,m,d) + Tabab(b,m,c,e)*Tabab(e,a,m,d) + KC_sta(ab,cd) = KC_sta(ab,cd) - num*dem/(dem**2 + eta**2) + + num = Taaaa(b,e,c,m)*Taaaa(m,a,e,d) + Tabab(b,e,c,m)*Tabab(m,a,e,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..a3ea11a 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,24 @@ 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) + 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 :: chi - double precision :: eps - integer :: i,j,k,l,ij,kl,ef,mn + double precision,external :: Kronecker_delta + double precision :: dem,num + 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 +50,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 = Tbaab(i,m,k,e)*Tbaab(e,j,m,l) + Tbaab(i,e,k,m)*Tbaab(m,j,e,l) + KD_sta(ij,kl) = KD_sta(ij,kl) - num*dem/(dem**2 + eta**2) + + num = Tbaab(j,m,k,e)*Tbaab(e,i,m,l) + Tbaab(j,e,k,m)*Tbaab(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 +99,25 @@ 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 +126,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..6138996 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 = 0.5d0*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 = 1.5d0*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 = 0.5d0*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 = 1.5d0*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 = 0.5d0*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 = 1.5d0*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 = 0.5d0*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 = 1.5d0*rho2t(a,b,kl)**2 EcGM = EcGM - num*eps/(eps**2 + eta**2) end do diff --git a/src/GT/print_GG0T0pp.f90 b/src/GT/print_GG0T0pp.f90 index 468c4c3..02787e1 100644 --- a/src/GT/print_GG0T0pp.f90 +++ b/src/GT/print_GG0T0pp.f90 @@ -1,4 +1,4 @@ -subroutine print_GG0T0pp(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eGT,EcGM,EcRPA) +subroutine print_GG0T0pp(nBas,nO,eHF,ENuc,EGHF,SigT,Z,eGT,EcGM,EcRPA) ! Print one-electron energies and other stuff for G0T0 @@ -8,7 +8,7 @@ subroutine print_GG0T0pp(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eGT,EcGM,EcRPA) integer,intent(in) :: nBas integer,intent(in) :: nO double precision,intent(in) :: ENuc - double precision,intent(in) :: ERHF + double precision,intent(in) :: EGHF double precision,intent(in) :: EcGM double precision,intent(in) :: EcRPA double precision,intent(in) :: eHF(nBas) @@ -48,14 +48,14 @@ subroutine print_GG0T0pp(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eGT,EcGM,EcRPA) end do write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A60,F15.6,A3)') 'G0T0pp@RHF HOMO energy = ',eGT(HOMO)*HaToeV,' eV' - write(*,'(2X,A60,F15.6,A3)') 'G0T0pp@RHF LUMO energy = ',eGT(LUMO)*HaToeV,' eV' - write(*,'(2X,A60,F15.6,A3)') 'G0T0pp@RHF HOMO-LUMO gap = ',Gap*HaToeV,' eV' + write(*,'(2X,A60,F15.6,A3)') 'G0T0pp@GHF HOMO energy = ',eGT(HOMO)*HaToeV,' eV' + write(*,'(2X,A60,F15.6,A3)') 'G0T0pp@GHF LUMO energy = ',eGT(LUMO)*HaToeV,' eV' + write(*,'(2X,A60,F15.6,A3)') 'G0T0pp@GHF HOMO-LUMO gap = ',Gap*HaToeV,' eV' write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A60,F15.6,A3)') 'ppRPA@G0T0pp@RHF correlation energy = ',EcRPA,' au' - write(*,'(2X,A60,F15.6,A3)') 'ppRPA@G0T0pp@RHF total energy = ',ENuc + ERHF + EcRPA,' au' - write(*,'(2X,A60,F15.6,A3)') ' GM@G0T0pp@RHF correlation energy = ',EcGM,' au' - write(*,'(2X,A60,F15.6,A3)') ' GM@G0T0pp@RHF total energy = ',ENuc + ERHF + EcGM,' au' + write(*,'(2X,A60,F15.6,A3)') 'ppRPA@G0T0pp@GHF correlation energy = ',EcRPA,' au' + write(*,'(2X,A60,F15.6,A3)') 'ppRPA@G0T0pp@GHF total energy = ',ENuc + EGHF + EcRPA,' au' + write(*,'(2X,A60,F15.6,A3)') ' GM@G0T0pp@GHF correlation energy = ',EcGM,' au' + write(*,'(2X,A60,F15.6,A3)') ' GM@G0T0pp@GHF total energy = ',ENuc + EGHF + EcGM,' au' write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/GT/ufRG0T0pp.f90 b/src/GT/ufRG0T0pp.f90 index 1e398f9..2652837 100644 --- a/src/GT/ufRG0T0pp.f90 +++ b/src/GT/ufRG0T0pp.f90 @@ -33,7 +33,6 @@ subroutine ufRG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) logical :: print_T = .false. logical :: dRPA integer :: ispin - integer :: iblock integer :: nOOs,nOOt integer :: nVVs,nVVt double precision :: EcRPA(nspin) @@ -74,18 +73,12 @@ subroutine ufRG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Dimensions of the ppRPA linear reponse matrices -! nOOs = nO*(nO + 1)/2 -! nVVs = nV*(nV + 1)/2 - - nOOs = nO*nO - nVVs = nV*nV + nOOs = nO*(nO + 1)/2 + nVVs = nV*(nV + 1)/2 nOOt = nO*(nO - 1)/2 nVVt = nV*(nV - 1)/2 -! nOO = nO*nO -! nVV = nV*nV - ! Dimension of the supermatrix n2h1p = (nOOs+nOOt)*nV @@ -117,16 +110,14 @@ subroutine ufRG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! alpha-beta block ispin = 1 -! iblock = 1 - iblock = 3 ! Compute linear response allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) - call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp) - call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,1d0,eHF,ERI,Cpp) - call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,1d0,eHF,ERI,Dpp) + 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,eHF,ERI,Cpp) + call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOOs,1d0,eHF,ERI,Dpp) call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin)) @@ -135,23 +126,21 @@ subroutine ufRG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Compute excitation densities - call RGTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s) + call RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s) deallocate(Bpp,Cpp,Dpp,X1s,Y1s,X2s,Y2s) ! alpha-alpha block ispin = 2 -! iblock = 2 - iblock = 4 ! Compute linear response allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) - call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp) - call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,1d0,eHF,ERI,Cpp) - call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,1d0,eHF,ERI,Dpp) + 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,eHF,ERI,Cpp) + call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOOt,1d0,eHF,ERI,Dpp) call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin)) @@ -160,10 +149,14 @@ subroutine ufRG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Compute excitation densities - call RGTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t) + call RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t) deallocate(Bpp,Cpp,Dpp,X1t,Y1t,X2t,Y2t) + else + + allocate(rho1s(0,0,0),rho1t(0,0,0),rho2s(0,0,0),rho2t(0,0,0)) + end if ! Memory allocation diff --git a/src/GW/GGW_ppBSE_dynamic_kernel_B.f90 b/src/GW/GGW_ppBSE_dynamic_kernel_B.f90 index 52dc5b4..dbdf24a 100644 --- a/src/GW/GGW_ppBSE_dynamic_kernel_B.f90 +++ b/src/GW/GGW_ppBSE_dynamic_kernel_B.f90 @@ -50,28 +50,23 @@ subroutine GGW_ppBSE_dynamic_kernel_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGW ij = ij + 1 do m=1,nS + + num = (rho(a,i,m)*rho(b,j,m) - rho(b,i,m)*rho(a,j,m))/2 + dem = - Om(m) - eGW(b) + eGW(j) + KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2) - dem = eGW(j) - Om(m) - eGW(b) - num = rho(a,i,m)*rho(b,j,m) + dem = - Om(m) - eGW(a) + eGW(i) + KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2) - KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2) + dem = - Om(m) - eGW(a) + eGW(j) + KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2) + + dem = - Om(m) - eGW(b) + eGW(i) + KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2) - dem = eGW(j) - Om(m) - eGW(a) - num = rho(b,i,m)*rho(a,j,m) + end do - KB_dyn(ab,ij) = KB_dyn(ab,ij) - num*dem/(dem**2 + eta**2) - - dem = eGW(i) - Om(m) - eGW(a) - num = rho(a,i,m)*rho(b,j,m) - - KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2) - - dem = eGW(i) - Om(m) - eGW(b) - num = rho(b,i,m)*rho(a,j,m) - - KB_dyn(ab,ij) = KB_dyn(ab,ij) - num*dem/(dem**2 + eta**2) - - end do + KB_dyn(ab,ij) = 0.5d0*KB_dyn(ab,ij) end do end do diff --git a/src/GW/GGW_ppBSE_dynamic_kernel_C.f90 b/src/GW/GGW_ppBSE_dynamic_kernel_C.f90 index 75c0050..c61a280 100644 --- a/src/GW/GGW_ppBSE_dynamic_kernel_C.f90 +++ b/src/GW/GGW_ppBSE_dynamic_kernel_C.f90 @@ -52,36 +52,28 @@ subroutine GGW_ppBSE_dynamic_kernel_C(eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,eGW,Om, cd = cd + 1 do m=1,nS + + num = (rho(a,c,m)*rho(b,d,m) - rho(b,c,m)*rho(a,d,m))/2 + dem = OmBSE - Om(m) - eGW(b) - eGW(d) + KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2) + ZC_dyn(ab,cd) = ZC_dyn(ab,cd) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 - dem = OmBSE - (eGW(a) + eGW(c) + Om(m)) -! num = 0.5d0*(rho(a,c,m)*rho(b,d,m) - rho(b,c,m)*rho(a,d,m)) - num = - rho(b,c,m)*rho(a,d,m) + dem = OmBSE - Om(m) - eGW(a) - eGW(c) + KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2) + ZC_dyn(ab,cd) = ZC_dyn(ab,cd) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 - KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2) - ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 + dem = OmBSE - Om(m) - eGW(a) - eGW(d) + KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2) + ZC_dyn(ab,cd) = ZC_dyn(ab,cd) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 + + dem = OmBSE - Om(m) - eGW(b) - eGW(c) + KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2) + ZC_dyn(ab,cd) = ZC_dyn(ab,cd) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 - dem = OmBSE - (eGW(b) + eGW(d) + Om(m)) -! num = 0.5d0*(rho(a,c,m)*rho(b,d,m) - rho(b,c,m)*rho(a,d,m)) - num = - rho(b,c,m)*rho(a,d,m) + end do - KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2) - ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 - - dem = OmBSE - (eGW(b) + eGW(c) + Om(m)) - num = rho(a,c,m)*rho(b,d,m) -! num = 0.5d0*(rho(a,c,m)*rho(b,d,m) - rho(b,c,m)*rho(a,d,m)) - - KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2) - ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 - - dem = OmBSE - (eGW(a) + eGW(d) + Om(m)) - num = rho(a,c,m)*rho(b,d,m) -! num = 0.5d0*(rho(a,c,m)*rho(b,d,m) - rho(b,c,m)*rho(a,d,m)) - - KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2) - ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 - - end do + KC_dyn(ab,cd) = 0.5d0*KC_dyn(ab,cd) + ZC_dyn(ab,cd) = 0.5d0*ZC_dyn(ab,cd) end do end do diff --git a/src/GW/RGW_phACFDT.f90 b/src/GW/RGW_phACFDT.f90 index e2b2254..154fa6f 100644 --- a/src/GW/RGW_phACFDT.f90 +++ b/src/GW/RGW_phACFDT.f90 @@ -99,8 +99,8 @@ subroutine RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nBas, call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA) call RGW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) - call RGW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA) - call RGW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB) + call RGW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA) + if(.not.TDA) call RGW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB) ! Singlet manifold @@ -129,8 +129,8 @@ subroutine RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nBas, call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA) call RGW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) - call RGW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KA) - call RGW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KB) + call RGW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KA) + if(.not.TDA) call RGW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KB) end if @@ -186,8 +186,8 @@ subroutine RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nBas, call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA) call RGW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) - call RGW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KA) - call RGW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KB) + call RGW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KA) + if(.not.TDA) call RGW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KB) end if diff --git a/src/GW/RGW_phBSE.f90 b/src/GW/RGW_phBSE.f90 index 32ddaae..4e1a926 100644 --- a/src/GW/RGW_phBSE.f90 +++ b/src/GW/RGW_phBSE.f90 @@ -89,8 +89,8 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA) call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) - call RGW_phBSE_static_kernel_A(eta,nOrb,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA_sta) - call RGW_phBSE_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB_sta) + call RGW_phBSE_static_kernel_A(eta,nOrb,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA_sta) + if(.not.TDA) call RGW_phBSE_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB_sta) !------------------- ! Singlet manifold @@ -116,8 +116,8 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple write(*,*) call RGW_phBSE_static_kernel(eta,nOrb,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,W) - call RGW_phBSE2_static_kernel_A(eta,nOrb,nC,nO,nV,nR,nS,1d0,eW,W,KA_sta) + call RGW_phBSE2_static_kernel_A(eta,nOrb,nC,nO,nV,nR,nS,1d0,eW,W,KA_sta) if(.not.TDA) call RGW_phBSE2_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nS,1d0,eW,W,KB_sta) deallocate(W) diff --git a/src/GW/RGW_ppBSE.f90 b/src/GW/RGW_ppBSE.f90 index baf1d8a..5e0f77c 100644 --- a/src/GW/RGW_ppBSE.f90 +++ b/src/GW/RGW_ppBSE.f90 @@ -137,55 +137,23 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS allocate(Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO)) allocate(KB_sta(nVV,nOO),KC_sta(nVV,nVV),KD_sta(nOO,nOO)) - !print*, 'RGW_ppBSE_static_kernel_C:' - !call wall_time(tt1) call RGW_ppBSE_static_kernel_C(ispin,eta,nOrb,nC,nO,nV,nR,nS,nVV,1d0,ERI,OmRPA,rho_RPA,KC_sta) - !call wall_time(tt2) - !write(*,'(A,1X,F10.3)'), 'wall time for RGW_ppBSE_static_kernel_C (sec)', tt2-tt1 - - !print*, 'RGW_ppBSE_static_kernel_D:' - !call wall_time(tt1) call RGW_ppBSE_static_kernel_D(ispin,eta,nOrb,nC,nO,nV,nR,nS,nOO,1d0,ERI,OmRPA,rho_RPA,KD_sta) - !call wall_time(tt2) - !write(*,'(A,1X,F10.3)'), 'wall time for RGW_ppBSE_static_kernel_D (sec)', tt2-tt1 - if(.not.TDA) then - !print*, 'RGW_ppBSE_static_kernel_B:' - !call wall_time(tt1) call RGW_ppBSE_static_kernel_B(ispin,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KB_sta) - !call wall_time(tt2) - !write(*,'(A,1X,F10.3)'), 'wall time for RGW_ppBSE_static_kernel_B (sec)', tt2-tt1 endif - !print*, 'ppLR_C:' - !call wall_time(tt1) call ppLR_C(ispin,nOrb,nC,nO,nV,nR,nVV,1d0,eGW,ERI,Cpp) - !call wall_time(tt2) - !write(*,'(A,1X,F10.3)'), 'wall time for ppLR_C (sec)', tt2-tt1 - - !print*, 'ppLR_D:' - !call wall_time(tt1) call ppLR_D(ispin,nOrb,nC,nO,nV,nR,nOO,1d0,eGW,ERI,Dpp) - !call wall_time(tt2) - !write(*,'(A,1X,F10.3)'), 'wall time for ppLR_D (sec)', tt2-tt1 - if(.not.TDA) then - !print*, 'ppLR_B:' - !call wall_time(tt1) call ppLR_B(ispin,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) - !call wall_time(tt2) - !write(*,'(A,1X,F10.3)'), 'wall time for ppLR_B (sec)', tt2-tt1 endif Bpp(:,:) = Bpp(:,:) + KB_sta(:,:) Cpp(:,:) = Cpp(:,:) + KC_sta(:,:) Dpp(:,:) = Dpp(:,:) + KD_sta(:,:) - !print*, 'ppLR:' - !call wall_time(tt1) call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin)) - !call wall_time(tt2) - !write(*,'(A,1X,F10.3)'), 'wall time for ppLR (sec)', tt2-tt1 deallocate(Bpp,Cpp,Dpp) ! @@ -197,50 +165,50 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS - ! --- - ! Davidson - ! --- + !! --- + !! Davidson + !! --- -! n_states = nOO + 5 -! n_states_diag = n_states + 4 -! allocate(Om(nOO+nVV), R(nOO+nVV,n_states_diag)) -! -! supp_data_int = 1 -! allocate(supp_data_int(supp_data_int_size)) -! supp_data_int(1) = nS -! -! supp_data_dbl_size = nS + nOrb*nOrb*nS + 1 -! allocate(supp_data_dbl(supp_data_dbl_size)) -! ! scalars -! supp_data_dbl(1) = eta -! i_data = 1 -! ! rho_RPA -! do q = 1, nOrb -! do p = 1, nOrb -! do m = 1, nS -! i_data = i_data + 1 -! supp_data_dbl(i_data) = rho_RPA(p,q,m) -! enddo -! enddo -! enddo -! ! OmRPA -! do m = 1, nS -! i_data = i_data + 1 -! supp_data_dbl(i_data) = OmRPA(m) -! enddo -! -! call ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, & -! 1.d0, & ! lambda -! eGW(1), & -! 0.d0, & ! eF -! ERI(1,1,1,1), & -! supp_data_int(1), supp_data_int_size, & -! supp_data_dbl(1), supp_data_dbl_size, & -! Om(1), R(1,1), n_states, n_states_diag, "GW") -! -! deallocate(Om, R) -! deallocate(supp_data_dbl, supp_data_int) -! stop + !n_states = nOO + 5 + !n_states_diag = n_states + 4 + !allocate(Om(nOO+nVV), R(nOO+nVV,n_states_diag)) + + !supp_data_int_size = 1 + !allocate(supp_data_int(supp_data_int_size)) + !supp_data_int(1) = nS + + !supp_data_dbl_size = nS + nOrb*nOrb*nS + 1 + !allocate(supp_data_dbl(supp_data_dbl_size)) + !! scalars + !supp_data_dbl(1) = eta + !i_data = 1 + !! rho_RPA + !do q = 1, nOrb + ! do p = 1, nOrb + ! do m = 1, nS + ! i_data = i_data + 1 + ! supp_data_dbl(i_data) = rho_RPA(p,q,m) + ! enddo + ! enddo + !enddo + !! OmRPA + !do m = 1, nS + ! i_data = i_data + 1 + ! supp_data_dbl(i_data) = OmRPA(m) + !enddo + + !call ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, & + ! 1.d0, & ! lambda + ! eGW(1), & + ! 0.d0, & ! eF + ! ERI(1,1,1,1), & + ! supp_data_int(1), supp_data_int_size, & + ! supp_data_dbl(1), supp_data_dbl_size, & + ! Om(1), R(1,1), n_states, n_states_diag, "GW", 1) + + !deallocate(Om, R) + !deallocate(supp_data_dbl, supp_data_int) + !stop ! --- @@ -308,10 +276,6 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin)) deallocate(Bpp,Cpp,Dpp) - !print*, 'LAPACK:' - !print*, Om2 - !print*, Om1 - ! --- ! Davidson ! --- @@ -320,7 +284,7 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS !n_states_diag = n_states + 4 !allocate(Om(nOO+nVV), R(nOO+nVV,n_states_diag)) - !supp_data_int = 1 + !supp_data_int_size = 1 !allocate(supp_data_int(supp_data_int_size)) !supp_data_int(1) = nS @@ -351,7 +315,7 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS ! ERI(1,1,1,1), & ! supp_data_int(1), supp_data_int_size, & ! supp_data_dbl(1), supp_data_dbl_size, & - ! Om(1), R(1,1), n_states, n_states_diag, "GW") + ! Om(1), R(1,1), n_states, n_states_diag, "GW", 1) !deallocate(Om, R) !deallocate(supp_data_dbl, supp_data_int) @@ -371,15 +335,9 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS call RGW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, & Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta) - - !----------------! - ! Upfolded ppBSE ! - !----------------! - - ! call RGW_ppBSE_upfolded(ispin,nOrb,nC,nO,nV,nR,nS,ERI,rho_RPA,OmRPA,eGW) - deallocate(KB_sta,KC_sta,KD_sta) deallocate(Om1,X1,Y1,Om2,X2,Y2) + end if end subroutine diff --git a/src/GW/RGW_ppBSE_dynamic_kernel_B.f90 b/src/GW/RGW_ppBSE_dynamic_kernel_B.f90 index 838b017..87c42da 100644 --- a/src/GW/RGW_ppBSE_dynamic_kernel_B.f90 +++ b/src/GW/RGW_ppBSE_dynamic_kernel_B.f90 @@ -53,30 +53,23 @@ subroutine RGW_ppBSE_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lamb ij = ij + 1 do m=1,nS + num = (rho(a,i,m)*rho(b,j,m) + rho(b,i,m)*rho(a,j,m))/2 - dem = eGW(j) - Om(m) - eGW(b) - num = rho(a,i,m)*rho(b,j,m) + dem = - Om(m) - eGW(b) + eGW(j) + KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2) - KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2) + dem = - Om(m) - eGW(a) + eGW(i) + KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2) - dem = eGW(j) - Om(m) - eGW(a) - num = rho(b,i,m)*rho(a,j,m) + dem = - Om(m) - eGW(a) + eGW(j) + KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2) - KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2) - - dem = eGW(i) - Om(m) - eGW(a) - num = rho(a,i,m)*rho(b,j,m) - - KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2) - - dem = eGW(i) - Om(m) - eGW(b) - num = rho(b,i,m)*rho(a,j,m) - - KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2) + dem = - Om(m) - eGW(b) + eGW(i) + KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2) end do - KB_dyn(ab,ij) = 2d0*KB_dyn(ab,ij)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j))) + KB_dyn(ab,ij) = KB_dyn(ab,ij)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j))) end do end do @@ -99,31 +92,22 @@ subroutine RGW_ppBSE_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lamb ij = ij + 1 do m=1,nS + num = (rho(a,i,m)*rho(b,j,m) - rho(b,i,m)*rho(a,j,m))/2 - dem = eGW(j) - Om(m) - eGW(b) - num = rho(a,i,m)*rho(b,j,m) + dem = - Om(m) - eGW(b) + eGW(j) + KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2) - KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2) + dem = - Om(m) - eGW(a) + eGW(i) + KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2) - dem = eGW(j) - Om(m) - eGW(a) - num = rho(b,i,m)*rho(a,j,m) + dem = - Om(m) - eGW(a) + eGW(j) + KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2) - KB_dyn(ab,ij) = KB_dyn(ab,ij) - num*dem/(dem**2 + eta**2) - - dem = eGW(i) - Om(m) - eGW(a) - num = rho(a,i,m)*rho(b,j,m) - - KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2) - - dem = eGW(i) - Om(m) - eGW(b) - num = rho(b,i,m)*rho(a,j,m) - - KB_dyn(ab,ij) = KB_dyn(ab,ij) - num*dem/(dem**2 + eta**2) + dem = - Om(m) - eGW(b) + eGW(i) + KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2) end do - KB_dyn(ab,ij) = 2d0*KB_dyn(ab,ij) - end do end do diff --git a/src/GW/RGW_ppBSE_dynamic_kernel_C.f90 b/src/GW/RGW_ppBSE_dynamic_kernel_C.f90 index 04ce120..a90217b 100644 --- a/src/GW/RGW_ppBSE_dynamic_kernel_C.f90 +++ b/src/GW/RGW_ppBSE_dynamic_kernel_C.f90 @@ -55,35 +55,28 @@ subroutine RGW_ppBSE_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,e cd = cd + 1 do m=1,nS + num = (rho(a,c,m)*rho(b,d,m) + rho(b,c,m)*rho(a,d,m))/2 - dem = OmBSE - eGW(c) - Om(m) - eGW(b) - num = rho(a,c,m)*rho(b,d,m) + dem = OmBSE - Om(m) - eGW(b) - eGW(d) + KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2) + ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 - KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2) - ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 + dem = OmBSE - Om(m) - eGW(a) - eGW(c) + KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2) + ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 - dem = OmBSE - eGW(c) - Om(m) - eGW(a) - num = rho(b,c,m)*rho(a,d,m) + dem = OmBSE - Om(m) - eGW(a) - eGW(d) + KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2) + ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 - KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2) - ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 - - dem = OmBSE - eGW(d) - Om(m) - eGW(a) - num = rho(a,c,m)*rho(b,d,m) - - KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2) - ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 - - dem = OmBSE - eGW(d) - Om(m) - eGW(b) - num = rho(b,c,m)*rho(a,d,m) - - KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2) - ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 + dem = OmBSE - Om(m) - eGW(b) - eGW(c) + KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2) + ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 end do - KC_dyn(ab,cd) = 2d0*KC_dyn(ab,cd)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) - ZC_dyn(ab,cd) = 2d0*ZC_dyn(ab,cd)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) + KC_dyn(ab,cd) = KC_dyn(ab,cd)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) + ZC_dyn(ab,cd) = ZC_dyn(ab,cd)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) end do end do @@ -106,35 +99,25 @@ subroutine RGW_ppBSE_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,e cd = cd + 1 do m=1,nS + num = (rho(a,c,m)*rho(b,d,m) - rho(b,c,m)*rho(a,d,m))/2 - dem = OmBSE - eGW(c) - Om(m) - eGW(b) - num = rho(a,c,m)*rho(b,d,m) + dem = OmBSE - Om(m) - eGW(b) - eGW(d) + KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2) + ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 - KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2) - ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 + dem = OmBSE - Om(m) - eGW(a) - eGW(c) + KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2) + ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 - dem = OmBSE - eGW(c) - Om(m) - eGW(a) - num = rho(b,c,m)*rho(a,d,m) + dem = OmBSE - Om(m) - eGW(a) - eGW(d) + KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2) + ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 - KC_dyn(ab,cd) = KC_dyn(ab,cd) - num*dem/(dem**2 + eta**2) - ZC_dyn(ab,cd) = ZC_dyn(ab,cd) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 - - dem = OmBSE - eGW(d) - Om(m) - eGW(a) - num = rho(a,c,m)*rho(b,d,m) - - KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2) - ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 - - dem = OmBSE - eGW(d) - Om(m) - eGW(b) - num = rho(b,c,m)*rho(a,d,m) - - KC_dyn(ab,cd) = KC_dyn(ab,cd) - num*dem/(dem**2 + eta**2) - ZC_dyn(ab,cd) = ZC_dyn(ab,cd) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 + dem = OmBSE - Om(m) - eGW(b) - eGW(c) + KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2) + ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 end do - - KC_dyn(ab,cd) = 2d0*KC_dyn(ab,cd) - ZC_dyn(ab,cd) = 2d0*ZC_dyn(ab,cd) end do end do diff --git a/src/GW/RGW_ppBSE_dynamic_kernel_D.f90 b/src/GW/RGW_ppBSE_dynamic_kernel_D.f90 index 6e5e332..e8d4ecf 100644 --- a/src/GW/RGW_ppBSE_dynamic_kernel_D.f90 +++ b/src/GW/RGW_ppBSE_dynamic_kernel_D.f90 @@ -57,7 +57,6 @@ subroutine RGW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,e do m=1,nS num = (rho(i,k,m)*rho(j,l,m) + rho(j,k,m)*rho(i,l,m))/2 -! dem = - Om(m) dem = - OmBSE - Om(m) + eGW(j) + eGW(l) KD_dyn(ij,kl) = KD_dyn(ij,kl) + num*dem/(dem**2 + eta**2) ZD_dyn(ij,kl) = ZD_dyn(ij,kl) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 @@ -101,6 +100,7 @@ subroutine RGW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,e do m=1,nS num = (rho(i,k,m)*rho(j,l,m) - rho(j,k,m)*rho(i,l,m))/2 + dem = - OmBSE - Om(m) + eGW(j) + eGW(l) KD_dyn(ij,kl) = KD_dyn(ij,kl) + num*dem/(dem**2 + eta**2) ZD_dyn(ij,kl) = ZD_dyn(ij,kl) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 @@ -119,9 +119,6 @@ subroutine RGW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,e end do - KD_dyn(ij,kl) = KD_dyn(ij,kl) - ZD_dyn(ij,kl) = ZD_dyn(ij,kl) - end do end do diff --git a/src/GW/RGW_ppBSE_upfolded.f90 b/src/GW/RGW_ppBSE_upfolded.f90 deleted file mode 100644 index d5e8deb..0000000 --- a/src/GW/RGW_ppBSE_upfolded.f90 +++ /dev/null @@ -1,263 +0,0 @@ -subroutine RGW_ppBSE_upfolded(ispin,nOrb,nC,nO,nV,nR,nS,ERI,rho,Om,eGW) - -! Upfolded ppBSE@GW (TDA only!) - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: ispin - integer,intent(in) :: nOrb - integer,intent(in) :: nC - integer,intent(in) :: nO - integer,intent(in) :: nV - integer,intent(in) :: nR - integer,intent(in) :: nS - double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) - double precision,intent(in) :: rho(nOrb,nOrb,nS) - double precision,intent(in) :: Om(nS) - double precision,intent(in) :: eGW(nOrb) - -! Local variables - - integer :: s - integer :: i,j,k,l - integer :: a,b,c,d - integer :: m,ij,kl,ijm - integer,parameter :: maxH = 100 - double precision :: tmp,tmp1,tmp2,tmp3,tmp4 - - integer :: n1h,n1p,n2h,n2p,n1h1p,n3h1p,n3p1h,n2h2p,nH - double precision,external :: Kronecker_delta - - integer,allocatable :: order(:) - double precision,allocatable :: H(:,:) - double precision,allocatable :: X(:,:) - double precision,allocatable :: OmBSE(:) - double precision,allocatable :: Z(:) - -! Output variables - -! Hello world - - write(*,*) - write(*,*)'*********************************' - write(*,*)'* Upfolded ppBSE@GW Calculation *' - write(*,*)'*********************************' - write(*,*) - -! TDA for W - - write(*,*) 'Tamm-Dancoff approximation by default!' - write(*,*) - -! Dimension of the supermatrix - - n1h = nO - n1p = nV - - if(ispin == 1) then - - n2h = nO*(nO+1)/2 - n2p = nV*(nV+1)/2 - - end if - - if(ispin == 2) then - - n2h = nO*(nO-1)/2 - n2p = nV*(nV-1)/2 - - end if - - n1h1p = n1h*n1p - - n3h1p = n2h*n1h1p - n3p1h = n2p*n1h1p - - nH = n2h + 4*n3h1p - -! Memory allocation - - allocate(order(nH),H(nH,nH),X(nH,nH),OmBSE(nH),Z(nH)) - -! Initialization - - H(:,:) = 0d0 - -!----------------------------------------! -! Compute BSE supermatrix ! -!----------------------------------------! -! ! -! | D -M1 -M1 -M2 -M2 | ! -! | | ! -! | +M1 E1 0 0 0 | ! -! | | ! -! H = | +M1 0 E2 0 0 | ! -! | | ! -! | +M2 0 0 E3 0 | ! -! | | ! -! | +M2 0 0 0 E4 | ! -! ! -!----------------------------------------! - - !----------------------! - ! Block D for singlets ! - !----------------------! - - if(ispin == 1) then - - ij = 0 - do i=nC+1,nO - do j=i,nO - ij = ij + 1 - - kl = 0 - do k=nC+1,nO - do l=k,nO - kl = kl + 1 - - H(ij,kl) = - (eGW(i) + eGW(j))*Kronecker_delta(i,k)*Kronecker_delta(j,l) & - + (ERI(i,j,k,l) + ERI(i,j,l,k))/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) - - end do - end do - - end do - end do - - end if - - !----------------------! - ! Block D for triplets ! - !----------------------! - - if(ispin == 2) then - - ij = 0 - do i=nC+1,nO - do j=i+1,nO - ij = ij + 1 - - kl = 0 - do k=nC+1,nO - do l=k+1,nO - kl = kl + 1 - - H(ij,kl) = - (eGW(i) + eGW(j))*Kronecker_delta(i,k)*Kronecker_delta(j,l) & - + (ERI(i,j,k,l) - ERI(i,j,l,k)) - - end do - end do - - end do - end do - - end if - - !----------------! - ! Blocks M1 & M2 ! - !----------------! - - ijm = 0 - do i=nC+1,nO - do j=i+1,nO - do m=1,nS - ijm = ijm + 1 - - kl = 0 - do k=nC+1,nO - do l=k+1,nO - kl = kl + 1 - - tmp1 = sqrt(2d0)*Kronecker_delta(j,l)*rho(i,k,m) - tmp2 = sqrt(2d0)*Kronecker_delta(j,k)*rho(i,l,m) - tmp3 = sqrt(2d0)*Kronecker_delta(i,l)*rho(j,k,m) - tmp4 = sqrt(2d0)*Kronecker_delta(i,k)*rho(j,l,m) - - H(n2h+0*n3h1p+ijm,kl ) = +tmp1 - H(kl ,n2h+0*n3h1p+ijm) = +tmp4 - - H(n2h+1*n3h1p+ijm,kl ) = +tmp1 - H(kl ,n2h+1*n3h1p+ijm) = +tmp4 - - H(n2h+2*n3h1p+ijm,kl ) = +tmp2 - H(kl ,n2h+2*n3h1p+ijm) = +tmp3 - - H(n2h+3*n3h1p+ijm,kl ) = +tmp2 - H(kl ,n2h+3*n3h1p+ijm) = +tmp4 - - end do - end do - - end do - end do - end do - - !------------! - ! Block 3h1p ! - !------------! - - ijm = 0 - do i=nC+1,nO - do j=i+1,nO - do m=1,nS - ijm = ijm + 1 - - tmp = - eGW(i) - eGW(j) + Om(m) - - H(n2h+0*n3h1p+ijm,n2h+0*n3h1p+ijm) = tmp - H(n2h+1*n3h1p+ijm,n2h+1*n3h1p+ijm) = tmp - H(n2h+2*n3h1p+ijm,n2h+2*n3h1p+ijm) = tmp - H(n2h+3*n3h1p+ijm,n2h+3*n3h1p+ijm) = tmp - - end do - end do - end do - -!-------------------------! -! Diagonalize supermatrix ! -!-------------------------! - - call diagonalize_general_matrix(nH,H,OmBSE,X) - - do s=1,nH - order(s) = s - end do - - call quick_sort(OmBSE,order,nH) - call set_order(X,order,nH,nH) - -!-----------------! -! Compute weights ! -!-----------------! - - Z(:) = 0d0 - do s=1,nH - do ij=1,n2h - Z(s) = Z(s) + X(ij,s)**2 - end do - end do - -!--------------! -! Dump results ! -!--------------! - - write(*,*)'-------------------------------------------' - write(*,*)' Upfolded ppBSE excitation energies (eV) ' - write(*,*)'-------------------------------------------' - write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X)') & - '|','#','|','OmBSE (eV)','|','Z','|' - write(*,*)'-------------------------------------------' - - do s=1,min(nH,maxH) - if(Z(s) > 1d-7) & - write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & - '|',s,'|',-OmBSE(s)*HaToeV,'|',Z(s),'|' - end do - - write(*,*)'-------------------------------------------' - write(*,*) - -end subroutine diff --git a/src/GW/UG0W0.f90 b/src/GW/UG0W0.f90 index 44773a2..493f464 100644 --- a/src/GW/UG0W0.f90 +++ b/src/GW/UG0W0.f90 @@ -194,10 +194,10 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@UHF correlation energy (spin-conserved) = ',EcBSE(1),' au' - write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@UHF correlation energy (spin-flip) = ',EcBSE(2),' au' - write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@UHF correlation energy = ',sum(EcBSE),' au' - write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@UHF total energy = ',ENuc + EUHF + sum(EcBSE),' au' + write(*,'(2X,A60,F20.10,A3)') 'Tr@BSE@G0W0@UHF correlation energy (spin-conserved) = ',EcBSE(1),' au' + write(*,'(2X,A60,F20.10,A3)') 'Tr@BSE@G0W0@UHF correlation energy (spin-flip) = ',EcBSE(2),' au' + write(*,'(2X,A60,F20.10,A3)') 'Tr@BSE@G0W0@UHF correlation energy = ',sum(EcBSE),' au' + write(*,'(2X,A60,F20.10,A3)') 'Tr@BSE@G0W0@UHF total energy = ',ENuc + EUHF + sum(EcBSE),' au' write(*,*)'-------------------------------------------------------------------------------' write(*,*) @@ -210,10 +210,10 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@G0W0@UHF correlation energy (spin-conserved) = ',EcBSE(1),' au' - write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@G0W0@UHF correlation energy (spin-flip) = ',EcBSE(2),' au' - write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@G0W0@UHF correlation energy = ',sum(EcBSE),' au' - write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@G0W0@UHF total energy = ',ENuc + EUHF + sum(EcBSE),' au' + write(*,'(2X,A60,F20.10,A3)') 'AC@BSE@G0W0@UHF correlation energy (spin-conserved) = ',EcBSE(1),' au' + write(*,'(2X,A60,F20.10,A3)') 'AC@BSE@G0W0@UHF correlation energy (spin-flip) = ',EcBSE(2),' au' + write(*,'(2X,A60,F20.10,A3)') 'AC@BSE@G0W0@UHF correlation energy = ',sum(EcBSE),' au' + write(*,'(2X,A60,F20.10,A3)') 'AC@BSE@G0W0@UHF total energy = ',ENuc + EUHF + sum(EcBSE),' au' write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/GW/UGW_phACFDT.f90 b/src/GW/UGW_phACFDT.f90 index 9801c0c..9693d8d 100644 --- a/src/GW/UGW_phACFDT.f90 +++ b/src/GW/UGW_phACFDT.f90 @@ -112,14 +112,14 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,spin_conserved,spin_flip, allocate(OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin)) allocate(Aph(nS_sc,nS_sc),Bph(nS_sc,nS_sc),KA(nS_sc,nS_sc),KB(nS_sc,nS_sc)) - call phULR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) - if(.not.TDA) call phULR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) + call phULR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) + if(.not.TDA_W) call phULR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) call phULR(TDA_W,nS_aa,nS_bb,nS_sc,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA) call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) - call UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,OmRPA,rho_RPA,KA) - call UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,OmRPA,rho_RPA,KB) + call UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,OmRPA,rho_RPA,KA) + if(.not.TDA) call UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,OmRPA,rho_RPA,KB) ! Spin-conserved manifold @@ -144,24 +144,24 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,spin_conserved,spin_flip, if(doXBS) then - call phULR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) - if(.not.TDA) call phULR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) + call phULR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) + if(.not.TDA_W) call phULR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) call phULR(TDA_W,nS_aa,nS_bb,nS_sc,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA) call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) - call UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,OmRPA,rho_RPA,KA) - call UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,OmRPA,rho_RPA,KB) + call UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,OmRPA,rho_RPA,KA) + if(.not.TDA) call UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,OmRPA,rho_RPA,KB) end if - call phULR_A(isp_W,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) - if(.not.TDA) call phULR_B(isp_W,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) + call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) + if(.not.TDA) call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) Aph(:,:) = Aph(:,:) + KA(:,:) if(.not.TDA) Bph(:,:) = Bph(:,:) + KB(:,:) - call phULR(TDA_W,nS_aa,nS_bb,nS_sc,Aph,Bph,EcAC(ispin),OmRPA,XpY_RPA,XmY_RPA) + call phULR(TDA,nS_aa,nS_bb,nS_sc,Aph,Bph,EcAC(ispin),Om,XpY,XmY) call phUACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc, & ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,XmY,Ec(iAC,ispin)) @@ -208,24 +208,24 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,spin_conserved,spin_flip, if(doXBS) then - call phULR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) - if(.not.TDA) call phULR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) + call phULR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) + if(.not.TDA_W) call phULR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) call phULR(TDA_W,nS_aa,nS_bb,nS_sc,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA) call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) - call UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,OmRPA,rho_RPA,KA) - call UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,OmRPA,rho_RPA,KB) + call UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,OmRPA,rho_RPA,KA) + if(.not.TDA) call UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,OmRPA,rho_RPA,KB) end if - call phULR_A(isp_W,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,lambda,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) - if(.not.TDA) call phULR_B(isp_W,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) + call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,lambda,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) + if(.not.TDA) call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) Aph(:,:) = Aph(:,:) + KA(:,:) if(.not.TDA) Bph(:,:) = Bph(:,:) + KB(:,:) - call phULR(TDA_W,nS_aa,nS_bb,nS_sf,Aph,Bph,EcAC(ispin),OmRPA,XpY_RPA,XmY_RPA) + call phULR(TDA,nS_aa,nS_bb,nS_sf,Aph,Bph,EcAC(ispin),Om,XpY,XmY) call phUACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf, & ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,XmY,Ec(iAC,ispin)) diff --git a/src/GW/UGW_phBSE.f90 b/src/GW/UGW_phBSE.f90 index 1bf4cdc..e876484 100644 --- a/src/GW/UGW_phBSE.f90 +++ b/src/GW/UGW_phBSE.f90 @@ -117,11 +117,11 @@ subroutine UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_fli call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) if(.not.TDA) call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) - call UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,OmRPA,rho_RPA,KA) + call UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,OmRPA,rho_RPA,KA) if(.not.TDA) call UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,OmRPA,rho_RPA,KB) Aph(:,:) = Aph(:,:) + KA(:,:) -! if(.not.TDA) Bph(:,:) = Bph(:,:) + KB(:,:) + if(.not.TDA) Bph(:,:) = Bph(:,:) + KB(:,:) call phULR(TDA,nS_aa,nS_bb,nS_sc,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE) @@ -155,10 +155,11 @@ subroutine UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_fli allocate(OmBSE(nS_sf),XpY_BSE(nS_sf,nS_sf),XmY_BSE(nS_sf,nS_sf)) ! Compute spin-conserved BSE excitation energies + call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,1d0,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) if(.not.TDA) call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) - call UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sf,1d0,OmRPA,rho_RPA,KA) + call UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sf,1d0,OmRPA,rho_RPA,KA) if(.not.TDA) call UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sf,1d0,OmRPA,rho_RPA,KB) Aph(:,:) = Aph(:,:) + KA(:,:) diff --git a/src/GW/UGW_phBSE_static_kernel_A.f90 b/src/GW/UGW_phBSE_static_kernel_A.f90 index ef08ef3..6b2f0b7 100644 --- a/src/GW/UGW_phBSE_static_kernel_A.f90 +++ b/src/GW/UGW_phBSE_static_kernel_A.f90 @@ -32,6 +32,10 @@ subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s double precision,intent(out) :: KA(nSt,nSt) +! Initialization + + KA(:,:) = 0d0 + !--------------------------------------------------! ! Build BSE matrix for spin-conserving transitions ! !--------------------------------------------------! @@ -55,7 +59,7 @@ subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s chi = chi + rho(i,j,kc,1)*rho(a,b,kc,1)*Om(kc)/eps end do - KA(ia,jb) = KA(ia,jb) + 2d0*lambda*chi + KA(ia,jb) = 2d0*lambda**2*chi end do end do @@ -79,7 +83,7 @@ subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s chi = chi + rho(i,j,kc,2)*rho(a,b,kc,2)*Om(kc)/eps end do - KA(nSa+ia,nSa+jb) = KA(nSa+ia,nSa+jb) + 2d0*lambda*chi + KA(nSa+ia,nSa+jb) = 2d0*lambda**2*chi end do end do @@ -111,7 +115,7 @@ subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s chi = chi + rho(i,j,kc,1)*rho(a,b,kc,2)*Om(kc)/eps end do - KA(ia,jb) = KA(ia,jb) + 2d0*lambda*chi + KA(ia,jb) = 2d0*lambda**2*chi end do end do @@ -135,7 +139,7 @@ subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s chi = chi + rho(i,j,kc,2)*rho(a,b,kc,1)*Om(kc)/eps end do - KA(nSa+ia,nSa+jb) = KA(nSa+ia,nSa+jb) + 2d0*lambda*chi + KA(nSa+ia,nSa+jb) = 2d0*lambda**2*chi end do end do diff --git a/src/GW/UGW_phBSE_static_kernel_B.f90 b/src/GW/UGW_phBSE_static_kernel_B.f90 index c374bf2..d6c6d21 100644 --- a/src/GW/UGW_phBSE_static_kernel_B.f90 +++ b/src/GW/UGW_phBSE_static_kernel_B.f90 @@ -32,6 +32,10 @@ subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s double precision,intent(out) :: KB(nSt,nSt) +! Initialization + + KB(:,:) = 0d0 + !--------------------------------------------------! ! Build BSE matrix for spin-conserving transitions ! !--------------------------------------------------! @@ -55,7 +59,7 @@ subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s chi = chi + rho(i,b,kc,1)*rho(a,j,kc,1)*Om(kc)/eps end do - KB(ia,jb) = KB(ia,jb) + 2d0*lambda*chi + KB(ia,jb) = 2d0*lambda**2*chi end do end do @@ -80,7 +84,7 @@ subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s chi = chi + rho(i,b,kc,2)*rho(a,j,kc,2)*Om(kc)/eps end do - KB(nSa+ia,nSa+jb) = KB(nSa+ia,nSa+jb) + 2d0*lambda*chi + KB(nSa+ia,nSa+jb) = 2d0*lambda**2*chi end do end do @@ -113,7 +117,7 @@ subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s chi = chi + rho(i,b,kc,1)*rho(a,j,kc,2)*Om(kc)/eps end do - KB(ia,nSa+jb) = KB(ia,nSa+jb) + 2d0*lambda*chi + KB(ia,nSa+jb) = 2d0*lambda**2*chi end do end do @@ -137,7 +141,7 @@ subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s chi = chi + rho(i,b,kc,2)*rho(a,j,kc,1)*Om(kc)/eps end do - KB(nSa+ia,jb) = KB(nSa+ia,jb) + 2d0*lambda*chi + KB(nSa+ia,jb) = 2d0*lambda**2*chi end do 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/GW/ufRG0W0.f90 b/src/GW/ufRG0W0.f90 index bf89929..010dd45 100644 --- a/src/GW/ufRG0W0.f90 +++ b/src/GW/ufRG0W0.f90 @@ -111,6 +111,10 @@ subroutine ufRG0W0(dotest,TDA_W,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) deallocate(Aph,Bph,XpY,XmY) + else + + allocate(rho(0,0,0)) + end if !-------------------------! diff --git a/src/GW/ufRGW.f90 b/src/GW/ufRGW.f90 index 1861ba0..970ef14 100644 --- a/src/GW/ufRGW.f90 +++ b/src/GW/ufRGW.f90 @@ -111,6 +111,10 @@ subroutine ufRGW(dotest,TDA_W,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) deallocate(Aph,Bph,XpY,XmY) + else + + allocate(rho(0,0,0)) + end if ! Initialization diff --git a/src/HF/print_UHF.f90 b/src/HF/print_UHF.f90 index 70dc834..82e744c 100644 --- a/src/HF/print_UHF.f90 +++ b/src/HF/print_UHF.f90 @@ -31,7 +31,6 @@ subroutine print_UHF(nBas,nO,S,eHF,c,P,ENuc,ET,EV,EJ,Ex,EUHF,dipole) double precision :: Sz double precision :: Sx2,Sy2,Sz2 integer :: mu,nu - double precision,allocatable :: qa(:),qb(:) logical :: dump_orb = .false. @@ -119,18 +118,4 @@ subroutine print_UHF(nBas,nO,S,eHF,c,P,ENuc,ET,EV,EJ,Ex,EUHF,dipole) call vecout(nBas,eHF(:,2)) write(*,*) - allocate(qa(nBas),qb(nBas)) - - qa(:) = 0d0 - qb(:) = 0d0 - do mu=1,nBas - do nu=1,nBas - qa(mu) = qa(mu) + P(mu,nu,1)*S(nu,mu) - qb(mu) = qb(mu) + P(mu,nu,2)*S(nu,mu) - end do - end do - - call vecout(nBas,qa) - call vecout(nBas,qb) - end subroutine diff --git a/src/LR/phULR_A.f90 b/src/LR/phULR_A.f90 index 19fed4e..e3fe0d6 100644 --- a/src/LR/phULR_A.f90 +++ b/src/LR/phULR_A.f90 @@ -39,6 +39,10 @@ subroutine phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,eHF,ERI_aaaa,E delta_dRPA = 0d0 if(dRPA) delta_dRPA = 1d0 +! Initialization + + Aph(:,:) = 0d0 + !---------------------------------------------- ! Build A matrix for spin-conserved transitions !---------------------------------------------- @@ -127,8 +131,6 @@ subroutine phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,eHF,ERI_aaaa,E if(ispin == 2) then - Aph(:,:) = 0d0 - ! abab block ia = 0 diff --git a/src/LR/phULR_B.f90 b/src/LR/phULR_B.f90 index a9a231e..69c0e8c 100644 --- a/src/LR/phULR_B.f90 +++ b/src/LR/phULR_B.f90 @@ -38,6 +38,10 @@ subroutine phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,ERI_aaaa,ERI_a delta_dRPA = 0d0 if(dRPA) delta_dRPA = 1d0 +! Initialization + + Bph(:,:) = 0d0 + !----------------------------------------------- ! Build B matrix for spin-conserving transitions !----------------------------------------------- @@ -124,8 +128,6 @@ subroutine phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,ERI_aaaa,ERI_a if(ispin == 2) then - Bph(:,:) = 0d0 - ! abba block ia = 0 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/LR/ppLR_GW_davidson.f90 b/src/LR/ppLR_GW_davidson.f90 index a8c1af2..2c8ac12 100644 --- a/src/LR/ppLR_GW_davidson.f90 +++ b/src/LR/ppLR_GW_davidson.f90 @@ -33,7 +33,7 @@ subroutine ppLR_GW_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, nS, lambda, e, eF, call wall_time(t1) - if((nOO+nVV) .le. 20000) then + if((nOO+nVV) .le. 30000) then call ppLR_GW_HR_calc_oneshot(ispin, nOrb, nC, nO, nR, nOO, nVV, nS, lambda, e(1), eF, n_states_diag, & ERI(1,1,1,1), eta, rho(1,1,1), Om(1), U(1,1), W(1,1)) diff --git a/src/LR/ppLR_davidson.f90 b/src/LR/ppLR_davidson.f90 index 6a4d7ed..c2f5aef 100644 --- a/src/LR/ppLR_davidson.f90 +++ b/src/LR/ppLR_davidson.f90 @@ -4,7 +4,7 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, ERI, & supp_data_int, supp_data_int_size, & supp_data_dbl, supp_data_dbl_size, & - Om, R, n_states, n_states_diag, kernel) + Om, R, n_states, n_states_diag, kernel, mode_dav) ! ! Extract the low n_states @@ -18,6 +18,54 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, ! (-B.T -D) ! + implicit none + + logical, intent(in) :: TDA + integer, intent(in) :: ispin + integer, intent(in) :: nC, nO, nR, nOrb, nOO, nVV + integer, intent(in) :: n_states ! nb of physical states + integer, intent(in) :: n_states_diag ! nb of states used to get n_states + integer, intent(in) :: supp_data_int_size + integer, intent(in) :: supp_data_dbl_size + integer, intent(in) :: mode_dav + character(len=*), intent(in) :: kernel + double precision, intent(in) :: lambda, eF + double precision, intent(in) :: e(nOrb) + double precision, intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + integer, intent(in) :: supp_data_int(supp_data_int_size) + double precision, intent(in) :: supp_data_dbl(supp_data_dbl_size) + double precision, intent(out) :: Om(n_states) + double precision, intent(out) :: R(nOO+nVV,n_states_diag) + + if(mode_dav .eq. 1) then + + call ppLR_davidson_1(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e(1), eF, ERI(1,1,1,1), & + supp_data_int(1), supp_data_int_size, supp_data_dbl(1), supp_data_dbl_size, & + Om(1), R(1,1), n_states, n_states_diag, kernel) + + elseif(mode_dav .eq. 2) then + + call ppLR_davidson_2(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e(1), eF, ERI(1,1,1,1), & + supp_data_int(1), supp_data_int_size, supp_data_dbl(1), supp_data_dbl_size, & + Om(1), R(1,1), n_states, n_states_diag, kernel) + + else + + print*, " unknown Davidson's variant" + stop + + endif + + return +end + +! --- + +subroutine ppLR_davidson_1(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, ERI, & + supp_data_int, supp_data_int_size, & + supp_data_dbl, supp_data_dbl_size, & + Om, R, n_states, n_states_diag, kernel) + use omp_lib implicit none @@ -43,16 +91,15 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, integer :: shift1, shift2 integer :: i, j, k, l, ab integer :: p, q, mm, i_data, nS - integer :: i_omax(n_states) logical :: converged - character*(16384) :: write_buffer + character(len=6+41*n_states) :: write_buffer double precision :: r1, r2, dtwo_pi double precision :: lambda_tmp - double precision :: to_print(2,n_states) double precision :: mem double precision :: eta double precision :: t1, t2, tt1, tt2 character(len=len(kernel)) :: kernel_name + integer, allocatable :: i_omax(:) double precision, allocatable :: H_diag(:) double precision, allocatable :: W(:,:) double precision, allocatable :: U(:,:) @@ -61,6 +108,7 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, double precision, allocatable :: overlap(:) double precision, allocatable :: S_check(:,:) double precision, allocatable :: rho_tmp(:,:,:), Om_tmp(:) + double precision, allocatable :: to_print(:,:) double precision, external :: u_dot_u @@ -92,7 +140,8 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, write(*,'(A40, A12)') 'Kernel: ', kernel_name - + allocate(i_omax(n_states)) + allocate(to_print(2,n_states)) allocate(H_diag(N)) allocate(U(N,M)) allocate(W(N,M)) @@ -101,7 +150,7 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, allocate(residual_norm(n_states_diag)) mem = 8.d0 * dble(nOrb + nOrb**4 + N*n_states) & - + 8.d0 * dble(supp_data_dbl_size) + 4.d0 * dble(supp_data_int_size) + + 8.d0 * dble(2*supp_data_dbl_size) + 4.d0 * dble(2*supp_data_int_size) write(*,'(A40, F12.4)') 'I/O mem (GB) = ', mem / (1024.d0*1024.d0*1024.d0) @@ -118,6 +167,10 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, call ppLR_RPA_H_diag(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e(1), eF, & ERI(1,1,1,1), H_diag(1)) + ! to avoid compiler warnings + allocate(rho_tmp(0,0,0)) + allocate(Om_tmp(0)) + elseif(kernel_name .eq. "gw") then nS = supp_data_int(1) @@ -214,6 +267,7 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, write(6,'(A)') write_buffer(1:6+41*n_states) + W = 0.d0 converged = .False. itertot = 0 @@ -325,7 +379,6 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, W(1,1), size(W, 1), h_vec(1,1), size(h_vec, 1), & 0.d0, W(1,shift2+1), size(W, 1)) - ! check if W1 = U1 h_val !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(i, k) & !$OMP SHARED(n_states, n_states_diag, N, shift2, U, h_val, W, H_diag, residual_norm, to_print) @@ -342,13 +395,11 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, enddo !$OMP END DO !$OMP END PARALLEL - !print*, " to_print", to_print if((itertot > 1) .and. (iter == 1)) then continue else - write(*,'(1X, I3, 1X, 100(1X, F16.10, 1X, F12.6))') iter-1, to_print(1:2,1:n_states) - !write(*, '(1X, I3, 1X, 100(1X, F16.10, 1X, F16.10, 1X, F16.10))') iter-1, to_print(1:2,1:n_states) + write(*,'(1X, I3, 1X, 10000(1X, F16.10, 1X, F12.6))') iter-1, to_print(1:2,1:n_states) endif !call wall_time(tt2) @@ -361,7 +412,7 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, endif do k = 1, n_states - if(residual_norm(k) > 1.d8) then + if(residual_norm(k) > 1.d10) then print *, 'Davidson failed' stop -1 endif @@ -377,7 +428,8 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, call dgemm('N', 'N', N, n_states_diag, shift2, 1.d0, & W(1,1), size(W, 1), h_vec(1,1), size(h_vec, 1), & - 0.d0, R, size(R, 1)) + 0.d0, R(1,1), size(R, 1)) + do k = 1, n_states_diag do i = 1, N W(i,k) = R(i,k) @@ -427,6 +479,8 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, print*, k, Om(k) enddo + deallocate(i_omax) + deallocate(to_print) deallocate(H_diag) deallocate(U) deallocate(W) @@ -449,3 +503,404 @@ end ! --- +subroutine ppLR_davidson_2(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, ERI, & + supp_data_int, supp_data_int_size, & + supp_data_dbl, supp_data_dbl_size, & + Om, R, n_states, n_states_diag, kernel) + + use omp_lib + + implicit none + + logical, intent(in) :: TDA + integer, intent(in) :: ispin + integer, intent(in) :: nC, nO, nR, nOrb, nOO, nVV + integer, intent(in) :: n_states ! nb of physical states + integer, intent(in) :: n_states_diag ! nb of states used to get n_states + integer, intent(in) :: supp_data_int_size + integer, intent(in) :: supp_data_dbl_size + character(len=*), intent(in) :: kernel + double precision, intent(in) :: lambda, eF + double precision, intent(in) :: e(nOrb) + double precision, intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + integer, intent(in) :: supp_data_int(supp_data_int_size) + double precision, intent(in) :: supp_data_dbl(supp_data_dbl_size) + double precision, intent(out) :: Om(n_states) + double precision, intent(out) :: R(nOO+nVV,n_states_diag) + + integer :: N, M, num_threads, n_states_delta + integer :: it_start, it_delta, it_size + integer :: iter, itermax, itertot + integer :: i, j, k, l, ab + integer :: p, q, mm, i_data, nS + logical :: converged + double precision :: r1, r2, dtwo_pi + double precision :: mem + double precision :: eta + double precision :: t1, t2, tt1, tt2 + character(len=len(kernel)) :: kernel_name + integer, allocatable :: i_omax(:) + character(len=:), allocatable :: write_buffer + double precision, allocatable :: to_print(:,:) + double precision, allocatable :: H_diag(:) + double precision, allocatable :: W0(:,:), W1(:,:) + double precision, allocatable :: U0(:,:), U1(:,:) + double precision, allocatable :: h(:,:), h_vec(:,:), h_val(:) + double precision, allocatable :: residual_norm(:) + double precision, allocatable :: rho_tmp(:,:,:), Om_tmp(:) + + double precision, external :: u_dot_u + + call wall_time(t1) + + dtwo_pi = 6.283185307179586d0 + + N = nOO + nVV + + n_states_delta = min(max(25, n_states_diag/2), n_states_diag) + itermax = 8 + M = n_states_diag + itermax * n_states_delta + + call lower_case(trim(kernel), kernel_name) + + if(M .ge. N) then + print*, 'N = ', N + print*, 'M = ', M + print*, ' use Lapack or decrease n_states and/or itermax ' + stop + endif + + write(6,'(A)') '' + write(6,'(A)') 'Davidson Diagonalization' + write(6,'(A)') '------------------------' + write(6,'(A)') '' + + write(*,'(A40, I12)') 'Number of states = ', n_states + write(*,'(A40, I12)') 'Number of states in diag = ', n_states_diag + write(*,'(A40, I12)') 'Number of states to add = ', n_states_delta + write(*,'(A40, I12)') 'Number of basis functions = ', N + write(*,'(A40, A12)') 'Kernel: ', kernel_name + + + + allocate(character(len=50*n_states) :: write_buffer) + allocate(i_omax(n_states)) + allocate(to_print(2,n_states)) + allocate(H_diag(N)) + allocate(U0(N,M), U1(N,n_states_diag)) + allocate(W0(N,M), W1(N,n_states_diag)) + allocate(h(M,M), h_vec(M,M), h_val(M)) + allocate(residual_norm(n_states_diag)) + + mem = 8.d0 * dble(nOrb) + 8.d0 * dble(nOrb)**4 + 8.d0 * dble(N*n_states) & + + 8.d0 * dble(2*supp_data_dbl_size) + 4.d0 * dble(2*supp_data_int_size) + + write(*,'(A40, F12.4)') 'I/O mem (GB) = ', mem / (1024.d0*1024.d0*1024.d0) + + mem = 8.d0 * dble(N) & + + 8.d0 * dble(N*M) & + + 8.d0 * dble(N*M) & + + 8.d0 * dble(N*n_states_diag) & + + 8.d0 * dble(N*n_states_diag) & + + 8.d0 * dble(M*M) & + + 8.d0 * dble(M*M) & + + 8.d0 * dble(M) & + + 8.d0 * dble(n_states_diag) & + + 1.d0 * dble(50*n_states) + + write(*,'(A40, F12.4)') 'tmp mem (GB) = ', mem / (1024.d0*1024.d0*1024.d0) + + num_threads = omp_get_max_threads() + write(*,'(A40, I12)') 'Number of threads = ', num_threads + + + if(kernel_name .eq. "rpa") then + + allocate(rho_tmp(0,0,0)) + allocate(Om_tmp(0)) + + call ppLR_RPA_H_diag(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e(1), eF, & + ERI(1,1,1,1), H_diag(1)) + + elseif(kernel_name .eq. "gw") then + + nS = supp_data_int(1) + + allocate(rho_tmp(nS,nOrb,nOrb)) + allocate(Om_tmp(nS)) + + eta = supp_data_dbl(1) + i_data = 1 + do q = 1, nOrb + do p = 1, nOrb + do mm = 1, nS + i_data = i_data + 1 + rho_tmp(mm,p,q) = supp_data_dbl(i_data) + enddo + enddo + enddo + do mm = 1, nS + i_data = i_data + 1 + Om_tmp(mm) = supp_data_dbl(i_data) + enddo + + call ppLR_GW_H_diag(ispin, nOrb, nC, nO, nR, nOO, nVV, nS, lambda, e(1), eF, & + ERI(1,1,1,1), eta, rho_tmp(1,1,1), Om_tmp(1), H_diag(1)) + + !! TODO + !elseif(kernel_name .eq. "gf2") then + + else + + print*, ' kernel not supported', kernel + stop + + endif + + U0 = 0.d0 + W0 = 0.d0 + U1 = 0.d0 + W1 = 0.d0 + + ! TODO: improve guess + ! initialize guess + R = 0.d0 + do k = 1, n_states + R(k,k) = 1.d0 + enddo + do k = n_states+1, n_states_diag + do i = 1, N + call random_number(r1) + call random_number(r2) + r1 = dsqrt(-2.d0*dlog(r1)) + r2 = dtwo_pi*r2 + R(i,k) = r1*dcos(r2) + enddo + R(k,k) = R(k,k) + 10.d0 + call normalize(R(1,k), N) + enddo + + do k = 1, n_states_diag + U0(:,k) = R(:,k) + enddo + + + write(6,'(A)') '' + write_buffer = '=====' + do i = 1, n_states + write_buffer = trim(write_buffer)//' ================ ===========' + enddo + write(6,'(A)') write_buffer(1:6+41*n_states) + write_buffer = 'Iter' + do i = 1, n_states + write_buffer = trim(write_buffer)//' Energy Residual ' + enddo + write(6,'(A)') write_buffer(1:6+41*n_states) + write_buffer = '=====' + do i = 1, n_states + write_buffer = trim(write_buffer)//' ================ ===========' + enddo + write(6,'(A)') write_buffer(1:6+41*n_states) + + + converged = .False. + itertot = 0 + + do while (.not.converged) + + itertot = itertot + 1 + if(itertot == itermax) then + print*, 'exit before convergence !' + print*, 'itertot == itermax', itertot + exit + endif + + do iter = 1, itermax-1 + + if(iter .eq. 1) then + it_start = 0 + it_delta = n_states_diag + else + it_start = n_states_diag + n_states_delta * (iter - 2) + it_delta = n_states_delta + endif + + it_size = it_start + it_delta + + if((iter > 1) .or. (itertot == 1)) then + + !call wall_time(tt1) + + call ortho_qr(U0(1,1), size(U0, 1), N, it_size) + + if(kernel_name .eq. "rpa") then + + call ppLR_RPA_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e(1), eF, it_delta, & + ERI(1,1,1,1), & + U0(1,it_start+1), W0(1,it_start+1)) + + elseif(kernel_name .eq. "gw") then + + call ppLR_GW_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, nS, lambda, e(1), eF, it_delta, & + ERI(1,1,1,1), eta, rho_tmp(1,1,1), Om_tmp(1), & + U0(1,it_start+1), W0(1,it_start+1)) + + !! TODO + !elseif(kernel_name .eq. "gf2") then + + endif + + else + + ! computed below + continue + endif + + ! h = U0.T H U0 + call dgemm('T', 'N', it_size, it_size, N, 1.d0, & + U0(1,1), size(U0, 1), W0(1,1), size(W0, 1), & + 0.d0, h(1,1), size(h, 1)) + + ! h h_vec = h_val h_vec + call diag_nonsym_right(it_size, h(1,1), size(h, 1), h_vec(1,1), size(h_vec, 1), h_val(1), size(h_val, 1)) + + ! U1(:,1:it_delta) = U0 h_vec(:,1:it_delta) + call dgemm('N', 'N', N, it_delta, it_size, 1.d0, & + U0(1,1), size(U0, 1), h_vec(1,1), size(h_vec, 1), & + 0.d0, U1(1,1), size(U1, 1)) + + do k = 1, it_delta + call normalize(U1(1,k), N) + enddo + + ! W1(:,1:it_delta) = W0 h_vec(:,1:it_delta) + call dgemm('N', 'N', N, it_delta, it_size, 1.d0, & + W0(1,1), size(W0, 1), h_vec(1,1), size(h_vec, 1), & + 0.d0, W1(1,1), size(W1, 1)) + + !$OMP PARALLEL & + !$OMP DEFAULT(NONE) & + !$OMP PRIVATE(i, k) & + !$OMP SHARED(n_states, it_size, it_delta, N, U0, U1, & + !$OMP h_val, W1, H_diag, residual_norm, to_print) + !$OMP DO + do k = 1, it_delta + do i = 1, N + U1(i,k) = (h_val(k) * U1(i,k) - W1(i,k)) / max(dabs(H_diag(i) - h_val(k)), 1.d-2) + U0(i,it_size+k) = U1(i,k) + enddo + if(k <= n_states) then + residual_norm(k) = u_dot_u(U1(1,k), N) + to_print(1,k) = h_val(k) + to_print(2,k) = residual_norm(k) + endif + enddo + !$OMP END DO + !$OMP END PARALLEL + + + if((itertot > 1) .and. (iter == 1)) then + continue + else + write(*,'(1X, I3, 1X, 10000(1X, F16.10, 1X, F12.6))') iter-1, to_print(1:2,1:n_states) + endif + + !call wall_time(tt2) + !write(*,'(A50, F12.4)') 'wall time for one Davidson iteration (sec): ', tt2-tt1 + !stop + + if(iter > 1) then + converged = dabs(maxval(residual_norm(1:n_states))) < 1d-15 + endif + + do k = 1, n_states + if(residual_norm(k) > 1.d10) then + print *, 'Davidson failed' + stop -1 + endif + enddo + + if(converged) exit + + enddo ! loop over iter + + + ! Re-contract U0 and update W0 + ! -------------------------------- + + call dgemm('N', 'N', N, n_states_diag, it_size, 1.d0, & + W0(1,1), size(W0, 1), h_vec(1,1), size(h_vec, 1), & + 0.d0, R(1,1), size(R, 1)) + + do k = 1, n_states_diag + do i = 1, N + W0(i,k) = R(i,k) + enddo + enddo + + call dgemm('N', 'N', N, n_states_diag, it_size, 1.d0, & + U0(1,1), size(U0, 1), h_vec(1,1), size(h_vec, 1), & + 0.d0, R(1,1), size(R, 1)) + + do k = 1, n_states_diag + do i = 1, N + U0(i,k) = R(i,k) + enddo + enddo + + call ortho_qr(U0(1,1), size(U0, 1), N, n_states_diag) + + do j = 1, n_states_diag + k = 1 + do while((k < N) .and. (U0(k,j) == 0.d0)) + k = k+1 + enddo + if(U0(k,j) * R(k,j) < 0.d0) then + do i = 1, N + W0(i,j) = -W0(i,j) + enddo + endif + enddo + + enddo ! loop over while + + ! --- + + write_buffer = '=====' + do i = 1, n_states + write_buffer = trim(write_buffer)//' ================ ===========' + enddo + write(6,'(A)') trim(write_buffer) + write(6,'(A)') '' + + + print*, " Davidson eigenvalues" + do k = 1, n_states + Om(k) = h_val(k) + print*, k, Om(k) + enddo + + deallocate(write_buffer) + deallocate(i_omax) + deallocate(to_print) + deallocate(H_diag) + deallocate(U0, U1) + deallocate(W0, W1) + deallocate(h) + deallocate(h_vec) + deallocate(h_val) + deallocate(residual_norm) + + if(kernel_name .eq. "gw") then + deallocate(rho_tmp) + deallocate(Om_tmp) + endif + + call wall_time(t2) + write(*,'(A50, F12.4)') 'total wall time for Davidson (sec): ', t2-t1 + + return +end + +! --- + 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 8c31986..b58660c 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -249,15 +249,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) !-----------!