fix bug in CCGW

input/methods

@@ -1,7 +1,7 @@
 # RHF UHF KS MOM 
   T   F   F  F   
 # MP2* MP3 
-  T   F   
+  F   F   
 # CCD pCCD DCD CCSD CCSD(T) 
   F   F    F   F    F
 # drCCD rCCD crCCD lCCD
@@ -13,7 +13,7 @@
 # G0F2* evGF2* qsGF2* G0F3 evGF3
   F     F      F      F    F
 # G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW
-  F     F     T     F        F      F
+  T     F     F     F        F      T
 # G0T0 evGT qsGT 
   F    F    F
 # * unrestricted version available

input/options

@@ -9,10 +9,10 @@
 # GF: maxSCF thresh  DIIS n_diis lin eta renorm reg
       256    0.00001 T    5      T   0.0 3      F
 # GW: maxSCF thresh  DIIS n_diis lin eta COHSEX SOSEX TDA_W reg
-      256     0.00001  T   5       T   0.01 F      F     F  F 
+      256     0.00001  T  5     T   0.01 F      F     F  F 
 # GT: maxSCF thresh  DIIS n_diis lin eta TDA_T reg
       10    0.00001  T    5      T   0.0 F     F  
 # ACFDT: AC Kx  XBS
          F  T   T
 # BSE: BSE dBSE dTDA evDyn ppBSE BSE2
-        F    F    T    F   F     T
+        F    F    T    F   F     F

int/ERI.Hu.dat

@@ -1,2 +1,2 @@
-1 1 1 1 1.
+1 1 1 1 5.
-2 2 2 2 1.
+2 2 2 2 5.

src/CC/CCGW.f90

@@ -29,9 +29,6 @@ subroutine CCGW(maxSCF,thresh,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,e)
   integer                       :: nSCF
   double precision              :: Conv
 
-  double precision,allocatable  :: eO(:)
-  double precision,allocatable  :: eV(:)
-
   double precision,allocatable  :: OVVO(:,:,:,:)
   double precision,allocatable  :: VOOV(:,:,:,:)
 
@@ -74,21 +71,16 @@ subroutine CCGW(maxSCF,thresh,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,e)
  
 ! Form energy denominator and guess amplitudes
 
-  allocate(eO(nO),eV(nV))
   allocate(delta_2h1p(nO,nO,nV,nBas),delta_2p1h(nO,nV,nV,nBas))
   allocate(V_2h1p(nBas,nO,nO,nV),V_2p1h(nBas,nO,nV,nV))
   allocate(t_2h1p(nO,nO,nV,nBas),t_2p1h(nO,nV,nV,nBas))
   allocate(x_2h1p(nBas,nBas),x_2p1h(nBas,nBas))
 
-  eO(:) = e(1:nO)
-  eV(:) = e(nO+1:nBas)
-
   do i=nC+1,nO
     do j=nC+1,nO
       do a=1,nV-nR
         do p=nC+1,nBas-nR
 
-          delta_2h1p(i,j,a,p) = eO(i) + eO(j) - eV(a) - e(p)
           V_2h1p(p,i,j,a) = sqrt(2d0)*ERI(p,nO+a,i,j)
 
         end do
@@ -101,7 +93,6 @@ subroutine CCGW(maxSCF,thresh,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,e)
       do b=1,nV-nR
         do p=nC+1,nBas-nR
 
-          delta_2p1h(i,a,b,p) = eV(a) + eV(b) - eO(i) - e(p)
           V_2p1h(p,i,a,b) = sqrt(2d0)*ERI(p,i,nO+b,nO+a)
 
         end do
@@ -117,6 +108,7 @@ subroutine CCGW(maxSCF,thresh,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,e)
 
   Conv = 1d0
   nSCF = 0
+  eGW(:) = e(:)
 
   t_2h1p(:,:,:,:) = 0d0
   t_2p1h(:,:,:,:) = 0d0
@@ -138,6 +130,32 @@ subroutine CCGW(maxSCF,thresh,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,e)
 
     nSCF = nSCF + 1
 
+! Compute energy differences
+
+    do i=nC+1,nO
+      do j=nC+1,nO
+        do a=1,nV-nR
+          do p=nC+1,nBas-nR
+ 
+            delta_2h1p(i,j,a,p) = eGW(i) + eGW(j) - eGW(nO+a) - e(p)
+ 
+          end do
+        end do
+      end do
+    end do
+
+    do i=nC+1,nO
+      do a=1,nV-nR
+        do b=1,nV-nR
+          do p=nC+1,nBas-nR
+ 
+            delta_2p1h(i,a,b,p) = eGW(nO+a) + eGW(nO+b) - eGW(i) - e(p)
+ 
+          end do
+        end do
+      end do
+    end do
+
 !   Compute intermediates
 
     x_2h1p(:,:) = 0d0 
@@ -278,6 +296,8 @@ subroutine CCGW(maxSCF,thresh,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,e)
       end do
     end do
 
+!   Diagonalize non-Hermitian matrix
+
     call diagonalize_general_matrix(nBas,SigGW,eGW,cGW)
 
     do p=1,nBas

src/QuAcK/QuAcK.f90

@@ -984,7 +984,7 @@ program QuAcK
       else
         call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE,singlet,triplet, &
                   linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0)
-!       call soG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE,singlet,triplet, &
+!       call ehTM(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE,singlet,triplet, &
 !                 linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0)
       end if