clean up pCCD

input/methods

@@ -1,11 +1,11 @@
 # RHF UHF MOM
-  F   T   F
-# MP2* MP3 MP2-F12
   T   F   F
+# MP2* MP3 MP2-F12
+  F   F   F
 # CCD CCSD CCSD(T) 
-  F   F    F
+  T   F    F
 # drCCD rCCD lCCD pCCD
-  F      F   F    F
+  F      F   F    T
 # CIS* CIS(D) CID CISD
   F    F      F   F    
 # RPA* RPAx* ppRPA 
@@ -13,7 +13,7 @@
 # G0F2 evGF2 G0F3 evGF3
   F    F     F    F
 # G0W0* evGW* qsGW 
-  T    F    F
+  F    F    F
 # G0T0 evGT qsGT 
   F    F    F
 # MCMP2

input/options

@@ -5,7 +5,7 @@
 # CC: maxSCF thresh   DIIS n_diis
       64     0.0000001  T    5
 # spin: TDA singlet triplet spin_conserved spin_flip 
-         F   T       T       T              T 
+         T   T       T       T              T 
 # GF: maxSCF thresh  DIIS n_diis lin eta renorm 
       256    0.00001 T    5      T   0.0  3      
 # GW/GT: maxSCF thresh  DIIS n_diis lin eta        COHSEX SOSEX TDA_W G0W GW0 
@@ -13,6 +13,6 @@
 # ACFDT: AC Kx  XBS
          T  F   T
 # BSE: BSE dBSE dTDA evDyn
-        T    F    T    F
+        T    T    T    F
 # MCMP2: nMC    nEq    nWalk dt  nPrint iSeed doDrift
          1000000 100000 10    0.3 10000 1234  T

src/CC/form_delta_OV.f90

@@ -7,8 +7,8 @@ subroutine form_delta_OV(nC,nO,nV,nR,eO,eV,delta)
 ! Input variables
 
   integer,intent(in)            :: nC,nO,nV,nR
-  double precision,intent(in)   :: eO(nO)
+  double precision,intent(in)   :: eO(nO-nC)
-  double precision,intent(in)   :: eV(nV)
+  double precision,intent(in)   :: eV(nV-nR)
 
 ! Local variables
 
@@ -18,7 +18,7 @@ subroutine form_delta_OV(nC,nO,nV,nR,eO,eV,delta)
 
   double precision,intent(out)  :: delta(nO,nV)
 
-    do i=nC+1,nO
+    do i=1,nO-nC
       do a=1,nV-nR
         delta(i,a) = eV(a) - eO(i)
       enddo

src/CC/pCCD.f90

@@ -23,7 +23,9 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
   double precision              :: Conv
   double precision              :: ECCD,EcCCD
 
-  double precision,allocatable  :: delta_OOVV(:,:)
+  double precision,allocatable  :: eO(:)
+  double precision,allocatable  :: eV(:)
+  double precision,allocatable  :: delta_OV(:,:)
 
   double precision,allocatable  :: OOOO(:,:)
   double precision,allocatable  :: OOVV(:,:)
@@ -40,6 +42,7 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
   double precision              :: rcond
   double precision,allocatable  :: error_diis(:,:)
   double precision,allocatable  :: t_diis(:,:)
+  double precision,external     :: trace_matrix
           
 ! Hello world
 
@@ -51,36 +54,28 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
 
 ! Form energy denominator
 
-  allocate(delta_OOVV(nO,nV))
+  allocate(eO(nO-nC),eV(nV-nR),delta_OV(nO-nC,nV-nR))
 
-  delta_OOVV(:,:) = 0d0
+  eO(:) = eHF(nC+1:nO)
+  eV(:) = eHF(nO+1:nBas-nR)
 
-  do i=nC+1,nO
+  call form_delta_OV(nC,nO,nV,nR,eO,eV,delta_OV)
-    do a=1,nV-nR
-      delta_OOVV(i,a) = 2d0*(eHF(nO+a) - eHF(i))
-    enddo
-  enddo
 
 ! Create integral batches
 
-  allocate(OOOO(nO,nO),OOVV(nO,nV),OVOV(nO,nV),OVVO(nO,nV),VVVV(nV,nV))
+  allocate(OOOO(nO-nC,nO-nC),OOVV(nO-nC,nV-nR),OVOV(nO-nC,nV-nR),OVVO(nO-nC,nV-nR),VVVV(nV-nR,nV-nR))
 
-  OOOO(:,:) = 0d0
+  do i=1,nO-nC
-  OOVV(:,:) = 0d0
+    do j=1,nO-nC
-  OVVO(:,:) = 0d0
+      OOOO(i,j) = ERI(nC+i,nC+i,nC+j,nC+j)
-  VVVV(:,:) = 0d0
-
-  do i=nC+1,nO
-    do j=nC+1,nO
-      OOOO(i,j) = ERI(i,i,j,j)
     end do
   end do
 
