From 08bf6632dfbdd65f5aa18b9c1e8d67cf4d7d8049 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Fri, 30 Aug 2024 16:26:59 +0200 Subject: [PATCH] saving work in pCCD --- src/CC/RCC.f90 | 22 +- src/CC/pCCD.f90 | 1141 ++++++++++++++++++++++-------------------- src/QuAcK/RQuAcK.f90 | 2 +- 3 files changed, 603 insertions(+), 562 deletions(-) diff --git a/src/CC/RCC.f90 b/src/CC/RCC.f90 index 21f6866..aa8923b 100644 --- a/src/CC/RCC.f90 +++ b/src/CC/RCC.f90 @@ -1,5 +1,5 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, & - maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,eHF,cHF) + maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI_AO,ERI_MO,ENuc,ERHF,eHF,cHF) ! Coupled-cluster module @@ -34,7 +34,8 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: cHF(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas) ! Local variables @@ -47,7 +48,7 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d if(doCCD) then call wall_time(start_CC) - call CCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) + call CCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC @@ -64,7 +65,7 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d call wall_time(start_CC) call DCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR, & - ERI,ENuc,ERHF,eHF) + ERI_MO,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC @@ -82,7 +83,7 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d if(doCCSD) then call wall_time(start_CC) - call CCSD(dotest,maxSCF,thresh,max_diis,doCCSDT,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) + call CCSD(dotest,maxSCF,thresh,max_diis,doCCSDT,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC @@ -98,7 +99,7 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d if(dodrCCD) then call wall_time(start_CC) - call drCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) + call drCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC @@ -114,7 +115,7 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d if(dorCCD) then call wall_time(start_CC) - call rCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) + call rCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC @@ -130,7 +131,7 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d if(docrCCD) then call wall_time(start_CC) - call crCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) + call crCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC @@ -146,7 +147,7 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d if(dolCCD) then call wall_time(start_CC) - call lCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) + call lCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC @@ -162,7 +163,8 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d if(dopCCD) then call wall_time(start_CC) - call pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,eHF,cHF) + call pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI_AO,ENuc,ERHF,eHF,cHF) + call wall_time(end_CC) t_CC = end_CC - start_CC diff --git a/src/CC/pCCD.f90 b/src/CC/pCCD.f90 index 570a679..bcfd50b 100644 --- a/src/CC/pCCD.f90 +++ b/src/CC/pCCD.f90 @@ -1,4 +1,4 @@ -subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,eHF,cHF) +subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI_AO,ENuc,ERHF,eHF,cHF) ! pair CCD module @@ -8,25 +8,32 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF, logical,intent(in) :: dotest - integer,intent(in) :: maxSCF + integer,intent(in) :: maxIt integer,intent(in) :: max_diis double precision,intent(in) :: thresh - integer,intent(in) :: nBas,nC,nO,nV,nR - double precision,intent(in) :: ENuc,ERHF + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + double precision,intent(in) :: ENuc + double precision,intent(in) :: ERHF double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: cHF(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) ! Local variables - integer :: p,q,r,s,t,u + integer :: p,q,r,s,t,u,w integer :: pq,rs integer :: i,j,a,b - integer :: nSCF - double precision :: Conv + integer :: nItAmp + integer :: nItOrb + double precision :: CvgAmp + double precision :: CvgOrb double precision :: ECC double precision :: EcCC @@ -56,14 +63,15 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF, double precision :: tr_2rdm double precision :: E1,E2 - double precision,allocatable :: h(:,:) double precision,allocatable :: c(:,:) + double precision,allocatable :: h(:,:) + double precision,allocatable :: ERI_MO(:,:,:,:) double precision,allocatable :: grad(:) double precision,allocatable :: tmp(:,:,:,:) double precision,allocatable :: hess(:,:) double precision,allocatable :: hessInv(:,:) - double precision,allocatable :: kappa(:,:) - double precision,allocatable :: ekappa(:,:) + double precision,allocatable :: Kap(:,:) + double precision,allocatable :: ExpKap(:,:) integer :: O,V,N integer :: n_diis @@ -77,9 +85,9 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF, ! Hello world write(*,*) - write(*,*)'**************************************' - write(*,*)'| pair CCD calculation |' - write(*,*)'**************************************' + write(*,*)'*******************************' + write(*,*)'* Restricted pCCD Calculation *' + write(*,*)'*******************************' write(*,*) ! Useful quantities @@ -88,612 +96,643 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF, V = nV - nR N = O + V -! Form energy denominator + !------------------------------------! + ! Star Loop for orbital optimization ! + !------------------------------------! + allocate(ERI_MO(N,N,N,N)) + allocate(c(N,N),h(N,N)) allocate(eO(O),eV(V),delta_OV(O,V)) - - eO(:) = eHF(nC+1:nO) - eV(:) = eHF(nO+1:nBas-nR) - - call form_delta_OV(nC,nO,nV,nR,eO,eV,delta_OV) - -! Create integral batches - allocate(OOOO(O,O),OOVV(O,V),OVOV(O,V),OVVO(O,V),VVVV(V,V)) - do i=1,O - do j=1,O - OOOO(i,j) = ERI(nC+i,nC+i,nC+j,nC+j) - end do - end do + c(:,:) = cHF(nC+1:nBas-nR,nC+1:nBas-nR) - do i=1,O - 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 + CvgOrb = 1d0 + nItOrb = 0 - do a=1,V - do b=1,V - VVVV(a,b) = ERI(nO+a,nO+a,nO+b,nO+b) - end do - end do - -! Initialization - - allocate(t2(O,V),r2(O,V),yO(O,O),yV(V,V)) - -! Memory allocation for DIIS - - allocate(err_diis(O*V,max_diis),t2_diis(O*V,max_diis)) - -!------------------------------------------------------------------------ -! Start orbital optimization -!------------------------------------------------------------------------ - - allocate(c(N,N)) - c(:,:) = cHF(:,:) - -!------------------------------------------------------------------------ -! Compute t ampltiudes -!------------------------------------------------------------------------ - - Conv = 1d0 - nSCF = 0 - ECC = ERHF - EcCC = 0d0 - - n_diis = 0 - t2(:,:) = 0d0 - t2_diis(:,:) = 0d0 - err_diis(:,:) = 0d0 - -!------------------------------------------------------------------------ -! Main SCF loop -!------------------------------------------------------------------------ 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(*,*)'| Orbital Optimization for pCCD |' write(*,*)'----------------------------------------------------' - do while(Conv > thresh .and. nSCF < maxSCF) + do while(CvgOrb > thresh .and. nItOrb < 1) - ! Increment + nItOrb = nItOrb + 1 - nSCF = nSCF + 1 + ! Transform integrals - ! Form intermediate array - - yO(:,:) = matmul(t2,transpose(OOVV)) - - ! Compute residual + h = matmul(transpose(c),matmul(Hc(nC+1:nBas-nR,nC+1:nBas-nR),c)) + call AOtoMO_ERI_RHF(N,c,ERI_AO(nC+1:nBas-nR,nC+1:nBas-nR,nC+1:nBas-nR,nC+1:nBas-nR),ERI_MO) + ! Form energy denominator - r2(:,:) = OOVV(:,:) + 2d0*delta_OV(:,:)*t2(:,:) & - - 2d0*(2d0*OVOV(:,:) - OVVO(:,:) - OOVV(:,:)*t2(:,:))*t2(:,:) + eO(:) = eHF(nC+1:nO) + eV(:) = eHF(nO+1:nBas-nR) 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 - + delta_OV(i,a) = eV(a) - eO(i) end do end do - ! Check convergence - - Conv = maxval(abs(r2(:,:))) - - ! Update amplitudes - - t2(:,:) = t2(:,:) - 0.5d0*r2(:,:)/delta_OV(:,:) - - ! Compute correlation energy - - EcCC = trace_matrix(V,matmul(transpose(OOVV),t2)) - - ! Dump results - - ECC = ERHF + EcCC - - ! DIIS extrapolation - - if(max_diis > 1) then - - n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nO*nV,nO*nV,n_diis,err_diis,t2_diis,-0.5d0*r2/delta_OV,t2) - - 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(*,*)'----------------------------------------------------' -!------------------------------------------------------------------------ -! 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,t2_diis) - -! Memory allocation - - allocate(z2(O,V)) - -! Memory allocation for DIIS - - allocate(err_diis(O*V,max_diis),z2_diis(O*V,max_diis)) - -!------------------------------------------------------------------------ -! Compute z ampltiudes -!------------------------------------------------------------------------ - - Conv = 1d0 - nSCF = 0 - - n_diis = 0 - z2_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(t2)) - yV(:,:) = matmul(transpose(OOVV),t2) - - ! Compute residual - - r2(:,:) = OOVV(:,:) + 2d0*delta_OV(:,:)*z2(:,:) & - - 2d0*(2d0*OVOV(:,:) - OVVO(:,:) - 2d0*OOVV(:,:)*t2(:,:))*z2(:,:) + ! Create integral batches + do i=1,O + do j=1,O + OOOO(i,j) = ERI_MO(i,i,j,j) + end do + end do + 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 - + OOVV(i,a) = ERI_MO(i,i,O+a,O+a) + OVOV(i,a) = ERI_MO(i,O+a,i,O+a) + OVVO(i,a) = ERI_MO(i,O+a,O+a,i) + end do + end do + + do a=1,V + do b=1,V + VVVV(a,b) = ERI_MO(O+a,O+a,O+b,O+b) end do end do - ! Check convergence + !----------------------------! + ! Star Loop for t amplitudes ! + !----------------------------! + + allocate(t2(O,V),r2(O,V),yO(O,O)) + allocate(err_diis(O*V,max_diis),t2_diis(O*V,max_diis)) - Conv = maxval(abs(r2(:,:))) + CvgAmp = 1d0 + nItAmp = 0 + ECC = ERHF + EcCC = 0d0 + + n_diis = 0 + t2(:,:) = 0d0 + t2_diis(:,:) = 0d0 + 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(*,*)'----------------------------------------------------' + + do while(CvgAmp > thresh .and. nItAmp < maxIt) + + ! 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 + + ! Check convergence + + CvgAmp = maxval(abs(r2(:,:))) + + ! Update amplitudes + + t2(:,:) = t2(:,:) - 0.5d0*r2(:,:)/delta_OV(:,:) + + ! Compute correlation energy + + EcCC = 0d0 + do i=1,O + do a=1,V + EcCC = EcCC + OOVV(i,a)*t2(i,a) + end do + end do + + ! Dump results + + ECC = ERHF + EcCC + + ! DIIS extrapolation + + if(max_diis > 1) then + + n_diis = min(n_diis+1,max_diis) + call DIIS_extrapolation(rcond,nO*nV,nO*nV,n_diis,err_diis,t2_diis,-0.5d0*r2/delta_OV,t2) + + 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,'|' + + end do + write(*,*)'----------------------------------------------------' + + !---------------------------! + ! End Loop for t amplitudes ! + !---------------------------! + + deallocate(r2,yO) + deallocate(err_diis,t2_diis) + + ! Did it actually converge? + + if(nItAmp == maxIt) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)'! Convergence failed for t ampitudes !' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + + stop - ! Update amplitudes - - z2(:,:) = z2(:,:) - 0.5d0*r2(:,:)/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,z2_diis,-0.5d0*r2/delta_OV,z2) - end if + + !-----------------------------! + ! Start Loop for z amplitudes ! + !-----------------------------! + + allocate(z2(O,V),r2(O,V),yO(O,O),yV(V,V)) + allocate(err_diis(O*V,max_diis),z2_diis(O*V,max_diis)) - 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 + CvgAmp = 1d0 + nItAmp = 0 + + n_diis = 0 + z2_diis(:,:) = 0d0 + err_diis(:,:) = 0d0 write(*,*) - write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' - write(*,*)' Convergence failed ' - write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' - - stop + 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(CvgAmp > thresh .and. nItAmp < maxIt) + + ! 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 + + ! Check convergence + + CvgAmp = maxval(abs(r2(:,:))) + + ! Update amplitudes + + z2(:,:) = z2(:,:) - 0.5d0*r2(:,:)/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,z2_diis,-0.5d0*r2/delta_OV,z2) + + 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,'|' - end if - -! Deallocate memory - - deallocate(err_diis,z2_diis,r2) - -!--------------------------! -! Compute density matrices ! -!--------------------------! - - allocate(xOO(O,O),xVV(V,V),xOV(O,V)) - - xOO(:,:) = matmul(t2,transpose(z2)) - xVV(:,:) = matmul(transpose(z2),t2) - xOV(:,:) = matmul(t2,matmul(transpose(z2),t2)) - -! Form 1RDM - - allocate(rdm1(N,N)) - - 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(N,rdm1) - write(*,'(A25,F16.10)') ' --> Trace of the 1RDM = ',tr_1rdm - - if( abs(dble(2*O) - tr_1rdm) > thresh ) & - write(*,*) ' !!! Your 1RDM seems broken !!! ' - write(*,*) - -! write(*,*) '1RDM is diagonal at the pCCD level:' -! call matout(N,N,rdm1) - -! Form 2RM - - allocate(rdm2(N,N,N,N)) - - rdm2(:,:,:,:) = 0d0 - - ! iijj - - do i=1,O - do j=1,O - rdm2(i,i,j,j) = 2d0*xOO(i,j) end do - end do + write(*,*)'----------------------------------------------------' + write(*,*) - ! iiaa + !---------------------------! + ! End Loop for z ampltiudes ! + !---------------------------! - do i=1,O + deallocate(r2,yO,yV) + deallocate(err_diis,z2_diis) + + ! Did it actually converge? + + if(nItAmp == maxIt) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)'! Convergence failed for z ampltiudes !' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + + stop + + end if + + !--------------------------! + ! Compute density matrices ! + !--------------------------! + + allocate(rdm1(N,N),rdm2(N,N,N,N)) + allocate(xOO(O,O),xVV(V,V),xOV(O,V)) + + xOO(:,:) = matmul(t2,transpose(z2)) + xVV(:,:) = matmul(transpose(z2),t2) + xOV(:,:) = matmul(t2,matmul(transpose(z2),t2)) + + ! Form 1RDM + + rdm1(:,:) = 0d0 + + do i=1,O + rdm1(i,i) = 2d0*(1d0 - xOO(i,i)) + end do + do a=1,V - rdm2(i,i,O+a,O+a) = 2d0*(t2(i,a) + xOV(i,a) - 2d0*t2(i,a)*(xVV(a,a) + xOO(i,i) - t2(i,a)*z2(i,a))) + rdm1(O+a,O+a) = 2d0*xVV(a,a) end do - end do - - ! aaii - - do i=1,O + + ! Check 1RDM + + tr_1rdm = trace_matrix(N,rdm1) + write(*,'(A25,F16.10)') ' --> Trace of the 1RDM = ',tr_1rdm + + if( abs(dble(2*O) - tr_1rdm) > thresh ) & + write(*,*) ' !!! Your 1RDM seems broken !!! ' + write(*,*) + + ! write(*,*) '1RDM is diagonal at the pCCD level:' + ! call matout(N,N,rdm1) + + ! Form 2RM + + 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*(t2(i,a) + xOV(i,a) - 2d0*t2(i,a)*(xVV(a,a) + xOO(i,i) - t2(i,a)*z2(i,a))) + end do + end do + + ! aaii + + do i=1,O + do a=1,V + rdm2(O+a,O+a,i,i) = 2d0*z2(i,a) + end do + end do + + ! aabb + do a=1,V - rdm2(O+a,O+a,i,i) = 2d0*z2(i,a) + do b=1,V + rdm2(O+a,O+a,O+b,O+b) = 2d0*xVV(a,b) + end do 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) + + ! 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 - end do - - ! ijij - - do i=1,O - do j=1,O - rdm2(i,j,i,j) = 4d0*(1d0 - xOO(i,i) - xOO(j,j)) + + ! 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 - end do - - ! ijji - - do i=1,O - do j=1,O - rdm2(i,j,j,i) = - 2d0*(1d0 - xOO(i,i) - xOO(j,j)) + + ! iiii + + do i=1,O + rdm2(i,i,i,i) = 2d0*(1d0 - xOO(i,i)) 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 + + ! iaia + + do i=1,O + do a=1,V + rdm2(i,O+a,i,O+a) = 4d0*(xVV(a,a) - t2(i,a)*z2(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) - t2(i,a)*z2(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) - t2(i,a)*z2(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) - t2(i,a)*z2(i,a)) + end do + end do + + ! abab + do a=1,V - rdm2(i,O+a,i,O+a) = 4d0*(xVV(a,a) - t2(i,a)*z2(i,a)) + rdm2(O+a,O+a,O+a,O+a) = 2d0*xVV(a,a) end do - end do + + ! Check 2RDM + + tr_2rdm = trace_matrix(N**2,rdm2) + write(*,'(A25,F16.10)') ' --> Trace of the 2RDM = ',tr_2rdm + + if( abs(dble(2*O*(2*O-1)) - tr_2rdm) > thresh ) & + write(*,*) ' !!! Your 2RDM seems broken !!! ' + write(*,*) + + ! write(*,*) '2RDM is not diagonal at the pCCD level:' + ! call matout(N**2,N**2,rdm2) - ! iaai - - do i=1,O - do a=1,V - rdm2(i,O+a,O+a,i) = - 2d0*(xVV(a,a) - t2(i,a)*z2(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) - t2(i,a)*z2(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) - t2(i,a)*z2(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(N**2,rdm2) - write(*,'(A25,F16.10)') ' --> Trace of the 2RDM = ',tr_2rdm - - if( abs(dble(2*O*(2*O-1)) - tr_2rdm) > thresh ) & - write(*,*) ' !!! Your 2RDM seems broken !!! ' - write(*,*) - -! write(*,*) '2RDM is not diagonal at the pCCD level:' -! call matout(N**2,N**2,rdm2) - -! Compute electronic energy - - allocate(h(N,N)) - h = matmul(transpose(cHF),matmul(Hc,cHF)) - - E1 = 0d0 - E2 = 0d0 - - do p=1,N - do q=1,N - E1 = E1 + rdm1(p,q)*h(p,q) - do r=1,N - do s=1,N - E2 = E2 + rdm2(p,q,r,s)*ERI(p,q,r,s) + deallocate(xOO,xVV,xOV) + deallocate(t2,z2) + + ! Compute electronic energy + + E1 = 0d0 + E2 = 0d0 + + do p=1,N + do q=1,N + E1 = E1 + rdm1(p,q)*h(p,q) + do r=1,N + do s=1,N + E2 = E2 + rdm2(p,q,r,s)*ERI_MO(p,q,r,s) + end do end do end do end do - end do - - E2 = 0.5d0*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(*,*) - -! Compute gradient - - allocate(grad(N**2)) - - grad(:) = 0d0 - - pq = 0 - do p=1,N - do q=1,N - - pq = pq + 1 - - do r=1,N - grad(pq) = grad(pq) + h(r,p)*rdm1(r,q) - h(q,r)*rdm1(p,r) - end do - - do r=1,N - do s=1,N - do t=1,N - grad(pq) = grad(pq) + (ERI(r,s,p,t)*rdm2(r,s,q,t) - ERI(q,t,r,s)*rdm2(p,t,r,s)) + + E2 = 0.5d0*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(*,*) + + !--------------------------! + ! Compute orbital gradient ! + !--------------------------! + + allocate(grad(N**2)) + + grad(:) = 0d0 + + pq = 0 + do p=1,N + do q=1,N + + pq = pq + 1 + + do r=1,N + grad(pq) = grad(pq) + h(r,p)*rdm1(r,q) - h(q,r)*rdm1(p,r) + end do + + do r=1,N + do s=1,N + do t=1,N + grad(pq) = grad(pq) + (ERI_MO(r,s,p,t)*rdm2(r,s,q,t) - ERI_MO(q,t,r,s)*rdm2(p,t,r,s)) + end do end do end do - end do - - end do - end do - - write(*,*) 'Orbital gradient at the pCCD level:' - call matout(N,N,grad) - write(*,*) - -! Convergence - - Conv = maxval(abs(grad)) - write(*,*) ' Convergence of orbtial gradient = ',Conv - write(*,*) -! Compute Hessian - - allocate(hess(N**2,N**2),tmp(N,N,N,N)) - - tmp(:,:,:,:) = 0d0 - - do p=1,N - do q=1,N - - do r=1,N - do s=1,N - - tmp(p,q,r,s) = - h(s,p)*rdm1(r,q) - h(q,r)*rdm1(p,s) - - do u=1,N - - tmp(p,q,r,s) = tmp(p,q,r,s) + 0.5d0*( & - Kronecker_delta(q,r)*(h(u,p)*rdm1(u,s) + h(s,u)*rdm1(p,u)) & - + Kronecker_delta(p,s)*(h(u,r)*rdm1(u,q) + h(q,u)*rdm1(r,u)) ) - - end do - - do u=1,N - do v=1,N - - tmp(p,q,r,s) = tmp(p,q,r,s) + ERI(u,v,p,r)*rdm2(u,v,q,s) + ERI(q,s,u,v)*rdm2(p,r,u,v) - - end do - end do - - do t=1,N + end do + end do + + write(*,*) 'Orbital gradient at the pCCD level:' + call matout(N,N,grad) + write(*,*) + + ! Check convergence of orbital optimization + + CvgOrb = maxval(abs(grad)) + write(*,*) ' Iteration',nItOrb,'for pCCD orbital optimization' + write(*,*) ' Convergence of orbital gradient = ',CvgOrb + write(*,*) + + !-------------------------! + ! Compute orbital Hessian ! + !-------------------------! + + allocate(hess(N**2,N**2),tmp(N,N,N,N)) + + tmp(:,:,:,:) = 0d0 + + do p=1,N + do q=1,N + + do r=1,N + do s=1,N + + tmp(p,q,r,s) = - h(s,p)*rdm1(r,q) - h(q,r)*rdm1(p,s) + do u=1,N - - tmp(p,q,r,s) = tmp(p,q,r,s) - ( & - ERI(s,t,p,u)*rdm2(r,t,q,u) + ERI(t,s,p,u)*rdm2(t,r,q,u) & - + ERI(q,u,r,t)*rdm2(p,u,s,t) + ERI(q,u,t,r)*rdm2(p,u,t,s) ) - + + tmp(p,q,r,s) = tmp(p,q,r,s) + 0.5d0*( & + Kronecker_delta(q,r)*(h(u,p)*rdm1(u,s) + h(s,u)*rdm1(p,u)) & + + Kronecker_delta(p,s)*(h(u,r)*rdm1(u,q) + h(q,u)*rdm1(r,u)) ) + end do - end do - - do t=1,N + do u=1,N - do v=1,N - - tmp(p,q,r,s) = tmp(p,q,r,s) + 0.5d0*( & - Kronecker_delta(q,r)*(ERI(u,v,p,t)*rdm2(u,v,s,t) + ERI(s,t,u,v)*rdm2(p,t,u,v)) & - + Kronecker_delta(p,s)*(ERI(q,t,u,v)*rdm2(r,t,u,v) + ERI(u,v,r,t)*rdm2(u,v,q,t)) ) - + do w=1,N + + tmp(p,q,r,s) = tmp(p,q,r,s) + ERI_MO(u,w,p,r)*rdm2(u,w,q,s) + ERI_MO(q,s,u,w)*rdm2(p,r,u,w) + end do end do + + do t=1,N + do u=1,N + + tmp(p,q,r,s) = tmp(p,q,r,s) - ( & + ERI_MO(s,t,p,u)*rdm2(r,t,q,u) + ERI_MO(t,s,p,u)*rdm2(t,r,q,u) & + + ERI_MO(q,u,r,t)*rdm2(p,u,s,t) + ERI_MO(q,u,t,r)*rdm2(p,u,t,s) ) + + end do + end do + + do t=1,N + do u=1,N + do w=1,N + + tmp(p,q,r,s) = tmp(p,q,r,s) + 0.5d0*( & + Kronecker_delta(q,r)*(ERI_MO(u,w,p,t)*rdm2(u,w,s,t) + ERI_MO(s,t,u,w)*rdm2(p,t,u,w)) & + + Kronecker_delta(p,s)*(ERI_MO(q,t,u,w)*rdm2(r,t,u,w) + ERI_MO(u,w,r,t)*rdm2(u,w,q,t)) ) + + end do + end do + end do + end do - end do + end do - end do + + ! Flatten Hessian matrix and add permutations + + pq = 0 + do p=1,N + do q=1,N + + pq = pq + 1 + + rs = 0 + do r=1,N + do s=1,N + + rs = rs + 1 + + hess(pq,rs) = tmp(p,r,q,s) - tmp(r,p,q,s) - tmp(p,r,s,q) + tmp(r,p,s,q) +!! hess(pq,rs) = tmp(p,q,r,s) - tmp(q,p,r,s) - tmp(p,q,s,r) + tmp(q,p,s,r) + + end do + end do + + end do + end do + + deallocate(rdm1,rdm2,tmp) + + allocate(hessInv(N**2,N**2)) + + call inverse_matrix(N**2,hess,hessInv) + + deallocate(hess) + + allocate(Kap(N,N)) + + Kap(:,:) = 0d0 + + pq = 0 + do p=1,N + do q=1,N + + pq = pq + 1 + + rs = 0 + do r=1,N + do s=1,N + + rs = rs + 1 + + Kap(p,q) = Kap(p,q) - hessInv(pq,rs)*grad(rs) + + end do + end do + + end do + end do + + deallocate(hessInv,grad) + + write(*,*) 'kappa' + call matout(N,N,Kap) + write(*,*) + + allocate(ExpKap(N,N)) + call matrix_exponential(N,Kap,ExpKap) + deallocate(Kap) + + write(*,*) 'e^kappa' + call matout(N,N,ExpKap) + write(*,*) + + write(*,*) 'Old orbitals' + call matout(N,N,c) + write(*,*) + + c = matmul(c,ExpKap) + deallocate(ExpKap) + + write(*,*) 'Rotated orbitals' + call matout(N,N,c) + write(*,*) + end do -! Flatten Hessian matrix and add permutations + !-----------------------------------! + ! End Loop for orbital optimization ! + !-----------------------------------! - pq = 0 - do p=1,N - do q=1,N + ! Did it actually converge? - pq = pq + 1 - - rs = 0 - do r=1,N - do s=1,N + if(nItOrb == maxIt) then - rs = rs + 1 + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)'! Convergence failed for orbital optimization !' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' - hess(pq,rs) = tmp(p,r,q,s) - tmp(r,p,q,s) - tmp(p,r,s,q) + tmp(r,p,s,q) -! hess(pq,rs) = tmp(p,q,r,s) - tmp(q,p,r,s) - tmp(p,q,s,r) + tmp(q,p,s,r) + stop - end do - end do - - end do - end do - - call matout(N**2,N**2,hess) - -! allocate(eig(N**2)) - -! call diagonalize_matrix(N**2,hess,eig) - -! call vecout(N**2,eig) - - allocate(hessInv(N**2,N**2)) - - call inverse_matrix(N**2,hess,hessInv) - - allocate(kappa(N,N),ekappa(N,N)) - - kappa(:,:) = 0d0 - - pq = 0 - do p=1,N - do q=1,N - - pq = pq + 1 - - rs = 0 - do r=1,N - do s=1,N - - rs = rs + 1 - - kappa(p,q) = kappa(p,q) + hessInv(pq,rs)*grad(rs) - - end do - end do - - end do - end do - - write(*,*) 'kappa' - call matout(N,N,kappa) - write(*,*) - - call matrix_exponential(N,kappa,ekappa) - - write(*,*) 'e^kappa' - call matout(N,N,ekappa) - write(*,*) - - write(*,*) 'Old orbitals' - call matout(N,N,c) - write(*,*) - - c = matmul(c,ekappa) - - write(*,*) 'Rotated orbitals' - call matout(N,N,c) - write(*,*) + end if ! Testing zone diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90 index 3897f37..9632c32 100644 --- a/src/QuAcK/RQuAcK.f90 +++ b/src/QuAcK/RQuAcK.f90 @@ -228,7 +228,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d call wall_time(start_CC) call RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, & - maxSCF_CC,thresh_CC,max_diis_CC,nBas,nC,nO,nV,nR,Hc,ERI_MO,ENuc,ERHF,eHF,cHF) + maxSCF_CC,thresh_CC,max_diis_CC,nBas,nC,nO,nV,nR,Hc,ERI_AO,ERI_MO,ENuc,ERHF,eHF,cHF) call wall_time(end_CC) t_CC = end_CC - start_CC