G0W0+SOSEX (to be tested)

src/MCQC/G0W0.f90

@@ -18,7 +18,7 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, &
 
   logical                       :: dRPA
   integer                       :: ispin
-  double precision              :: EcRPA
+  double precision              :: EcRPA,EcGM
   double precision,allocatable  :: H(:,:),SigmaC(:),Z(:)
   double precision,allocatable  :: Omega(:,:),XpY(:,:,:),rho(:,:,:,:),rhox(:,:,:,:)
 
@@ -64,12 +64,12 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, &
 
 ! Compute correlation part of the self-energy 
 
-  call excitation_density(nBas,nC,nO,nR,nS,cHF,ERI_AO_basis,XpY(:,:,ispin),rho(:,:,:,ispin))
+  call excitation_density_from_MO(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin))
 
-  if(SOSEX) call excitation_density_SOSEX(nBas,nC,nO,nR,nS,cHF,ERI_AO_basis,XpY(:,:,ispin),rhox(:,:,:,ispin))
+  if(SOSEX) call excitation_density_SOSEX_from_MO(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rhox(:,:,:,ispin))
 
   call self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,eHF, & 
-                                    Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),SigmaC)
+                                    Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigmaC)
 
 ! COHSEX static approximation
 
@@ -90,7 +90,7 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, &
 ! Dump results
 
   call print_excitation('RPA  ',ispin,nS,Omega(:,ispin))
-  call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigmaC,Z,eG0W0,EcRPA)
+  call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigmaC,Z,eG0W0,EcRPA,EcGM)
 
 ! Plot stuff
 

src/MCQC/excitation_density_SOSEX_from_MO.f90

@@ -0,0 +1,35 @@
+subroutine excitation_density_SOSEX_from_MO(nBas,nC,nO,nR,nS,G,XpY,rho)
+
+! Compute excitation densities
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: nBas,nC,nO,nR,nS
+  double precision,intent(in)   :: G(nBas,nBas,nBas,nBas),XpY(nS,nS)
+
+! Local variables
+
+  integer                       :: ia,jb,x,y,j,b
+
+! Output variables
+
+  double precision,intent(out)  :: rho(nBas,nBas,nS)
+
+  rho(:,:,:) = 0d0
+  do ia=1,nS
+    do x=nC+1,nBas-nR
+      do y=nC+1,nBas-nR
+        jb = 0
+        do j=nC+1,nO
+          do b=nO+1,nBas-nR
+            jb = jb + 1
+            rho(x,y,ia) = rho(x,y,ia) + G(x,y,b,j)*XpY(ia,jb)
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+
+end subroutine excitation_density_SOSEX_from_MO

src/MCQC/excitation_density_from_MO.f90

@@ -0,0 +1,35 @@
+subroutine excitation_density_from_MO(nBas,nC,nO,nR,nS,G,XpY,rho)
+
+! Compute excitation densities
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: nBas,nC,nO,nR,nS
+  double precision,intent(in)   :: G(nBas,nBas,nBas,nBas),XpY(nS,nS)
+
+! Local variables
+
+  integer                       :: ia,jb,x,y,j,b
+
+! Output variables
+
+  double precision,intent(out)  :: rho(nBas,nBas,nS)
+
+  rho(:,:,:) = 0d0
+  do ia=1,nS
+    do x=nC+1,nBas-nR
+      do y=nC+1,nBas-nR
+        jb = 0
+        do j=nC+1,nO
+          do b=nO+1,nBas-nR
+            jb = jb + 1
+            rho(x,y,ia) = rho(x,y,ia) + G(x,y,j,b)*XpY(ia,jb)
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+
+end subroutine excitation_density_from_MO

src/MCQC/self_energy_correlation_diag.f90

@@ -64,6 +64,8 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega
       enddo
     enddo
 
+    ! GM correlation energy
+
     EcGM=0d0
     do i=nC+1,nO
       EcGM = EcGM + SigC(i)
@@ -105,7 +107,7 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega
       enddo
     enddo
 
-! GM correlation energy
+    ! GM correlation energy
 
     EcGM=0d0
     do i=nC+1,nO
@@ -153,6 +155,21 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega
         enddo
       enddo
 
+      ! GM correlation energy
+
+      do i=nC+1,nO
+        do a=nO+1,nBas-nR
+          jb = 0
+          do j=nC+1,nO
+            do b=nO+1,nBas-nR
+              jb = jb + 1
+              eps = e(a) - e(i) + Omega(jb)
+              EcGM = EcGM + 2d0*rho(a,i,jb)*rhox(a,i,jb)/eps
+            enddo
+          enddo
+        enddo
+      enddo
+
     endif
 
   endif