-  do i=nC+1,nO
+  do i=1,nO-nC
     do a=1,nV-nR
-      OOVV(i,a) = ERI(i,i,nO+a,nO+a)
+      OOVV(i,a) = ERI(nC+i,nC+i,nO+a,nO+a)
-      OVOV(i,a) = ERI(i,nO+a,i,nO+a)
+      OVOV(i,a) = ERI(nC+i,nO+a,nC+i,nO+a)
-      OVVO(i,a) = ERI(i,nO+a,nO+a,i)
+      OVVO(i,a) = ERI(nC+i,nO+a,nO+a,nC+i)
     end do
   end do
 
@@ -92,17 +87,17 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
 
 ! MP2 guess amplitudes
 
-  allocate(t(nO,nV))
+  allocate(t(nO-nC,nV-nR))
 
-  t(:,:) = - OOVV(:,:)/delta_OOVV(:,:)
+  t(:,:) = -0.5d0*OOVV(:,:)/delta_OV(:,:)
 
 ! Memory allocation for DIIS
 
-  allocate(error_diis(nO*nV,max_diis),t_diis(nO*nV,max_diis))
+  allocate(error_diis((nO-nC)*(nV-nR),max_diis),t_diis((nO-nC)*(nV-nR),max_diis))
 
 ! Initialization
 
-  allocate(r(nO,nV),y(nO,nO))
+  allocate(r(nO-nC,nV-nR),y(nO-nC,nO-nC))
 
   Conv = 1d0
   nSCF = 0
@@ -130,24 +125,17 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
 
   ! Form intermediate array
     
-   y(:,:) = 0d0
+   y(:,:) = matmul(t,transpose(OOVV))
-   do i=nC+1,nO
-     do j=nC+1,nO
-       do b=1,nV-nR
-         y(i,j) = y(i,j) + OOVV(j,b)*t(i,b)
-       end do
-     end do
-   end do
     
    ! Compute residual
 
-    do i=nC+1,nO
+    do i=1,nO-nC
       do a=1,nV-nR
 
-        r(i,a) = OOVV(i,a) + delta_OOVV(i,a)*t(i,a) & 
+        r(i,a) = OOVV(i,a) + 2d0*delta_OV(i,a)*t(i,a) & 
                - 2d0*(2d0*OVOV(i,a) - OVVO(i,a) - OOVV(i,a)*t(i,a))*t(i,a)
 
-        do j=nC+1,nO
+        do j=1,nO-nC
           r(i,a) = r(i,a) - 2d0*OOVV(j,a)*t(j,a)*t(i,a) + OOOO(j,i)*t(j,a) + y(i,j)*t(j,a) 
         end do 
 
@@ -164,16 +152,11 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
   
    ! Update amplitudes
 
-    t(:,:) = t(:,:) - r(:,:)/delta_OOVV(:,:)
+   t(:,:) = t(:,:) - 0.5d0*r(:,:)/delta_OV(:,:)
 
    ! Compute correlation energy
 
-    EcCCD = 0d0
+    EcCCD = trace_matrix(nO,matmul(t,transpose(OOVV)))
-    do i=nC+1,nO
-      do a=1,nV-nR
-        EcCCD = EcCCD + t(i,a)*OOVV(i,a)
-      end do
-    end do
 
    ! Dump results
 
@@ -182,7 +165,7 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
     ! DIIS extrapolation
 
     n_diis = min(n_diis+1,max_diis)
-    call DIIS_extrapolation(rcond,nO*nV,nO*nV,n_diis,error_diis,t_diis,-r/delta_OOVV,t)
+    call DIIS_extrapolation(rcond,nO*nV,nO*nV,n_diis,error_diis,t_diis,-0.5d0*r/delta_OV,t)
 
     !  Reset DIIS if required
 

src/QuAcK/QuAcK.f90

@@ -268,6 +268,14 @@ program QuAcK
 
   if(doRHF) then
 
+   ! Check that RHF calculation is worth doing...
+
+   if(nO(1) /= nO(2)) then
+     write(*,*) ' !!! The system does not appear to be closed shell !!!'
+     write(*,*) 
+     stop
+   end if
+
     call cpu_time(start_HF)
     call RHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,nNuc,ZNuc,rNuc,ENuc, &
              nBas,nO,S,T,V,Hc,ERI_AO,dipole_int,X,ERHF,eHF,cHF,PHF)

src/modules/.gitignore

@@ -0,0 +1 @@
+*.mod