GM in GT

input/dft

@@ -13,7 +13,7 @@
 # GGA          = 2: LYP,PBE
 # MGGA         = 3:
 # Hybrid       = 4: HF,B3LYP,PBE 
-0 H
+1 VWN5
 # quadrature grid SG-n
   1
 # Number of states in ensemble (nEns)
@@ -31,7 +31,7 @@
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 # Ensemble weights: wEns(1),...,wEns(nEns-1)
-  1.0 0.0 0.0 
+  0.6 0.0 0.0 
 # Ncentered ? 
 F
 # Parameters for CC weight-dependent exchange functional

input/methods

@@ -1,5 +1,5 @@
 # RHF UHF KS MOM 
-   F   F   T  F   
+  T   F   F  F   
 # MP2* MP3 MP2-F12
   F   F   F
 # CCD pCCD DCD CCSD CCSD(T) 
@@ -15,7 +15,7 @@
 # G0W0* evGW* qsGW* ufG0W0 ufGW
   F    F    F       F      F
 # G0T0 evGT qsGT 
-  F    F    F
+  F    F    T
 # MCMP2
   F
 # * unrestricted version available

input/options

@@ -1,5 +1,5 @@
 # HF: maxSCF thresh   DIIS n_diis guess_type ortho_type mix_guess stability
-      1024   0.00001  T    5      1          1          T         F
+      1024   0.00001  T    5      1          1          F         F
 # MP: 
       
 # CC: maxSCF thresh   DIIS n_diis
@@ -15,6 +15,6 @@
 # ACFDT: AC Kx  XBS
          F  F   F
 # BSE: BSE dBSE dTDA evDyn
-        T    T    T    F
+        F    T    T    F
 # MCMP2: nMC    nEq    nWalk dt  nPrint iSeed doDrift
          1000000 100000 10    0.3 10000 1234  T

mol/h2.xyz

@@ -1,4 +1,4 @@
 2
  
 H 0. 0. 0.
-H 0. 0. 2.000000
+H 0. 0. 0.7

src/GT/G0T0.f90

@@ -48,6 +48,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
   double precision              :: EcRPA(nspin)
   double precision              :: EcBSE(nspin)
   double precision              :: EcAC(nspin)
+  double precision              :: EcGM
   double precision,allocatable  :: Omega1s(:),Omega1t(:)
   double precision,allocatable  :: X1s(:,:),X1t(:,:)
   double precision,allocatable  :: Y1s(:,:),Y1t(:,:)
@@ -134,6 +135,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
 ! Compute T-matrix version of the self-energy 
 !----------------------------------------------
 
+  EcGM    = 0d0
   SigT(:) = 0d0
   Z(:)    = 0d0
 
@@ -141,7 +143,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
 
   call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI_MO,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
 
-  call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF,Omega1s,rho1s,Omega2s,rho2s,SigT)
+  call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF,Omega1s,rho1s,Omega2s,rho2s,EcGM,SigT)
 
   call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF,Omega1s,rho1s,Omega2s,rho2s,Z)
 
@@ -149,7 +151,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
 
   call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI_MO,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
 
-  call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF,Omega1t,rho1t,Omega2t,rho2t,SigT)
+  call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF,Omega1t,rho1t,Omega2t,rho2t,EcGM,SigT)
 
   call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF,Omega1t,rho1t,Omega2t,rho2t,Z)
 
@@ -194,7 +196,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
   EcRPA(1) = EcRPA(1) - EcRPA(2)
   EcRPA(2) = 3d0*EcRPA(2)
 
-  call print_G0T0(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eG0T0,EcRPA)
+  call print_G0T0(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eG0T0,EcGM,EcRPA)
 
 ! Perform BSE calculation
 

src/GT/evGT.f90

@@ -55,6 +55,7 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
   integer                       :: iblock
   integer                       :: nOOs,nOOt
   integer                       :: nVVs,nVVt
+  double precision              :: EcGM
   double precision              :: EcRPA(nspin)
   double precision              :: EcBSE(nspin)
   double precision              :: EcAC(nspin)
@@ -148,10 +149,14 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
     call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT,ERI_MO,  & 
                           Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin))
 
+    EcRPA(1) = EcRPA(1) - EcRPA(2)
+    EcRPA(2) = 3d0*EcRPA(2)
+
   !----------------------------------------------
   ! Compute T-matrix version of the self-energy 
   !----------------------------------------------
 
+    EcGM    = 0d0
     SigT(:) = 0d0
     Z(:)    = 0d0
  
@@ -161,7 +166,7 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
                                     X1s,Y1s,rho1s,X2s,Y2s,rho2s)
  
     call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, & 
-                                  Omega1s,rho1s,Omega2s,rho2s,SigT)
+                                  Omega1s,rho1s,Omega2s,rho2s,EcGM,SigT)
  
     call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, & 
                                         Omega1s,rho1s,Omega2s,rho2s,Z)
