mirror of
https://github.com/pfloos/quack
synced 2024-12-23 04:43:53 +01:00
Merge branch 'master' of https://github.com/pfloos/QuAcK
This commit is contained in:
Clotilde Marut
commit
f8b7cc125d
@ -1,15 +1,15 @@
|
|||||||
# RHF UHF KS MOM
|
# RHF UHF KS MOM
|
||||||
F F T F
|
T F F F
|
||||||
# MP2* MP3 MP2-F12
|
# MP2* MP3 MP2-F12
|
||||||
F F F
|
F F F
|
||||||
# CCD DCD CCSD CCSD(T)
|
# CCD DCD CCSD CCSD(T)
|
||||||
F F F F
|
F F F F
|
||||||
# drCCD rCCD lCCD pCCD
|
# drCCD rCCD lCCD pCCD
|
||||||
F F F F
|
T T F F
|
||||||
# CIS* CIS(D) CID CISD FCI
|
# CIS* CIS(D) CID CISD FCI
|
||||||
F F F F F
|
F F F F F
|
||||||
# RPA* RPAx* ppRPA
|
# RPA* RPAx* ppRPA
|
||||||
F F F
|
T T F
|
||||||
# G0F2* evGF2* qsGF2* G0F3 evGF3
|
# G0F2* evGF2* qsGF2* G0F3 evGF3
|
||||||
F F F F F
|
F F F F F
|
||||||
# G0W0* evGW* qsGW* ufG0W0 ufGW
|
# G0W0* evGW* qsGW* ufG0W0 ufGW
|
||||||
|
|||||||
@ -5,7 +5,7 @@
|
|||||||
# CC: maxSCF thresh DIIS n_diis
|
# CC: maxSCF thresh DIIS n_diis
|
||||||
64 0.0000000001 T 5
|
64 0.0000000001 T 5
|
||||||
# spin: TDA singlet triplet spin_conserved spin_flip
|
# spin: TDA singlet triplet spin_conserved spin_flip
|
||||||
F T F T T
|
F T T T T
|
||||||
# GF: maxSCF thresh DIIS n_diis lin eta renorm
|
# GF: maxSCF thresh DIIS n_diis lin eta renorm
|
||||||
256 0.00001 T 5 T 0.00367493 3
|
256 0.00001 T 5 T 0.00367493 3
|
||||||
# GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0
|
# GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0
|
||||||
|
|||||||
@ -1,4 +1,4 @@
|
|||||||
2
|
2
|
||||||
|
|
||||||
H 0. 0. 0.
|
H 0. 0. 0.
|
||||||
H 0. -0.740848 0.740848
|
H 0. 0. 2.000000
|
||||||
|
|||||||
@ -220,6 +220,15 @@ subroutine CCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,e
|
|||||||
|
|
||||||
endif
|
endif
|
||||||
|
|
||||||
|
write(*,*)
|
||||||
|
write(*,*)'----------------------------------------------------'
|
||||||
|
write(*,*)' CCD energy '
|
||||||
|
write(*,*)'----------------------------------------------------'
|
||||||
|
write(*,'(1X,A30,1X,F15.10)')' E(CCD) = ',ECCD
|
||||||
|
write(*,'(1X,A30,1X,F15.10)')' Ec(CCD) = ',EcCCD
|
||||||
|
write(*,*)'----------------------------------------------------'
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
! Moller-Plesset energies
|
! Moller-Plesset energies
|
||||||
|
|
||||||
write(*,*)
|
write(*,*)
|
||||||
|
|||||||
237
src/CC/crCCD.f90
Normal file
237
src/CC/crCCD.f90
Normal file
@ -0,0 +1,237 @@
|
|||||||
|
subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF)
|
||||||
|
|
||||||
|
! Crossed-ring CCD module
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
! Input variables
|
||||||
|
|
||||||
|
integer,intent(in) :: maxSCF
|
||||||
|
integer,intent(in) :: max_diis
|
||||||
|
double precision,intent(in) :: thresh
|
||||||
|
|
||||||
|
integer,intent(in) :: nBasin
|
||||||
|
integer,intent(in) :: nCin
|
||||||
|
integer,intent(in) :: nOin
|
||||||
|
integer,intent(in) :: nVin
|
||||||
|
integer,intent(in) :: nRin
|
||||||
|
double precision,intent(in) :: ENuc,ERHF
|
||||||
|
double precision,intent(in) :: eHF(nBasin)
|
||||||
|
double precision,intent(in) :: ERI(nBasin,nBasin,nBasin,nBasin)
|
||||||
|
|
||||||
|
! Local variables
|
||||||
|
|
||||||
|
integer :: nBas
|
||||||
|
integer :: nC
|
||||||
|
integer :: nO
|
||||||
|
integer :: nV
|
||||||
|
integer :: nR
|
||||||
|
integer :: nSCF
|
||||||
|
double precision :: Conv
|
||||||
|
double precision :: EcMP2
|
||||||
|
double precision :: ECCD,EcCCD
|
||||||
|
double precision,allocatable :: seHF(:)
|
||||||
|
double precision,allocatable :: sERI(:,:,:,:)
|
||||||
|
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 :: OVVO(:,:,:,:)
|
||||||
|
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 :: r2r(:,:,:,:)
|
||||||
|
double precision,allocatable :: r2l(:,:,:,:)
|
||||||
|
|
||||||
|
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(*,*)'| crossed-ring CCD calculation |'
|
||||||
|
write(*,*)'**************************************'
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
! Spatial to spin orbitals
|
||||||
|
|
||||||
|
nBas = 2*nBasin
|
||||||
|
nC = 2*nCin
|
||||||
|
nO = 2*nOin
|
||||||
|
nV = 2*nVin
|
||||||
|
nR = 2*nRin
|
||||||
|
|
||||||
|
allocate(seHF(nBas),sERI(nBas,nBas,nBas,nBas))
|
||||||
|
|
||||||
|
call spatial_to_spin_MO_energy(nBasin,eHF,nBas,seHF)
|
||||||
|
call spatial_to_spin_ERI(nBasin,ERI,nBas,sERI)
|
||||||
|
|
||||||
|
! Antysymmetrize ERIs
|
||||||
|
|
||||||
|
allocate(dbERI(nBas,nBas,nBas,nBas))
|
||||||
|
|
||||||
|
call antisymmetrize_ERI(2,nBas,sERI,dbERI)
|
||||||
|
|
||||||
|
deallocate(sERI)
|
||||||
|
|
||||||
|
! Form energy denominator
|
||||||
|
|
||||||
|
allocate(eO(nO),eV(nV))
|
||||||
|
allocate(delta_OOVV(nO,nO,nV,nV))
|
||||||
|
|
||||||
|
eO(:) = seHF(1:nO)
|
||||||
|
eV(:) = seHF(nO+1:nBas)
|
||||||
|
|
||||||
|
call form_delta_OOVV(nC,nO,nV,nR,eO,eV,delta_OOVV)
|
||||||
|
|
||||||
|
deallocate(seHF)
|
||||||
|
|
||||||
|
! Create integral batches
|
||||||
|
|
||||||
|
allocate(OOOO(nO,nO,nO,nO),OOVV(nO,nO,nV,nV),OVOV(nO,nV,nO,nV),VVVV(nV,nV,nV,nV),OVVO(nO,nV,nV,nO))
|
||||||
|
|
||||||
|
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)
|
||||||
|
OVVO(:,:,:,:) = dbERI( 1:nO ,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(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(r2r(nO,nO,nV,nV),r2l(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(*,*)'| 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
|
||||||
|
|
||||||
|
! Form linear array
|
||||||
|
|
||||||
|
call form_u(nC,nO,nV,nR,OOOO,VVVV,OVOV,t2,u)
|
||||||
|
|
||||||
|
! Form interemediate arrays
|
||||||
|
|
||||||
|
call form_X(nC,nO,nV,nR,OOVV,t2,X1,X2,X3,X4)
|
||||||
|
|
||||||
|
! Form quadratic array
|
||||||
|
|
||||||
|
call form_v(nC,nO,nV,nR,X1,X2,X3,X4,t2,v)
|
||||||
|
|
||||||
|
call form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2r)
|
||||||
|
|
||||||
|
call form_ladder_r(nC,nO,nV,nR,OOOO,OOVV,VVVV,t2,r2l)
|
||||||
|
|
||||||
|
r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + u(:,:,:,:) + v(:,:,:,:) - r2r(:,:,:,:) - r2l(:,:,:,:)
|
||||||
|
|
||||||
|
! 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,EcCCD)
|
||||||
|
|
||||||
|
! Dump results
|
||||||
|
|
||||||
|
ECCD = ERHF + EcCCD
|
||||||
|
|
||||||
|
! 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,'|',ECCD+ENuc,'|',EcCCD,'|',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) = ',ECCD
|
||||||
|
write(*,'(1X,A30,1X,F15.10)')' Ec(crCCD) = ',EcCCD
|
||||||
|
write(*,*)'----------------------------------------------------'
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
end subroutine crCCD
|
||||||
@ -39,6 +39,7 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
|
|||||||
|
|
||||||
double precision,allocatable :: OOVV(:,:,:,:)
|
double precision,allocatable :: OOVV(:,:,:,:)
|
||||||
double precision,allocatable :: OVVO(:,:,:,:)
|
double precision,allocatable :: OVVO(:,:,:,:)
|
||||||
|
double precision,allocatable :: VOOV(:,:,:,:)
|
||||||
|
|
||||||
double precision,allocatable :: r2(:,:,:,:)
|
double precision,allocatable :: r2(:,:,:,:)
|
||||||
double precision,allocatable :: t2(:,:,:,:)
|
double precision,allocatable :: t2(:,:,:,:)
|
||||||
@ -84,10 +85,11 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
|
|||||||
|
|
||||||
! Create integral batches
|
! Create integral batches
|
||||||
|
|
||||||
allocate(OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO))
|
allocate(OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO),VOOV(nV,nO,nO,nV))
|
||||||
|
|
||||||
OOVV(:,:,:,:) = sERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas)
|
OOVV(:,:,:,:) = sERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas)
|
||||||
OVVO(:,:,:,:) = sERI( 1:nO ,nO+1:nBas,nO+1:nBas, 1:nO )
|
OVVO(:,:,:,:) = sERI( 1:nO ,nO+1:nBas,nO+1:nBas, 1:nO )
|
||||||
|
VOOV(:,:,:,:) = sERI(nO+1:nBas, 1:nO , 1:nO ,nO+1:nBas)
|
||||||
|
|
||||||
deallocate(sERI)
|
deallocate(sERI)
|
||||||
|
|
||||||
@ -133,7 +135,7 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
|
|||||||
|
|
||||||
! Compute residual
|
! Compute residual
|
||||||
|
|
||||||
call form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2)
|
call form_ring_r(nC,nO,nV,nR,OVVO,VOOV,OOVV,t2,r2)
|
||||||
|
|
||||||
r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:)
|
r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:)
|
||||||
|
|
||||||
|
|||||||
@ -1,4 +1,4 @@
|
|||||||
subroutine form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2)
|
subroutine form_ring_r(nC,nO,nV,nR,OVVO,VOOV,OOVV,t2,r2)
|
||||||
|
|
||||||
! Form residuals for ring CCD
|
! Form residuals for ring CCD
|
||||||
|
|
||||||
@ -9,6 +9,7 @@ subroutine form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2)
|
|||||||
integer,intent(in) :: nC,nO,nV,nR
|
integer,intent(in) :: nC,nO,nV,nR
|
||||||
double precision,intent(in) :: t2(nO,nO,nV,nV)
|
double precision,intent(in) :: t2(nO,nO,nV,nV)
|
||||||
double precision,intent(in) :: OVVO(nO,nV,nV,nO)
|
double precision,intent(in) :: OVVO(nO,nV,nV,nO)
|
||||||
|
double precision,intent(in) :: VOOV(nV,nO,nO,nV)
|
||||||
double precision,intent(in) :: OOVV(nO,nO,nV,nV)
|
double precision,intent(in) :: OOVV(nO,nO,nV,nV)
|
||||||
|
|
||||||
! Local variables
|
! Local variables
|
||||||
@ -49,7 +50,7 @@ subroutine form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2)
|
|||||||
do k=nC+1,nO
|
do k=nC+1,nO
|
||||||
do c=1,nV-nR
|
do c=1,nV-nR
|
||||||
|
|
||||||
r2(i,j,a,b) = r2(i,j,a,b) + OVVO(i,c,a,k)*t2(k,j,c,b) + OVVO(j,c,b,k)*t2(i,k,a,c)
|
r2(i,j,a,b) = r2(i,j,a,b) + VOOV(a,k,i,c)*t2(k,j,c,b) + OVVO(k,b,c,j)*t2(i,k,a,c)
|
||||||
|
|
||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
|
|||||||
@ -40,6 +40,7 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,
|
|||||||
|
|
||||||
double precision,allocatable :: OOVV(:,:,:,:)
|
double precision,allocatable :: OOVV(:,:,:,:)
|
||||||
double precision,allocatable :: OVVO(:,:,:,:)
|
double precision,allocatable :: OVVO(:,:,:,:)
|
||||||
|
double precision,allocatable :: VOOV(:,:,:,:)
|
||||||
|
|
||||||
double precision,allocatable :: r2(:,:,:,:)
|
double precision,allocatable :: r2(:,:,:,:)
|
||||||
double precision,allocatable :: t2(:,:,:,:)
|
double precision,allocatable :: t2(:,:,:,:)
|
||||||
@ -92,10 +93,11 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,
|
|||||||
|
|
||||||
! Create integral batches
|
! Create integral batches
|
||||||
|
|
||||||
allocate(OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO))
|
allocate(OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO),VOOV(nV,nO,nO,nV))
|
||||||
|
|
||||||
OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas)
|
OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas)
|
||||||
OVVO(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas,nO+1:nBas, 1:nO )
|
OVVO(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas,nO+1:nBas, 1:nO )
|
||||||
|
VOOV(:,:,:,:) = dbERI(nO+1:nBas, 1:nO , 1:nO ,nO+1:nBas)
|
||||||
|
|
||||||
deallocate(dbERI)
|
deallocate(dbERI)
|
||||||
|
|
||||||
@ -141,7 +143,7 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,
|
|||||||
|
|
||||||
! Compute residual
|
! Compute residual
|
||||||
|
|
||||||
call form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2)
|
call form_ring_r(nC,nO,nV,nR,OVVO,VOOV,OOVV,t2,r2)
|
||||||
|
|
||||||
r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:)
|
r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:)
|
||||||
|
|
||||||
|
|||||||
@ -1,4 +1,4 @@
|
|||||||
subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rho_RPA,OmBSE,A_dyn)
|
subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rho_RPA,OmBSE,A_dyn,ZA_dyn)
|
||||||
|
|
||||||
! Compute the dynamic part of the Bethe-Salpeter equation matrices
|
! Compute the dynamic part of the Bethe-Salpeter equation matrices
|
||||||
|
|
||||||
@ -25,10 +25,12 @@ subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,Om
|
|||||||
! Output variables
|
! Output variables
|
||||||
|
|
||||||
double precision,intent(out) :: A_dyn(nS,nS)
|
double precision,intent(out) :: A_dyn(nS,nS)
|
||||||
|
double precision,intent(out) :: ZA_dyn(nS,nS)
|
||||||
|
|
||||||
! Initialization
|
! Initialization
|
||||||
|
|
||||||
A_dyn(:,:) = 0d0
|
A_dyn(:,:) = 0d0
|
||||||
|
ZA_dyn(:,:) = 0d0
|
||||||
|
|
||||||
! Number of poles taken into account
|
! Number of poles taken into account
|
||||||
|
|
||||||
@ -67,6 +69,19 @@ subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,Om
|
|||||||
|
|
||||||
A_dyn(ia,jb) = A_dyn(ia,jb) - 2d0*lambda*chi
|
A_dyn(ia,jb) = A_dyn(ia,jb) - 2d0*lambda*chi
|
||||||
|
|
||||||
|
chi = 0d0
|
||||||
|
do kc=1,maxS
|
||||||
|
|
||||||
|
eps = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j))
|
||||||
|
chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
||||||
|
|
||||||
|
eps = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i))
|
||||||
|
chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
||||||
|
|
||||||
|
enddo
|
||||||
|
|
||||||
|
ZA_dyn(ia,jb) = ZA_dyn(ia,jb) + 2d0*lambda*chi
|
||||||
|
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
|||||||
@ -88,11 +88,7 @@ subroutine Bethe_Salpeter_Tmatrix_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR
|
|||||||
|
|
||||||
! Resonant part of the BSE correction for dynamical TDA
|
! Resonant part of the BSE correction for dynamical TDA
|
||||||
|
|
||||||
call dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,eGT,Omega1,Omega2,rho1,rho2,OmBSE(ia),Ap_dyn)
|
call dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,eGT,Omega1,Omega2,rho1,rho2,OmBSE(ia),Ap_dyn,Zap_dyn)
|
||||||
|
|
||||||
! Renormalization factor of the resonant parts for dynamical TDA
|
|
||||||
|
|
||||||
call dynamic_Tmatrix_ZA(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,eGT,Omega1,Omega2,rho1,rho2,OmBSE(ia),ZAp_dyn)
|
|
||||||
|
|
||||||
ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X))
|
ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X))
|
||||||
OmDyn(ia) = dot_product(X,matmul( Ap_dyn,X))
|
OmDyn(ia) = dot_product(X,matmul( Ap_dyn,X))
|
||||||
|
|||||||
@ -82,11 +82,7 @@ subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,e
|
|||||||
|
|
||||||
! Resonant part of the BSE correction for dynamical TDA
|
! Resonant part of the BSE correction for dynamical TDA
|
||||||
|
|
||||||
call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia),Ap_dyn)
|
call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia),Ap_dyn,ZAp_dyn)
|
||||||
|
|
||||||
! Renormalization factor of the resonant parts for dynamical TDA
|
|
||||||
|
|
||||||
call Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia),ZAp_dyn)
|
|
||||||
|
|
||||||
ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X))
|
ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X))
|
||||||
OmDyn(ia) = dot_product(X,matmul( Ap_dyn,X))
|
OmDyn(ia) = dot_product(X,matmul( Ap_dyn,X))
|
||||||
|
|||||||
@ -1,4 +1,4 @@
|
|||||||
subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,Omega2,rho1,rho2,OmBSE,A_dyn)
|
subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,Omega2,rho1,rho2,OmBSE,A_dyn,ZA_dyn)
|
||||||
|
|
||||||
! Compute the dynamic part of the Bethe-Salpeter equation matrices for GT
|
! Compute the dynamic part of the Bethe-Salpeter equation matrices for GT
|
||||||
|
|
||||||
@ -36,11 +36,13 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O
|
|||||||
|
|
||||||
! Output variables
|
! Output variables
|
||||||
|
|
||||||
double precision,intent(out) :: A_dyn(nS,nS)
|
double precision,intent(out) :: A_dyn(nS,nS)
|
||||||
|
double precision,intent(out) :: ZA_dyn(nS,nS)
|
||||||
|
|
||||||
! Initialization
|
! Initialization
|
||||||
|
|
||||||
A_dyn(:,:) = 0d0
|
A_dyn(:,:) = 0d0
|
||||||
|
ZA_dyn(:,:) = 0d0
|
||||||
|
|
||||||
! Build dynamic A matrix
|
! Build dynamic A matrix
|
||||||
|
|
||||||
@ -60,10 +62,10 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O
|
|||||||
end do
|
end do
|
||||||
|
|
||||||
do kl=1,nOO
|
do kl=1,nOO
|
||||||
chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/(Omega2(kl)**2 + eta**2)
|
chi = chi - rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/(Omega2(kl)**2 + eta**2)
|
||||||
end do
|
end do
|
||||||
|
|
||||||
A_dyn(ia,jb) = A_dyn(ia,jb) - 2d0*lambda*chi
|
A_dyn(ia,jb) = A_dyn(ia,jb) - 1d0*lambda*chi
|
||||||
|
|
||||||
chi = 0d0
|
chi = 0d0
|
||||||
|
|
||||||
@ -77,7 +79,21 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O
|
|||||||
chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2)
|
chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2)
|
||||||
end do
|
end do
|
||||||
|
|
||||||
A_dyn(ia,jb) = A_dyn(ia,jb) - 2d0*lambda*chi
|
A_dyn(ia,jb) = A_dyn(ia,jb) + 1d0*lambda*chi
|
||||||
|
|
||||||
|
chi = 0d0
|
||||||
|
|
||||||
|
do cd=1,nVV
|
||||||
|
eps = + OmBSE - Omega1(cd) + (eGT(i) + eGT(j))
|
||||||
|
chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
||||||
|
end do
|
||||||
|
|
||||||
|
do kl=1,nOO
|
||||||
|
eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b))
|
||||||
|
chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
||||||
|
end do
|
||||||
|
|
||||||
|
ZA_dyn(ia,jb) = ZA_dyn(ia,jb) - 1d0*lambda*chi
|
||||||
|
|
||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
|
|||||||
@ -1,73 +0,0 @@
|
|||||||
subroutine dynamic_Tmatrix_ZA(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,Omega2,rho1,rho2,OmBSE,ZA_dyn)
|
|
||||||
|
|
||||||
! Compute the dynamic part of the Bethe-Salpeter equation matrices
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
include 'parameters.h'
|
|
||||||
|
|
||||||
! Input variables
|
|
||||||
|
|
||||||
double precision,intent(in) :: eta
|
|
||||||
integer,intent(in) :: nBas
|
|
||||||
integer,intent(in) :: nC
|
|
||||||
integer,intent(in) :: nO
|
|
||||||
integer,intent(in) :: nV
|
|
||||||
integer,intent(in) :: nR
|
|
||||||
integer,intent(in) :: nS
|
|
||||||
|
|
||||||
integer,intent(in) :: nOO
|
|
||||||
integer,intent(in) :: nVV
|
|
||||||
|
|
||||||
double precision,intent(in) :: lambda
|
|
||||||
double precision,intent(in) :: eGT(nBas)
|
|
||||||
double precision,intent(in) :: OmBSE
|
|
||||||
|
|
||||||
double precision,intent(in) :: Omega1(nVV)
|
|
||||||
double precision,intent(in) :: Omega2(nOO)
|
|
||||||
double precision,intent(in) :: rho1(nBas,nBas,nVV)
|
|
||||||
double precision,intent(in) :: rho2(nBas,nBas,nOO)
|
|
||||||
|
|
||||||
! Local variables
|
|
||||||
|
|
||||||
double precision :: chi
|
|
||||||
double precision :: eps
|
|
||||||
integer :: i,j,a,b,ia,jb,cd,kl
|
|
||||||
|
|
||||||
! Output variables
|
|
||||||
|
|
||||||
double precision,intent(out) :: ZA_dyn(nS,nS)
|
|
||||||
|
|
||||||
! Initialization
|
|
||||||
|
|
||||||
ZA_dyn(:,:) = 0d0
|
|
||||||
|
|
||||||
! Build dynamic A matrix
|
|
||||||
|
|
||||||
ia = 0
|
|
||||||
do i=nC+1,nO
|
|
||||||
do a=nO+1,nBas-nR
|
|
||||||
ia = ia + 1
|
|
||||||
jb = 0
|
|
||||||
do j=nC+1,nO
|
|
||||||
do b=nO+1,nBas-nR
|
|
||||||
jb = jb + 1
|
|
||||||
|
|
||||||
chi = 0d0
|
|
||||||
do cd=1,nVV
|
|
||||||
eps = + OmBSE - Omega1(cd) + (eGT(i) + eGT(j))
|
|
||||||
chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
|
||||||
end do
|
|
||||||
|
|
||||||
do kl=1,nOO
|
|
||||||
eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b))
|
|
||||||
chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
|
||||||
end do
|
|
||||||
|
|
||||||
ZA_dyn(ia,jb) = ZA_dyn(ia,jb) + 2d0*lambda*chi
|
|
||||||
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
|
|
||||||
end subroutine dynamic_Tmatrix_ZA
|
|
||||||
@ -26,7 +26,6 @@ subroutine static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r
|
|||||||
! Local variables
|
! Local variables
|
||||||
|
|
||||||
double precision :: chi
|
double precision :: chi
|
||||||
double precision :: eps
|
|
||||||
integer :: i,j,a,b,ia,jb,kl,cd
|
integer :: i,j,a,b,ia,jb,kl,cd
|
||||||
|
|
||||||
! Output variables
|
! Output variables
|
||||||
@ -45,18 +44,16 @@ subroutine static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r
|
|||||||
chi = 0d0
|
chi = 0d0
|
||||||
|
|
||||||
do cd=1,nVV
|
do cd=1,nVV
|
||||||
eps = Omega1(cd)**2 + eta**2
|
! chi = chi + lambda*rho1(i,j,cd)*rho1(a,b,cd)*Omega1(cd)/(Omega1(cd)**2 + eta**2)
|
||||||
! chi = chi + lambda*rho1(i,j,cd)*rho1(a,b,cd)*Omega1(cd)/eps
|
chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*Omega1(cd)/(Omega1(cd)**2 + eta**2)
|
||||||
chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*Omega1(cd)/eps
|
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
do kl=1,nOO
|
do kl=1,nOO
|
||||||
eps = Omega2(kl)**2 + eta**2
|
! chi = chi - lambda*rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/(Omega2(kl)**2 + eta**2)
|
||||||
! chi = chi + lambda*rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/eps
|
chi = chi - rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/(Omega2(kl)**2 + eta**2)
|
||||||
chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/eps
|
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
TA(ia,jb) = TA(ia,jb) + 2d0*lambda*chi
|
TA(ia,jb) = TA(ia,jb) + 1d0*lambda*chi
|
||||||
|
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
|||||||
@ -26,7 +26,6 @@ subroutine static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r
|
|||||||
! Local variables
|
! Local variables
|
||||||
|
|
||||||
double precision :: chi
|
double precision :: chi
|
||||||
double precision :: eps
|
|
||||||
integer :: i,j,a,b,ia,jb,kl,cd
|
integer :: i,j,a,b,ia,jb,kl,cd
|
||||||
|
|
||||||
! Output variables
|
! Output variables
|
||||||
@ -45,18 +44,16 @@ subroutine static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r
|
|||||||
chi = 0d0
|
chi = 0d0
|
||||||
|
|
||||||
do cd=1,nVV
|
do cd=1,nVV
|
||||||
eps = Omega1(cd)**2 + eta**2
|
! chi = chi + lambda*rho1(i,b,cd)*rho1(a,j,cd)*Omega1(cd)/Omega1(cd)**2 + eta**2
|
||||||
! chi = chi + lambda*rho1(i,b,cd)*rho1(a,j,cd)*Omega1(cd)/eps
|
chi = chi + rho1(i,b,cd)*rho1(a,j,cd)*Omega1(cd)/Omega1(cd)**2 + eta**2
|
||||||
chi = chi + rho1(i,b,cd)*rho1(a,j,cd)*Omega1(cd)/eps
|
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
do kl=1,nOO
|
do kl=1,nOO
|
||||||
eps = Omega2(kl)**2 + eta**2
|
! chi = chi - lambda*rho2(i,b,kl)*rho2(a,j,kl)*Omega2(kl)/Omega2(kl)**2 + eta**2
|
||||||
! chi = chi + lambda*rho2(i,b,kl)*rho2(a,j,kl)*Omega2(kl)/eps
|
chi = chi - rho2(i,b,kl)*rho2(a,j,kl)*Omega2(kl)/Omega2(kl)**2 + eta**2
|
||||||
chi = chi + rho2(i,b,kl)*rho2(a,j,kl)*Omega2(kl)/eps
|
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
TB(ia,jb) = TB(ia,jb) + 2d0*lambda*chi
|
TB(ia,jb) = TB(ia,jb) + 1d0*lambda*chi
|
||||||
|
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
|||||||
@ -10,7 +10,7 @@ program QuAcK
|
|||||||
logical :: doKS
|
logical :: doKS
|
||||||
logical :: doMP2,doMP3,doMP2F12
|
logical :: doMP2,doMP3,doMP2F12
|
||||||
logical :: doCCD,doDCD,doCCSD,doCCSDT
|
logical :: doCCD,doDCD,doCCSD,doCCSDT
|
||||||
logical :: do_drCCD,do_rCCD,do_lCCD,do_pCCD
|
logical :: do_drCCD,do_rCCD,do_lCCD,do_crCCD,do_pCCD
|
||||||
logical :: doCIS,doCIS_D,doCID,doCISD,doFCI
|
logical :: doCIS,doCIS_D,doCID,doCISD,doFCI
|
||||||
logical :: doRPA,doRPAx,doppRPA
|
logical :: doRPA,doRPAx,doppRPA
|
||||||
logical :: doADC
|
logical :: doADC
|
||||||
@ -662,6 +662,27 @@ program QuAcK
|
|||||||
|
|
||||||
end if
|
end if
|
||||||
|
|
||||||
|
!------------------------------------------------------------------------
|
||||||
|
! Perform crossed-ring CCD calculation
|
||||||
|
!------------------------------------------------------------------------
|
||||||
|
|
||||||
|
do_crCCD = .false.
|
||||||
|
|
||||||
|
if(do_crCCD) then
|
||||||
|
|
||||||
|
call cpu_time(start_CCD)
|
||||||
|
call ehRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
|
||||||
|
call crCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR, &
|
||||||
|
ERI_MO,ENuc,ERHF,eHF)
|
||||||
|
|
||||||
|
call cpu_time(end_CCD)
|
||||||
|
|
||||||
|
t_CCD = end_CCD - start_CCD
|
||||||
|
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for crossed-ring CCD = ',t_CCD,' seconds'
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
!------------------------------------------------------------------------
|
!------------------------------------------------------------------------
|
||||||
! Perform pair CCD calculation
|
! Perform pair CCD calculation
|
||||||
!------------------------------------------------------------------------
|
!------------------------------------------------------------------------
|
||||||
|
|||||||
130
src/RPA/ehRPA.f90
Normal file
130
src/RPA/ehRPA.f90
Normal file
@ -0,0 +1,130 @@
|
|||||||
|
subroutine ehRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, &
|
||||||
|
ERI,dipole_int,eHF)
|
||||||
|
|
||||||
|
! Perform random phase approximation calculation with exchange (aka TDHF)
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
include 'parameters.h'
|
||||||
|
include 'quadrature.h'
|
||||||
|
|
||||||
|
! Input variables
|
||||||
|
|
||||||
|
logical,intent(in) :: TDA
|
||||||
|
logical,intent(in) :: doACFDT
|
||||||
|
logical,intent(in) :: exchange_kernel
|
||||||
|
logical,intent(in) :: singlet
|
||||||
|
double precision,intent(in) :: eta
|
||||||
|
logical,intent(in) :: triplet
|
||||||
|
integer,intent(in) :: nBas
|
||||||
|
integer,intent(in) :: nC
|
||||||
|
integer,intent(in) :: nO
|
||||||
|
integer,intent(in) :: nV
|
||||||
|
integer,intent(in) :: nR
|
||||||
|
integer,intent(in) :: nS
|
||||||
|
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)
|
||||||
|
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
|
||||||
|
|
||||||
|
! Local variables
|
||||||
|
|
||||||
|
integer :: ispin
|
||||||
|
double precision,allocatable :: Omega(:,:)
|
||||||
|
double precision,allocatable :: XpY(:,:,:)
|
||||||
|
double precision,allocatable :: XmY(:,:,:)
|
||||||
|
|
||||||
|
double precision :: rho
|
||||||
|
double precision :: EcRPAx(nspin)
|
||||||
|
double precision :: EcAC(nspin)
|
||||||
|
|
||||||
|
! Hello world
|
||||||
|
|
||||||
|
write(*,*)
|
||||||
|
write(*,*)'***********************************************************'
|
||||||
|
write(*,*)'| Random phase approximation calculation: eh channel |'
|
||||||
|
write(*,*)'***********************************************************'
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
! TDA
|
||||||
|
|
||||||
|
if(TDA) then
|
||||||
|
write(*,*) 'Tamm-Dancoff approximation activated!'
|
||||||
|
write(*,*)
|
||||||
|
end if
|
||||||
|
|
||||||
|
! Initialization
|
||||||
|
|
||||||
|
EcRPAx(:) = 0d0
|
||||||
|
EcAC(:) = 0d0
|
||||||
|
|
||||||
|
! Memory allocation
|
||||||
|
|
||||||
|
allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin))
|
||||||
|
|
||||||
|
! Singlet manifold
|
||||||
|
|
||||||
|
if(singlet) then
|
||||||
|
|
||||||
|
ispin = 1
|
||||||
|
|
||||||
|
call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI,Omega(:,ispin),rho, &
|
||||||
|
EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
|
||||||
|
call print_excitation('ehRPA@HF ',ispin,nS,Omega(:,ispin))
|
||||||
|
call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
|
||||||
|
|
||||||
|
endif
|
||||||
|
|
||||||
|
! Triplet manifold
|
||||||
|
|
||||||
|
if(triplet) then
|
||||||
|
|
||||||
|
ispin = 2
|
||||||
|
|
||||||
|
call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI,Omega(:,ispin),rho, &
|
||||||
|
EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
|
||||||
|
call print_excitation('ehRPA@HF ',ispin,nS,Omega(:,ispin))
|
||||||
|
call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
|
||||||
|
|
||||||
|
endif
|
||||||
|
|
||||||
|
if(exchange_kernel) then
|
||||||
|
|
||||||
|
EcRPAx(1) = 0.5d0*EcRPAx(1)
|
||||||
|
EcRPAx(2) = 1.5d0*EcRPAx(2)
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
write(*,*)
|
||||||
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
write(*,'(2X,A50,F20.10)') 'Tr@ehRPA correlation energy (singlet) =',EcRPAx(1)
|
||||||
|
write(*,'(2X,A50,F20.10)') 'Tr@ehRPA correlation energy (triplet) =',EcRPAx(2)
|
||||||
|
write(*,'(2X,A50,F20.10)') 'Tr@ehRPA correlation energy =',EcRPAx(1) + EcRPAx(2)
|
||||||
|
write(*,'(2X,A50,F20.10)') 'Tr@ehRPA total energy =',ENuc + ERHF + EcRPAx(1) + EcRPAx(2)
|
||||||
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
! Compute the correlation energy via the adiabatic connection
|
||||||
|
|
||||||
|
! if(doACFDT) then
|
||||||
|
|
||||||
|
! write(*,*) '-------------------------------------------------------'
|
||||||
|
! write(*,*) 'Adiabatic connection version of ehRPA correlation energy'
|
||||||
|
! write(*,*) '-------------------------------------------------------'
|
||||||
|
! write(*,*)
|
||||||
|
|
||||||
|
! call ACFDT(exchange_kernel,.false.,.false.,.false.,TDA,.false.,singlet,triplet,eta, &
|
||||||
|
! nBas,nC,nO,nV,nR,nS,ERI,eHF,eHF,EcAC)
|
||||||
|
|
||||||
|
! write(*,*)
|
||||||
|
! write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
! write(*,'(2X,A50,F20.10)') 'AC@RPAx correlation energy (singlet) =',EcAC(1)
|
||||||
|
! write(*,'(2X,A50,F20.10)') 'AC@RPAx correlation energy (triplet) =',EcAC(2)
|
||||||
|
! write(*,'(2X,A50,F20.10)') 'AC@RPAx correlation energy =',EcAC(1) + EcAC(2)
|
||||||
|
! write(*,'(2X,A50,F20.10)') 'AC@RPAx total energy =',ENuc + ERHF + EcAC(1) + EcAC(2)
|
||||||
|
! write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
! write(*,*)
|
||||||
|
|
||||||
|
! end if
|
||||||
|
|
||||||
|
end subroutine ehRPA
|
||||||
Loading…
Reference in New Issue
Block a user