FC in pCCD

examples/molecule.F2

@@ -2,4 +2,4 @@
  2   9     9   0     0
 # Znuc   x            y           z
   F     0.           0.         0.
-  F     0.           0.         2
+  F     0.           0.         3.3

src/QuAcK/G0T0.f90

@@ -194,8 +194,6 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,*)
 
-!   Compute the BSE correlation energy via the adiabatic connection 
-
 !   Compute the BSE correlation energy via the adiabatic connection 
 
     if(doACFDT) then

src/QuAcK/QuAcK.f90

@@ -447,7 +447,7 @@ program QuAcK
   if(do_pCCD) then
 
     call cpu_time(start_CCD)
-    call pCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nO(1),nV(1),ERI_MO_basis,ENuc,ERHF,eHF)
+    call pCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO_basis,ENuc,ERHF,eHF)
     call cpu_time(end_CCD)
 
     t_CCD = end_CCD - start_CCD
@@ -513,7 +513,8 @@ program QuAcK
   if(doppRPA) then
 
     call cpu_time(start_ppRPA)
-    call ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO_basis,eHF)
+    call ppRPA(singlet_manifold,triplet_manifold, & 
+               nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF)
     call cpu_time(end_ppRPA)
 
     t_ppRPA = end_ppRPA - start_ppRPA
@@ -529,7 +530,8 @@ program QuAcK
   if(doADC) then
 
     call cpu_time(start_ADC)
-    call ADC(singlet_manifold,triplet_manifold,maxSCF_GF,thresh_GF,n_diis_GF,nBas,nC(1),nO(1),nV(1),nR(1),eHF,ERI_MO_basis)
+    call ADC(singlet_manifold,triplet_manifold,maxSCF_GF,thresh_GF,n_diis_GF, & 
+             nBas,nC(1),nO(1),nV(1),nR(1),eHF,ERI_MO_basis)
     call cpu_time(end_ADC)
 
     t_ADC = end_ADC - start_ADC

src/QuAcK/pCCD.f90

@@ -1,4 +1,4 @@
-subroutine pCCD(maxSCF,thresh,max_diis,nBas,nO,nV,ERI,ENuc,ERHF,eHF)
+subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
 
 ! pair CCD module
 
@@ -10,7 +10,7 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nO,nV,ERI,ENuc,ERHF,eHF)
   integer,intent(in)            :: max_diis
   double precision,intent(in)   :: thresh
 
-  integer,intent(in)            :: nBas,nO,nV
+  integer,intent(in)            :: nBas,nC,nO,nV,nR
   double precision,intent(in)   :: ENuc,ERHF
   double precision,intent(in)   :: eHF(nBas)
   double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
@@ -31,8 +31,7 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nO,nV,ERI,ENuc,ERHF,eHF)
   double precision,allocatable  :: OVVO(:,:)
   double precision,allocatable  :: VVVV(:,:)
 
-  double precision,allocatable  :: X(:,:)
+  double precision,allocatable  :: y(:,:)
-  double precision,allocatable  :: Y(:,:)
 
   double precision,allocatable  :: r(:,:)
   double precision,allocatable  :: t(:,:)
@@ -51,8 +50,8 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nO,nV,ERI,ENuc,ERHF,eHF)
 
   allocate(delta_OOVV(nO,nV))
 
-  do i=1,nO
+  do i=nC+1,nO
-    do a=1,nV
+    do a=1,nV-nR
       delta_OOVV(i,a) = 2d0*(eHF(nO+a) - eHF(i))
     enddo
   enddo
@@ -61,22 +60,22 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nO,nV,ERI,ENuc,ERHF,eHF)
 
   allocate(OOOO(nO,nO),OOVV(nO,nV),OVOV(nO,nV),OVVO(nO,nV),VVVV(nV,nV))
 
-  do i=1,nO
+  do i=nC+1,nO
-    do j=1,nO
+    do j=nC+1,nO
       OOOO(i,j) = ERI(i,i,j,j)
     end do
   end do
 
-  do i=1,nO
+  do i=nC+1,nO
-    do a=1,nV
+    do a=1,nV-nR
       OOVV(i,a) = ERI(i,i,nO+a,nO+a)
       OVOV(i,a) = ERI(i,nO+a,i,nO+a)
       OVVO(i,a) = ERI(i,nO+a,nO+a,i)
     end do
   end do
 
-  do a=1,nV
+  do a=1,nV-nR
-    do b=1,nV
+    do b=1,nV-nR
       VVVV(a,b) = ERI(nO+a,nO+a,nO+b,nO+b)
     end do
   end do
@@ -89,7 +88,7 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nO,nV,ERI,ENuc,ERHF,eHF)
 
 ! Initialization
 
-  allocate(r(nO,nV),X(nV,nV),Y(nO,nO))
+  allocate(r(nO,nV),y(nO,nO))
 
   Conv = 1d0
   nSCF = 0
@@ -113,22 +112,34 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nO,nV,ERI,ENuc,ERHF,eHF)
 
   ! Form intermediate array
     
-    X(:,:) = matmul(transpose(OOVV(:,:)),t(:,:))
+   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=1,nO
+    do i=nC+1,nO
-      do a=1,nV
+      do a=1,nV-nR
-        r(i,a) = - 2d0*(X(a,a) + Y(i,i))*t(i,a)
+
+        r(i,a) = OOVV(i,a) + delta_OOVV(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
+          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 
+
+        do b=1,nV-nR
+          r(i,a) = r(i,a) - 2d0*OOVV(i,b)*t(i,b)*t(i,a) + VVVV(a,b)*t(i,b)
+        end do 
+
       end do
     end do
 
-    r(:,:) = r(:,:) + OOVV(:,:) + delta_OOVV(:,:)*t(:,:)               &
-           - 2d0*(2d0*OVOV(:,:) - OVVO(:,:) - OOVV(:,:)*t(:,:))*t(:,:) &
-           + matmul(t(:,:),transpose(VVVV(:,:)))                       &
-           + matmul(transpose(OOOO(:,:)),t(:,:))                       &
-           + matmul(Y(:,:),t)
 
    ! Check convergence 
 
@@ -140,7 +151,12 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nO,nV,ERI,ENuc,ERHF,eHF)
 
    ! Compute correlation energy
 
-    EcCCD = trace_matrix(nO,matmul(t(:,:),transpose(OOVV(:,:))))
+    EcCCD = 0d0
+    do i=nC+1,nO
+      do a=1,nV-nR
+        EcCCD = EcCCD + t(i,a)*OOVV(i,a)
+      end do
+    end do
 
 !   Dump results