diff --git a/src/GT/GG0T0pp.f90 b/src/GT/GG0T0pp.f90 new file mode 100644 index 0000000..f7ced74 --- /dev/null +++ b/src/GT/GG0T0pp.f90 @@ -0,0 +1,255 @@ +subroutine GG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,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) + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: dotest + + logical,intent(in) :: doACFDT + logical,intent(in) :: exchange_kernel + logical,intent(in) :: doXBS + logical,intent(in) :: dophBSE + logical,intent(in) :: doppBSE + logical,intent(in) :: TDA_T + logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA + logical,intent(in) :: linearize + double precision,intent(in) :: eta + logical,intent(in) :: regularize + + 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) :: ENuc + double precision,intent(in) :: EGHF + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart) + +! Local variables + + logical :: print_T = .false. + + integer :: nOO + integer :: nVV + double precision :: EcRPA + double precision :: EcBSE + double precision :: EcGM + double precision,allocatable :: Bpp(:,:) + double precision,allocatable :: Cpp(:,:) + double precision,allocatable :: Dpp(:,:) + double precision,allocatable :: Om1(:) + double precision,allocatable :: X1(:,:) + double precision,allocatable :: Y1(:,:) + double precision,allocatable :: rho1(:,:,:) + double precision,allocatable :: Om2(:) + double precision,allocatable :: X2(:,:) + double precision,allocatable :: Y2(:,:) + double precision,allocatable :: rho2(:,:,:) + double precision,allocatable :: Sig(:) + double precision,allocatable :: Z(:) + double precision,allocatable :: eGT(:) + double precision,allocatable :: eGTlin(:) + double precision, allocatable :: Om(:), R(:,:) + + +! Output variables + +! Hello world + + write(*,*) + write(*,*)'**********************************' + write(*,*)'* Generalized G0T0pp Calculation *' + write(*,*)'**********************************' + write(*,*) + + +! TDA for T + + if(TDA_T) then + write(*,*) 'Tamm-Dancoff approximation activated for pp T-matrix!' + write(*,*) + end if + +! TDA + + if(TDA) then + write(*,*) 'Tamm-Dancoff approximation activated!' + write(*,*) + end if + +! Dimensions of the pp-RPA linear reponse matrices + + nOO = nO*(nO - 1)/2 + nVV = nV*(nV - 1)/2 + +! Memory allocation + + allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO), & + rho1(nOrb,nOrb,nVV),rho2(nOrb,nOrb,nOO),Sig(nOrb),Z(nOrb),eGT(nOrb),eGTlin(nOrb)) + +!---------------------------------------------- +! Compute T-matrix +!---------------------------------------------- + +! Compute linear response + + allocate(Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO)) + + 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) + +!---------------------------------------------- +! Compute excitation densities +!---------------------------------------------- + + call GGTpp_excitation_density(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2) + +!---------------------------------------------- +! Compute T-matrix version of the self-energy +!---------------------------------------------- + + if(regularize) call GTpp_regularization(nOrb,nC,nO,nV,nR,nOO,nVV,eHF,Om1,rho1,Om2,rho2) + + call GGTpp_self_energy_diag(eta,nOrb,nC,nO,nV,nR,nOO,nVV,eHF,Om1,rho1,Om2,rho2,EcGM,Sig,Z) + +!---------------------------------------------- +! Solve the quasi-particle equation +!---------------------------------------------- + + eGTlin(:) = eHF(:) + Z(:)*Sig(:) + + if(linearize) then + + write(*,*) ' *** Quasiparticle energies obtained by linearization *** ' + write(*,*) + + eGT(:) = eGTlin(:) + + else + + write(*,*) ' *** Quasiparticle energies obtained by root search *** ' + write(*,*) + + call GGTpp_QP_graph(eta,nOrb,nC,nO,nV,nR,nOO,nVV,eHF,Om1,rho1,Om2,rho2,eGTlin,eHF,eGT,Z) + + end if + +! call GGTpp_plot_self_energy(nOrb,nC,nO,nV,nR,nOO,nVV,eHF,eGT,Om1,rho1,Om2,rho2) + +!---------------------------------------------- +! Dump results +!---------------------------------------------- + +! Compute the ppRPA correlation energy + + allocate(Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO)) + + call ppGLR_C(nOrb,nC,nO,nV,nR,nVV,1d0,eGT,ERI,Cpp) + call ppGLR_D(nOrb,nC,nO,nV,nR,nOO,1d0,eGT,ERI,Dpp) + if(.not.TDA_T) call ppGLR_B(nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) + + call ppGLR(TDA_T,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA) + + deallocate(Bpp,Cpp,Dpp) + + call print_GG0T0pp(nOrb,nO,eHF,ENuc,EGHF,Sig,Z,eGT,EcGM,EcRPA) + +! Perform BSE calculation + +! if(dophBSE) then + +! call GGTpp_phBSE(TDA_T,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV, & +! Om1,X1,Y1,Om2,X2,Y2,rho1,rho2,ERI,dipole_int,eHF,eGT,EcBSE) + +! if(exchange_kernel) then + +! EcBSE = 0.5d0*EcBSE + +! end if + +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp@GHF correlation energy = ',EcBSE,' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp@GHF total energy = ',ENuc + EGHF + EcBSE,' au' +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) + +! Compute the BSE correlation energy via the adiabatic connection + +! if(doACFDT) then + +! write(*,*) '--------------------------------------------------------' +! write(*,*) 'Adiabatic connection version of phBSE correlation energy' +! write(*,*) '--------------------------------------------------------' +! write(*,*) + +! if(doXBS) then + +! write(*,*) '*** scaled screening version (XBS) ***' +! write(*,*) + +! end if + +! call RGTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,dophBSE,eta,nOrb,nC,nO,nV,nR,nS, & +! nOO,nVV,Om1,X1,Y1,Om2,X2,Y2,rho1,rho2,ERI,eHF,eGT,EcBSE) + +! if(exchange_kernel) then + +! EcBSE = 0.5d0*EcBSE + +! end if + +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp@GHF correlation energy = ',EcBSE,' au' +! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp@GHF total energy = ',ENuc + EGHF + EcBSE,' au' +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) + +! end if + +! end if + +! if(doppBSE) then + +! call GGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nOO,nVV, & +! Om1,X1,Y1,Om2,X2,Y2,rho1,rho2,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 + + call dump_test_value('G','G0T0pp correlation energy',EcRPA) + call dump_test_value('G','G0T0pp HOMO energy',eGT(nO)) + call dump_test_value('G','G0T0pp LUMO energy',eGT(nO+1)) + + end if + +end subroutine diff --git a/src/GT/GGTpp_QP_graph.f90 b/src/GT/GGTpp_QP_graph.f90 new file mode 100644 index 0000000..bd23bf4 --- /dev/null +++ b/src/GT/GGTpp_QP_graph.f90 @@ -0,0 +1,87 @@ +subroutine GGTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Om1,rho1,Om2,rho2,eGTlin,eOld,eGT,Z) + +! Compute the graphical solution of the QP equation + + 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) :: nOO + integer,intent(in) :: nVV + + double precision,intent(in) :: eta + double precision,intent(in) :: eHF(nBas) + 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) :: eGTlin(nBas) + double precision,intent(in) :: eOld(nBas) + +! Local variables + + integer :: p + integer :: nIt + integer,parameter :: maxIt = 64 + double precision,parameter :: thresh = 1d-6 + double precision,external :: GGTpp_SigC,GGTpp_dSigC + double precision :: SigC,dSigC + double precision :: f,df + double precision :: w + +! Output variables + + double precision,intent(out) :: eGT(nBas) + double precision,intent(out) :: Z(nBas) + +! Run Newton's algorithm to find the root + + write(*,*)'-----------------------------------------------------' + write(*,'(A5,1X,A3,1X,A15,1X,A15,1X,A10)') 'Orb.','It.','e_GTpplin (eV)','e_GTpplin (eV)','Z' + write(*,*)'-----------------------------------------------------' + + do p=nC+1,nBas-nR + + w = eGTlin(p) + nIt = 0 + f = 1d0 + + do while (abs(f) > thresh .and. nIt < maxIt) + + nIt = nIt + 1 + + SigC = GGTpp_SigC(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,eOld,Om1,rho1,Om2,rho2) + dSigC = GGTpp_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,eOld,Om1,rho1,Om2,rho2) + f = w - eHF(p) - SigC + df = 1d0/(1d0 - dSigC) + w = w - df*f + + end do + + if(nIt == maxIt) then + + eGT(p) = eGTlin(p) + write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6,1X,A12)') p,nIt,eGTlin(p)*HaToeV,eGT(p)*HaToeV,Z(p),'Cvg Failed!' + + else + + eGT(p) = w + Z(p) = df + + write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6)') p,nIt,eGTlin(p)*HaToeV,eGT(p)*HaToeV,Z(p) + + end if + + end do + + write(*,*)'-----------------------------------------------------' + write(*,*) + +end subroutine diff --git a/src/GT/GGTpp_SigC.f90 b/src/GT/GGTpp_SigC.f90 new file mode 100644 index 0000000..3914806 --- /dev/null +++ b/src/GT/GGTpp_SigC.f90 @@ -0,0 +1,62 @@ +double precision function GGTpp_SigC(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho2) + +! Compute diagonal of the correlation part of the self-energy + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: p + double precision,intent(in) :: w + double precision,intent(in) :: eta + integer,intent(in) :: nBas + 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) :: e(nBas) + double precision,intent(in) :: Om1(nVV) + double precision,intent(in) :: rho1(nBas,nBas,nVV) + double precision,intent(in) :: Om2(nOO) + double precision,intent(in) :: rho2(nBas,nBas,nOO) + +! Local variables + + integer :: i,a,cd,kl + double precision :: eps + +! Initialize + + GGTpp_SigC = 0d0 + +!-------------------------------------------! +! Occupied part of the T-matrix self-energy ! +!-------------------------------------------! + + do i=nC+1,nO + + do cd=1,nVV + eps = w + e(i) - Om1(cd) + GGTpp_SigC = GGTpp_SigC + rho1(p,i,cd)**2*eps/(eps**2 + eta**2) + end do + + end do + +!------------------------------------------! +! Virtual part of the T-matrix self-energy ! +!------------------------------------------! + + do a=nO+1,nBas-nR + + do kl=1,nOO + eps = w + e(a) - Om2(kl) + GGTpp_SigC = GGTpp_SigC + rho2(p,a,kl)**2*eps/(eps**2 + eta**2) + end do + + end do + +end function diff --git a/src/GT/GGTpp_dSigC.f90 b/src/GT/GGTpp_dSigC.f90 new file mode 100644 index 0000000..22aa82e --- /dev/null +++ b/src/GT/GGTpp_dSigC.f90 @@ -0,0 +1,62 @@ +double precision function GGTpp_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho2) + +! Compute diagonal of the correlation part of the self-energy + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: p + double precision,intent(in) :: w + double precision,intent(in) :: eta + integer,intent(in) :: nBas + 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) :: e(nBas) + double precision,intent(in) :: Om1(nVV) + double precision,intent(in) :: rho1(nBas,nBas,nVV) + double precision,intent(in) :: Om2(nOO) + double precision,intent(in) :: rho2(nBas,nBas,nOO) + +! Local variables + + integer :: i,a,cd,kl + double precision :: eps + +! Initialize + + GGTpp_dSigC = 0d0 + +!-------------------------------------------! +! Occupied part of the T-matrix self-energy ! +!-------------------------------------------! + + do i=nC+1,nO + + do cd=1,nVV + eps = w + e(i) - Om1(cd) + GGTpp_dSigC = GGTpp_dSigC - rho1(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,nOO + eps = w + e(a) - Om2(kl) + GGTpp_dSigC = GGTpp_dSigC - rho2(p,a,kl)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + end do + + end do + +end function diff --git a/src/GT/GGTpp_excitation_density.f90 b/src/GT/GGTpp_excitation_density.f90 new file mode 100644 index 0000000..0550f25 --- /dev/null +++ b/src/GT/GGTpp_excitation_density.f90 @@ -0,0 +1,180 @@ +subroutine GGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2) + +! Compute excitation densities for T-matrix self-energy + + ! TODO + ! debug DGEMM for nC != 0 + ! and nR != 0 + + + + implicit none + +! Input variables + + integer,intent(in) :: ispin + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + integer,intent(in) :: nOO + integer,intent(in) :: nVV + double precision,intent(in) :: X1(nVV,nVV) + double precision,intent(in) :: Y1(nOO,nVV) + double precision,intent(in) :: X2(nVV,nOO) + double precision,intent(in) :: Y2(nOO,nOO) + +! Local variables + + integer :: i,j,k,l + integer :: a,b,c,d + integer :: p,q + integer :: ab,cd,ij,kl + double precision,external :: Kronecker_delta + +! Output variables + + double precision,intent(out) :: rho1(nBas,nBas,nVV) + double precision,intent(out) :: rho2(nBas,nBas,nOO) + + integer :: dim_1, dim_2 + double precision, allocatable :: ERI_1(:,:,:) + double precision, allocatable :: ERI_2(:,:,:) + +! Initialization + + rho1(:,:,:) = 0d0 + rho2(:,:,:) = 0d0 + + dim_1 = (nBas - nO) * (nBas - nO - 1) / 2 + dim_2 = nO * (nO - 1) / 2 + + if((dim_1 .eq. 0) .or. (dim_2 .eq. 0)) then + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(p, q, a, b, ab, c, d, cd, i, j, ij, k, l, kl) & + !$OMP SHARED(nC, nBas, nR, nO, rho1, rho2, ERI, X1, Y1, X2, Y2) + !$OMP DO COLLAPSE(2) + do q = nC+1, nBas-nR + do p = nC+1, nBas-nR + + ab = 0 + + do a = nO+1, nBas-nR + do b = a+1, nBas-nR + + ab = ab + 1 + + cd = 0 + do c = nO+1, nBas-nR + do d = c+1, nBas-nR + + cd = cd + 1 + + rho1(p,q,ab) = rho1(p,q,ab) & + + (ERI(p,q,c,d) - ERI(p,q,d,c))*X1(cd,ab) + end do ! d + end do ! c + + kl = 0 + do k = nC+1, nO + do l = k+1, nO + + kl = kl + 1 + + rho1(p,q,ab) = rho1(p,q,ab) & + + (ERI(p,q,k,l) - ERI(p,q,l,k))*Y1(kl,ab) + end do ! l + end do ! k + end do ! b + end do ! a + + ij = 0 + do i = nC+1, nO + do j = i+1, nO + + ij = ij + 1 + + cd = 0 + + do c = nO+1, nBas-nR + do d = c+1, nBas-nR + + cd = cd + 1 + + rho2(p,q,ij) = rho2(p,q,ij) & + + (ERI(p,q,c,d) - ERI(p,q,d,c))*X2(cd,ij) + end do ! d + end do ! c + + kl = 0 + do k = nC+1, nO + do l = k+1, nO + + kl = kl + 1 + + rho2(p,q,ij) = rho2(p,q,ij) & + + (ERI(p,q,k,l) - ERI(p,q,l,k))*Y2(kl,ij) + end do ! l + end do ! k + end do ! j + end do ! i + end do ! p + end do ! q + !$OMP END DO + !$OMP END PARALLEL + + else + + allocate(ERI_1(nBas,nBas,dim_1), ERI_2(nBas,nBas,dim_2)) + ERI_1 = 0.d0 + ERI_2 = 0.d0 + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(p, q, c, d, cd, k, l, kl) & + !$OMP SHARED(nC, nBas, nR, nO, ERI_1, ERI_2, ERI) + !$OMP DO COLLAPSE(2) + do q = nC+1, nBas-nR + do p = nC+1, nBas-nR + cd = 0 + do c = nO+1, nBas-nR + do d = c+1, nBas-nR + cd = cd + 1 + ERI_1(p,q,cd) = ERI(p,q,c,d) - ERI(p,q,d,c) + enddo + enddo + kl = 0 + do k = nC+1, nO + do l = k+1, nO + kl = kl + 1 + ERI_2(p,q,kl) = ERI(p,q,k,l) - ERI(p,q,l,k) + end do + end do + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + call dgemm("N", "N", nBas*nBas, dim_1, dim_1, 1.d0, & + ERI_1(1,1,1), nBas*nBas, X1(1,1), dim_1, & + 0.d0, rho1(1,1,1), nBas*nBas) + + call dgemm("N", "N", nBas*nBas, dim_1, dim_2, 1.d0, & + ERI_2(1,1,1), nBas*nBas, Y1(1,1), dim_2, & + 1.d0, rho1(1,1,1), nBas*nBas) + + call dgemm("N", "N", nBas*nBas, dim_2, dim_1, 1.d0, & + ERI_1(1,1,1), nBas*nBas, X2(1,1), dim_1, & + 0.d0, rho2(1,1,1), nBas*nBas) + + call dgemm("N", "N", nBas*nBas, dim_2, dim_2, 1.d0, & + ERI_2(1,1,1), nBas*nBas, Y2(1,1), dim_2, & + 1.d0, rho2(1,1,1), nBas*nBas) + + deallocate(ERI_1, ERI_2) + + endif + +end subroutine diff --git a/src/GT/GGTpp_self_energy_diag.f90 b/src/GT/GGTpp_self_energy_diag.f90 new file mode 100644 index 0000000..c34bb73 --- /dev/null +++ b/src/GT/GGTpp_self_energy_diag.f90 @@ -0,0 +1,105 @@ +subroutine GGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho2,EcGM,Sig,Z) + +! Compute diagonal of the correlation part of the T-matrix self-energy + + implicit none + include 'parameters.h' + +! Input variables + + double precision,intent(in) :: eta + integer,intent(in) :: nBas + 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) :: e(nBas) + double precision,intent(in) :: Om1(nVV) + double precision,intent(in) :: rho1(nBas,nBas,nVV) + double precision,intent(in) :: Om2(nOO) + double precision,intent(in) :: rho2(nBas,nBas,nOO) + +! Local variables + + integer :: i,j,a,b,p,cd,kl + double precision :: num,eps + +! Output variables + + double precision,intent(inout) :: EcGM + double precision,intent(inout) :: Sig(nBas) + double precision,intent(inout) :: Z(nBas) + +! Initialization + + Sig(:) = 0d0 + Z(:) = 0d0 + EcGM = 0d0 + +!--------------------------------------! +! Occupied part of the Tpp self-energy ! +!--------------------------------------! + + do p=nC+1,nBas-nR + do i=nC+1,nO + + do cd=1,nVV + eps = e(p) + e(i) - Om1(cd) + num = rho1(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 + + end do + end do + +!------------------------------------------! +! Virtual part of the T-matrix self-energy ! +!------------------------------------------! + + do p=nC+1,nBas-nR + do a=nO+1,nBas-nR + + do kl=1,nOO + eps = e(p) + e(a) - Om2(kl) + num = rho2(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 + + end do + end do + +!-------------------------------------! +! Galitskii-Migdal correlation energy ! +!-------------------------------------! + + do i=nC+1,nO + do j=nC+1,nO + + do cd=1,nVV + eps = e(i) + e(j) - Om1(cd) + num = rho1(i,j,cd)**2 + EcGM = EcGM + num*eps/(eps**2 + eta**2) + end do + + end do + end do + + do a=nO+1,nBas-nR + do b=nO+1,nBas-nR + + do kl=1,nOO + eps = e(a) + e(b) - Om2(kl) + num = rho2(a,b,kl)**2 + EcGM = EcGM - num*eps/(eps**2 + eta**2) + end do + + end do + end do + + Z(:) = 1d0/(1d0 - Z(:)) + +end subroutine diff --git a/src/GW/GG0W0.f90 b/src/GW/GG0W0.f90 index 9b1e30f..20451d5 100644 --- a/src/GW/GG0W0.f90 +++ b/src/GW/GG0W0.f90 @@ -207,8 +207,8 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@GG0W0 correlation energy =',EcBSE,' au' - write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@GG0W0 total energy =',ENuc + EGHF + EcBSE,' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@GHF correlation energy = ',EcBSE,' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@GHF total energy = ',ENuc + EGHF + EcBSE,' au' write(*,*)'-------------------------------------------------------------------------------' write(*,*)