implement 1RDM and 2RDM for pCCD

int/CAP.dat

@@ -1,16 +0,0 @@
-   1  1  9.1642021581097924E-03  6.2961947849362709E-02 9.1642021581097941E-03
-   1  2  2.9798815568270971E-02  1.0031339688416364E-01 2.9798815568270971E-02
-   1  3  4.8078353659559226E-03  5.1255302523161485E-03 4.8078353659559234E-03
-   1  4  2.3003539814844435E-02  4.1290024754715535E-02 2.3003539814844435E-02
-   2  1  2.9798815568270971E-02  1.0031339688416364E-01 2.9798815568270971E-02
-   2  2  3.5629639141443131E-01  5.7428563627799001E-01 3.5629639141443131E-01
-   2  3  2.3003539814844435E-02  4.1290024754715576E-02 2.3003539814844435E-02
-   2  4  3.0301481386007040E-01  3.0301481386007040E-01 3.0301481386007040E-01
-   3  1  4.8078353659559226E-03  5.1255302523161485E-03 4.8078353659559234E-03
-   3  2  2.3003539814844435E-02  4.1290024754715576E-02 2.3003539814844435E-02
-   3  3  9.1642021581097924E-03  6.2961947849362682E-02 9.1642021581097941E-03
-   3  4  2.9798815568270971E-02  1.0031339688416376E-01 2.9798815568270971E-02
-   4  1  2.3003539814844435E-02  4.1290024754715535E-02 2.3003539814844435E-02
-   4  2  3.0301481386007040E-01  3.0301481386007040E-01 3.0301481386007040E-01
-   4  3  2.9798815568270971E-02  1.0031339688416376E-01 2.9798815568270971E-02
-   4  4  3.5629639141443131E-01  5.7428563627799034E-01 3.5629639141443131E-01

mol/H2.xyz

@@ -1,4 +1,4 @@
 2
 
-H  0.0000  0.0000  -0.37500000
+H     0.00000000      -0.37500000       0.00000000
-H  0.0000  0.0000  0.37500000
+H     0.00000000       0.37500000       0.00000000

mol/N2.xyz

@@ -1,4 +1,4 @@
 2
 
-N  0.0000  0.0000  0.0000
+N  0.0000  0.0000  -0.5475132
-N  0.0000  0.0000  1.1007
+N  0.0000  0.0000  0.5475132

src/CC/pCCD.f90

@@ -35,15 +35,26 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF
   double precision,allocatable  :: OVVO(:,:)
   double precision,allocatable  :: VVVV(:,:)
 
-  double precision,allocatable  :: y(:,:)
+  double precision,allocatable  :: yO(:,:),yV(:,:)
 
   double precision,allocatable  :: r(:,:)
   double precision,allocatable  :: t(:,:)
+  double precision,allocatable  :: z(:,:)
 
+  double precision,allocatable  :: rdm1(:,:)
+  double precision,allocatable  :: rdm2(:,:,:,:)
+  double precision,allocatable  :: xOO(:,:)
+  double precision,allocatable  :: xVV(:,:)
+  double precision,allocatable  :: xOV(:,:)
+  double precision              :: tr_1rdm
+  double precision              :: tr_2rdm
+
+  integer                       :: O,V
   integer                       :: n_diis
   double precision              :: rcond
-  double precision,allocatable  :: error_diis(:,:)
+  double precision,allocatable  :: err_diis(:,:)
   double precision,allocatable  :: t_diis(:,:)
+  double precision,allocatable  :: z_diis(:,:)
   double precision,external     :: trace_matrix
           
 ! Hello world
@@ -54,9 +65,14 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF
   write(*,*)'**************************************'
   write(*,*)
 
+! Useful quantities
+
+  O = nO - nC
+  V = nV - NR
+
 ! Form energy denominator
 
-  allocate(eO(nO-nC),eV(nV-nR),delta_OV(nO-nC,nV-nR))
+  allocate(eO(O),eV(V),delta_OV(O,V))
 
   eO(:) = eHF(nC+1:nO)
   eV(:) = eHF(nO+1:nBas-nR)
@@ -65,62 +81,56 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF
 
 ! Create integral batches
 
-  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))
+  allocate(OOOO(O,O),OOVV(O,V),OVOV(O,V),OVVO(O,V),VVVV(V,V))
 
