mirror of
https://github.com/pfloos/quack
synced 2024-11-03 20:53:53 +01:00
plenty of CC for GHF reference
This commit is contained in:
parent
ce57b4b494
commit
3318daf241
@ -1,20 +1,20 @@
|
||||
# RHF UHF GHF ROHF
|
||||
T F F F
|
||||
F F T F
|
||||
# MP2 MP3
|
||||
T T
|
||||
# CCD pCCD DCD CCSD CCSD(T)
|
||||
T T T T F
|
||||
F F F F F
|
||||
# drCCD rCCD crCCD lCCD
|
||||
T T T T
|
||||
F F F F
|
||||
# CIS CIS(D) CID CISD FCI
|
||||
T F F F F
|
||||
F F F F F
|
||||
# phRPA phRPAx crRPA ppRPA
|
||||
T T T T
|
||||
F F F F
|
||||
# G0F2 evGF2 qsGF2 G0F3 evGF3
|
||||
T T F F F
|
||||
F F F F F
|
||||
# G0W0 evGW qsGW SRG-qsGW ufG0W0 ufGW
|
||||
T T F F F F
|
||||
F F F F F F
|
||||
# G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh
|
||||
T T F F F F
|
||||
F F F F F F
|
||||
# Rtest Utest Gtest
|
||||
T F F
|
||||
F F F
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
# HF: maxSCF thresh DIIS guess mix shift stab search
|
||||
10000 0.0000001 5 1 0.0 0.0 F F
|
||||
10000 0.0000001 5 2 0.0 0.0 F T
|
||||
# MP: reg
|
||||
F
|
||||
# CC: maxSCF thresh DIIS
|
||||
|
||||
172
src/CC/GCC.f90
Normal file
172
src/CC/GCC.f90
Normal file
@ -0,0 +1,172 @@
|
||||
subroutine GCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, &
|
||||
maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EGHF,eHF)
|
||||
|
||||
! Generalized Coupled-cluster module
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
logical :: dotest
|
||||
|
||||
logical :: doCCD
|
||||
logical :: dopCCD
|
||||
logical :: doDCD
|
||||
logical :: doCCSD
|
||||
logical :: doCCSDT
|
||||
logical :: dodrCCD
|
||||
logical :: dorCCD
|
||||
logical :: docrCCD
|
||||
logical :: dolCCD
|
||||
|
||||
integer,intent(in) :: maxSCF
|
||||
integer,intent(in) :: max_diis
|
||||
double precision,intent(in) :: thresh
|
||||
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC(nspin)
|
||||
integer,intent(in) :: nO(nspin)
|
||||
integer,intent(in) :: nV(nspin)
|
||||
integer,intent(in) :: nR(nspin)
|
||||
double precision,intent(in) :: ENuc
|
||||
double precision,intent(in) :: EGHF
|
||||
double precision,intent(in) :: eHF(nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
double precision :: start_CC ,end_CC ,t_CC
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Perform CCD calculation
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
if(doCCD) then
|
||||
|
||||
call wall_time(start_CC)
|
||||
call GCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EGHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCD = ',t_CC,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Perform DCD calculation
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
if(doDCD) then
|
||||
|
||||
call wall_time(start_CC)
|
||||
! call DCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR, &
|
||||
! ERI,ENuc,EGHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for DCD = ',t_CC,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Perform CCSD or CCSD(T) calculation
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
if(doCCSDT) doCCSD = .true.
|
||||
|
||||
if(doCCSD) then
|
||||
|
||||
call wall_time(start_CC)
|
||||
call GCCSD(dotest,maxSCF,thresh,max_diis,doCCSDT,nBas,nC,nO,nV,nR,ERI,ENuc,EGHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCSD or CCSD(T)= ',t_CC,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Perform direct ring CCD calculation
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
if(dodrCCD) then
|
||||
|
||||
call wall_time(start_CC)
|
||||
call drGCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EGHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for direct ring CCD = ',t_CC,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Perform ring CCD calculation
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
if(dorCCD) then
|
||||
|
||||
call wall_time(start_CC)
|
||||
call rGCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EGHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for rCCD = ',t_CC,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Perform crossed-ring CCD calculation
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
if(docrCCD) then
|
||||
|
||||
call wall_time(start_CC)
|
||||
call crGCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EGHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for crossed-ring CCD = ',t_CC,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Perform ladder CCD calculation
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
if(dolCCD) then
|
||||
|
||||
call wall_time(start_CC)
|
||||
call lGCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EGHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ladder CCD = ',t_CC,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Perform pair CCD calculation
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
if(dopCCD) then
|
||||
|
||||
call wall_time(start_CC)
|
||||
! call pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EGHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pair CCD = ',t_CC,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
end subroutine
|
||||
262
src/CC/GCCD.f90
Normal file
262
src/CC/GCCD.f90
Normal file
@ -0,0 +1,262 @@
|
||||
subroutine GCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EGHF,eHF)
|
||||
|
||||
! Generalized CCD module
|
||||
|
||||
implicit none
|
||||
|
||||
! Input variables
|
||||
|
||||
logical,intent(in) :: dotest
|
||||
|
||||
integer,intent(in) :: maxSCF
|
||||
integer,intent(in) :: max_diis
|
||||
double precision,intent(in) :: thresh
|
||||
|
||||
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) :: EGHF
|
||||
double precision,intent(in) :: eHF(nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: nSCF
|
||||
double precision :: Conv
|
||||
double precision :: EcMP2,EcMP3,EcMP4
|
||||
double precision :: ECC,EcCC
|
||||
double precision,allocatable :: dbERI(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: eO(:)
|
||||
double precision,allocatable :: eV(:)
|
||||
double precision,allocatable :: delta_OOVV(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: OOOO(:,:,:,:)
|
||||
double precision,allocatable :: OOVV(:,:,:,:)
|
||||
double precision,allocatable :: OVOV(:,:,:,:)
|
||||
double precision,allocatable :: VVVV(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: OVVO(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: X1(:,:,:,:)
|
||||
double precision,allocatable :: X2(:,:)
|
||||
double precision,allocatable :: X3(:,:)
|
||||
double precision,allocatable :: X4(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: u(:,:,:,:)
|
||||
double precision,allocatable :: v(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: r(:,:,:,:)
|
||||
double precision,allocatable :: t(:,:,:,:)
|
||||
|
||||
integer :: n_diis
|
||||
double precision :: rcond
|
||||
double precision,allocatable :: error_diis(:,:)
|
||||
double precision,allocatable :: t_diis(:,:)
|
||||
|
||||
logical :: do_EE_EOM_CC_1h1p = .false.
|
||||
logical :: do_EA_EOM_CC_1p = .false.
|
||||
logical :: do_IP_EOM_CC_1h = .false.
|
||||
logical :: do_DEA_EOM_CC_2p = .false.
|
||||
logical :: do_DIP_EOM_CC_2h = .false.
|
||||
|
||||
! Hello world
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'*******************************'
|
||||
write(*,*)'* Generalized CCD Calculation *'
|
||||
write(*,*)'*******************************'
|
||||
write(*,*)
|
||||
|
||||
! Antysymmetrize ERIs
|
||||
|
||||
allocate(dbERI(nBas,nBas,nBas,nBas))
|
||||
|
||||
call antisymmetrize_ERI(2,nBas,ERI,dbERI)
|
||||
|
||||
! Form energy denominator
|
||||
|
||||
allocate(eO(nO-nC),eV(nV-nR))
|
||||
allocate(delta_OOVV(nO-nC,nO-nC,nV-nR,nV-nR))
|
||||
|
||||
eO(:) = eHF(nC+1:nO)
|
||||
eV(:) = eHF(nO+1:nBas-nR)
|
||||
|
||||
call form_delta_OOVV(nC,nO,nV,nR,eO,eV,delta_OOVV)
|
||||
|
||||
! Create integral batches
|
||||
|
||||
allocate(OOOO(nO-nC,nO-nC,nO-nC,nO-nC),OOVV(nO-nC,nO-nC,nV-nR,nV-nR), &
|
||||
OVOV(nO-nC,nV-nR,nO-nC,nV-nR),VVVV(nV-nR,nV-nR,nV-nR,nV-nR))
|
||||
|
||||
OOOO(:,:,:,:) = dbERI(nC+1:nO ,nC+1:nO ,nC+1:nO ,nC+1:nO )
|
||||
OOVV(:,:,:,:) = dbERI(nC+1:nO ,nC+1:nO ,nO+1:nBas-nR,nO+1:nBas-nR)
|
||||
OVOV(:,:,:,:) = dbERI(nC+1:nO ,nO+1:nBas-nR,nC+1:nO ,nO+1:nBas-nR)
|
||||
VVVV(:,:,:,:) = dbERI(nO+1:nBas-nR,nO+1:nBas-nR,nO+1:nBas-nR,nO+1:nBas-nR)
|
||||
|
||||
allocate(OVVO(nO-nC,nV-nR,nV-nR,nO-nC))
|
||||
OVVO(:,:,:,:) = dbERI(nC+1:nO,nO+1:nBas-nR,nO+1:nBas-nR,nC+1:nO)
|
||||
|
||||
deallocate(dbERI)
|
||||
|
||||
! MP2 guess amplitudes
|
||||
|
||||
allocate(t(nO-nC,nO-nC,nV-nR,nV-nR))
|
||||
|
||||
t(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:)
|
||||
|
||||
call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t,EcMP2)
|
||||
EcMP4 = 0d0
|
||||
|
||||
! Memory allocation for DIIS
|
||||
|
||||
allocate(error_diis((nO-nR)**2*(nV-nR)**2,max_diis),t_diis((nO-nR)**2*(nV-nR)**2,max_diis))
|
||||
|
||||
! Initialization
|
||||
|
||||
allocate(r(nO-nC,nO-nC,nV-nR,nV-nR),u(nO-nC,nO-nC,nV-nR,nV-nR),v(nO-nC,nO-nC,nV-nR,nV-nR))
|
||||
allocate(X1(nO-nC,nO-nC,nO-nC,nO-nC),X2(nV-nR,nV-nR),X3(nO-nC,nO-nC),X4(nO-nC,nO-nC,nV-nR,nV-nR))
|
||||
|
||||
Conv = 1d0
|
||||
nSCF = 0
|
||||
|
||||
n_diis = 0
|
||||
t_diis(:,:) = 0d0
|
||||
error_diis(:,:) = 0d0
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Main SCF loop
|
||||
!------------------------------------------------------------------------
|
||||
write(*,*)
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,*)'| GCCD calculation |'
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') &
|
||||
'|','#','|','E(GCCD)','|','Ec(GCCD)','|','Conv','|'
|
||||
write(*,*)'----------------------------------------------------'
|
||||
|
||||
do while(Conv > thresh .and. nSCF < maxSCF)
|
||||
|
||||
! Increment
|
||||
|
||||
nSCF = nSCF + 1
|
||||
|
||||
! Form linear array
|
||||
|
||||
call form_u(nC,nO,nV,nR,OOOO,VVVV,OVOV,t,u)
|
||||
|
||||
! Form interemediate arrays
|
||||
|
||||
call form_X(nC,nO,nV,nR,OOVV,t,X1,X2,X3,X4)
|
||||
|
||||
! Form quadratic array
|
||||
|
||||
call form_v(nC,nO,nV,nR,X1,X2,X3,X4,t,v)
|
||||
|
||||
! Compute residual
|
||||
|
||||
r(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t(:,:,:,:) + u(:,:,:,:) + v(:,:,:,:)
|
||||
|
||||
! Check convergence
|
||||
|
||||
Conv = maxval(abs(r(:,:,:,:)))
|
||||
|
||||
! Update amplitudes
|
||||
|
||||
t(:,:,:,:) = t(:,:,:,:) - r(:,:,:,:)/delta_OOVV(:,:,:,:)
|
||||
|
||||
! Compute correlation energy
|
||||
|
||||
call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t,EcCC)
|
||||
|
||||
if(nSCF == 1) call MP3_correlation_energy(nC,nO,nV,nR,OOVV,t,v,delta_OOVV,EcMP3)
|
||||
|
||||
! Dump results
|
||||
|
||||
ECC = EGHF + EcCC
|
||||
|
||||
! DIIS extrapolation
|
||||
|
||||
n_diis = min(n_diis+1,max_diis)
|
||||
call DIIS_extrapolation(rcond,(nO-nC)**2*(nV-nR)**2,(nO-nC)**2*(nV-nR)**2,n_diis,error_diis,t_diis,-r/delta_OOVV,t)
|
||||
|
||||
! Reset DIIS if required
|
||||
|
||||
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,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|'
|
||||
|
||||
enddo
|
||||
write(*,*)'----------------------------------------------------'
|
||||
!------------------------------------------------------------------------
|
||||
! End of SCF loop
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
! Did it actually converge?
|
||||
|
||||
if(nSCF == maxSCF) then
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)' Convergence failed '
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)
|
||||
|
||||
stop
|
||||
|
||||
endif
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,*)' GCCD energy '
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,'(1X,A30,1X,F15.10)')' E(GCCD) = ',ECC
|
||||
write(*,'(1X,A30,1X,F15.10)')' Ec(GCCD) = ',EcCC
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
! Moller-Plesset energies
|
||||
|
||||
write(*,*)
|
||||
write(*,'(1X,A15,1X,F10.6)') 'Ec(GMP2) = ',EcMP2
|
||||
write(*,'(1X,A15,1X,F10.6)') 'Ec(GMP3) = ',EcMP3
|
||||
write(*,'(1X,A15,1X,F10.6)') 'Ec(GMP4-SDQ) = ',EcMP4
|
||||
write(*,*)
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! EOM section
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
! EE-EOM-CCD (1h1p)
|
||||
|
||||
if(do_EE_EOM_CC_1h1p) call EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t)
|
||||
|
||||
! EA-EOM (1p)
|
||||
|
||||
! if(do_EA-EOM-CC_1p) call EA-EOM-CCD_1p()
|
||||
|
||||
! IP-EOM-CCD(1h)
|
||||
|
||||
! if(do_IP-EOM-CC_1h) call IP-EOM-CCD_1h()
|
||||
|
||||
! DEA-EOM (2p)
|
||||
|
||||
if(do_DEA_EOM_CC_2p) call DEA_EOM_CCD_2p(nC,nO,nV,nR,eV,OOVV,VVVV,t)
|
||||
|
||||
! DIP-EOM-CCD(2h)
|
||||
|
||||
if(do_DIP_EOM_CC_2h) call DIP_EOM_CCD_2h(nC,nO,nV,nR,eO,OOVV,OOOO,t)
|
||||
|
||||
! Testing zone
|
||||
|
||||
if(dotest) then
|
||||
|
||||
call dump_test_value('G','GCCD correlation energy',EcCC)
|
||||
|
||||
end if
|
||||
|
||||
end subroutine
|
||||
291
src/CC/GCCSD.f90
Normal file
291
src/CC/GCCSD.f90
Normal file
@ -0,0 +1,291 @@
|
||||
subroutine GCCSD(dotest,maxSCF,thresh,max_diis,doCCSDT,nBas,nC,nO,nV,nR,ERI,ENuc,EGHF,eHF)
|
||||
|
||||
! Generalized CCSD module
|
||||
|
||||
implicit none
|
||||
|
||||
! Input variables
|
||||
|
||||
logical,intent(in) :: dotest
|
||||
|
||||
integer,intent(in) :: maxSCF
|
||||
integer,intent(in) :: max_diis
|
||||
double precision,intent(in) :: thresh
|
||||
|
||||
logical,intent(in) :: doCCSDT
|
||||
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) :: EGHF
|
||||
double precision,intent(in) :: eHF(nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
double precision :: start_CCSDT,end_CCSDT,t_CCSDT
|
||||
integer :: nSCF
|
||||
double precision :: Conv
|
||||
double precision :: EcMP2
|
||||
double precision :: ECC
|
||||
double precision :: EcCC
|
||||
double precision :: EcCCT
|
||||
|
||||
double precision,allocatable :: dbERI(:,:,:,:)
|
||||
double precision,allocatable :: delta_OV(:,:)
|
||||
double precision,allocatable :: delta_OOVV(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: OOOO(:,:,:,:)
|
||||
double precision,allocatable :: OOOV(:,:,:,:)
|
||||
double precision,allocatable :: OVOO(:,:,:,:)
|
||||
double precision,allocatable :: VOOO(:,:,:,:)
|
||||
double precision,allocatable :: OOVV(:,:,:,:)
|
||||
double precision,allocatable :: OVVO(:,:,:,:)
|
||||
double precision,allocatable :: OVVV(:,:,:,:)
|
||||
double precision,allocatable :: VOVV(:,:,:,:)
|
||||
double precision,allocatable :: VVVO(:,:,:,:)
|
||||
double precision,allocatable :: VVVV(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: eO(:)
|
||||
double precision,allocatable :: eV(:)
|
||||
double precision,allocatable :: hvv(:,:)
|
||||
double precision,allocatable :: hoo(:,:)
|
||||
double precision,allocatable :: hvo(:,:)
|
||||
double precision,allocatable :: gvv(:,:)
|
||||
double precision,allocatable :: goo(:,:)
|
||||
double precision,allocatable :: aoooo(:,:,:,:)
|
||||
double precision,allocatable :: bvvvv(:,:,:,:)
|
||||
double precision,allocatable :: hovvo(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: r1(:,:)
|
||||
double precision,allocatable :: r2(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: t1(:,:)
|
||||
double precision,allocatable :: t2(:,:,:,:)
|
||||
double precision,allocatable :: tau(:,:,:,:)
|
||||
|
||||
integer :: n_diis
|
||||
double precision :: rcond1
|
||||
double precision :: rcond2
|
||||
double precision,allocatable :: err1_diis(:,:)
|
||||
double precision,allocatable :: err2_diis(:,:)
|
||||
double precision,allocatable :: t1_diis(:,:)
|
||||
double precision,allocatable :: t2_diis(:,:)
|
||||
|
||||
! Hello world
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'********************************'
|
||||
write(*,*)'* Generalized CCSD Calculation *'
|
||||
write(*,*)'********************************'
|
||||
write(*,*)
|
||||
|
||||
! Antysymmetrize ERIs
|
||||
|
||||
allocate(dbERI(nBas,nBas,nBas,nBas))
|
||||
|
||||
call antisymmetrize_ERI(2,nBas,ERI,dbERI)
|
||||
|
||||
! Form energy denominator
|
||||
|
||||
allocate(eO(nO),eV(nV))
|
||||
allocate(delta_OV(nO,nV),delta_OOVV(nO,nO,nV,nV))
|
||||
|
||||
eO(:) = eHF(1:nO)
|
||||
eV(:) = eHF(nO+1:nBas)
|
||||
|
||||
call form_delta_OV(nC,nO,nV,nR,eO,eV,delta_OV)
|
||||
call form_delta_OOVV(nC,nO,nV,nR,eO,eV,delta_OOVV)
|
||||
|
||||
! Create integral batches
|
||||
|
||||
allocate(OOOO(nO,nO,nO,nO), &
|
||||
OOOV(nO,nO,nO,nV),OVOO(nO,nV,nO,nO),VOOO(nV,nO,nO,nO), &
|
||||
OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO), &
|
||||
OVVV(nO,nV,nV,nV),VOVV(nV,nO,nV,nV),VVVO(nV,nV,nV,nO), &
|
||||
VVVV(nV,nV,nV,nV))
|
||||
|
||||
OOOO(:,:,:,:) = dbERI( 1:nO , 1:nO , 1:nO , 1:nO )
|
||||
OOOV(:,:,:,:) = dbERI( 1:nO , 1:nO , 1:nO ,nO+1:nBas)
|
||||
OVOO(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas, 1:nO , 1:nO )
|
||||
VOOO(:,:,:,:) = dbERI(nO+1:nBas, 1:nO , 1:nO , 1:nO )
|
||||
OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas)
|
||||
OVVO(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas,nO+1:nBas, 1:nO )
|
||||
OVVV(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas,nO+1:nBas,nO+1:nBas)
|
||||
VOVV(:,:,:,:) = dbERI(nO+1:nBas, 1:nO ,nO+1:nBas,nO+1:nBas)
|
||||
VVVO(:,:,:,:) = dbERI(nO+1:nBas,nO+1:nBas,nO+1:nBas, 1:nO )
|
||||
VVVV(:,:,:,:) = dbERI(nO+1:nBas,nO+1:nBas,nO+1:nBas,nO+1:nBas)
|
||||
|
||||
deallocate(dbERI)
|
||||
|
||||
! MP2 guess amplitudes
|
||||
|
||||
allocate(t1(nO,nV),t2(nO,nO,nV,nV),tau(nO,nO,nV,nV))
|
||||
|
||||
t1(:,:) = 0d0
|
||||
t2(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:)
|
||||
call form_tau(nC,nO,nV,nR,t1,t2,tau)
|
||||
|
||||
EcMP2 = 0.5d0*dot_product(pack(OOVV,.true.),pack(tau,.true.))
|
||||
write(*,'(1X,A20,1X,F10.6)') 'Ec(MP2) = ',EcMP2
|
||||
|
||||
! Initialization
|
||||
|
||||
allocate(hvv(nV,nV),hoo(nO,nO),hvo(nV,nO), &
|
||||
gvv(nV,nV),goo(nO,nO), &
|
||||
aoooo(nO,nO,nO,nO),bvvvv(nV,nV,nV,nV),hovvo(nO,nV,nV,nO), &
|
||||
r1(nO,nV),r2(nO,nO,nV,nV))
|
||||
|
||||
! Memory allocation for DIIS
|
||||
|
||||
allocate(err1_diis(nO*nV ,max_diis),t1_diis(nO*nV ,max_diis), &
|
||||
err2_diis(nO*nO*nV*nV,max_diis),t2_diis(nO*nO*nV*nV,max_diis))
|
||||
|
||||
Conv = 1d0
|
||||
nSCF = 0
|
||||
ECC = EGHF
|
||||
|
||||
n_diis = 0
|
||||
t1_diis(:,:) = 0d0
|
||||
t2_diis(:,:) = 0d0
|
||||
err1_diis(:,:) = 0d0
|
||||
err2_diis(:,:) = 0d0
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Main SCF loop
|
||||
!------------------------------------------------------------------------
|
||||
write(*,*)
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,*)'| GCCSD calculation |'
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') &
|
||||
'|','#','|','E(GCCSD)','|','Ec(GCCSD)','|','Conv','|'
|
||||
write(*,*)'----------------------------------------------------'
|
||||
|
||||
do while(Conv > thresh .and. nSCF < maxSCF)
|
||||
|
||||
! Increment
|
||||
|
||||
nSCF = nSCF + 1
|
||||
|
||||
! Scuseria Eqs. (5), (6) and (7)
|
||||
|
||||
call form_h(nC,nO,nV,nR,eO,eV,OOVV,t1,tau,hvv,hoo,hvo)
|
||||
|
||||
! Scuseria Eqs. (9), (10), (11), (12) and (13)
|
||||
|
||||
call form_g(nC,nO,nV,nR,hvv,hoo,VOVV,OOOV,t1,gvv,goo)
|
||||
|
||||
call form_abh(nC,nO,nV,nR,OOOO,OVOO,OOVV,VVVV,VOVV,OVVO,OVVV,t1,tau,aoooo,bvvvv,hovvo)
|
||||
|
||||
! Compute residuals
|
||||
|
||||
call form_r1(nC,nO,nV,nR,OVVO,OVVV,OOOV,hvv,hoo,hvo,t1,t2,tau,r1)
|
||||
|
||||
call form_r2(nC,nO,nV,nR,OOVV,OVOO,OVVV,OVVO,gvv,goo,aoooo,bvvvv,hovvo,t1,t2,tau,r2)
|
||||
|
||||
! Check convergence
|
||||
|
||||
Conv = max(maxval(abs(r1(nC+1:nO,1:nV-nR))),maxval(abs(r2(nC+1:nO,nC+1:nO,1:nV-nR,1:nV-nR))))
|
||||
|
||||
! Update
|
||||
|
||||
t1(:,:) = t1(:,:) - r1(:,:) /delta_OV (:,:)
|
||||
t2(:,:,:,:) = t2(:,:,:,:) - r2(:,:,:,:)/delta_OOVV(:,:,:,:)
|
||||
|
||||
call form_tau(nC,nO,nV,nR,t1,t2,tau)
|
||||
|
||||
! Compute correlation energy
|
||||
|
||||
call CCSD_correlation_energy(nC,nO,nV,nR,OOVV,tau,EcCC)
|
||||
|
||||
! Dump results
|
||||
|
||||
ECC = EGHF + EcCC
|
||||
|
||||
! DIIS extrapolation
|
||||
|
||||
n_diis = min(n_diis+1,max_diis)
|
||||
! call DIIS_extrapolation(rcond1,nO*nV ,nO*nV ,n_diis,err1_diis,t1_diis,-r1/delta_OV ,t1)
|
||||
! call DIIS_extrapolation(rcond2,nO*nO*nV*nV,nO*nO*nV*nV,n_diis,err2_diis,t2_diis,-r2/delta_OOVV,t2)
|
||||
|
||||
! Reset DIIS if required
|
||||
|
||||
! if(min(abs(rcond1),abs(rcond2)) < 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,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|'
|
||||
|
||||
end do
|
||||
write(*,*)'----------------------------------------------------'
|
||||
!------------------------------------------------------------------------
|
||||
! End of SCF loop
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
! Did it actually converge?
|
||||
|
||||
if(nSCF == maxSCF) then
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)' Convergence failed '
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)
|
||||
|
||||
stop
|
||||
|
||||
end if
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,*)' CCSD energy '
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,'(1X,A20,1X,F15.10)')' E(CCSD) = ',ENuc+ECC
|
||||
write(*,'(1X,A20,1X,F10.6)')' Ec(CCSD) = ',EcCC
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
! Deallocate memory
|
||||
|
||||
deallocate(hvv,hoo,hvo, &
|
||||
delta_OV,delta_OOVV, &
|
||||
gvv,goo, &
|
||||
aoooo,bvvvv,hovvo, &
|
||||
tau, &
|
||||
r1,r2)
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! (T) correction
|
||||
!------------------------------------------------------------------------
|
||||
if(doCCSDT) then
|
||||
|
||||
call wall_time(start_CCSDT)
|
||||
call CCSDT(nC,nO,nV,nR,eO,eV,OOVV,VVVO,VOOO,t1,t2,EcCCT)
|
||||
call wall_time(end_CCSDT)
|
||||
|
||||
t_CCSDT = end_CCSDT - start_CCSDT
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for (T) = ',t_CCSDT,' seconds'
|
||||
write(*,*)
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,*)' GCCSD(T) energy '
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,'(1X,A20,1X,F15.10)')' E(GCCSD(T)) = ',ENuc + ECC + EcCCT
|
||||
write(*,'(1X,A20,1X,F10.6)')' Ec(GCCSD(T)) = ',EcCC + EcCCT
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
! Testing zone
|
||||
|
||||
if(dotest) then
|
||||
|
||||
call dump_test_value('R','GCCSD correlation energy',EcCC)
|
||||
|
||||
end if
|
||||
|
||||
end subroutine
|
||||
@ -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,ERI,ENuc,EHF,epsHF)
|
||||
maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
! Coupled-cluster module
|
||||
|
||||
@ -30,8 +30,8 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d
|
||||
integer,intent(in) :: nV(nspin)
|
||||
integer,intent(in) :: nR(nspin)
|
||||
double precision,intent(in) :: ENuc
|
||||
double precision,intent(in) :: EHF
|
||||
double precision,intent(in) :: epsHF(nBas)
|
||||
double precision,intent(in) :: ERHF
|
||||
double precision,intent(in) :: eHF(nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
@ -45,7 +45,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,EHF,epsHF)
|
||||
call CCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
@ -62,7 +62,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,EHF,epsHF)
|
||||
ERI,ENuc,ERHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
@ -80,7 +80,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,EHF,epsHF)
|
||||
call CCSD(dotest,maxSCF,thresh,max_diis,doCCSDT,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
@ -96,7 +96,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,EHF,epsHF)
|
||||
call drCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
@ -112,7 +112,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,EHF,epsHF,epsHF)
|
||||
call rCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
@ -128,7 +128,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,EHF,epsHF)
|
||||
call crCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
@ -144,7 +144,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,EHF,epsHF)
|
||||
call lCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
@ -160,8 +160,7 @@ 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,ERI,ENuc,EHF,epsHF)
|
||||
|
||||
call pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
|
||||
189
src/CC/crGCCD.f90
Normal file
189
src/CC/crGCCD.f90
Normal file
@ -0,0 +1,189 @@
|
||||
subroutine crGCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
! Generalized crossed-ring CCD module
|
||||
|
||||
implicit none
|
||||
|
||||
! Input variables
|
||||
|
||||
logical,intent(in) :: dotest
|
||||
|
||||
integer,intent(in) :: maxSCF
|
||||
integer,intent(in) :: max_diis
|
||||
double precision,intent(in) :: thresh
|
||||
|
||||
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,ERHF
|
||||
double precision,intent(in) :: eHF(nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: nSCF
|
||||
double precision :: Conv
|
||||
double precision :: EcMP2
|
||||
double precision :: ECC,EcCC
|
||||
double precision,allocatable :: dbERI(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: eO(:)
|
||||
double precision,allocatable :: eV(:)
|
||||
double precision,allocatable :: delta_OOVV(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: OOVV(:,:,:,:)
|
||||
double precision,allocatable :: OVOV(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: r2(:,:,:,:)
|
||||
double precision,allocatable :: t2(:,:,:,:)
|
||||
|
||||
integer :: n_diis
|
||||
double precision :: rcond
|
||||
double precision,allocatable :: error_diis(:,:)
|
||||
double precision,allocatable :: t_diis(:,:)
|
||||
|
||||
! Hello world
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'*********************************'
|
||||
write(*,*)'* Generalized crCCD Calculation *'
|
||||
write(*,*)'*********************************'
|
||||
write(*,*)
|
||||
|
||||
! Antysymmetrize ERIs
|
||||
|
||||
allocate(dbERI(nBas,nBas,nBas,nBas))
|
||||
|
||||
call antisymmetrize_ERI(2,nBas,ERI,dbERI)
|
||||
|
||||
! Form energy denominator
|
||||
|
||||
allocate(eO(nO),eV(nV))
|
||||
allocate(delta_OOVV(nO,nO,nV,nV))
|
||||
|
||||
eO(:) = eHF(1:nO)
|
||||
eV(:) = eHF(nO+1:nBas)
|
||||
|
||||
call form_delta_OOVV(nC,nO,nV,nR,eO,eV,delta_OOVV)
|
||||
|
||||
! Create integral batches
|
||||
|
||||
allocate(OOVV(nO,nO,nV,nV),OVOV(nO,nV,nO,nV))
|
||||
|
||||
OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas)
|
||||
OVOV(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas, 1:nO ,nO+1:nBas)
|
||||
|
||||
deallocate(dbERI)
|
||||
|
||||
! MP2 guess amplitudes
|
||||
|
||||
allocate(t2(nO,nO,nV,nV))
|
||||
|
||||
t2(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:)
|
||||
|
||||
call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcMP2)
|
||||
|
||||
! Memory allocation for DIIS
|
||||
|
||||
allocate(error_diis(nO*nO*nV*nV,max_diis),t_diis(nO*nO*nV*nV,max_diis))
|
||||
|
||||
! Initialization
|
||||
|
||||
allocate(r2(nO,nO,nV,nV))
|
||||
|
||||
Conv = 1d0
|
||||
nSCF = 0
|
||||
|
||||
n_diis = 0
|
||||
t_diis(:,:) = 0d0
|
||||
error_diis(:,:) = 0d0
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Main SCF loop
|
||||
!------------------------------------------------------------------------
|
||||
write(*,*)
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,*)'| crossed-ring CCD calculation |'
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') &
|
||||
'|','#','|','E(CCD)','|','Ec(CCD)','|','Conv','|'
|
||||
write(*,*)'----------------------------------------------------'
|
||||
|
||||
do while(Conv > thresh .and. nSCF < maxSCF)
|
||||
|
||||
! Increment
|
||||
|
||||
nSCF = nSCF + 1
|
||||
|
||||
! Compute residual
|
||||
|
||||
call form_crossed_ring_r(nC,nO,nV,nR,OVOV,OOVV,t2,r2)
|
||||
|
||||
r2(:,:,:,:) = OOVV(:,:,:,:) - delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:)
|
||||
|
||||
! Check convergence
|
||||
|
||||
Conv = maxval(abs(r2(nC+1:nO,nC+1:nO,1:nV-nR,1:nV-nR)))
|
||||
|
||||
! Update amplitudes
|
||||
|
||||
t2(:,:,:,:) = t2(:,:,:,:) + r2(:,:,:,:)/delta_OOVV(:,:,:,:)
|
||||
|
||||
! Compute correlation energy
|
||||
|
||||
call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcCC)
|
||||
|
||||
! Dump results
|
||||
|
||||
ECC = ERHF + EcCC
|
||||
|
||||
! DIIS extrapolation
|
||||
|
||||
n_diis = min(n_diis+1,max_diis)
|
||||
call DIIS_extrapolation(rcond,nO*nO*nV*nV,nO*nO*nV*nV,n_diis,error_diis,t_diis,-r2/delta_OOVV,t2)
|
||||
|
||||
! Reset DIIS if required
|
||||
|
||||
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,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|'
|
||||
|
||||
enddo
|
||||
write(*,*)'----------------------------------------------------'
|
||||
!------------------------------------------------------------------------
|
||||
! End of SCF loop
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
! Did it actually converge?
|
||||
|
||||
if(nSCF == maxSCF) then
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)' Convergence failed '
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)
|
||||
|
||||
stop
|
||||
|
||||
endif
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,*)' crossed-ring CCD energy '
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,'(1X,A30,1X,F15.10)')' E(crCCD) = ',ECC
|
||||
write(*,'(1X,A30,1X,F15.10)')' Ec(crCCD) = ',EcCC
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
if(dotest) then
|
||||
|
||||
call dump_test_value('R','crCCD correlation energy',EcCC)
|
||||
|
||||
end if
|
||||
|
||||
end subroutine
|
||||
181
src/CC/drGCCD.f90
Normal file
181
src/CC/drGCCD.f90
Normal file
@ -0,0 +1,181 @@
|
||||
subroutine drGCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
! Generalized Direct ring CCD module
|
||||
|
||||
implicit none
|
||||
|
||||
! Input variables
|
||||
|
||||
logical,intent(in) :: dotest
|
||||
|
||||
integer,intent(in) :: maxSCF
|
||||
integer,intent(in) :: max_diis
|
||||
double precision,intent(in) :: thresh
|
||||
|
||||
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,ERHF
|
||||
double precision,intent(in) :: eHF(nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: nSCF
|
||||
double precision :: Conv
|
||||
double precision :: EcMP2
|
||||
double precision :: ECC,EcCC
|
||||
|
||||
double precision,allocatable :: eO(:)
|
||||
double precision,allocatable :: eV(:)
|
||||
double precision,allocatable :: delta_OOVV(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: OOVV(:,:,:,:)
|
||||
double precision,allocatable :: OVVO(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: r2(:,:,:,:)
|
||||
double precision,allocatable :: t2(:,:,:,:)
|
||||
|
||||
|
||||
integer :: n_diis
|
||||
double precision :: rcond
|
||||
double precision,allocatable :: error_diis(:,:)
|
||||
double precision,allocatable :: t_diis(:,:)
|
||||
|
||||
! Hello world
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'*********************************'
|
||||
write(*,*)'* Generalized drCCD Calculation *'
|
||||
write(*,*)'*********************************'
|
||||
write(*,*)
|
||||
! Form energy denominator
|
||||
|
||||
allocate(eO(nO),eV(nV))
|
||||
allocate(delta_OOVV(nO,nO,nV,nV))
|
||||
|
||||
eO(:) = eHF(1:nO)
|
||||
eV(:) = eHF(nO+1:nBas)
|
||||
|
||||
call form_delta_OOVV(nC,nO,nV,nR,eO,eV,delta_OOVV)
|
||||
|
||||
! Create integral batches
|
||||
|
||||
allocate(OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO))
|
||||
|
||||
OOVV(:,:,:,:) = ERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas)
|
||||
OVVO(:,:,:,:) = ERI( 1:nO ,nO+1:nBas,nO+1:nBas, 1:nO )
|
||||
|
||||
! MP2 guess amplitudes
|
||||
|
||||
allocate(t2(nO,nO,nV,nV))
|
||||
|
||||
t2(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:)
|
||||
|
||||
call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcMP2)
|
||||
|
||||
! Memory allocation for DIIS
|
||||
|
||||
allocate(error_diis(nO*nO*nV*nV,max_diis),t_diis(nO*nO*nV*nV,max_diis))
|
||||
|
||||
! Initialization
|
||||
|
||||
allocate(r2(nO,nO,nV,nV))
|
||||
|
||||
Conv = 1d0
|
||||
nSCF = 0
|
||||
|
||||
n_diis = 0
|
||||
t_diis(:,:) = 0d0
|
||||
error_diis(:,:) = 0d0
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Main SCF loop
|
||||
!------------------------------------------------------------------------
|
||||
write(*,*)
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,*)'| direct ring GCCD calculation |'
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') &
|
||||
'|','#','|','E(GCCD)','|','Ec(GCCD)','|','Conv','|'
|
||||
write(*,*)'----------------------------------------------------'
|
||||
|
||||
do while(Conv > thresh .and. nSCF < maxSCF)
|
||||
|
||||
! Increment
|
||||
|
||||
nSCF = nSCF + 1
|
||||
|
||||
! Compute residual
|
||||
|
||||
call form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2)
|
||||
|
||||
r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:)
|
||||
|
||||
! Check convergence
|
||||
|
||||
Conv = maxval(abs(r2(nC+1:nO,nC+1:nO,1:nV-nR,1:nV-nR)))
|
||||
|
||||
! Update amplitudes
|
||||
|
||||
t2(:,:,:,:) = t2(:,:,:,:) - r2(:,:,:,:)/delta_OOVV(:,:,:,:)
|
||||
|
||||
! Compute correlation energy
|
||||
|
||||
call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcCC)
|
||||
EcCC = 2d0*EcCC
|
||||
|
||||
! Dump results
|
||||
|
||||
ECC = ERHF + EcCC
|
||||
|
||||
! DIIS extrapolation
|
||||
|
||||
n_diis = min(n_diis+1,max_diis)
|
||||
call DIIS_extrapolation(rcond,nO*nO*nV*nV,nO*nO*nV*nV,n_diis,error_diis,t_diis,-r2/delta_OOVV,t2)
|
||||
|
||||
! Reset DIIS if required
|
||||
|
||||
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,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|'
|
||||
|
||||
enddo
|
||||
write(*,*)'----------------------------------------------------'
|
||||
!------------------------------------------------------------------------
|
||||
! End of SCF loop
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
! Did it actually converge?
|
||||
|
||||
if(nSCF == maxSCF) then
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)' Convergence failed '
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)
|
||||
|
||||
stop
|
||||
|
||||
endif
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,*)' direct ring GCCD energy '
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,'(1X,A30,1X,F15.10)')' E(drGCCD) = ',ECC
|
||||
write(*,'(1X,A30,1X,F15.10)')' Ec(drGCCD) = ',EcCC
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
if(dotest) then
|
||||
|
||||
call dump_test_value('G','drCCD correlation energy',EcCC)
|
||||
|
||||
end if
|
||||
|
||||
end subroutine
|
||||
202
src/CC/lGCCD.f90
Normal file
202
src/CC/lGCCD.f90
Normal file
@ -0,0 +1,202 @@
|
||||
subroutine lGCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
! Generalized Ladder CCD module
|
||||
|
||||
implicit none
|
||||
|
||||
! Input variables
|
||||
|
||||
logical,intent(in) :: dotest
|
||||
|
||||
integer,intent(in) :: maxSCF
|
||||
integer,intent(in) :: max_diis
|
||||
double precision,intent(in) :: thresh
|
||||
|
||||
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,ERHF
|
||||
double precision,intent(in) :: eHF(nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: nSCF
|
||||
double precision :: Conv
|
||||
double precision :: EcMP2
|
||||
double precision :: ECC,EcCC
|
||||
double precision,allocatable :: dbERI(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: eO(:)
|
||||
double precision,allocatable :: eV(:)
|
||||
double precision,allocatable :: delta_OOVV(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: OOOO(:,:,:,:)
|
||||
double precision,allocatable :: OOVV(:,:,:,:)
|
||||
double precision,allocatable :: OVOV(:,:,:,:)
|
||||
double precision,allocatable :: VVVV(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: X1(:,:,:,:)
|
||||
double precision,allocatable :: X2(:,:)
|
||||
double precision,allocatable :: X3(:,:)
|
||||
double precision,allocatable :: X4(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: u(:,:,:,:)
|
||||
double precision,allocatable :: v(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: r2(:,:,:,:)
|
||||
double precision,allocatable :: t2(:,:,:,:)
|
||||
|
||||
integer :: n_diis
|
||||
double precision :: rcond
|
||||
double precision,allocatable :: error_diis(:,:)
|
||||
double precision,allocatable :: t_diis(:,:)
|
||||
|
||||
! Hello world
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'********************************'
|
||||
write(*,*)'* Generalized lCCD Calculation *'
|
||||
write(*,*)'********************************'
|
||||
write(*,*)
|
||||
|
||||
! Antysymmetrize ERIs
|
||||
|
||||
allocate(dbERI(nBas,nBas,nBas,nBas))
|
||||
|
||||
call antisymmetrize_ERI(2,nBas,ERI,dbERI)
|
||||
|
||||
! Form energy denominator
|
||||
|
||||
allocate(eO(nO),eV(nV))
|
||||
allocate(delta_OOVV(nO,nO,nV,nV))
|
||||
|
||||
eO(:) = eHF(1:nO)
|
||||
eV(:) = eHF(nO+1:nBas)
|
||||
|
||||
call form_delta_OOVV(nC,nO,nV,nR,eO,eV,delta_OOVV)
|
||||
|
||||
! Create integral batches
|
||||
|
||||
allocate(OOOO(nO,nO,nO,nO),OOVV(nO,nO,nV,nV),OVOV(nO,nV,nO,nV),VVVV(nV,nV,nV,nV))
|
||||
|
||||
OOOO(:,:,:,:) = dbERI( 1:nO , 1:nO , 1:nO , 1:nO )
|
||||
OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas)
|
||||
OVOV(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas, 1:nO ,nO+1:nBas)
|
||||
VVVV(:,:,:,:) = dbERI(nO+1:nBas,nO+1:nBas,nO+1:nBas,nO+1:nBas)
|
||||
|
||||
deallocate(dbERI)
|
||||
|
||||
! MP2 guess amplitudes
|
||||
|
||||
allocate(t2(nO,nO,nV,nV))
|
||||
|
||||
t2(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:)
|
||||
|
||||
call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcMP2)
|
||||
|
||||
! Memory allocation for DIIS
|
||||
|
||||
allocate(error_diis(nO*nO*nV*nV,max_diis),t_diis(nO*nO*nV*nV,max_diis))
|
||||
|
||||
! Initialization
|
||||
|
||||
allocate(r2(nO,nO,nV,nV),u(nO,nO,nV,nV),v(nO,nO,nV,nV))
|
||||
allocate(X1(nO,nO,nO,nO),X2(nV,nV),X3(nO,nO),X4(nO,nO,nV,nV))
|
||||
|
||||
Conv = 1d0
|
||||
nSCF = 0
|
||||
|
||||
n_diis = 0
|
||||
t_diis(:,:) = 0d0
|
||||
error_diis(:,:) = 0d0
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Main SCF loop
|
||||
!------------------------------------------------------------------------
|
||||
write(*,*)
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,*)'| ladder CCD calculation |'
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') &
|
||||
'|','#','|','E(CCD)','|','Ec(CCD)','|','Conv','|'
|
||||
write(*,*)'----------------------------------------------------'
|
||||
|
||||
do while(Conv > thresh .and. nSCF < maxSCF)
|
||||
|
||||
! Increment
|
||||
|
||||
nSCF = nSCF + 1
|
||||
|
||||
! Compute residual
|
||||
|
||||
call form_ladder_r(nC,nO,nV,nR,OOOO,OOVV,VVVV,t2,r2)
|
||||
|
||||
r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:)
|
||||
|
||||
! Check convergence
|
||||
|
||||
Conv = maxval(abs(r2(nC+1:nO,nC+1:nO,1:nV-nR,1:nV-nR)))
|
||||
|
||||
! Update amplitudes
|
||||
|
||||
t2(:,:,:,:) = t2(:,:,:,:) - r2(:,:,:,:)/delta_OOVV(:,:,:,:)
|
||||
|
||||
! Compute correlation energy
|
||||
|
||||
call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcCC)
|
||||
|
||||
! Dump results
|
||||
|
||||
ECC = ERHF + EcCC
|
||||
|
||||
! DIIS extrapolation
|
||||
|
||||
n_diis = min(n_diis+1,max_diis)
|
||||
call DIIS_extrapolation(rcond,nO*nO*nV*nV,nO*nO*nV*nV,n_diis,error_diis,t_diis,-r2/delta_OOVV,t2)
|
||||
|
||||
! Reset DIIS if required
|
||||
|
||||
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,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|'
|
||||
|
||||
enddo
|
||||
write(*,*)'----------------------------------------------------'
|
||||
!------------------------------------------------------------------------
|
||||
! End of SCF loop
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
! Did it actually converge?
|
||||
|
||||
if(nSCF == maxSCF) then
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)' Convergence failed '
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)
|
||||
|
||||
stop
|
||||
|
||||
endif
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,*)' ladder CCD energy '
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,'(1X,A30,1X,F15.10)')' E(lCCD) = ',ECC
|
||||
write(*,'(1X,A30,1X,F15.10)')' Ec(lCCD) = ',EcCC
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
if(dotest) then
|
||||
|
||||
call dump_test_value('R','lCCD correlation energy',EcCC)
|
||||
|
||||
end if
|
||||
|
||||
end subroutine
|
||||
@ -1,4 +1,4 @@
|
||||
subroutine rCCD(dotest,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF,eGW)
|
||||
subroutine rCCD(dotest,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
! Ring CCD module
|
||||
|
||||
@ -20,7 +20,6 @@ subroutine rCCD(dotest,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENu
|
||||
double precision,intent(in) :: ENuc
|
||||
double precision,intent(in) :: ERHF
|
||||
double precision,intent(in) :: eHF(nBasin)
|
||||
double precision,intent(in) :: eGW(nBasin)
|
||||
double precision,intent(in) :: ERI(nBasin,nBasin,nBasin,nBasin)
|
||||
|
||||
! Local variables
|
||||
@ -35,7 +34,6 @@ subroutine rCCD(dotest,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENu
|
||||
double precision :: EcMP2
|
||||
double precision :: ECC,EcCC
|
||||
double precision,allocatable :: seHF(:)
|
||||
double precision,allocatable :: seGW(:)
|
||||
double precision,allocatable :: sERI(:,:,:,:)
|
||||
double precision,allocatable :: dbERI(:,:,:,:)
|
||||
|
||||
@ -72,10 +70,9 @@ subroutine rCCD(dotest,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENu
|
||||
nV = 2*nVin
|
||||
nR = 2*nRin
|
||||
|
||||
allocate(seHF(nBas),seGW(nBas),sERI(nBas,nBas,nBas,nBas))
|
||||
allocate(seHF(nBas),sERI(nBas,nBas,nBas,nBas))
|
||||
|
||||
call spatial_to_spin_MO_energy(nBasin,eHF,nBas,seHF)
|
||||
call spatial_to_spin_MO_energy(nBasin,eGW,nBas,seGW)
|
||||
call spatial_to_spin_ERI(nBasin,ERI,nBas,sERI)
|
||||
|
||||
! Antysymmetrize ERIs
|
||||
@ -91,12 +88,12 @@ subroutine rCCD(dotest,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENu
|
||||
allocate(eO(nO),eV(nV))
|
||||
allocate(delta_OOVV(nO,nO,nV,nV))
|
||||
|
||||
eO(:) = seGW(1:nO)
|
||||
eV(:) = seGW(nO+1:nBas)
|
||||
eO(:) = seHF(1:nO)
|
||||
eV(:) = seHF(nO+1:nBas)
|
||||
|
||||
call form_delta_OOVV(nC,nO,nV,nR,eO,eV,delta_OOVV)
|
||||
|
||||
! deallocate(seHF,seGW)
|
||||
deallocate(seHF)
|
||||
|
||||
! Create integral batches
|
||||
|
||||
|
||||
208
src/CC/rGCCD.f90
Normal file
208
src/CC/rGCCD.f90
Normal file
@ -0,0 +1,208 @@
|
||||
subroutine rGCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
! Generalized ring CCD module
|
||||
|
||||
implicit none
|
||||
|
||||
! Input variables
|
||||
|
||||
logical,intent(in) :: dotest
|
||||
|
||||
integer,intent(in) :: maxSCF
|
||||
integer,intent(in) :: max_diis
|
||||
double precision,intent(in) :: thresh
|
||||
|
||||
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) :: ERI(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: nSCF
|
||||
double precision :: Conv
|
||||
double precision :: EcMP2
|
||||
double precision :: ECC,EcCC
|
||||
double precision,allocatable :: dbERI(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: eO(:)
|
||||
double precision,allocatable :: eV(:)
|
||||
double precision,allocatable :: delta_OOVV(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: OOVV(:,:,:,:)
|
||||
double precision,allocatable :: OVVO(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: r(:,:,:,:)
|
||||
double precision,allocatable :: t(:,:,:,:)
|
||||
|
||||
integer :: n_diis
|
||||
double precision :: rcond
|
||||
double precision,allocatable :: error_diis(:,:)
|
||||
double precision,allocatable :: t_diis(:,:)
|
||||
|
||||
logical :: do_EE_EOM_CC_1h1p = .true.
|
||||
|
||||
! Hello world
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'********************************'
|
||||
write(*,*)'* Generalized rCCD Calculation *'
|
||||
write(*,*)'********************************'
|
||||
write(*,*)
|
||||
|
||||
! Antysymmetrize ERIs
|
||||
|
||||
allocate(dbERI(nBas,nBas,nBas,nBas))
|
||||
|
||||
call antisymmetrize_ERI(2,nBas,ERI,dbERI)
|
||||
|
||||
! Form energy denominator
|
||||
|
||||
allocate(eO(nO),eV(nV))
|
||||
allocate(delta_OOVV(nO,nO,nV,nV))
|
||||
|
||||
eO(:) = eHF(1:nO)
|
||||
eV(:) = eHF(nO+1:nBas)
|
||||
|
||||
call form_delta_OOVV(nC,nO,nV,nR,eO,eV,delta_OOVV)
|
||||
|
||||
! Create integral batches
|
||||
|
||||
allocate(OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO))
|
||||
|
||||
OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas)
|
||||
OVVO(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas,nO+1:nBas, 1:nO )
|
||||
|
||||
deallocate(dbERI)
|
||||
|
||||
! MP2 guess amplitudes
|
||||
|
||||
allocate(t(nO,nO,nV,nV))
|
||||
|
||||
t(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:)
|
||||
|
||||
call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t,EcMP2)
|
||||
|
||||
! Memory allocation for DIIS
|
||||
|
||||
allocate(error_diis(nO*nO*nV*nV,max_diis),t_diis(nO*nO*nV*nV,max_diis))
|
||||
|
||||
! Initialization
|
||||
|
||||
allocate(r(nO,nO,nV,nV))
|
||||
|
||||
Conv = 1d0
|
||||
nSCF = 0
|
||||
|
||||
n_diis = 0
|
||||
t_diis(:,:) = 0d0
|
||||
error_diis(:,:) = 0d0
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Main SCF loop
|
||||
!------------------------------------------------------------------------
|
||||
write(*,*)
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,*)'| ring CCD calculation |'
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') &
|
||||
'|','#','|','E(CCD)','|','Ec(CCD)','|','Conv','|'
|
||||
write(*,*)'----------------------------------------------------'
|
||||
|
||||
do while(Conv > thresh .and. nSCF < maxSCF)
|
||||
|
||||
! Increment
|
||||
|
||||
nSCF = nSCF + 1
|
||||
|
||||
! Compute residual
|
||||
|
||||
call form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t,r)
|
||||
|
||||
r(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t(:,:,:,:) + r(:,:,:,:)
|
||||
|
||||
! Check convergence
|
||||
|
||||
Conv = maxval(abs(r(nC+1:nO,nC+1:nO,1:nV-nR,1:nV-nR)))
|
||||
|
||||
! Update amplitudes
|
||||
|
||||
t(:,:,:,:) = t(:,:,:,:) - r(:,:,:,:)/delta_OOVV(:,:,:,:)
|
||||
|
||||
! Compute correlation energy
|
||||
|
||||
call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t,EcCC)
|
||||
|
||||
! Dump results
|
||||
|
||||
ECC = ERHF + EcCC
|
||||
|
||||
! DIIS extrapolation
|
||||
|
||||
n_diis = min(n_diis+1,max_diis)
|
||||
call DIIS_extrapolation(rcond,nO*nO*nV*nV,nO*nO*nV*nV,n_diis,error_diis,t_diis,-r/delta_OOVV,t)
|
||||
|
||||
! Reset DIIS if required
|
||||
|
||||
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,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|'
|
||||
|
||||
enddo
|
||||
write(*,*)'----------------------------------------------------'
|
||||
!------------------------------------------------------------------------
|
||||
! End of SCF loop
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
! Did it actually converge?
|
||||
|
||||
if(nSCF == maxSCF) then
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)' Convergence failed '
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)
|
||||
|
||||
stop
|
||||
|
||||
endif
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,*)' ring CCD energy '
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,'(1X,A30,1X,F15.10)')' E(rCCD) = ',ECC
|
||||
write(*,'(1X,A30,1X,F15.10)')' Ec(rCCD) = ',EcCC
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
! write(*,*)
|
||||
! write(*,*)'----------------------------------------------------'
|
||||
! write(*,*)' ring CCD amplitudes '
|
||||
! write(*,*)'----------------------------------------------------'
|
||||
! call matout(nO*nO,nV*nV,t)
|
||||
! write(*,*)
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! EOM section
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
! EE-EOM-CCD (1h1p)
|
||||
|
||||
if(do_EE_EOM_CC_1h1p) call EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t)
|
||||
|
||||
|
||||
if(dotest) then
|
||||
|
||||
call dump_test_value('R','rCCD correlation energy',EcCC)
|
||||
|
||||
end if
|
||||
|
||||
end subroutine
|
||||
@ -98,7 +98,7 @@ subroutine GHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
|
||||
!---------------------!
|
||||
|
||||
call wall_time(start_HF)
|
||||
call GHF(maxSCF,thresh,max_diis,guess,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
|
||||
call GHF(.false.,maxSCF,thresh,max_diis,guess,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
|
||||
nBas,nBas2,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,e,c,P)
|
||||
call wall_time(end_HF)
|
||||
|
||||
|
||||
@ -91,7 +91,7 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
|
||||
!---------------------!
|
||||
|
||||
call wall_time(start_HF)
|
||||
call RHF(maxSCF,thresh,max_diis,guess,level_shift,nNuc,ZNuc,rNuc,ENuc, &
|
||||
call RHF(.false.,maxSCF,thresh,max_diis,guess,level_shift,nNuc,ZNuc,rNuc,ENuc, &
|
||||
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,e,c,P)
|
||||
call wall_time(end_HF)
|
||||
|
||||
|
||||
@ -102,7 +102,7 @@ subroutine UHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
|
||||
!---------------------!
|
||||
|
||||
call wall_time(start_HF)
|
||||
call UHF(maxSCF,thresh,max_diis,guess,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
|
||||
call UHF(.false.,maxSCF,thresh,max_diis,guess,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
|
||||
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,e,c,P)
|
||||
call wall_time(end_HF)
|
||||
|
||||
|
||||
@ -1,7 +1,9 @@
|
||||
subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,doppRPA, &
|
||||
subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
|
||||
dodrCCD,dorCCD,docrCCD,dolCCD,dophRPA,dophRPAx,doppRPA, &
|
||||
doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, &
|
||||
nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO, &
|
||||
maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,reg_MP, &
|
||||
maxSCF_CC,max_diis_CC,thresh_CC, &
|
||||
TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, &
|
||||
maxSCF_GW,max_diis_GW,thresh_GW,TDA_W,lin_GW,reg_GW,eta_GW, &
|
||||
dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS)
|
||||
@ -16,6 +18,8 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dopp
|
||||
logical,intent(in) :: dosearch
|
||||
logical,intent(in) :: doMP2
|
||||
logical,intent(in) :: doMP3
|
||||
logical,intent(in) :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT
|
||||
logical,intent(in) :: dodrCCD,dorCCD,docrCCD,dolCCD
|
||||
logical,intent(in) :: dophRPA,dophRPAx,doppRPA
|
||||
logical,intent(in) :: doG0F2,doevGF2,doqsGF2
|
||||
logical,intent(in) :: doG0W0,doevGW,doqsGW
|
||||
@ -43,6 +47,9 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dopp
|
||||
|
||||
logical,intent(in) :: reg_MP
|
||||
|
||||
integer,intent(in) :: maxSCF_CC,max_diis_CC
|
||||
double precision,intent(in) :: thresh_CC
|
||||
|
||||
logical,intent(in) :: TDA
|
||||
|
||||
integer,intent(in) :: maxSCF_GF,max_diis_GF
|
||||
@ -60,18 +67,19 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dopp
|
||||
|
||||
! Local variables
|
||||
|
||||
logical :: doMP,doRPA,doGF,doGW
|
||||
logical :: doMP,doCC,doRPA,doGF,doGW
|
||||
|
||||
double precision :: start_HF ,end_HF ,t_HF
|
||||
double precision :: start_stab ,end_stab ,t_stab
|
||||
double precision :: start_AOtoMO ,end_AOtoMO ,t_AOtoMO
|
||||
double precision :: start_MP ,end_MP ,t_MP
|
||||
double precision :: start_CC ,end_CC ,t_CC
|
||||
double precision :: start_RPA ,end_RPA ,t_RPA
|
||||
double precision :: start_GF ,end_GF ,t_GF
|
||||
double precision :: start_GW ,end_GW ,t_GW
|
||||
|
||||
double precision,allocatable :: cHF(:,:),epsHF(:),PHF(:,:)
|
||||
double precision :: EHF
|
||||
double precision,allocatable :: cHF(:,:),eHF(:),PHF(:,:)
|
||||
double precision :: EGHF
|
||||
double precision,allocatable :: dipole_int_MO(:,:,:)
|
||||
double precision,allocatable :: ERI_MO(:,:,:,:)
|
||||
double precision,allocatable :: ERI_tmp(:,:,:,:)
|
||||
@ -93,7 +101,7 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dopp
|
||||
|
||||
nBas2 = 2*nBas
|
||||
|
||||
allocate(cHF(nBas2,nBas2),epsHF(nBas2),PHF(nBas2,nBas2), &
|
||||
allocate(cHF(nBas2,nBas2),eHF(nBas2),PHF(nBas2,nBas2), &
|
||||
dipole_int_MO(nBas2,nBas2,ncart),ERI_MO(nBas2,nBas2,nBas2,nBas2))
|
||||
|
||||
!---------------------!
|
||||
@ -104,7 +112,7 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dopp
|
||||
|
||||
call wall_time(start_HF)
|
||||
call GHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
|
||||
nBas,nBas2,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF)
|
||||
nBas,nBas2,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EGHF,eHF,cHF,PHF)
|
||||
call wall_time(end_HF)
|
||||
|
||||
t_HF = end_HF - start_HF
|
||||
@ -165,7 +173,7 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dopp
|
||||
if(dostab) then
|
||||
|
||||
call wall_time(start_stab)
|
||||
call GHF_stability(nBas2,nC,nO,nV,nR,nS,epsHF,ERI_MO)
|
||||
call GHF_stability(nBas2,nC,nO,nV,nR,nS,eHF,ERI_MO)
|
||||
call wall_time(end_stab)
|
||||
|
||||
t_stab = end_stab - start_stab
|
||||
@ -178,7 +186,7 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dopp
|
||||
|
||||
call wall_time(start_stab)
|
||||
call GHF_search(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
|
||||
nBas,nBas2,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF)
|
||||
nBas,nBas2,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EGHF,eHF,cHF,PHF)
|
||||
call wall_time(end_stab)
|
||||
|
||||
t_stab = end_stab - start_stab
|
||||
@ -196,7 +204,7 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dopp
|
||||
if(doMP) then
|
||||
|
||||
call wall_time(start_MP)
|
||||
call GMP(dotest,doMP2,doMP3,reg_MP,nBas2,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF)
|
||||
call GMP(dotest,doMP2,doMP3,reg_MP,nBas2,nC,nO,nV,nR,ERI_MO,ENuc,EGHF,eHF)
|
||||
call wall_time(end_MP)
|
||||
|
||||
t_MP = end_MP - start_MP
|
||||
@ -205,6 +213,25 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dopp
|
||||
|
||||
end if
|
||||
|
||||
!------------------------!
|
||||
! Coupled-cluster module !
|
||||
!------------------------!
|
||||
|
||||
doCC = doCCD .or. doCCSD .or. doCCSDT .or. dodrCCD .or. dorCCD .or. docrCCD .or. dolCCD
|
||||
|
||||
if(doCC) then
|
||||
|
||||
call wall_time(start_CC)
|
||||
call GCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, &
|
||||
maxSCF_CC,thresh_CC,max_diis_CC,nBas2,nC,nO,nV,nR,ERI_MO,ENuc,EGHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CC = ',t_CC,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
!-----------------------------------!
|
||||
! Random-phase approximation module !
|
||||
!-----------------------------------!
|
||||
@ -214,8 +241,8 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dopp
|
||||
if(doRPA) then
|
||||
|
||||
call wall_time(start_RPA)
|
||||
call GRPA(dotest,dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas2,nC,nO,nV,nR,nS,ENuc,EHF, &
|
||||
ERI_MO,dipole_int_MO,epsHF,cHF,S)
|
||||
call GRPA(dotest,dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas2,nC,nO,nV,nR,nS,ENuc,EGHF, &
|
||||
ERI_MO,dipole_int_MO,eHF,cHF,S)
|
||||
call wall_time(end_RPA)
|
||||
|
||||
t_RPA = end_RPA - start_RPA
|
||||
@ -234,7 +261,7 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dopp
|
||||
|
||||
call wall_time(start_GF)
|
||||
call GGF(dotest,doG0F2,doevGF2,doqsGF2,maxSCF_GF,thresh_GF,max_diis_GF,dophBSE,doppBSE,TDA,dBSE,dTDA,lin_GF,eta_GF,reg_GF, &
|
||||
nNuc,ZNuc,rNuc,ENuc,nBas2,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF)
|
||||
nNuc,ZNuc,rNuc,ENuc,nBas2,nC,nO,nV,nR,nS,EGHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
|
||||
call wall_time(end_GF)
|
||||
|
||||
t_GF = end_GF - start_GF
|
||||
@ -254,7 +281,7 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dopp
|
||||
call wall_time(start_GW)
|
||||
call GGW(dotest,doG0W0,doevGW,doqsGW,maxSCF_GW,thresh_GW,max_diis_GW,doACFDT,exchange_kernel,doXBS, &
|
||||
dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc, &
|
||||
nBas,nBas2,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF)
|
||||
nBas,nBas2,nC,nO,nV,nR,nS,EGHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
|
||||
call wall_time(end_GW)
|
||||
|
||||
t_GW = end_GW - start_GW
|
||||
|
||||
@ -226,11 +226,12 @@ program QuAcK
|
||||
!--------------------------!
|
||||
|
||||
if(doGQuAcK) &
|
||||
call GQuAcK(doGtest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,doppRPA, &
|
||||
call GQuAcK(doGtest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
|
||||
dodrCCD,dorCCD,docrCCD,dolCCD,dophRPA,dophRPAx,doppRPA, &
|
||||
doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, &
|
||||
nNuc,nBas,sum(nC),sum(nO),sum(nV),sum(nR),ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO, &
|
||||
maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,reg_MP, &
|
||||
TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, &
|
||||
maxSCF_CC,max_diis_CC,thresh_CC,TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, &
|
||||
maxSCF_GW,max_diis_GW,thresh_GW,TDA_W,lin_GW,reg_GW,eta_GW, &
|
||||
dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS)
|
||||
|
||||
|
||||
@ -90,8 +90,8 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
|
||||
double precision :: start_GW ,end_GW ,t_GW
|
||||
double precision :: start_GT ,end_GT ,t_GT
|
||||
|
||||
double precision,allocatable :: cHF(:,:),epsHF(:),PHF(:,:)
|
||||
double precision :: EHF
|
||||
double precision,allocatable :: cHF(:,:),eHF(:),PHF(:,:)
|
||||
double precision :: ERHF
|
||||
double precision,allocatable :: dipole_int_MO(:,:,:)
|
||||
double precision,allocatable :: ERI_MO(:,:,:,:)
|
||||
integer :: ixyz
|
||||
@ -107,7 +107,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
|
||||
! Memory allocation !
|
||||
!-------------------!
|
||||
|
||||
allocate(cHF(nBas,nBas),epsHF(nBas),PHF(nBas,nBas), &
|
||||
allocate(cHF(nBas,nBas),eHF(nBas),PHF(nBas,nBas), &
|
||||
dipole_int_MO(nBas,nBas,ncart),ERI_MO(nBas,nBas,nBas,nBas))
|
||||
|
||||
!---------------------!
|
||||
@ -118,7 +118,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
|
||||
|
||||
call wall_time(start_HF)
|
||||
call RHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
|
||||
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF)
|
||||
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF)
|
||||
call wall_time(end_HF)
|
||||
|
||||
t_HF = end_HF - start_HF
|
||||
@ -131,7 +131,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
|
||||
|
||||
call wall_time(start_HF)
|
||||
call ROHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
|
||||
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF)
|
||||
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF)
|
||||
call wall_time(end_HF)
|
||||
|
||||
t_HF = end_HF - start_HF
|
||||
@ -175,7 +175,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
|
||||
if(dostab) then
|
||||
|
||||
call wall_time(start_stab)
|
||||
call RHF_stability(nBas,nC,nO,nV,nR,nS,epsHF,ERI_MO)
|
||||
call RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI_MO)
|
||||
call wall_time(end_stab)
|
||||
|
||||
t_stab = end_stab - start_stab
|
||||
@ -188,7 +188,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
|
||||
|
||||
call wall_time(start_stab)
|
||||
call RHF_search(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
|
||||
nBas,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF)
|
||||
nBas,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF)
|
||||
call wall_time(end_stab)
|
||||
|
||||
t_stab = end_stab - start_stab
|
||||
@ -206,7 +206,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
|
||||
if(doMP) then
|
||||
|
||||
call wall_time(start_MP)
|
||||
call RMP(dotest,doMP2,doMP3,reg_MP,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF)
|
||||
call RMP(dotest,doMP2,doMP3,reg_MP,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
|
||||
call wall_time(end_MP)
|
||||
|
||||
t_MP = end_MP - start_MP
|
||||
@ -226,7 +226,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,ERI_MO,ENuc,EHF,epsHF)
|
||||
maxSCF_CC,thresh_CC,max_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
@ -245,7 +245,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
|
||||
|
||||
call wall_time(start_CI)
|
||||
call RCI(dotest,doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO, &
|
||||
epsHF,EHF,cHF,S)
|
||||
eHF,ERHF,cHF,S)
|
||||
call wall_time(end_CI)
|
||||
|
||||
t_CI = end_CI - start_CI
|
||||
@ -264,7 +264,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
|
||||
|
||||
call wall_time(start_RPA)
|
||||
call RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,singlet,triplet, &
|
||||
nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,dipole_int_MO,epsHF,cHF,S)
|
||||
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF,cHF,S)
|
||||
call wall_time(end_RPA)
|
||||
|
||||
t_RPA = end_RPA - start_RPA
|
||||
@ -284,8 +284,8 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
|
||||
call wall_time(start_GF)
|
||||
call RGF(dotest,doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,renorm_GF,maxSCF_GF,thresh_GF,max_diis_GF, &
|
||||
dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,lin_GF,eta_GF,reg_GF, &
|
||||
nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI_MO, &
|
||||
dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF)
|
||||
nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO, &
|
||||
dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
|
||||
call wall_time(end_GF)
|
||||
|
||||
t_GF = end_GF - start_GF
|
||||
@ -305,8 +305,8 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
|
||||
call wall_time(start_GW)
|
||||
call RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF_GW,thresh_GW,max_diis_GW,doACFDT, &
|
||||
exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, &
|
||||
lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, &
|
||||
ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF)
|
||||
lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc, &
|
||||
ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
|
||||
call wall_time(end_GW)
|
||||
|
||||
t_GW = end_GW - start_GW
|
||||
@ -326,8 +326,8 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
|
||||
call wall_time(start_GT)
|
||||
call RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF_GT,thresh_GT,max_diis_GT,doACFDT, &
|
||||
exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, &
|
||||
lin_GT,eta_GT,reg_GT,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, &
|
||||
ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF)
|
||||
lin_GT,eta_GT,reg_GT,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc, &
|
||||
ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
|
||||
call wall_time(end_GT)
|
||||
|
||||
t_GT = end_GT - start_GT
|
||||
|
||||
@ -88,8 +88,8 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
|
||||
double precision :: start_GW ,end_GW ,t_GW
|
||||
double precision :: start_GT ,end_GT ,t_GT
|
||||
|
||||
double precision,allocatable :: cHF(:,:,:),epsHF(:,:),PHF(:,:,:)
|
||||
double precision :: EHF
|
||||
double precision,allocatable :: cHF(:,:,:),eHF(:,:),PHF(:,:,:)
|
||||
double precision :: EUHF
|
||||
double precision,allocatable :: dipole_int_aa(:,:,:),dipole_int_bb(:,:,:)
|
||||
double precision,allocatable :: ERI_aaaa(:,:,:,:),ERI_aabb(:,:,:,:),ERI_bbbb(:,:,:,:)
|
||||
integer :: ixyz
|
||||
@ -105,7 +105,7 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
|
||||
! Memory allocation !
|
||||
!-------------------!
|
||||
|
||||
allocate(cHF(nBas,nBas,nspin),epsHF(nBas,nspin),PHF(nBas,nBas,nspin), &
|
||||
allocate(cHF(nBas,nBas,nspin),eHF(nBas,nspin),PHF(nBas,nBas,nspin), &
|
||||
dipole_int_aa(nBas,nBas,ncart),dipole_int_bb(nBas,nBas,ncart), &
|
||||
ERI_aaaa(nBas,nBas,nBas,nBas),ERI_aabb(nBas,nBas,nBas,nBas),ERI_bbbb(nBas,nBas,nBas,nBas))
|
||||
|
||||
@ -117,7 +117,7 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
|
||||
|
||||
call wall_time(start_HF)
|
||||
call UHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
|
||||
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF)
|
||||
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EUHF,eHF,cHF,PHF)
|
||||
call wall_time(end_HF)
|
||||
|
||||
t_HF = end_HF - start_HF
|
||||
@ -138,7 +138,7 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
|
||||
! write(*,*)
|
||||
! call eDFT(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc,nBas,nC, &
|
||||
! nO,nV,nR,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell, &
|
||||
! max_ang_mom,min_exponent,max_exponent,S,T,V,Hc,X,ERI_AO,dipole_int_AO,EHF,epsHF,cHF,PHF,Vxc)
|
||||
! max_ang_mom,min_exponent,max_exponent,S,T,V,Hc,X,ERI_AO,dipole_int_AO,EUHF,eHF,cHF,PHF,Vxc)
|
||||
! call wall_time(end_KS)
|
||||
|
||||
! t_KS = end_KS - start_KS
|
||||
@ -191,7 +191,7 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
|
||||
if(dostab) then
|
||||
|
||||
call wall_time(start_stab)
|
||||
call UHF_stability(nBas,nC,nO,nV,nR,nS,epsHF,ERI_aaaa,ERI_aabb,ERI_bbbb)
|
||||
call UHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI_aaaa,ERI_aabb,ERI_bbbb)
|
||||
call wall_time(end_stab)
|
||||
|
||||
t_stab = end_stab - start_stab
|
||||
@ -204,7 +204,7 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
|
||||
|
||||
call wall_time(start_stab)
|
||||
call UHF_search(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
|
||||
nBas,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF)
|
||||
nBas,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EUHF,eHF,cHF,PHF)
|
||||
|
||||
call wall_time(end_stab)
|
||||
|
||||
@ -223,7 +223,7 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
|
||||
if(doMP) then
|
||||
|
||||
call wall_time(start_MP)
|
||||
call UMP(dotest,doMP2,doMP3,reg_MP,nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EHF,epsHF)
|
||||
call UMP(dotest,doMP2,doMP3,reg_MP,nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EUHF,eHF)
|
||||
call wall_time(end_MP)
|
||||
|
||||
t_MP = end_MP - start_MP
|
||||
@ -260,7 +260,7 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
|
||||
|
||||
call wall_time(start_CI)
|
||||
call UCI(dotest,doCIS,doCIS_D,doCID,doCISD,doFCI,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS, &
|
||||
ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,EHF,cHF,S)
|
||||
ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,eHF,EUHF,cHF,S)
|
||||
call wall_time(end_CI)
|
||||
|
||||
t_CI = end_CI - start_CI
|
||||
@ -279,7 +279,7 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
|
||||
|
||||
call wall_time(start_RPA)
|
||||
call URPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip, &
|
||||
nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,cHF,S)
|
||||
nBas,nC,nO,nV,nR,nS,ENuc,EUHF,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,eHF,cHF,S)
|
||||
call wall_time(end_RPA)
|
||||
|
||||
t_RPA = end_RPA - start_RPA
|
||||
@ -299,8 +299,8 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
|
||||
call wall_time(start_GF)
|
||||
call UGF(dotest,doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,renorm_GF,maxSCF_GF,thresh_GF,max_diis_GF, &
|
||||
dophBSE,doppBSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,lin_GF,eta_GF,reg_GF, &
|
||||
nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb, &
|
||||
dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF)
|
||||
nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EUHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb, &
|
||||
dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF)
|
||||
call wall_time(end_GF)
|
||||
|
||||
t_GF = end_GF - start_GF
|
||||
@ -320,8 +320,8 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
|
||||
call wall_time(start_GW)
|
||||
call UGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF_GW,thresh_GW,max_diis_GW,doACFDT, &
|
||||
exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
|
||||
lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, &
|
||||
ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF)
|
||||
lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EUHF,S,X,T,V,Hc, &
|
||||
ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF)
|
||||
call wall_time(end_GW)
|
||||
|
||||
t_GW = end_GW - start_GW
|
||||
@ -341,8 +341,8 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
|
||||
call wall_time(start_GT)
|
||||
call UGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF_GT,thresh_GT,max_diis_GT,doACFDT, &
|
||||
exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
|
||||
lin_GT,eta_GT,reg_GT,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, &
|
||||
ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF)
|
||||
lin_GT,eta_GT,reg_GT,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EUHF,S,X,T,V,Hc, &
|
||||
ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF)
|
||||
call wall_time(end_GT)
|
||||
|
||||
t_GT = end_GT - start_GT
|
||||
|
||||
Loading…
Reference in New Issue
Block a user