From 42565f95f1c6fc6cd508cd420b95c34d8b76ecb3 Mon Sep 17 00:00:00 2001 From: pfloos Date: Fri, 16 Aug 2024 19:05:11 +0200 Subject: [PATCH] implement 1RDM and 2RDM for pCCD --- int/CAP.dat | 16 --- mol/H2.xyz | 4 +- mol/N2.xyz | 4 +- src/CC/pCCD.f90 | 367 ++++++++++++++++++++++++++++++++++++++++-------- 4 files changed, 315 insertions(+), 76 deletions(-) delete mode 100644 int/CAP.dat diff --git a/int/CAP.dat b/int/CAP.dat deleted file mode 100644 index e8a0eda..0000000 --- a/int/CAP.dat +++ /dev/null @@ -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 diff --git a/mol/H2.xyz b/mol/H2.xyz index 8c244ab..48f4296 100644 --- a/mol/H2.xyz +++ b/mol/H2.xyz @@ -1,4 +1,4 @@ 2 -H 0.0000 0.0000 -0.37500000 -H 0.0000 0.0000 0.37500000 +H 0.00000000 -0.37500000 0.00000000 +H 0.00000000 0.37500000 0.00000000 diff --git a/mol/N2.xyz b/mol/N2.xyz index ecc0432..fb38cec 100644 --- a/mol/N2.xyz +++ b/mol/N2.xyz @@ -1,4 +1,4 @@ 2 -N 0.0000 0.0000 0.0000 -N 0.0000 0.0000 1.1007 +N 0.0000 0.0000 -0.5475132 +N 0.0000 0.0000 0.5475132 diff --git a/src/CC/pCCD.f90 b/src/CC/pCCD.f90 index 5306944..54c3c17 100644 --- a/src/CC/pCCD.f90 +++ b/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 j=1,nO-nC + do i=1,O + 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 a=1,nV-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 - do a=1,nV-nR - do b=1,nV-nR + do a=1,V + 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)) - - 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 + allocate(t(O,V),r(O,V),yO(O,O),yV(V,V)) ! 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 - - allocate(r(nO-nC,nV-nR),y(nO-nC,nO-nC)) +!------------------------------------------------------------------------ +! Compute t ampltiudes +!------------------------------------------------------------------------ Conv = 1d0 nSCF = 0 + ECC = ERHF + EcCC = 0d0 - n_diis = 0 - t_diis(:,:) = 0d0 - error_diis(:,:) = 0d0 + n_diis = 0 + t(:,:) = 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 - 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 + yO(:,:) = matmul(t,transpose(OOVV)) ! 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) & - - 2d0*(2d0*OVOV(i,a) - OVVO(i,a) - OOVV(i,a)*t(i,a))*t(i,a) + r(:,:) = OOVV(:,:) + 2d0*delta_OV(:,:)*t(:,:) & + - 2d0*(2d0*OVOV(:,:) - OVVO(:,:) - OOVV(:,:)*t(:,:))*t(:,:) - 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) + 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)*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 - do i=1,nO-nC - do a=1,nV-nR - EcCC = EcCC + OOVV(i,a)*t(i,a) - end do - end do + EcCC = trace_matrix(V,matmul(transpose(OOVV),t)) ! Dump results ECC = ERHF + EcCC - ! DIIS extrapolation + ! DIIS extrapolation -! n_diis = min(n_diis+1,max_diis) -! call DIIS_extrapolation(rcond,nO*nV,nO*nV,n_diis,error_diis,t_diis,-0.5d0*r/delta_OV,t) + if(max_diis > 1) then - ! 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 @@ -211,9 +320,155 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' stop - + 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)