@@ -172,7 +177,7 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
                                     X1t,Y1t,rho1t,X2t,Y2t,rho2t)
  
     call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, & 
-                                  Omega1t,rho1t,Omega2t,rho2t,SigT)
+                                  Omega1t,rho1t,Omega2t,rho2t,EcGM,SigT)
  
     call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, & 
                                         Omega1t,rho1t,Omega2t,rho2t,Z)
@@ -195,7 +200,7 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
   ! Dump results
   !----------------------------------------------
 
-    call print_evGT(nBas,nO,nSCF,Conv,eHF,SigT,Z,eGT)
+    call print_evGT(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigT,Z,eGT,EcGM,EcRPA)
 
     ! DIIS extrapolation
 
@@ -219,29 +224,6 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
 ! End main loop
 !------------------------------------------------------------------------
 
-! Compute the ppRPA correlation energy
-
-  ispin  = 1
-  iblock = 3
-  call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT,ERI_MO,  &
-                          Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin))
-  ispin  = 2
-  iblock = 4
-  call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT,ERI_MO,  &
-                          Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin))
-  EcRPA(1) = EcRPA(1) - EcRPA(2)
-  EcRPA(2) = 3d0*EcRPA(2)
-
-  write(*,*)
-  write(*,*)'-------------------------------------------------------------------------------'
-  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@evGT correlation energy (singlet) =',EcRPA(1)
-  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@evGT correlation energy (triplet) =',EcRPA(2)
-  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@evGT correlation energy           =',EcRPA(1) + EcRPA(2)
-  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@evGT total energy                 =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
-  write(*,*)'-------------------------------------------------------------------------------'
-  write(*,*)
-
-
 ! Perform BSE calculation
 
   if(BSE) then

src/GT/print_G0T0.f90

@@ -1,13 +1,15 @@
-subroutine print_G0T0(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eGT,EcRPA)
+subroutine print_G0T0(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eGT,EcGM,EcRPA)
 
 ! Print one-electron energies and other stuff for G0T0
 
   implicit none
   include 'parameters.h'
 
-  integer,intent(in)                 :: nBas,nO
+  integer,intent(in)                 :: nBas
+  integer,intent(in)                 :: nO
   double precision,intent(in)        :: ENuc
   double precision,intent(in)        :: ERHF
+  double precision,intent(in)        :: EcGM
   double precision,intent(in)        :: EcRPA(nspin)
   double precision,intent(in)        :: eHF(nBas)
   double precision,intent(in)        :: SigT(nBas)
@@ -50,10 +52,12 @@ subroutine print_G0T0(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eGT,EcRPA)
   write(*,'(2X,A50,F15.6,A3)') 'G0T0 LUMO      energy (eV)            =',eGT(LUMO)*HaToeV,' eV'
   write(*,'(2X,A50,F15.6,A3)') 'G0T0 HOMO-LUMO gap    (eV)            =',Gap*HaToeV,' eV'
   write(*,*)'-------------------------------------------------------------------------------'
-  write(*,'(2X,A50,F20.10,A3)') 'Tr@RPA@G0T0 correlation energy (singlet) =',EcRPA(1),' au'
+  write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@G0T0 correlation energy (singlet) =',EcRPA(1),' au'
-  write(*,'(2X,A50,F20.10,A3)') 'Tr@RPA@G0T0 correlation energy (triplet) =',EcRPA(2),' au'
+  write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@G0T0 correlation energy (triplet) =',EcRPA(2),' au'
-  write(*,'(2X,A50,F20.10,A3)') 'Tr@RPA@G0T0 correlation energy           =',EcRPA(1) + EcRPA(2),' au'
+  write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@G0T0 correlation energy           =',EcRPA(1) + EcRPA(2),' au'
-  write(*,'(2X,A50,F20.10,A3)') 'Tr@RPA@G0T0 total energy                 =',ENuc + ERHF + EcRPA(1) + EcRPA(2),' au'
+  write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@G0T0 total energy                 =',ENuc + ERHF + EcRPA(1) + EcRPA(2),' au'
+  write(*,'(2X,A50,F20.10,A3)') '       GM@G0T0 correlation energy           =',EcGM,' au'
+  write(*,'(2X,A50,F20.10,A3)') '       GM@G0T0 total energy                 =',ENuc + ERHF + EcGM,' au'
   write(*,*)'-------------------------------------------------------------------------------'
   write(*,*)
 

src/GT/print_evGT.f90

@@ -1,4 +1,4 @@
-subroutine print_evGT(nBas,nO,nSCF,Conv,eHF,SigT,Z,eGT)
+subroutine print_evGT(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigT,Z,eGT,EcGM,EcRPA)
 
 ! Print one-electron energies and other stuff for evGT
 
