bug fix evGF2

input/methods

@@ -11,7 +11,7 @@
 # RPA RPAx ppRPA 
   F   F    F    
 # G0F2 evGF2 G0F3 evGF3
-  T    F     F    F
+  F    T     F    F
 # G0W0 evGW qsGW 
   F    F    F
 # G0T0 evGT qsGT 

src/QuAcK/G0F2.f90

@@ -66,6 +66,7 @@ subroutine G0F2(BSE,TDA,singlet_manifold,triplet_manifold,linearize,eta,nBas,nC,
 
           eps = eHF(p) + eHF(a) - eHF(i) - eHF(j)
           V  = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)
+
           Sig(p) = Sig(p) + V*eps/(eps**2 + eta**2)
           Z(p)   = Z(p)   - V*(eps**2 - eta**2)/(eps**2 + eta**2)**2
 
@@ -81,6 +82,7 @@ subroutine G0F2(BSE,TDA,singlet_manifold,triplet_manifold,linearize,eta,nBas,nC,
 
           eps = eHF(p) + eHF(i) - eHF(a) - eHF(b)
           V  = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)
+
           Sig(p) = Sig(p) + V*eps/(eps**2 + eta**2)
           Z(p)   = Z(p)   - V*(eps**2 - eta**2)/(eps**2 + eta**2)**2
 
@@ -111,51 +113,6 @@ subroutine G0F2(BSE,TDA,singlet_manifold,triplet_manifold,linearize,eta,nBas,nC,
 
     call BSE2(TDA,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGF2,EcBSE)
 
-    write(*,*)
-    write(*,*)'-------------------------------------------------------------------------------'
-    write(*,'(2X,A50,F20.10)') 'Tr@BSE2@G0F correlation energy (singlet) =',EcBSE(1)
-    write(*,'(2X,A50,F20.10)') 'Tr@BSE2@G0F correlation energy (triplet) =',EcBSE(2)
-    write(*,'(2X,A50,F20.10)') 'Tr@BSE2@G0F correlation energy           =',EcBSE(1) + EcBSE(2)
-    write(*,'(2X,A50,F20.10)') 'Tr@BSE2@G0F total energy                 =',ENuc + ERHF + EcBSE(1) + EcBSE(2)
-    write(*,*)'-------------------------------------------------------------------------------'
-    write(*,*)
-
-!   Compute the BSE correlation energy via the adiabatic connection 
-
-!   if(doACFDT) then
-
-!     write(*,*) '------------------------------------------------------'
-!     write(*,*) 'Adiabatic connection version of BSE correlation energy'
-!     write(*,*) '------------------------------------------------------'
-!     write(*,*) 
-
-!     if(doXBS) then 
-
-!       write(*,*) '*** scaled screening version (XBS) ***'
-!       write(*,*)
-
-!     end if
-
-!     call ACFDT(exchange_kernel,doXBS,.true.,TDA,BSE,singlet_manifold,triplet_manifold,eta, & 
-!                nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,Omega,XpY,XmY,rho,EcAC)
-
-!     if(exchange_kernel) then
-! 
-!       EcAC(1) = 0.5d0*EcAC(1)
-!       EcAC(2) = 1.5d0*EcAC(1)
-! 
-!     end if
-
-!     write(*,*)
-!     write(*,*)'-------------------------------------------------------------------------------'
-!     write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 correlation energy (singlet) =',EcAC(1)
-!     write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 correlation energy (triplet) =',EcAC(2)
-!     write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 correlation energy           =',EcAC(1) + EcAC(2)
-!     write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 total energy                 =',ENuc + ERHF + EcAC(1) + EcAC(2)
-!     write(*,*)'-------------------------------------------------------------------------------'
-!     write(*,*)
-
-!   end if
-
   end if
+
 end subroutine G0F2

src/QuAcK/evGF2.f90

@@ -1,5 +1,5 @@
 subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifold,linearize, & 
-                 eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,V,eHF)
+                 eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
 
 ! Perform eigenvalue self-consistent second-order Green function calculation
 
@@ -26,12 +26,14 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifol
   double precision,intent(in)   :: ENuc
   double precision,intent(in)   :: ERHF
   double precision,intent(in)   :: eHF(nBas)
-  double precision,intent(in)   :: V(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
 
 ! Local variables
 
   integer                       :: nSCF
   integer                       :: n_diis
+  double precision              :: EcBSE(nspin)
+  double precision              :: num
   double precision              :: eps
   double precision              :: Conv
   double precision              :: rcond
@@ -75,6 +77,7 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifol
     ! Frequency-dependent second-order contribution
 
     Sig(:) = 0d0
+    Z(:)   = 0d0
 
     do p=nC+1,nBas-nR
       do i=nC+1,nO
@@ -82,9 +85,10 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifol
           do a=nO+1,nBas-nR
 
             eps = eGF2(p) + eHF(a) - eHF(i) - eHF(j)
+            num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)
 
-            Sig(p) = Sig(p) &
+            Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
-                     + (2d0*V(p,a,i,j) - V(p,a,j,i))*V(p,a,i,j)*eps/(eps**2 + eta**2)
+            Z(p)   = Z(p)   - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
 
           end do
         end do
@@ -97,43 +101,10 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifol
           do b=nO+1,nBas-nR
 
             eps = eGF2(p) + eHF(i) - eHF(a) - eHF(b)
+            num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)
 
-            Sig(p) = Sig(p) &
+            Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
-                     + (2d0*V(p,i,a,b) - V(p,i,b,a))*V(p,i,a,b)*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
-    end do
-
-    ! Compute the renormalization factor
-
-    Z(:) = 0d0
-
-    do p=nC+1,nBas-nR
-      do i=nC+1,nO
-        do j=nC+1,nO
-          do a=nO+1,nBas-nR
-
-            eps = eGF2(p) + eHF(a) - eHF(i) - eHF(j)
-
-            Z(p) = Z(p) & 
-                 - (2d0*V(p,a,i,j) - V(p,a,j,i))*V(p,a,i,j)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
-
-          end do
-        end do
-      end do
-    end do
-
-    do p=nC+1,nBas-nR
-      do i=nC+1,nO
-        do a=nO+1,nBas-nR
-          do b=nO+1,nBas-nR
-
-            eps = eGF2(p) + eHF(i) - eHF(a) - eHF(b)
-
-            Z(p) = Z(p) & 
-                 - (2d0*V(p,i,a,b) - V(p,i,b,a))*V(p,i,a,b)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
 
           end do
         end do
@@ -156,7 +127,7 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifol
 
     ! Print results
 
-    call print_evGF2(nBas,nO,nSCF,Conv,eHF,eGF2)
+    call print_evGF2(nBas,nO,nSCF,Conv,eHF,Sig,Z,eGF2)
 
     ! DIIS extrapolation
 
@@ -188,4 +159,12 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifol
 
   end if
 
+! Perform BSE2 calculation
+
+  if(BSE) then
+
+    call BSE2(TDA,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGF2,EcBSE)
+
+  end if
+
 end subroutine evGF2