From 5b4dd9982ced04055d6efd10530aa31c8391c7df Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Thu, 5 Nov 2020 15:51:33 +0100 Subject: [PATCH] fix bug in pCC --- input/methods | 6 +++--- mol/h2.xyz | 2 +- src/CC/form_delta_OV.f90 | 7 +++++-- src/CC/pCCD.f90 | 31 +++++++++++++++++++++++++------ src/QuAcK/QuAcK.f90 | 3 +-- 5 files changed, 35 insertions(+), 14 deletions(-) diff --git a/input/methods b/input/methods index af2bd81..db30272 100644 --- a/input/methods +++ b/input/methods @@ -1,11 +1,11 @@ # RHF UHF KS MOM - F T F F + T F F F # MP2* MP3 MP2-F12 F F F # CCD CCSD CCSD(T) F F F # drCCD rCCD lCCD pCCD - F F F F + F F F T # CIS* CIS(D) CID CISD F F F F # RPA* RPAx* ppRPA @@ -13,7 +13,7 @@ # G0F2 evGF2 G0F3 evGF3 F F F F # G0W0* evGW* qsGW* - F F T + F F F # G0T0 evGT qsGT F F F # MCMP2 diff --git a/mol/h2.xyz b/mol/h2.xyz index dcaf1a5..6a4e902 100644 --- a/mol/h2.xyz +++ b/mol/h2.xyz @@ -1,4 +1,4 @@ 2 H 0.0 0.0 0.0 -H 0.0 0.0 3.0 +H 0.0 0.0 0.7 diff --git a/src/CC/form_delta_OV.f90 b/src/CC/form_delta_OV.f90 index a4a671a..565534f 100644 --- a/src/CC/form_delta_OV.f90 +++ b/src/CC/form_delta_OV.f90 @@ -6,7 +6,10 @@ subroutine form_delta_OV(nC,nO,nV,nR,eO,eV,delta) ! Input variables - integer,intent(in) :: nC,nO,nV,nR + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR double precision,intent(in) :: eO(nO-nC) double precision,intent(in) :: eV(nV-nR) @@ -16,7 +19,7 @@ subroutine form_delta_OV(nC,nO,nV,nR,eO,eV,delta) ! Output variables - double precision,intent(out) :: delta(nO,nV) + double precision,intent(out) :: delta(nO-nC,nV-nR) do i=1,nO-nC do a=1,nV-nR diff --git a/src/CC/pCCD.f90 b/src/CC/pCCD.f90 index 29e4629..03b6e2d 100644 --- a/src/CC/pCCD.f90 +++ b/src/CC/pCCD.f90 @@ -91,6 +91,14 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) t(:,:) = -0.5d0*OOVV(:,:)/delta_OV(:,:) + EcCCD = 0d0 + do i=1,nO-nC + do a=1,nV-nR + EcCCD = EcCCD + OOVV(i,a)*t(i,a) + end do + end do + print*,'Ec = ',EcCCD + ! Memory allocation for DIIS allocate(error_diis((nO-nC)*(nV-nR),max_diis),t_diis((nO-nC)*(nV-nR),max_diis)) @@ -125,7 +133,14 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) ! Form intermediate array - y(:,:) = matmul(t,transpose(OOVV)) + 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 ! Compute residual @@ -156,7 +171,12 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) ! Compute correlation energy - EcCCD = trace_matrix(nO,matmul(t,transpose(OOVV))) + EcCCD = 0d0 + do i=1,nO-nC + do a=1,nV-nR + EcCCD = EcCCD + OOVV(i,a)*t(i,a) + end do + end do ! Dump results @@ -164,12 +184,12 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) ! 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) +! 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) ! Reset DIIS if required - if(abs(rcond) < 1d-15) n_diis = 0 +! if(abs(rcond) < 1d-15) n_diis = 0 write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & '|',nSCF,'|',ECCD+ENuc,'|',EcCCD,'|',Conv,'|' @@ -188,7 +208,6 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*)' Convergence failed ' write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' - write(*,*) stop diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 4ce5392..395f3e9 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -605,8 +605,7 @@ program QuAcK if(do_pCCD) then call cpu_time(start_CCD) - call pCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR, & - ERI_MO,ENuc,ERHF,eHF) + call pCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call cpu_time(end_CCD) t_CCD = end_CCD - start_CCD