-  do i=1,nO-nC
+  do i=1,O
-    do j=1,nO-nC
+    do j=1,O
       OOOO(i,j) = ERI(nC+i,nC+i,nC+j,nC+j)
     end do
   end do
 
-  do i=1,nO-nC
+  do i=1,O
-    do a=1,nV-nR
+    do a=1,V
       OOVV(i,a) = ERI(nC+i,nC+i,nO+a,nO+a)
       OVOV(i,a) = ERI(nC+i,nO+a,nC+i,nO+a)
       OVVO(i,a) = ERI(nC+i,nO+a,nO+a,nC+i)
     end do
   end do
 
-  do a=1,nV-nR
+  do a=1,V
-    do b=1,nV-nR
+    do b=1,V
       VVVV(a,b) = ERI(nO+a,nO+a,nO+b,nO+b)
     end do
   end do
 
-! MP2 guess amplitudes
+! Initialization
 
-  allocate(t(nO-nC,nV-nR))
+  allocate(t(O,V),r(O,V),yO(O,O),yV(V,V))
-
-  t(:,:) = -0.5d0*OOVV(:,:)/delta_OV(:,:)
-
-  EcCC = 0d0
-  do i=1,nO-nC
-    do a=1,nV-nR
-      EcCC = EcCC + OOVV(i,a)*t(i,a)
-    end do
-  end do
 
 ! Memory allocation for DIIS
 
-  allocate(error_diis((nO-nC)*(nV-nR),max_diis),t_diis((nO-nC)*(nV-nR),max_diis))
+  allocate(err_diis(O*V,max_diis),t_diis(O*V,max_diis))
 
-! Initialization
+!------------------------------------------------------------------------
-
+! Compute t ampltiudes
-  allocate(r(nO-nC,nV-nR),y(nO-nC,nO-nC))
+!------------------------------------------------------------------------
 
   Conv = 1d0
   nSCF = 0
+  ECC  = ERHF
+  EcCC = 0d0
 
-  n_diis          = 0
+  n_diis        = 0
-  t_diis(:,:)     = 0d0
+  t(:,:)        = 0d0
-  error_diis(:,:) = 0d0
+  t_diis(:,:)   = 0d0
+  err_diis(:,:) = 0d0
 
 !------------------------------------------------------------------------
 ! Main SCF loop
 !------------------------------------------------------------------------
   write(*,*)
   write(*,*)'----------------------------------------------------'
-  write(*,*)'| pair CCD calculation                             |'
+  write(*,*)'| pCCD calculation: t amplitudes                   |'
   write(*,*)'----------------------------------------------------'
   write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') &
             '|','#','|','E(pCCD)','|','Ec(pCCD)','|','Conv','|'
@@ -134,28 +144,22 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF
 
   ! Form intermediate array
     
-   y(:,:) = 0d0
+   yO(:,:) = matmul(t,transpose(OOVV))
-   do i=1,nO-nC
-     do j=1,nO-nC
-       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-nC
-      do a=1,nV-nR
 
-        r(i,a) = OOVV(i,a) + 2d0*delta_OV(i,a)*t(i,a) & 
+    r(:,:) = OOVV(:,:) + 2d0*delta_OV(:,:)*t(:,:) & 
-               - 2d0*(2d0*OVOV(i,a) - OVVO(i,a) - OOVV(i,a)*t(i,a))*t(i,a)
+           - 2d0*(2d0*OVOV(:,:) - OVVO(:,:) - OOVV(:,:)*t(:,:))*t(:,:)
 
-        do j=1,nO-nC
+    do i=1,O
-          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) 
+      do a=1,V
+
+        do j=1,O
+          r(i,a) = r(i,a) - 2d0*OOVV(j,a)*t(j,a)*t(i,a) + OOOO(j,i)*t(j,a) + yO(i,j)*t(j,a) 
         end do 
 
-        do b=1,nV-nR
+        do b=1,V
           r(i,a) = r(i,a) - 2d0*OOVV(i,b)*t(i,b)*t(i,a) + VVVV(a,b)*t(i,b)
         end do 
 
@@ -172,25 +176,20 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF
 
    ! Compute correlation energy
 
-    EcCC = 0d0
+    EcCC = trace_matrix(V,matmul(transpose(OOVV),t))
-    do i=1,nO-nC    
-      do a=1,nV-nR
-        EcCC = EcCC + OOVV(i,a)*t(i,a)
-     end do
-   end do
 
    ! Dump results
 
     ECC = ERHF + EcCC
 
