From 348577f72a4667d5cd404a0b5e0a9946c723a2e4 Mon Sep 17 00:00:00 2001 From: pfloos Date: Tue, 3 Sep 2024 10:59:54 +0200 Subject: [PATCH] debug OO-pCCD --- src/CC/pCCD.f90 | 152 +++++++++++-------------------- src/CC/pCCD_orbital_gradient.f90 | 9 +- src/CC/pCCD_orbital_hessian.f90 | 11 ++- src/CC/pCCD_rdm.f90 | 9 +- src/CC/pCCD_t_residual.f90 | 57 ++++++++++++ src/CC/pCCD_z_residual.f90 | 62 +++++++++++++ 6 files changed, 193 insertions(+), 107 deletions(-) create mode 100644 src/CC/pCCD_t_residual.f90 create mode 100644 src/CC/pCCD_z_residual.f90 diff --git a/src/CC/pCCD.f90 b/src/CC/pCCD.f90 index d16e117..d71e89e 100644 --- a/src/CC/pCCD.f90 +++ b/src/CC/pCCD.f90 @@ -19,7 +19,8 @@ subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR, & integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR - double precision,intent(in) :: ENuc,ERHF + double precision,intent(in) :: ENuc + double precision,intent(in) :: ERHF double precision,intent(in) :: eHF(nOrb) double precision,intent(in) :: cHF(nBas,nOrb) double precision,intent(in) :: PHF(nBas,nBas) @@ -40,7 +41,7 @@ subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR, & double precision :: CvgOrb double precision :: ECC double precision :: EcCC - double precision :: dECC + double precision :: EOld double precision,allocatable :: eO(:) double precision,allocatable :: eV(:) @@ -52,9 +53,6 @@ subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR, & double precision,allocatable :: OVVO(:,:) double precision,allocatable :: VVVV(:,:) - double precision,allocatable :: yO(:,:) - double precision,allocatable :: yV(:,:) - double precision,allocatable :: r2(:,:) double precision,allocatable :: t2(:,:) double precision,allocatable :: z2(:,:) @@ -73,6 +71,7 @@ subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR, & double precision,allocatable :: ExpKap(:,:) integer :: O,V,N + integer :: Np integer :: n_diis double precision :: rcond double precision,allocatable :: err_diis(:,:) @@ -93,6 +92,8 @@ subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR, & V = nV - nR N = O + V + Np = N*N + !------------------------------------! ! Star Loop for orbital optimization ! !------------------------------------! @@ -102,17 +103,18 @@ subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR, & allocate(eO(O),eV(V),delta_OV(O,V)) allocate(OOOO(O,O),OOVV(O,V),OVOV(O,V),OVVO(O,V),VVVV(V,V)) - do i=1,N - c(:,i) = cHF(:,nC+i) + do p=1,N + c(:,p) = cHF(:,nC+p) enddo CvgOrb = 1d0 nItOrb = 0 + EOld = ECC write(*,*) - write(*,*)'----------------------------------------------------' - write(*,*)'| Orbital Optimization for pCCD |' - write(*,*)'----------------------------------------------------' + write(*,*)'---------------------------------------' + write(*,*)'| Orbital Optimization for pCCD |' + write(*,*)'---------------------------------------' do while(CvgOrb > thresh .and. nItOrb < maxIt) @@ -122,7 +124,7 @@ subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR, & h = matmul(transpose(c),matmul(Hc,c)) - call AOtoMO_ERI_RHF(nBas,N,c(1,1),ERI_AO(1,1,1,1),ERI_MO(1,1,1,1)) + call AOtoMO_ERI_RHF(nBas,N,c,ERI_AO,ERI_MO) ! Form energy denominator @@ -175,14 +177,13 @@ subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR, & ! Star Loop for t amplitudes ! !----------------------------! - allocate(t2(O,V),r2(O,V),yO(O,O)) + allocate(t2(O,V),r2(O,V)) allocate(err_diis(O*V,max_diis),t2_diis(O*V,max_diis)) CvgAmp = 1d0 nItAmp = 0 ECC = ERHF EcCC = 0d0 - dECC = ECC n_diis = 0 t2(:,:) = 0d0 @@ -190,41 +191,22 @@ subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR, & err_diis(:,:) = 0d0 write(*,*) - write(*,*)'----------------------------------------------------' - 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','|' - write(*,*)'----------------------------------------------------' + write(*,*)'---------------------------------------' + write(*,*)'| pCCD calculation: t amplitudes |' + write(*,*)'---------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X)') & + '|','#','|','Ec(pCCD)','|','Conv','|' + write(*,*)'---------------------------------------' do while(CvgAmp > thresh .and. nItAmp < maxIt) - ! Increment + ! Increment nItAmp = nItAmp + 1 - ! Form intermediate array - - yO(:,:) = matmul(t2,transpose(OOVV)) - - ! Compute residual - - r2(:,:) = OOVV(:,:) + 2d0*delta_OV(:,:)*t2(:,:) & - - 2d0*(2d0*OVOV(:,:) - OVVO(:,:) - OOVV(:,:)*t2(:,:))*t2(:,:) - - do i=1,O - do a=1,V - - do j=1,O - r2(i,a) = r2(i,a) - 2d0*OOVV(j,a)*t2(j,a)*t2(i,a) + OOOO(j,i)*t2(j,a) + yO(i,j)*t2(j,a) - end do - - do b=1,V - r2(i,a) = r2(i,a) - 2d0*OOVV(i,b)*t2(i,b)*t2(i,a) + VVVV(a,b)*t2(i,b) - end do - - end do - end do + ! Compute residual for t amplitudes + + call pCCD_t_residual(O,V,N,OOVV,OVOV,OVVO,OOOO,VVVV,delta_OV,t2,r2) ! Check convergence @@ -243,10 +225,6 @@ subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR, & end do end do - ! Dump results - - ECC = ERHF + EcCC - ! DIIS extrapolation if(max_diis > 1) then @@ -256,17 +234,17 @@ subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR, & 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)') & - '|',nItAmp,'|',ECC+ENuc,'|',EcCC,'|',CvgAmp,'|' + write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X)') & + '|',nItAmp,'|',EcCC,'|',CvgAmp,'|' end do - write(*,*)'----------------------------------------------------' + write(*,*)'---------------------------------------' !---------------------------! ! End Loop for t amplitudes ! !---------------------------! - deallocate(r2,yO) + deallocate(r2) deallocate(err_diis,t2_diis) ! Did it actually converge? @@ -286,7 +264,7 @@ subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR, & ! Start Loop for z amplitudes ! !-----------------------------! - allocate(z2(O,V),r2(O,V),yO(O,O),yV(V,V)) + allocate(z2(O,V),r2(O,V)) allocate(err_diis(O*V,max_diis),z2_diis(O*V,max_diis)) CvgAmp = 1d0 @@ -297,44 +275,22 @@ subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR, & err_diis(:,:) = 0d0 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(*,*)'----------------------------------------------------' + write(*,*)'---------------------------------------' + write(*,*)'| pCCD calculation: z amplitudes |' + write(*,*)'---------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X)') & + '|','#','|','Ec(pCCD)','|','Conv','|' + write(*,*)'---------------------------------------' do while(CvgAmp > thresh .and. nItAmp < maxIt) - ! Increment + ! Increment nItAmp = nItAmp + 1 - ! Form intermediate array - - yO(:,:) = matmul(OOVV,transpose(t2)) - yV(:,:) = matmul(transpose(OOVV),t2) - - ! Compute residual - - r2(:,:) = OOVV(:,:) + 2d0*delta_OV(:,:)*z2(:,:) & - - 2d0*(2d0*OVOV(:,:) - OVVO(:,:) - 2d0*OOVV(:,:)*t2(:,:))*z2(:,:) - - do i=1,O - do a=1,V - - do j=1,O - r2(i,a) = r2(i,a) - 2d0*OOVV(j,a)*t2(j,a)*z2(i,a) - 2d0*OOVV(i,a)*z2(j,a)*t2(j,a) & - + OOOO(i,j)*z2(j,a) + yO(i,j)*z2(j,a) - end do - - do b=1,V - r2(i,a) = r2(i,a) - 2d0*OOVV(i,b)*t2(i,b)*z2(i,a) - 2d0*OOVV(i,a)*z2(i,b)*t2(i,b) & - + VVVV(b,a)*z2(i,b) + yV(a,b)*z2(i,b) - end do - - end do - end do + ! Compute residual for the z amplitudes + + call pCCD_z_residual(O,V,N,OOVV,OVOV,OVVO,OOOO,VVVV,delta_OV,t2,z2,r2) ! Check convergence @@ -353,18 +309,18 @@ subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR, & 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)') & - '|',nItAmp,'|',ECC+ENuc,'|',EcCC,'|',CvgAmp,'|' + write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X)') & + '|',nItAmp,'|',EcCC,'|',CvgAmp,'|' end do - write(*,*)'----------------------------------------------------' + write(*,*)'---------------------------------------' write(*,*) !---------------------------! ! End Loop for z ampltiudes ! !---------------------------! - deallocate(r2,yO,yV) + deallocate(r2) deallocate(err_diis,z2_diis) ! Did it actually converge? @@ -386,7 +342,7 @@ subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR, & allocate(rdm1(N,N),rdm2(N,N,N,N)) - call pCCD_rdm(O,V,N,ENuc,h,ERI_MO,t2,z2,rdm1,rdm2) + call pCCD_rdm(O,V,N,ENuc,h,ERI_MO,t2,z2,rdm1,rdm2,ECC) deallocate(t2,z2) @@ -394,34 +350,36 @@ subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR, & ! Compute orbital gradient ! !--------------------------! - allocate(grad(N**2)) + allocate(grad(Np)) - call pCCD_orbital_gradient(O,V,N,h,ERI_MO,rdm1,rdm2,grad) + call pCCD_orbital_gradient(O,V,N,Np,h,ERI_MO,rdm1,rdm2,grad) ! Check convergence of orbital optimization CvgOrb = maxval(abs(grad)) + write(*,*) '----------------------------------------------------------' write(*,'(A10,I4,A30)') ' Iteration',nItOrb,'for pCCD orbital optimization' write(*,*) '----------------------------------------------------------' write(*,'(A40,F16.10,A3)') ' Convergence of orbital gradient = ',CvgOrb,' au' - write(*,'(A40,F16.10,A3)') ' Energy difference = ',ECC-dECC,' au' + write(*,'(A40,F16.10,A3)') ' Energy difference = ',ECC-EOld,' au' + write(*,*) '----------------------------------------------------------' write(*,*) - dECC = ECC - + EOld = ECC + !-------------------------! ! Compute orbital Hessian ! !-------------------------! - allocate(hess(N**2,N**2)) + allocate(hess(Np,Np)) - call pCCD_orbital_hessian(O,V,N,h,ERI_MO,rdm1,rdm2,hess) + call pCCD_orbital_hessian(O,V,N,Np,h,ERI_MO,rdm1,rdm2,hess) deallocate(rdm1,rdm2) - allocate(hessInv(N**2,N**2)) + allocate(hessInv(Np,Np)) - call inverse_matrix(N**2,hess,hessInv) + call inverse_matrix(Np,hess,hessInv) deallocate(hess) diff --git a/src/CC/pCCD_orbital_gradient.f90 b/src/CC/pCCD_orbital_gradient.f90 index a991a42..64260a6 100644 --- a/src/CC/pCCD_orbital_gradient.f90 +++ b/src/CC/pCCD_orbital_gradient.f90 @@ -1,4 +1,4 @@ -subroutine pCCD_orbital_gradient(O,V,N,h,ERI_MO,rdm1,rdm2,grad) +subroutine pCCD_orbital_gradient(O,V,N,Np,h,ERI_MO,rdm1,rdm2,grad) ! Compute the orbital gradient at the pCCD level @@ -9,6 +9,7 @@ subroutine pCCD_orbital_gradient(O,V,N,h,ERI_MO,rdm1,rdm2,grad) integer,intent(in) :: O integer,intent(in) :: V integer,intent(in) :: N + integer,intent(in) :: Np double precision,intent(in) :: h(N,N) double precision,intent(in) :: ERI_MO(N,N,N,N) double precision,intent(in) :: rdm1(N,N) @@ -19,11 +20,11 @@ subroutine pCCD_orbital_gradient(O,V,N,h,ERI_MO,rdm1,rdm2,grad) integer :: p,q,r,s,t integer :: pq - logical,parameter :: debug = .false. + logical,parameter :: debug = .true. ! Output variables - double precision,intent(out) :: grad(N**2) + double precision,intent(out) :: grad(Np) ! Compute gradient @@ -50,6 +51,8 @@ subroutine pCCD_orbital_gradient(O,V,N,h,ERI_MO,rdm1,rdm2,grad) end do end do +! Dump gradient + if(debug) then write(*,*) 'Orbital gradient at the pCCD level:' diff --git a/src/CC/pCCD_orbital_hessian.f90 b/src/CC/pCCD_orbital_hessian.f90 index b8987ac..49e77b5 100644 --- a/src/CC/pCCD_orbital_hessian.f90 +++ b/src/CC/pCCD_orbital_hessian.f90 @@ -1,4 +1,4 @@ -subroutine pCCD_orbital_hessian(O,V,N,h,ERI_MO,rdm1,rdm2,hess) +subroutine pCCD_orbital_hessian(O,V,N,Np,h,ERI_MO,rdm1,rdm2,hess) ! Compute the orbital hessian at the pCCD level @@ -9,6 +9,7 @@ subroutine pCCD_orbital_hessian(O,V,N,h,ERI_MO,rdm1,rdm2,hess) integer,intent(in) :: O integer,intent(in) :: V integer,intent(in) :: N + integer,intent(in) :: Np double precision,intent(in) :: h(N,N) double precision,intent(in) :: ERI_MO(N,N,N,N) double precision,intent(in) :: rdm1(N,N) @@ -19,7 +20,7 @@ subroutine pCCD_orbital_hessian(O,V,N,h,ERI_MO,rdm1,rdm2,hess) integer :: p,q,r,s,t,u,w integer :: pq,rs - logical,parameter :: debug = .false. + logical,parameter :: debug = .true. double precision,allocatable :: tmp(:,:,:,:) @@ -27,7 +28,7 @@ subroutine pCCD_orbital_hessian(O,V,N,h,ERI_MO,rdm1,rdm2,hess) ! Output variables - double precision,intent(out) :: hess(N**2,N**2) + double precision,intent(out) :: hess(Np,Np) ! Compute intermediate array @@ -110,10 +111,12 @@ subroutine pCCD_orbital_hessian(O,V,N,h,ERI_MO,rdm1,rdm2,hess) end do end do +! Dump Hessian + if(debug) then write(*,*) 'Orbital Hessian at the pCCD level:' - call matout(N**2,N**2,hess) + call matout(Np,Np,hess) write(*,*) end if diff --git a/src/CC/pCCD_rdm.f90 b/src/CC/pCCD_rdm.f90 index ac8b4a6..4257623 100644 --- a/src/CC/pCCD_rdm.f90 +++ b/src/CC/pCCD_rdm.f90 @@ -1,4 +1,4 @@ -subroutine pCCD_rdm(O,V,N,ENuc,h,ERI_MO,t2,z2,rdm1,rdm2) +subroutine pCCD_rdm(O,V,N,ENuc,h,ERI_MO,t2,z2,rdm1,rdm2,ECC) ! Compute the 1RDM and 2RDM at the pCCD level @@ -38,6 +38,7 @@ subroutine pCCD_rdm(O,V,N,ENuc,h,ERI_MO,t2,z2,rdm1,rdm2) double precision,intent(out) :: rdm1(N,N) double precision,intent(out) :: rdm2(N,N,N,N) + double precision,intent(out) :: ECC ! Allocate memory @@ -209,10 +210,12 @@ subroutine pCCD_rdm(O,V,N,ENuc,h,ERI_MO,t2,z2,rdm1,rdm2) E2 = 0.5d0*E2 + ECC = E1 + E2 + write(*,'(A25,F16.10)') ' One-electron energy = ',E1 write(*,'(A25,F16.10)') ' Two-electron energy = ',E2 - write(*,'(A25,F16.10)') ' Electronic energy = ',E1 + E2 - write(*,'(A25,F16.10)') ' Total energy = ',E1 + E2 + ENuc + write(*,'(A25,F16.10)') ' Electronic energy = ',ECC + write(*,'(A25,F16.10)') ' Total pCCD energy = ',ECC + ENuc write(*,*) end diff --git a/src/CC/pCCD_t_residual.f90 b/src/CC/pCCD_t_residual.f90 new file mode 100644 index 0000000..36fe12c --- /dev/null +++ b/src/CC/pCCD_t_residual.f90 @@ -0,0 +1,57 @@ +subroutine pCCD_t_residual(O,V,N,OOVV,OVOV,OVVO,OOOO,VVVV,delta_OV,t2,r2) + +! Compute the residual for the t amplitudes at the pCCD level + + implicit none + +! Input variables + + integer,intent(in) :: O + integer,intent(in) :: V + integer,intent(in) :: N + double precision,intent(in) :: OOVV(O,V) + double precision,intent(in) :: OVOV(O,V) + double precision,intent(in) :: OVVO(O,V) + double precision,intent(in) :: OOOO(O,O) + double precision,intent(in) :: VVVV(V,V) + double precision,intent(in) :: delta_OV(O,V) + double precision,intent(in) :: t2(O,V) + +! Local variables + + integer :: i,j,a,b + + double precision,allocatable :: yO(:,:) + +! Output variables + + double precision,intent(out) :: r2(O,V) + +! Allocate memory + + allocate(yO(O,O)) + +! Form intermediate array + + yO(:,:) = matmul(t2,transpose(OOVV)) + +! Compute residual + + r2(:,:) = OOVV(:,:) + 2d0*delta_OV(:,:)*t2(:,:) & + - 2d0*(2d0*OVOV(:,:) - OVVO(:,:) - OOVV(:,:)*t2(:,:))*t2(:,:) + + do i=1,O + do a=1,V + + do j=1,O + r2(i,a) = r2(i,a) - 2d0*OOVV(j,a)*t2(j,a)*t2(i,a) + OOOO(j,i)*t2(j,a) + yO(i,j)*t2(j,a) + end do + + do b=1,V + r2(i,a) = r2(i,a) - 2d0*OOVV(i,b)*t2(i,b)*t2(i,a) + VVVV(a,b)*t2(i,b) + end do + + end do + end do + +end diff --git a/src/CC/pCCD_z_residual.f90 b/src/CC/pCCD_z_residual.f90 new file mode 100644 index 0000000..32d0188 --- /dev/null +++ b/src/CC/pCCD_z_residual.f90 @@ -0,0 +1,62 @@ +subroutine pCCD_z_residual(O,V,N,OOVV,OVOV,OVVO,OOOO,VVVV,delta_OV,t2,z2,r2) + +! Compute the residual for the z amplitudes at the pCCD level + + implicit none + +! Input variables + + integer,intent(in) :: O + integer,intent(in) :: V + integer,intent(in) :: N + double precision,intent(in) :: OOVV(O,V) + double precision,intent(in) :: OVOV(O,V) + double precision,intent(in) :: OVVO(O,V) + double precision,intent(in) :: OOOO(O,O) + double precision,intent(in) :: VVVV(V,V) + double precision,intent(in) :: delta_OV(O,V) + double precision,intent(in) :: t2(O,V) + double precision,intent(in) :: z2(O,V) + +! Local variables + + integer :: i,j,a,b + + double precision,allocatable :: yO(:,:) + double precision,allocatable :: yV(:,:) + +! Output variables + + double precision,intent(out) :: r2(O,V) + +! Allocate memory + + allocate(yO(O,O),yV(V,V)) + +! Form intermediate array + + yO(:,:) = matmul(OOVV,transpose(t2)) + yV(:,:) = matmul(transpose(OOVV),t2) + +! Compute residual + + r2(:,:) = OOVV(:,:) + 2d0*delta_OV(:,:)*z2(:,:) & + - 2d0*(2d0*OVOV(:,:) - OVVO(:,:) - 2d0*OOVV(:,:)*t2(:,:))*z2(:,:) + + do i=1,O + do a=1,V + + do j=1,O + r2(i,a) = r2(i,a) - 2d0*OOVV(j,a)*t2(j,a)*z2(i,a) - 2d0*OOVV(i,a)*z2(j,a)*t2(j,a) & + + OOOO(i,j)*z2(j,a) + yO(i,j)*z2(j,a) + end do + + do b=1,V + r2(i,a) = r2(i,a) - 2d0*OOVV(i,b)*t2(i,b)*z2(i,a) - 2d0*OOVV(i,a)*z2(i,b)*t2(i,b) & + + VVVV(b,a)*z2(i,b) + yV(a,b)*z2(i,b) + end do + + end do + end do + +end