@@ -10,6 +10,10 @@ subroutine print_evGT(nBas,nO,nSCF,Conv,eHF,SigT,Z,eGT)
   integer,intent(in)                 :: nSCF
   double precision,intent(in)        :: Conv
   double precision,intent(in)        :: eHF(nBas)
+  double precision,intent(in)        :: ENuc
+  double precision,intent(in)        :: ERHF
+  double precision,intent(in)        :: EcGM
+  double precision,intent(in)        :: EcRPA(nspin)
   double precision,intent(in)        :: SigT(nBas)
   double precision,intent(in)        :: Z(nBas)
   double precision,intent(in)        :: eGT(nBas)
@@ -49,6 +53,13 @@ subroutine print_evGT(nBas,nO,nSCF,Conv,eHF,SigT,Z,eGT)
   write(*,'(2X,A30,F15.6)') 'evGT LUMO      energy (eV):',eGT(LUMO)*HaToeV
   write(*,'(2X,A30,F15.6)') 'evGT HOMO-LUMO gap    (eV):',Gap*HaToeV
   write(*,*)'-------------------------------------------------------------------------------'
+  write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@G0T0 correlation energy (singlet) =',EcRPA(1),' au'
+  write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@G0T0 correlation energy (triplet) =',EcRPA(2),' au'
+  write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@G0T0 correlation energy           =',EcRPA(1) + EcRPA(2),' au'
+  write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@G0T0 total energy                 =',ENuc + ERHF + EcRPA(1) + EcRPA(2),' au'
+  write(*,'(2X,A50,F20.10,A3)') '       GM@G0T0 correlation energy           =',EcGM,' au'
+  write(*,'(2X,A50,F20.10,A3)') '       GM@G0T0 total energy                 =',ENuc + ERHF + EcGM,' au'
+  write(*,*)'-------------------------------------------------------------------------------'
   write(*,*)
 
 end subroutine print_evGT

src/GT/print_qsGT.f90

@@ -48,9 +48,9 @@ subroutine print_qsGT(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ,Ex
 
   write(*,*)'-------------------------------------------------------------------------------'
   if(nSCF < 10) then
-    write(*,'(1X,A21,I1,A1,I1,A12)')'  Self-consistent qsG',nSCF,'W',nSCF,' calculation'
+    write(*,'(1X,A21,I1,A1,I1,A12)')'  Self-consistent qsG',nSCF,'T',nSCF,' calculation'
   else
-    write(*,'(1X,A21,I2,A1,I2,A12)')'  Self-consistent qsG',nSCF,'W',nSCF,' calculation'
+    write(*,'(1X,A21,I2,A1,I2,A12)')'  Self-consistent qsG',nSCF,'T',nSCF,' calculation'
   endif
   write(*,*)'-------------------------------------------------------------------------------'
   write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
@@ -73,7 +73,7 @@ subroutine print_qsGT(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ,Ex
   write(*,'(2X,A30,F15.6,A3)') '      qsGT total       energy:',ENuc + EqsGT,' au'
   write(*,'(2X,A30,F15.6,A3)') '      qsGT exchange    energy:',Ex,' au'
   write(*,'(2X,A30,F15.6,A3)') '   GM@qsGT correlation energy:',EcGM,' au'
-  write(*,'(2X,A30,F15.6,A3)') 'RPA@qsGT correlation energy:',sum(EcRPA(:)),' au'
+  write(*,'(2X,A30,F15.6,A3)') 'ppRPA@qsGT correlation energy:',sum(EcRPA(:)),' au'
   write(*,*)'-------------------------------------------'
   write(*,*)
 

src/GT/qsGT.f90

@@ -198,8 +198,12 @@ subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,T
     call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT,ERI_MO,  &
                           Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin))
 
+    EcRPA(1) = EcRPA(1) - EcRPA(2)
+    EcRPA(2) = 3d0*EcRPA(2)
+
     ! Compute correlation part of the self-energy 
 
+    EcGM      = 0d0
     SigT(:,:) = 0d0
     Z(:)      = 0d0
 
@@ -209,7 +213,7 @@ subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,T
                                     X1s,Y1s,rho1s,X2s,Y2s,rho2s)
 
     call self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, &
-                             Omega1s,rho1s,Omega2s,rho2s,SigT)
+                             Omega1s,rho1s,Omega2s,rho2s,EcGM,SigT)
 
     call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, &
                                         Omega1s,rho1s,Omega2s,rho2s,Z)
@@ -220,7 +224,7 @@ subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,T
                                     X1t,Y1t,rho1t,X2t,Y2t,rho2t)
 
     call self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, &