-    ! DIIS extrapolation
+   ! DIIS extrapolation
 
-!   n_diis = min(n_diis+1,max_diis)
+   if(max_diis > 1) then
-!   call DIIS_extrapolation(rcond,nO*nV,nO*nV,n_diis,error_diis,t_diis,-0.5d0*r/delta_OV,t)
 
-    !  Reset DIIS if required
+     n_diis = min(n_diis+1,max_diis)
+     call DIIS_extrapolation(rcond,nO*nV,nO*nV,n_diis,err_diis,t_diis,-0.5d0*r/delta_OV,t)
 
-!   if(abs(rcond) < 1d-15) n_diis = 0
+    end if
 
     write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') &
       '|',nSCF,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|'
@@ -201,6 +200,116 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF
 ! End of SCF loop
 !------------------------------------------------------------------------
 
+! Did it actually converge?
+
+  if(nSCF == maxSCF) then
+
+    write(*,*)
+    write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
+    write(*,*)'        Convergence failed for t ampitudes          '
+    write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
+
+    stop
+ 
+  end if
+
+! Deallocate memory
+
+  deallocate(err_diis,t_diis)
+
+! Memory allocation 
+
+  allocate(z(O,V))
+
+! Memory allocation for DIIS
+
+  allocate(err_diis(O*V,max_diis),z_diis(O*V,max_diis))
+
+!------------------------------------------------------------------------
+! Compute z ampltiudes
+!------------------------------------------------------------------------
+
+  Conv = 1d0
+  nSCF = 0
+
+  n_diis          = 0
+  z_diis(:,:)     = 0d0
+  err_diis(:,:) = 0d0
+
+!------------------------------------------------------------------------
+! Main SCF loop
+!------------------------------------------------------------------------
+  write(*,*)
+  write(*,*)'----------------------------------------------------'
+  write(*,*)'| pCCD calculation: z amplitudes                   |'
+  write(*,*)'----------------------------------------------------'
+  write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') &
+            '|','#','|','E(pCCD)','|','Ec(pCCD)','|','Conv','|'
+  write(*,*)'----------------------------------------------------'
+
+  do while(Conv > thresh .and. nSCF < maxSCF)
+
+  ! Increment 
+
+    nSCF = nSCF + 1
+
+  ! Form intermediate array
+  
+   yO(:,:) = matmul(OOVV,transpose(t))
+   yV(:,:) = matmul(transpose(OOVV),t)
+
+   ! Compute residual
+
+    r(:,:) = OOVV(:,:) + 2d0*delta_OV(:,:)*z(:,:) &
+           - 2d0*(2d0*OVOV(:,:) - OVVO(:,:) - 2d0*OOVV(:,:)*t(:,:))*z(:,:)
+
+    do i=1,O
+      do a=1,V
+
+        do j=1,O
+          r(i,a) = r(i,a) - 2d0*OOVV(j,a)*t(j,a)*z(i,a) &
+                          - 2d0*OOVV(i,a)*z(j,a)*t(j,a) & 
+                          + OOOO(i,j)*z(j,a)            &
+                          + yO(i,j)*z(j,a)
+        end do
+
+        do b=1,V
+          r(i,a) = r(i,a) - 2d0*OOVV(i,b)*t(i,b)*z(i,a) & 
+                          - 2d0*OOVV(i,a)*z(i,b)*t(i,b) & 
+                          + VVVV(b,a)*z(i,b)            & 
+                          + yV(a,b)*z(i,b)
+        end do
+
+      end do
+    end do
+
+   ! Check convergence 
+
+    Conv = maxval(abs(r(:,:)))
+  
+   ! Update amplitudes
+
+   z(:,:) = z(:,:) - 0.5d0*r(:,:)/delta_OV(:,:)
+
+   ! DIIS extrapolation
+
+   if(max_diis > 1) then
+
+     n_diis = min(n_diis+1,max_diis)
+     call DIIS_extrapolation(rcond,O*V,O*V,n_diis,err_diis,z_diis,-0.5d0*r/delta_OV,z)
+
+    end if
+
+    write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') &
+      '|',nSCF,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|'
+
+  end do
+  write(*,*)'----------------------------------------------------'
+  write(*,*)
+!------------------------------------------------------------------------
+! End of SCF loop
+!------------------------------------------------------------------------
+
 ! Did it actually converge?
 
   if(nSCF == maxSCF) then
@@ -214,6 +323,152 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF
  
   end if
 
