PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-09-27 20:11:05 +02:00
Code Releases Activity

saving work in pCCD

Browse Source
This commit is contained in:
Pierre-Francois Loos 2024-08-30 16:26:59 +02:00
parent 9910e6e889
commit 08bf6632df
3 changed files with 603 additions and 562 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
src/CC/RCC.f90
View File

@ -1,5 +1,5 @@
subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, & 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 ! 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) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas) double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: Hc(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 ! Local variables
@ -47,7 +48,7 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d
if(doCCD) then if(doCCD) then
call wall_time(start_CC) 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) call wall_time(end_CC)
t_CC = end_CC - start_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 wall_time(start_CC)
call DCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR, & 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) call wall_time(end_CC)
t_CC = end_CC - start_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 if(doCCSD) then
call wall_time(start_CC) 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) call wall_time(end_CC)
t_CC = end_CC - start_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 if(dodrCCD) then
call wall_time(start_CC) 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) call wall_time(end_CC)
t_CC = end_CC - start_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 if(dorCCD) then
call wall_time(start_CC) 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) call wall_time(end_CC)
t_CC = end_CC - start_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 if(docrCCD) then
call wall_time(start_CC) 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) call wall_time(end_CC)
t_CC = end_CC - start_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 if(dolCCD) then
call wall_time(start_CC) 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) call wall_time(end_CC)
t_CC = end_CC - start_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 if(dopCCD) then
call wall_time(start_CC) 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) call wall_time(end_CC)
t_CC = end_CC - start_CC t_CC = end_CC - start_CC

321
src/CC/pCCD.f90
View File

