4
1
mirror of https://github.com/pfloos/quack synced 2024-11-07 06:33:55 +01:00

working on RPA third channel

This commit is contained in:
Pierre-Francois Loos 2021-11-10 09:42:30 +01:00
parent e13c9866fb
commit cec0fe9205
17 changed files with 468 additions and 122 deletions

View File

@ -5,17 +5,17 @@
# 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
T F F F F F F F F F
# G0W0* evGW* qsGW* ufG0W0 ufGW # G0W0* evGW* qsGW* ufG0W0 ufGW
T F F F F F F F F F
# G0T0 evGT qsGT # G0T0 evGT qsGT
T F F F F F
# MCMP2 # MCMP2
F F
# * unrestricted version available # * unrestricted version available

View File

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

View File

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

View File

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

View File

@ -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(:,:,:,:)

View File

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

View File

@ -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(:,:,:,:)

View File

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

View File

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

View File

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

View File

@ -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
@ -63,7 +65,7 @@ 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)*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

View File

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

View File

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

View File

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

View File

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