+! Deallocate memory
+
+  deallocate(err_diis,z_diis,r)
+
+!--------------------------!
+! Compute density matrices !
+!--------------------------!
+
+  allocate(xOO(O,O),xVV(V,V),xOV(O,V))
+
+  xOO(:,:) = matmul(t,transpose(z))
+  xVV(:,:) = matmul(transpose(z),t)
+  xOV(:,:) = matmul(t,matmul(transpose(z),t))
+
+! Form 1RDM
+
+  allocate(rdm1(O+V,O+V))
+
+  rdm1(:,:) = 0d0
+
+  do i=1,O
+    rdm1(i,i) = 2d0*(1d0 - xOO(i,i))
+  end do
+
+  do a=1,V
+    rdm1(O+a,O+a) = 2d0*xVV(a,a)
+  end do
+
+! Check 1RDM
+
+  tr_1rdm = trace_matrix(O+V,rdm1)
+  write(*,*) ' --> Trace of the 1RDM = ',tr_1rdm
+
+  if( abs(dble(2*O) - tr_1rdm) > thresh ) & 
+  write(*,*) ' !!! Your 1RDM seems broken !!! '
+  write(*,*)
+
+! Form 2RM
+
+  allocate(rdm2(O+V,O+V,O+V,O+V))
+
+  rdm2(:,:,:,:) = 0d0
+
+  ! iijj
+
+  do i=1,O
+    do j=1,O
+      rdm2(i,i,j,j) = 2d0*xOO(i,j)
+    end do
+  end do
+
+  ! iiaa
+
+  do i=1,O
+    do a=1,V
+      rdm2(i,i,O+a,O+a) = 2d0*(t(i,a) + xOV(i,a) - 2d0*t(i,a)*(xVV(a,a) + xOO(i,i) - t(i,a)*z(i,a)))
+    end do
+  end do
+
+  ! aaii
+
+  do i=1,O
+    do a=1,V
+      rdm2(O+a,O+a,i,i) = 2d0*z(a,i)
+    end do
+  end do
+
+  ! aabb
+
+  do a=1,V
+    do b=1,V
+      rdm2(O+a,O+a,O+b,O+b) = 2d0*xVV(a,b)
+    end do
+  end do
+
+  ! ijij
+
+  do i=1,O
+    do j=1,O
+      rdm2(i,j,i,j) = 4d0*(1d0 - xOO(i,i) - xOO(j,j))
+    end do
+  end do
+
+  ! ijji
+
+  do i=1,O
+    do j=1,O
+      rdm2(i,j,j,i) = - 2d0*(1d0 - xOO(i,i) - xOO(j,j))
+    end do
+  end do
+
+  ! iiii
+
+  do i=1,O
+    rdm2(i,i,i,i) = 2d0*(1d0 - xOO(i,i))
+  end do
+
+  ! iaia
+
+  do i=1,O
+    do a=1,V
+      rdm2(i,O+a,i,O+a) = 4d0*(xVV(a,a) - t(i,a)*z(i,a))
+    end do
+  end do
+
+  ! iaai
+
+  do i=1,O
+    do a=1,V
+      rdm2(i,O+a,O+a,i) = - 2d0*(xVV(a,a) - t(i,a)*z(i,a))
+    end do
+  end do
+
+  ! aiai
+
+  do i=1,O
+    do a=1,V
+      rdm2(O+a,i,O+a,i) = 4d0*(xVV(a,a) - t(i,a)*z(i,a))
+    end do
+  end do
+
+  ! aiia
+
+  do i=1,O
+    do a=1,V
+      rdm2(O+a,i,i,O+a) = - 2d0*(xVV(a,a) - t(i,a)*z(i,a))
+    end do
+  end do
+
+  ! abab
+
+  do a=1,V
+    rdm2(O+a,O+a,O+a,O+a) = 2d0*xVV(a,a)
+  end do
+
+! Check 2RDM
+
+  tr_2rdm = trace_matrix((O+V)**2,rdm2)
+  write(*,*) ' --> Trace of the 2RDM = ',tr_2rdm
+
+  if( abs(dble(2*O*(2*O-1)) - tr_2rdm) > thresh ) & 
+  write(*,*) ' !!! Your 2RDM seems broken !!! '
+  write(*,*)
+
+! Testing zone
+
   if(dotest) then
 
     call dump_test_value('R','pCCD correlation energy',EcCC)