4
1
mirror of https://github.com/pfloos/quack synced 2024-07-25 04:07:35 +02:00

fix bug in pCC

This commit is contained in:
Pierre-Francois Loos 2020-11-05 15:51:33 +01:00
parent fedc3a5d23
commit 5b4dd9982c
5 changed files with 35 additions and 14 deletions

View File

@ -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

View File

@ -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

View File

@ -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

View File

@ -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

View File

@ -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