@ -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 ! 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 logical,intent(in) :: dotest
integer,intent(in) :: maxSCF integer,intent(in) :: maxIt
integer,intent(in) :: max_diis integer,intent(in) :: max_diis
double precision,intent(in) :: thresh double precision,intent(in) :: thresh
integer,intent(in) :: nBas,nC,nO,nV,nR integer,intent(in) :: nBas
double precision,intent(in) :: ENuc,ERHF 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) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas) double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: Hc(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 ! Local variables
integer :: p,q,r,s,t,u integer :: p,q,r,s,t,u,w
integer :: pq,rs integer :: pq,rs
integer :: i,j,a,b integer :: i,j,a,b
integer :: nSCF integer :: nItAmp
double precision :: Conv integer :: nItOrb
double precision :: CvgAmp
double precision :: CvgOrb
double precision :: ECC double precision :: ECC
double precision :: EcCC 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 :: tr_2rdm
double precision :: E1,E2 double precision :: E1,E2
double precision,allocatable :: h(:,:)
double precision,allocatable :: c(:,:) double precision,allocatable :: c(:,:)
double precision,allocatable :: h(:,:)
double precision,allocatable :: ERI_MO(:,:,:,:)
double precision,allocatable :: grad(:) double precision,allocatable :: grad(:)
double precision,allocatable :: tmp(:,:,:,:) double precision,allocatable :: tmp(:,:,:,:)
double precision,allocatable :: hess(:,:) double precision,allocatable :: hess(:,:)
double precision,allocatable :: hessInv(:,:) double precision,allocatable :: hessInv(:,:)
double precision,allocatable :: kappa(:,:) double precision,allocatable :: Kap(:,:)
double precision,allocatable :: ekappa(:,:) double precision,allocatable :: ExpKap(:,:)
integer :: O,V,N integer :: O,V,N
integer :: n_diis 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 ! Hello world
write(*,*) write(*,*)
write(*,*)'**************************************' write(*,*)'*******************************'
write(*,*)'| pair CCD calculation |' write(*,*)'* Restricted pCCD Calculation *'
write(*,*)'**************************************' write(*,*)'*******************************'
write(*,*) write(*,*)
! Useful quantities ! Useful quantities
@ -88,60 +96,76 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
V = nV - nR V = nV - nR
N = O + V 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)) 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))
c(:,:) = cHF(nC+1:nBas-nR,nC+1:nBas-nR)
CvgOrb = 1d0
nItOrb = 0
write(*,*)
write(*,*)'----------------------------------------------------'
write(*,*)'| Orbital Optimization for pCCD |'
write(*,*)'----------------------------------------------------'
do while(CvgOrb > thresh .and. nItOrb < 1)
nItOrb = nItOrb + 1
! Transform integrals
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
eO(:) = eHF(nC+1:nO) eO(:) = eHF(nC+1:nO)
eV(:) = eHF(nO+1:nBas-nR) eV(:) = eHF(nO+1:nBas-nR)
call form_delta_OV(nC,nO,nV,nR,eO,eV,delta_OV) do i=1,O
do a=1,V
delta_OV(i,a) = eV(a) - eO(i)
end do
end do
! Create integral batches ! Create integral batches
allocate(OOOO(O,O),OOVV(O,V),OVOV(O,V),OVVO(O,V),VVVV(V,V))
do i=1,O do i=1,O
do j=1,O do j=1,O
OOOO(i,j) = ERI(nC+i,nC+i,nC+j,nC+j) OOOO(i,j) = ERI_MO(i,i,j,j)
end do end do
end do end do
do i=1,O do i=1,O
do a=1,V do a=1,V
OOVV(i,a) = ERI(nC+i,nC+i,nO+a,nO+a) OOVV(i,a) = ERI_MO(i,i,O+a,O+a)
OVOV(i,a) = ERI(nC+i,nO+a,nC+i,nO+a) OVOV(i,a) = ERI_MO(i,O+a,i,O+a)
OVVO(i,a) = ERI(nC+i,nO+a,nO+a,nC+i) OVVO(i,a) = ERI_MO(i,O+a,O+a,i)
end do end do
end do end do
do a=1,V do a=1,V
do b=1,V do b=1,V
VVVV(a,b) = ERI(nO+a,nO+a,nO+b,nO+b) VVVV(a,b) = ERI_MO(O+a,O+a,O+b,O+b)
end do end do
end do end do
! Initialization !----------------------------!
! Star Loop for t amplitudes !
allocate(t2(O,V),r2(O,V),yO(O,O),yV(V,V)) !----------------------------!
! Memory allocation for DIIS
allocate(t2(O,V),r2(O,V),yO(O,O))
allocate(err_diis(O*V,max_diis),t2_diis(O*V,max_diis)) allocate(err_diis(O*V,max_diis),t2_diis(O*V,max_diis))
!------------------------------------------------------------------------ CvgAmp = 1d0
! Start orbital optimization nItAmp = 0
!------------------------------------------------------------------------
allocate(c(N,N))
c(:,:) = cHF(:,:)
!------------------------------------------------------------------------
! Compute t ampltiudes
!------------------------------------------------------------------------
Conv = 1d0
nSCF = 0
ECC = ERHF ECC = ERHF
EcCC = 0d0 EcCC = 0d0
@ -150,9 +174,6 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
t2_diis(:,:) = 0d0 t2_diis(:,:) = 0d0
err_diis(:,:) = 0d0 err_diis(:,:) = 0d0
!------------------------------------------------------------------------
! Main SCF loop
!------------------------------------------------------------------------
write(*,*) write(*,*)
write(*,*)'----------------------------------------------------' write(*,*)'----------------------------------------------------'
write(*,*)'| pCCD calculation: t amplitudes |' write(*,*)'| pCCD calculation: t amplitudes |'
@ -161,11 +182,11 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
'|','#','|','E(pCCD)','|','Ec(pCCD)','|','Conv','|' '|','#','|','E(pCCD)','|','Ec(pCCD)','|','Conv','|'
write(*,*)'----------------------------------------------------' write(*,*)'----------------------------------------------------'
do while(Conv > thresh .and. nSCF < maxSCF) do while(CvgAmp > thresh .and. nItAmp < maxIt)
! Increment ! Increment
nSCF = nSCF + 1 nItAmp = nItAmp + 1
! Form intermediate array ! Form intermediate array
@ -173,7 +194,6 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
! Compute residual ! Compute residual
r2(:,:) = OOVV(:,:) + 2d0*delta_OV(:,:)*t2(:,:) & r2(:,:) = OOVV(:,:) + 2d0*delta_OV(:,:)*t2(:,:) &
- 2d0*(2d0*OVOV(:,:) - OVVO(:,:) - OOVV(:,:)*t2(:,:))*t2(:,:) - 2d0*(2d0*OVOV(:,:) - OVVO(:,:) - OOVV(:,:)*t2(:,:))*t2(:,:)
@ -193,7 +213,7 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
! Check convergence ! Check convergence
Conv = maxval(abs(r2(:,:))) CvgAmp = maxval(abs(r2(:,:)))
! Update amplitudes ! Update amplitudes
@ -201,7 +221,12 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
! Compute correlation energy ! Compute correlation energy
EcCC = trace_matrix(V,matmul(transpose(OOVV),t2)) EcCC = 0d0
do i=1,O
do a=1,V
EcCC = EcCC + OOVV(i,a)*t2(i,a)
end do
end do
! Dump results ! Dump results
@ -217,53 +242,45 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
end if 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)') & 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,'|' '|',nItAmp,'|',ECC+ENuc,'|',EcCC,'|',CvgAmp,'|'
end do end do
write(*,*)'----------------------------------------------------' write(*,*)'----------------------------------------------------'
!------------------------------------------------------------------------
! End of SCF loop
!------------------------------------------------------------------------
! Did it actually converge? !---------------------------!
! End Loop for t amplitudes !
!---------------------------!
if(nSCF == maxSCF) then deallocate(r2,yO)
deallocate(err_diis,t2_diis)
! Did it actually converge?
if(nItAmp == maxIt) then
write(*,*) write(*,*)
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)' Convergence failed for t ampitudes ' write(*,*)'! Convergence failed for t ampitudes !'
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
stop stop
end if end if
! Deallocate memory !-----------------------------!
! Start Loop for z amplitudes !
deallocate(err_diis,t2_diis) !-----------------------------!
! Memory allocation
allocate(z2(O,V))
! Memory allocation for DIIS
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)) allocate(err_diis(O*V,max_diis),z2_diis(O*V,max_diis))
!------------------------------------------------------------------------ CvgAmp = 1d0
! Compute z ampltiudes nItAmp = 0
!------------------------------------------------------------------------
Conv = 1d0
nSCF = 0
n_diis = 0 n_diis = 0
z2_diis(:,:) = 0d0 z2_diis(:,:) = 0d0
err_diis(:,:) = 0d0 err_diis(:,:) = 0d0
!------------------------------------------------------------------------
! Main SCF loop
!------------------------------------------------------------------------
write(*,*) write(*,*)
write(*,*)'----------------------------------------------------' write(*,*)'----------------------------------------------------'
write(*,*)'| pCCD calculation: z amplitudes |' write(*,*)'| pCCD calculation: z amplitudes |'
@ -272,11 +289,11 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
'|','#','|','E(pCCD)','|','Ec(pCCD)','|','Conv','|' '|','#','|','E(pCCD)','|','Ec(pCCD)','|','Conv','|'
write(*,*)'----------------------------------------------------' write(*,*)'----------------------------------------------------'
do while(Conv > thresh .and. nSCF < maxSCF) do while(CvgAmp > thresh .and. nItAmp < maxIt)
! Increment ! Increment
nSCF = nSCF + 1 nItAmp = nItAmp + 1
! Form intermediate array ! Form intermediate array
@ -306,7 +323,7 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
! Check convergence ! Check convergence
Conv = maxval(abs(r2(:,:))) CvgAmp = maxval(abs(r2(:,:)))
! Update amplitudes ! Update amplitudes
@ -322,45 +339,44 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
end if 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)') & 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,'|' '|',nItAmp,'|',ECC+ENuc,'|',EcCC,'|',CvgAmp,'|'
end do end do
write(*,*)'----------------------------------------------------' write(*,*)'----------------------------------------------------'
write(*,*) write(*,*)
!------------------------------------------------------------------------
! End of SCF loop
!------------------------------------------------------------------------
! Did it actually converge? !---------------------------!
! End Loop for z ampltiudes !
!---------------------------!
if(nSCF == maxSCF) then deallocate(r2,yO,yV)
deallocate(err_diis,z2_diis)
! Did it actually converge?
if(nItAmp == maxIt) then
write(*,*) write(*,*)
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)' Convergence failed ' write(*,*)'! Convergence failed for z ampltiudes !'
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
stop stop
end if end if
! Deallocate memory !--------------------------!
! Compute density matrices !
deallocate(err_diis,z2_diis,r2) !--------------------------!
!--------------------------!
! Compute density matrices !
!--------------------------!
allocate(rdm1(N,N),rdm2(N,N,N,N))
allocate(xOO(O,O),xVV(V,V),xOV(O,V)) allocate(xOO(O,O),xVV(V,V),xOV(O,V))
xOO(:,:) = matmul(t2,transpose(z2)) xOO(:,:) = matmul(t2,transpose(z2))
xVV(:,:) = matmul(transpose(z2),t2) xVV(:,:) = matmul(transpose(z2),t2)
xOV(:,:) = matmul(t2,matmul(transpose(z2),t2)) xOV(:,:) = matmul(t2,matmul(transpose(z2),t2))
! Form 1RDM ! Form 1RDM
allocate(rdm1(N,N))
rdm1(:,:) = 0d0 rdm1(:,:) = 0d0
@ -372,7 +388,7 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
rdm1(O+a,O+a) = 2d0*xVV(a,a) rdm1(O+a,O+a) = 2d0*xVV(a,a)
end do end do
! Check 1RDM ! Check 1RDM
tr_1rdm = trace_matrix(N,rdm1) tr_1rdm = trace_matrix(N,rdm1)
write(*,'(A25,F16.10)') ' --> Trace of the 1RDM = ',tr_1rdm write(*,'(A25,F16.10)') ' --> Trace of the 1RDM = ',tr_1rdm
@ -381,12 +397,10 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
write(*,*) ' !!! Your 1RDM seems broken !!! ' write(*,*) ' !!! Your 1RDM seems broken !!! '
write(*,*) write(*,*)
! write(*,*) '1RDM is diagonal at the pCCD level:' ! write(*,*) '1RDM is diagonal at the pCCD level:'
! call matout(N,N,rdm1) ! call matout(N,N,rdm1)
! Form 2RM ! Form 2RM
allocate(rdm2(N,N,N,N))
rdm2(:,:,:,:) = 0d0 rdm2(:,:,:,:) = 0d0
@ -482,7 +496,7 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
rdm2(O+a,O+a,O+a,O+a) = 2d0*xVV(a,a) rdm2(O+a,O+a,O+a,O+a) = 2d0*xVV(a,a)
end do end do
! Check 2RDM ! Check 2RDM
tr_2rdm = trace_matrix(N**2,rdm2) tr_2rdm = trace_matrix(N**2,rdm2)
write(*,'(A25,F16.10)') ' --> Trace of the 2RDM = ',tr_2rdm write(*,'(A25,F16.10)') ' --> Trace of the 2RDM = ',tr_2rdm
@ -491,13 +505,13 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
write(*,*) ' !!! Your 2RDM seems broken !!! ' write(*,*) ' !!! Your 2RDM seems broken !!! '
write(*,*) write(*,*)
! write(*,*) '2RDM is not diagonal at the pCCD level:' ! write(*,*) '2RDM is not diagonal at the pCCD level:'
! call matout(N**2,N**2,rdm2) ! call matout(N**2,N**2,rdm2)
! Compute electronic energy deallocate(xOO,xVV,xOV)
deallocate(t2,z2)
allocate(h(N,N)) ! Compute electronic energy
h = matmul(transpose(cHF),matmul(Hc,cHF))
E1 = 0d0 E1 = 0d0
E2 = 0d0 E2 = 0d0
@ -507,7 +521,7 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
E1 = E1 + rdm1(p,q)*h(p,q) E1 = E1 + rdm1(p,q)*h(p,q)
do r=1,N do r=1,N
do s=1,N do s=1,N
E2 = E2 + rdm2(p,q,r,s)*ERI(p,q,r,s) E2 = E2 + rdm2(p,q,r,s)*ERI_MO(p,q,r,s)
end do end do
end do end do
end do end do
@ -521,7 +535,9 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
write(*,'(A25,F16.10)') ' Total energy = ',E1 + E2 + ENuc write(*,'(A25,F16.10)') ' Total energy = ',E1 + E2 + ENuc
write(*,*) write(*,*)
! Compute gradient !--------------------------!
! Compute orbital gradient !
!--------------------------!
allocate(grad(N**2)) allocate(grad(N**2))
@ -540,7 +556,7 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
do r=1,N do r=1,N
do s=1,N do s=1,N
do t=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)) 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
@ -552,13 +568,16 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
call matout(N,N,grad) call matout(N,N,grad)
write(*,*) write(*,*)
! Convergence ! Check convergence of orbital optimization
Conv = maxval(abs(grad)) CvgOrb = maxval(abs(grad))
write(*,*) ' Convergence of orbtial gradient = ',Conv write(*,*) ' Iteration',nItOrb,'for pCCD orbital optimization'
write(*,*) ' Convergence of orbital gradient = ',CvgOrb
write(*,*) write(*,*)
! Compute Hessian !-------------------------!
! Compute orbital Hessian !
!-------------------------!
allocate(hess(N**2,N**2),tmp(N,N,N,N)) allocate(hess(N**2,N**2),tmp(N,N,N,N))
@ -581,9 +600,9 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
end do end do
do u=1,N do u=1,N
do v=1,N do w=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) 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
end do end do
@ -592,19 +611,19 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
do u=1,N do u=1,N
tmp(p,q,r,s) = tmp(p,q,r,s) - ( & 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_MO(s,t,p,u)*rdm2(r,t,q,u) + ERI_MO(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) ) + 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
end do end do
do t=1,N do t=1,N
do u=1,N do u=1,N
do v=1,N do w=1,N
tmp(p,q,r,s) = tmp(p,q,r,s) + 0.5d0*( & 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(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(q,t,u,v)*rdm2(r,t,u,v) + ERI(u,v,r,t)*rdm2(u,v,q,t)) ) + 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
@ -616,7 +635,7 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
end do end do
end do end do
! Flatten Hessian matrix and add permutations ! Flatten Hessian matrix and add permutations
pq = 0 pq = 0
do p=1,N do p=1,N
@ -631,7 +650,7 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
rs = rs + 1 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,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) !! 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
@ -639,21 +658,17 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
end do end do
end do end do
call matout(N**2,N**2,hess) deallocate(rdm1,rdm2,tmp)
! allocate(eig(N**2))
! call diagonalize_matrix(N**2,hess,eig)
! call vecout(N**2,eig)
allocate(hessInv(N**2,N**2)) allocate(hessInv(N**2,N**2))
call inverse_matrix(N**2,hess,hessInv) call inverse_matrix(N**2,hess,hessInv)
allocate(kappa(N,N),ekappa(N,N)) deallocate(hess)
kappa(:,:) = 0d0 allocate(Kap(N,N))
Kap(:,:) = 0d0
pq = 0 pq = 0
do p=1,N do p=1,N
@ -667,7 +682,7 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
rs = rs + 1 rs = rs + 1
kappa(p,q) = kappa(p,q) + hessInv(pq,rs)*grad(rs) Kap(p,q) = Kap(p,q) - hessInv(pq,rs)*grad(rs)
end do end do
end do end do
@ -675,26 +690,50 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,
end do end do
end do end do
deallocate(hessInv,grad)
write(*,*) 'kappa' write(*,*) 'kappa'
call matout(N,N,kappa) call matout(N,N,Kap)
write(*,*) write(*,*)
call matrix_exponential(N,kappa,ekappa) allocate(ExpKap(N,N))
call matrix_exponential(N,Kap,ExpKap)
deallocate(Kap)
write(*,*) 'e^kappa' write(*,*) 'e^kappa'
call matout(N,N,ekappa) call matout(N,N,ExpKap)
write(*,*) write(*,*)
write(*,*) 'Old orbitals' write(*,*) 'Old orbitals'
call matout(N,N,c) call matout(N,N,c)
write(*,*) write(*,*)
c = matmul(c,ekappa) c = matmul(c,ExpKap)
deallocate(ExpKap)
write(*,*) 'Rotated orbitals' write(*,*) 'Rotated orbitals'
call matout(N,N,c) call matout(N,N,c)
write(*,*) write(*,*)
end do
!-----------------------------------!
! End Loop for orbital optimization !
!-----------------------------------!
! Did it actually converge?
if(nItOrb == maxIt) then
write(*,*)
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)'! Convergence failed for orbital optimization !'
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
stop
end if
! Testing zone ! Testing zone
if(dotest) then if(dotest) then

2
src/QuAcK/RQuAcK.f90
View File

@ -228,7 +228,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
call wall_time(start_CC) call wall_time(start_CC)
call RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, & 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) call wall_time(end_CC)
t_CC = end_CC - start_CC t_CC = end_CC - start_CC