-                             Omega1t,rho1t,Omega2t,rho2t,SigT)
+                             Omega1t,rho1t,Omega2t,rho2t,EcGM,SigT)
 
     call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, &
                                         Omega1t,rho1t,Omega2t,rho2t,Z)
@@ -302,28 +306,6 @@ subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,T
 ! End main loop
 !------------------------------------------------------------------------
 
-! Compute the ppRPA correlation energy
-
-  ispin  = 1
-  iblock = 3
-  call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT,ERI_MO,  &
-                          Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin))
-  ispin  = 2
-  iblock = 4
-  call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT,ERI_MO,  &
-                          Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin))
-  EcRPA(1) = EcRPA(1) - EcRPA(2)
-  EcRPA(2) = 3d0*EcRPA(2)
-
-  write(*,*)
-  write(*,*)'-------------------------------------------------------------------------------'
-  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@qsGT correlation energy (singlet) =',EcRPA(1)
-  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@qsGT correlation energy (triplet) =',EcRPA(2)
-  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@qsGT correlation energy           =',EcRPA(1) + EcRPA(2)
-  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@qsGT total energy                 =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
-  write(*,*)'-------------------------------------------------------------------------------'
-  write(*,*)
-
 ! Did it actually converge?
 
   if(nSCF == maxSCF+1) then

src/GT/self_energy_Tmatrix.f90

@@ -1,4 +1,4 @@
-subroutine self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,SigT)
+subroutine self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,EcGM,SigT)
 
 ! Compute the correlation part of the T-matrix self-energy
 
@@ -23,11 +23,12 @@ subroutine self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2
 
 ! Local variables
 
-  integer                       :: i,a,p,q,cd,kl
+  integer                       :: i,j,a,b,p,q,cd,kl
   double precision              :: eps
 
 ! Output variables
 
+  double precision,intent(inout):: EcGM
   double precision,intent(inout):: SigT(nBas,nBas)
 
 !----------------------------------------------
@@ -60,4 +61,26 @@ subroutine self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2
     enddo
   enddo
 
+!----------------------------------------------
+! Galitskii-Migdal correlation energy
+!----------------------------------------------
+
+  do i=nC+1,nO
+    do j=nC+1,nO
+      do cd=1,nVV
+        eps = e(i) + e(j) - Omega2(cd)
+        EcGM = EcGM - rho1(i,j,cd)*rho1(i,j,cd)*eps/(eps**2 + eta**2)
+      enddo
+    enddo
+  enddo
+
+  do a=nO+1,nBas-nR
+    do b=nO+1,nBas-nR
+      do kl=1,nOO
+        eps = e(a) + e(b) - Omega2(kl)
+        EcGM = EcGM + rho2(a,b,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2)
+      enddo
+    enddo
+  enddo
+
 end subroutine self_energy_Tmatrix

src/GT/self_energy_Tmatrix_diag.f90

@@ -1,4 +1,4 @@
-subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,SigT)
+subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,EcGM,SigT)
 
 ! Compute diagonal of the correlation part of the T-matrix self-energy
 
@@ -23,11 +23,12 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,O
 
 ! Local variables
 
-  integer                       :: i,a,p,cd,kl
+  integer                       :: i,j,a,b,p,cd,kl
   double precision              :: eps
 
 ! Output variables
 
+  double precision,intent(inout)  :: EcGM
   double precision,intent(inout)  :: SigT(nBas)
 
 !----------------------------------------------
@@ -56,4 +57,26 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,O
     enddo
   enddo
 
+!----------------------------------------------
+! Galitskii-Migdal correlation energy
+!----------------------------------------------
+
+  do i=nC+1,nO
+    do j=nC+1,nO
+      do cd=1,nVV
+        eps = e(i) + e(j) - Omega2(cd)
+        EcGM = EcGM - rho1(i,j,cd)*rho1(i,j,cd)*eps/(eps**2 + eta**2)
+      enddo
+    enddo
+  enddo
+
+  do a=nO+1,nBas-nR
+    do b=nO+1,nBas-nR
+      do kl=1,nOO
+        eps = e(a) + e(b) - Omega2(kl)
+        EcGM = EcGM + rho2(a,b,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2)
+      enddo
+    enddo
+  enddo
+
 end subroutine self_energy_Tmatrix_diag

src/GW/self_energy_correlation.f90

@@ -28,8 +28,8 @@ subroutine self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Ec
 
 ! Output variables
 
-  double precision,intent(out)  :: SigC(nBas,nBas)
   double precision,intent(out)  :: EcGM
+  double precision,intent(out)  :: SigC(nBas,nBas)
 
 ! Initialize 
 
@@ -102,7 +102,7 @@ subroutine self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Ec
       end do
     end do
 
-    ! GM correlation energy
+    ! Galitskii-Migdal correlation energy
 
     EcGM = 0d0
     do i=nC+1,nO