LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-06-18 11:15:30 +02:00
Code Releases Activity

Merge branch 'master' of https://github.com/pfloos/QuAcK

Browse Source
This commit is contained in:
Clotilde Marut 2021-11-22 22:55:10 +01:00
commit 210e4b5e7c
30 changed files with 940 additions and 284 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
input/methods
View File

@ -2,14 +2,14 @@
T F F F
# MP2* MP3 MP2-F12
F F F
# CCD DCD CCSD CCSD(T)
F F F F
# drCCD rCCD lCCD pCCD
T T F F
# CCD pCCD DCD CCSD CCSD(T)
F F F F F
# drCCD rCCD crCCD lCCD
F F F F
# CIS* CIS(D) CID CISD FCI
F F F F F
# RPA* RPAx* ppRPA
T T F
# RPA* RPAx* crRPA ppRPA
F F F T
# G0F2* evGF2* qsGF2* G0F3 evGF3
F F F F F
# G0W0* evGW* qsGW* ufG0W0 ufGW

4
input/options
View File

@ -11,8 +11,8 @@
# GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0
256 0.00001 T 5 T 0.00367493 F F F F F
# ACFDT: AC Kx XBS
F F T
T T F
# BSE: BSE dBSE dTDA evDyn
T T T F
T F T F
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift
1000000 100000 10 0.3 10000 1234 T

15
mol/phenol.xyz Normal file
View File

@ -0,0 +1,15 @@
13
Phenol; structure form HCP92 revised structure; l
C 4.555420 5.661760 4.489060
C 4.584960 2.843420 4.498910
C 3.378760 3.548270 4.498890
C 3.359200 4.943960 4.500070
C 5.758650 4.948670 4.502800
C 5.789840 3.550210 4.500330
H 4.622340 1.760920 4.495600
H 3.631910 7.321310 4.523190
H 2.474610 2.963630 4.498170
H 2.384290 5.419790 4.496670
H 6.695040 5.492670 4.498160
H 6.727200 3.021210 4.497160
O 4.537700 7.024010 4.500450

2
src/AOtoMO/AOtoMO_integral_transform.f90
View File

@ -79,7 +79,7 @@ subroutine AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,c,ERI_AO_basis,ERI
do nu=1,nBas
ERI_MO_basis(i,j,k,l) = ERI_MO_basis(i,j,k,l) + c(nu,j,ket1)*scr(i,nu,k,l)
enddo
! print*,i,k,j,l,ERI_MO_basis(i,j,k,l)
! write(11,'(I5,I5,I5,I5,F16.10)') i,j,k,l,ERI_MO_basis(i,j,k,l)
enddo
enddo
enddo

9
src/AOtoMO/AOtoMO_transform.f90
View File

@ -9,10 +9,19 @@ subroutine AOtoMO_transform(nBas,c,A)
integer,intent(in) :: nBas
double precision,intent(in) :: c(nBas,nBas)
integer :: i,j
! Output variables
double precision,intent(inout):: A(nBas,nBas)
A = matmul(transpose(c),matmul(A,c))
! do j=1,nBas
! do i=1,nBas
! write(10,'(I5,I5,F16.10)') i,j,A(i,j)
! enddo
! enddo
end subroutine AOtoMO_transform

54
src/CC/crCCD.f90
View File

@ -38,21 +38,8 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
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(:,:,:,:)
@ -66,7 +53,7 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
write(*,*)
write(*,*)'**************************************'
write(*,*)'| crossed-ring CCD calculation |'
write(*,*)'| Crossed-ring CCD calculation |'
write(*,*)'**************************************'
write(*,*)
@ -105,13 +92,10 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
! 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))
allocate(OOVV(nO,nO,nV,nV),OVOV(nO,nV,nO,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)
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)
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)
@ -129,9 +113,7 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
! 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))
allocate(r2(nO,nO,nV,nV))
Conv = 1d0
nSCF = 0
@ -159,23 +141,9 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
! Compute residual
! Form linear array
call form_crossed_ring_r(nC,nO,nV,nR,OVOV,OOVV,t2,r2)
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(:,:,:,:)
r2(:,:,:,:) = OOVV(:,:,:,:) - delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:)
! Check convergence
@ -183,7 +151,7 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
! Update amplitudes
t2(:,:,:,:) = t2(:,:,:,:) - r2(:,:,:,:)/delta_OOVV(:,:,:,:)
t2(:,:,:,:) = t2(:,:,:,:) + r2(:,:,:,:)/delta_OOVV(:,:,:,:)
! Compute correlation energy
@ -227,10 +195,10 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
write(*,*)
write(*,*)'----------------------------------------------------'
write(*,*)' crossed-ring CCD energy '
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(*,'(1X,A30,1X,F15.10)')' E(crCCD) = ',ECCD
write(*,'(1X,A30,1X,F15.10)')' Ec(crCCD) = ',EcCCD
write(*,*)'----------------------------------------------------'
write(*,*)

6
src/CC/drCCD.f90
View File

@ -39,7 +39,6 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
double precision,allocatable :: OOVV(:,:,:,:)
double precision,allocatable :: OVVO(:,:,:,:)
double precision,allocatable :: VOOV(:,:,:,:)
double precision,allocatable :: r2(:,:,:,:)
double precision,allocatable :: t2(:,:,:,:)
@ -85,11 +84,10 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
! Create integral batches
allocate(OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO),VOOV(nV,nO,nO,nV))
allocate(OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO))
OOVV(:,:,:,:) = sERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas)
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)
@ -135,7 +133,7 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
! Compute residual
call form_ring_r(nC,nO,nV,nR,OVVO,VOOV,OOVV,t2,r2)
call form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2)
r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:)

70
src/CC/form_crossed_ring_r.f90 Normal file
View File

@ -0,0 +1,70 @@
subroutine form_crossed_ring_r(nC,nO,nV,nR,OVOV,OOVV,t2,r2)
! Form residuals for crossed-ring CCD
implicit none
! Input variables
integer,intent(in) :: nC,nO,nV,nR
double precision,intent(in) :: t2(nO,nO,nV,nV)
double precision,intent(in) :: OVOV(nO,nV,nO,nV)
double precision,intent(in) :: OOVV(nO,nO,nV,nV)
! Local variables
integer :: i,j,k,l
integer :: a,b,c,d
double precision,allocatable :: y(:,:,:,:)
! Output variables
double precision,intent(out) :: r2(nO,nO,nV,nV)
r2(:,:,:,:) = 0d0
allocate(y(nV,nO,nO,nV))
y(:,:,:,:) = 0d0
do i=nC+1,nO
do b=1,nV-nR
do l=nC+1,nO
do d=1,nV-nR
do k=nC+1,nO
do c=1,nV-nR
y(b,i,l,d) = y(b,i,l,d) + t2(i,k,c,b)*OOVV(k,l,c,d)
end do
end do
end do
end do
end do
end do
do i=nC+1,nO
do j=nC+1,nO
do a=1,nV-nR
do b=1,nV-nR
do k=nC+1,nO
do c=1,nV-nR
r2(i,j,a,b) = r2(i,j,a,b) - OVOV(k,b,i,c)*t2(k,j,a,c) - OVOV(k,a,j,c)*t2(i,k,c,b)
end do
end do
do l=nC+1,nO
do d=1,nV-nR
r2(i,j,a,b) = r2(i,j,a,b) - y(b,i,l,d)*t2(l,j,a,d)
end do
end do
end do
end do
end do
end do
end subroutine form_crossed_ring_r

5
src/CC/form_ring_r.f90
View File

@ -1,4 +1,4 @@
subroutine form_ring_r(nC,nO,nV,nR,OVVO,VOOV,OOVV,t2,r2)
subroutine form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2)
! Form residuals for ring CCD
@ -9,7 +9,6 @@ subroutine form_ring_r(nC,nO,nV,nR,OVVO,VOOV,OOVV,t2,r2)
integer,intent(in) :: nC,nO,nV,nR
double precision,intent(in) :: t2(nO,nO,nV,nV)
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)
! Local variables
@ -50,7 +49,7 @@ subroutine form_ring_r(nC,nO,nV,nR,OVVO,VOOV,OOVV,t2,r2)
do k=nC+1,nO
do c=1,nV-nR
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)
r2(i,j,a,b) = r2(i,j,a,b) + OVVO(k,a,c,i)*t2(k,j,c,b) + OVVO(k,b,c,j)*t2(i,k,a,c)
end do
end do

6
src/CC/rCCD.f90
View File

@ -40,7 +40,6 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,
double precision,allocatable :: OOVV(:,:,:,:)
double precision,allocatable :: OVVO(:,:,:,:)
double precision,allocatable :: VOOV(:,:,:,:)
double precision,allocatable :: r2(:,:,:,:)
double precision,allocatable :: t2(:,:,:,:)
@ -93,11 +92,10 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,
! Create integral batches
allocate(OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO),VOOV(nV,nO,nO,nV))
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 )
VOOV(:,:,:,:) = dbERI(nO+1:nBas, 1:nO , 1:nO ,nO+1:nBas)
deallocate(dbERI)
@ -143,7 +141,7 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,
! Compute residual
call form_ring_r(nC,nO,nV,nR,OVVO,VOOV,OOVV,t2,r2)
call form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2)
r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:)

4
src/LR/linear_response_B_pp.f90
View File

@ -1,4 +1,4 @@
subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,B_pp)
subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,B_pp)
! Compute the B matrix of the pp channel
@ -10,7 +10,7 @@ subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,B_pp
integer,intent(in) :: ispin
integer,intent(in) :: nBas,nC,nO,nV,nR,nOO,nVV
double precision,intent(in) :: lambda
double precision,intent(in) :: e(nBas),ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables

91
src/LR/linear_response_C_pp_od.f90 Normal file
View File

@ -0,0 +1,91 @@
subroutine linear_response_C_pp_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,C_pp)
! Compute the C matrix of the pp channel (without diagonal term)
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: ispin
integer,intent(in) :: nBas,nC,nO,nV,nR,nOO,nVV
double precision,intent(in) :: lambda
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
double precision,external :: Kronecker_delta
integer :: a,b,c,d,ab,cd
! Output variables
double precision,intent(out) :: C_pp(nVV,nVV)
! Build C matrix for the singlet manifold
if(ispin == 1) then
ab = 0
do a=nO+1,nBas-nR
do b=a,nBas-nR
ab = ab + 1
cd = 0
do c=nO+1,nBas-nR
do d=c,nBas-nR
cd = cd + 1
C_pp(ab,cd) = lambda*(ERI(a,b,c,d) + ERI(a,b,d,c))/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d)))
end do
end do
end do
end do
end if
! Build C matrix for the triplet manifold, or alpha-alpha block, or in the spin-orbital basis
if(ispin == 2 .or. ispin == 4) then
ab = 0
do a=nO+1,nBas-nR
do b=a+1,nBas-nR
ab = ab + 1
cd = 0
do c=nO+1,nBas-nR
do d=c+1,nBas-nR
cd = cd + 1
C_pp(ab,cd) = lambda*(ERI(a,b,c,d) - ERI(a,b,d,c))
end do
end do
end do
end do
end if
! Build the alpha-beta block of the C matrix
if(ispin == 3) then
ab = 0
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
ab = ab + 1
cd = 0
do c=nO+1,nBas-nR
do d=nO+1,nBas-nR
cd = cd + 1
C_pp(ab,cd) = lambda*ERI(a,b,c,d)
end do
end do
end do
end do
end if
end subroutine linear_response_C_pp_od

91
src/LR/linear_response_D_pp_od.f90 Normal file
View File

@ -0,0 +1,91 @@
subroutine linear_response_D_pp_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,D_pp)
! Compute the D matrix of the pp channel (without the diagonal term)
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: ispin
integer,intent(in) :: nBas,nC,nO,nV,nR,nOO,nVV
double precision,intent(in) :: lambda
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
double precision,external :: Kronecker_delta
integer :: i,j,k,l,ij,kl
! Output variables
double precision,intent(out) :: D_pp(nOO,nOO)
! Build the D matrix for the singlet manifold
if(ispin == 1) then
ij = 0
do i=nC+1,nO
do j=i,nO
ij = ij + 1
kl = 0
do k=nC+1,nO
do l=k,nO
kl = kl + 1
D_pp(ij,kl) = lambda*(ERI(i,j,k,l) + ERI(i,j,l,k))/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l)))
end do
end do
end do
end do
end if
! Build the D matrix for the triplet manifold, the alpha-alpha block, or in the spin-orbital basis
if(ispin == 2 .or. ispin == 4) then
ij = 0
do i=nC+1,nO
do j=i+1,nO
ij = ij + 1
kl = 0
do k=nC+1,nO
do l=k+1,nO
kl = kl + 1
D_pp(ij,kl) = lambda*(ERI(i,j,k,l) - ERI(i,j,l,k))
end do
end do
end do
end do
end if
! Build the alpha-beta block of the D matrix
if(ispin == 3) then
ij = 0
do i=nC+1,nO
do j=nC+1,nO
ij = ij + 1
kl = 0
do k=nC+1,nO
do l=nC+1,nO
kl = kl + 1
D_pp(ij,kl) = lambda*ERI(i,j,k,l)
end do
end do
end do
end do
end if
end subroutine linear_response_D_pp_od

2
src/LR/linear_response_pp.f90
View File

@ -75,7 +75,7 @@ subroutine linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,Om
else
call linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,B)
call linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,B)
! Diagonal blocks

11
src/MBPT/G0T0.f90
View File

@ -218,6 +218,11 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
! Free memory
deallocate(Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s, &
Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,rho1t,rho2t)
! Compute the BSE correlation energy via the adiabatic connection
if(doACFDT) then
@ -234,13 +239,13 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
end if
call ACFDT(exchange_kernel,doXBS,.true.,TDA_T,TDA,BSE,singlet,triplet,eta, &
nBas,nC,nO,nV,nR,nS,ERI_MO,eHF,eG0T0,EcAC)
call ACFDT_Tmatrix(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, &
ERI_MO,eHF,eG0T0,EcAC)
if(exchange_kernel) then
EcAC(1) = 0.5d0*EcAC(1)
EcAC(2) = 1.5d0*EcAC(1)
EcAC(2) = 1.5d0*EcAC(2)
end if

4
src/MBPT/dynamic_Tmatrix_A.f90
View File

@ -76,7 +76,7 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O
do kl=1,nOO
eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b))
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
A_dyn(ia,jb) = A_dyn(ia,jb) + 1d0*lambda*chi
@ -89,7 +89,7 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O
end do
do kl=1,nOO
eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b))
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

6
src/MBPT/evGT.f90
View File

@ -281,13 +281,13 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
end if
call ACFDT(exchange_kernel,doXBS,.true.,TDA_T,TDA,BSE,singlet,triplet,eta, &
nBas,nC,nO,nV,nR,nS,ERI_MO,eGT,eGT,EcAC)
call ACFDT_Tmatrix(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, &
ERI_MO,eGT,eGT,EcAC)
if(exchange_kernel) then
EcAC(1) = 0.5d0*EcAC(1)
EcAC(2) = 1.5d0*EcAC(1)
EcAC(2) = 1.5d0*EcAC(2)
end if

3
src/MBPT/qsGT.f90
View File

@ -381,7 +381,8 @@ subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,T
end if
call ACFDT(exchange_kernel,doXBS,.true.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGT,eGT,EcAC)
call ACFDT_Tmatrix(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, &
ERI_MO,eGT,eGT,EcAC)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'

2
src/MBPT/static_Tmatrix_TA.f90
View File

@ -50,7 +50,7 @@ subroutine static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r
do kl=1,nOO
! chi = chi - lambda*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)
chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/(Omega2(kl)**2 + eta**2)
enddo
TA(ia,jb) = TA(ia,jb) + 1d0*lambda*chi

4
src/MBPT/static_Tmatrix_TB.f90
View File

@ -49,8 +49,8 @@ subroutine static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r
enddo
do kl=1,nOO
! chi = chi - lambda*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)/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 + rho2(i,b,kl)*rho2(a,j,kl)*Omega2(kl)/Omega2(kl)**2 + eta**2
enddo
TB(ia,jb) = TB(ia,jb) + 1d0*lambda*chi

107
src/QuAcK/QuAcK.f90
View File

@ -9,10 +9,10 @@ program QuAcK
logical :: dostab
logical :: doKS
logical :: doMP2,doMP3,doMP2F12
logical :: doCCD,doDCD,doCCSD,doCCSDT
logical :: do_drCCD,do_rCCD,do_lCCD,do_crCCD,do_pCCD
logical :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT
logical :: do_drCCD,do_rCCD,do_crCCD,do_lCCD
logical :: doCIS,doCIS_D,doCID,doCISD,doFCI
logical :: doRPA,doRPAx,doppRPA
logical :: doRPA,doRPAx,docrRPA,doppRPA
logical :: doADC
logical :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3
logical :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW
@ -91,8 +91,6 @@ program QuAcK
double precision :: start_CISD ,end_CISD ,t_CISD
double precision :: start_FCI ,end_FCI ,t_FCI
double precision :: start_RPA ,end_RPA ,t_RPA
double precision :: start_RPAx ,end_RPAx ,t_RPAx
double precision :: start_ppRPA ,end_ppRPA ,t_ppRPA
double precision :: start_ADC ,end_ADC ,t_ADC
double precision :: start_GF2 ,end_GF2 ,t_GF2
double precision :: start_GF3 ,end_GF3 ,t_GF3
@ -160,17 +158,17 @@ program QuAcK
! Which calculations do you want to do?
call read_methods(doRHF,doUHF,doKS,doMOM, &
doMP2,doMP3,doMP2F12, &
doCCD,doDCD,doCCSD,doCCSDT, &
do_drCCD,do_rCCD,do_lCCD,do_pCCD, &
doCIS,doCIS_D,doCID,doCISD,doFCI, &
doRPA,doRPAx,doppRPA, &
doG0F2,doevGF2,doqsGF2, &
doG0F3,doevGF3, &
doG0W0,doevGW,doqsGW, &
doufG0W0,doufGW, &
doG0T0,doevGT,doqsGT, &
call read_methods(doRHF,doUHF,doKS,doMOM, &
doMP2,doMP3,doMP2F12, &
doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
do_drCCD,do_rCCD,do_crCCD,do_lCCD, &
doCIS,doCIS_D,doCID,doCISD,doFCI, &
doRPA,doRPAx,docrRPA,doppRPA, &
doG0F2,doevGF2,doqsGF2, &
doG0F3,doevGF3, &
doG0W0,doevGW,doqsGW, &
doufG0W0,doufGW, &
doG0T0,doevGT,doqsGT, &
doMCMP2)
! Read options for methods
@ -444,7 +442,7 @@ program QuAcK
ket1 = 1
ket2 = 1
call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO)
! call AOtoMO_transform(nBas,cHF,T+V)
end if
end if
@ -645,6 +643,24 @@ program QuAcK
end if
!------------------------------------------------------------------------
! Perform crossed-ring CCD calculation
!------------------------------------------------------------------------
if(do_crCCD) then
call cpu_time(start_CCD)
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 ladder CCD calculation
!------------------------------------------------------------------------
@ -662,32 +678,11 @@ program QuAcK
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
!------------------------------------------------------------------------
if(do_pCCD) then
if(dopCCD) then
call cpu_time(start_CCD)
call pCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
@ -787,7 +782,7 @@ program QuAcK
if(doRPAx) then
call cpu_time(start_RPAx)
call cpu_time(start_RPA)
if(unrestricted) then
call URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, &
@ -798,10 +793,26 @@ program QuAcK
call RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
end if
call cpu_time(end_RPAx)
call cpu_time(end_RPA)
t_RPAx = end_RPAx - start_RPAx
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPAx = ',t_RPAx,' seconds'
t_RPA = end_RPA - start_RPA
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPAx = ',t_RPA,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Compute cr-RPA excitations
!------------------------------------------------------------------------
if(docrRPA) then
call cpu_time(start_RPA)
call crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
call cpu_time(end_RPA)
t_RPA = end_RPA - start_RPA
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pp-RPA = ',t_RPA,' seconds'
write(*,*)
end if
@ -812,12 +823,12 @@ program QuAcK
if(doppRPA) then
call cpu_time(start_ppRPA)
call ppRPA(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF)
call cpu_time(end_ppRPA)
call cpu_time(start_RPA)
call ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF)
call cpu_time(end_RPA)
t_ppRPA = end_ppRPA - start_ppRPA
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pp-RPA = ',t_ppRPA,' seconds'
t_RPA = end_RPA - start_RPA
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pp-RPA = ',t_RPA,' seconds'
write(*,*)
end if

52
src/QuAcK/read_methods.f90
View File

@ -1,14 +1,14 @@
subroutine read_methods(doRHF,doUHF,doKS,doMOM, &
doMP2,doMP3,doMP2F12, &
doCCD,doDCD,doCCSD,doCCSDT, &
do_drCCD,do_rCCD,do_lCCD,do_pCCD, &
doCIS,doCIS_D,doCID,doCISD,doFCI, &
doRPA,doRPAx,doppRPA, &
doG0F2,doevGF2,doqsGF2, &
doG0F3,doevGF3, &
doG0W0,doevGW,doqsGW, &
doufG0W0,doufGW, &
doG0T0,doevGT,doqsGT, &
subroutine read_methods(doRHF,doUHF,doKS,doMOM, &
doMP2,doMP3,doMP2F12, &
doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
do_drCCD,do_rCCD,do_crCCD,do_lCCD, &
doCIS,doCIS_D,doCID,doCISD,doFCI, &
doRPA,doRPAx,docrRPA,doppRPA, &
doG0F2,doevGF2,doqsGF2, &
doG0F3,doevGF3, &
doG0W0,doevGW,doqsGW, &
doufG0W0,doufGW, &
doG0T0,doevGT,doqsGT, &
doMCMP2)
! Read desired methods
@ -19,10 +19,10 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, &
logical,intent(out) :: doRHF,doUHF,doKS,doMOM
logical,intent(out) :: doMP2,doMP3,doMP2F12
logical,intent(out) :: doCCD,doDCD,doCCSD,doCCSDT
logical,intent(out) :: do_drCCD,do_rCCD,do_lCCD,do_pCCD
logical,intent(out) :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT
logical,intent(out) :: do_drCCD,do_rCCD,do_crCCD,do_lCCD
logical,intent(out) :: doCIS,doCIS_D,doCID,doCISD,doFCI
logical,intent(out) :: doRPA,doRPAx,doppRPA
logical,intent(out) :: doRPA,doRPAx,docrRPA,doppRPA
logical,intent(out) :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3
logical,intent(out) :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW
logical,intent(out) :: doG0T0,doevGT,doqsGT
@ -48,14 +48,15 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, &
doMP2F12 = .false.
doCCD = .false.
dopCCD = .false.
doDCD = .false.
doCCSD = .false.
doCCSDT = .false.
do_drCCD = .false.
do_rCCD = .false.
do_crCCD = .false.
do_lCCD = .false.
do_pCCD = .false.
doCIS = .false.
doCIS_D = .false.
@ -65,6 +66,7 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, &
doRPA = .false.
doRPAx = .false.
docrRPA = .false.
doppRPA = .false.
doG0F2 = .false.
@ -105,19 +107,20 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, &
! Read CC methods
read(1,*)
read(1,*) answer1,answer2,answer3,answer4
if(answer1 == 'T') doCCD = .true.
if(answer2 == 'T') doDCD = .true.
if(answer3 == 'T') doCCSD = .true.
if(answer4 == 'T') doCCSDT = .true.
read(1,*) answer1,answer2,answer3,answer4,answer5
if(answer1 == 'T') doCCD = .true.
if(answer2 == 'T') dopCCD = .true.
if(answer3 == 'T') doDCD = .true.
if(answer4 == 'T') doCCSD = .true.
if(answer5 == 'T') doCCSDT = .true.
! Read weird CC methods
read(1,*)
read(1,*) answer1,answer2,answer3,answer4
if(answer1 == 'T') do_drCCD = .true.
if(answer2 == 'T') do_rCCD = .true.
if(answer3 == 'T') do_lCCD = .true.
if(answer4 == 'T') do_pCCD = .true.
if(answer3 == 'T') do_crCCD = .true.
if(answer4 == 'T') do_lCCD = .true.
! Read excited state methods
@ -131,10 +134,11 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, &
if(doCIS_D) doCIS = .true.
read(1,*)
read(1,*) answer1,answer2,answer3
read(1,*) answer1,answer2,answer3,answer4
if(answer1 == 'T') doRPA = .true.
if(answer2 == 'T') doRPAx = .true.
if(answer3 == 'T') doppRPA = .true.
if(answer3 == 'T') docrRPA = .true.
if(answer4 == 'T') doppRPA = .true.
! Read Green function methods

117
src/RPA/ACFDT_Tmatrix.f90
View File

@ -1,5 +1,4 @@
subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s,Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,rho1t,rho2t, &
subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, &
ERI,eT,eGT,EcAC)
! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem for the T-matrix
@ -27,32 +26,10 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
integer,intent(in) :: nR
integer,intent(in) :: nS
integer,intent(in) :: nOOs
integer,intent(in) :: nOOt
integer,intent(in) :: nVVs
integer,intent(in) :: nVVt
double precision,intent(in) :: eT(nBas)
double precision,intent(in) :: eGT(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: Omega1s(nVVs)
double precision,intent(in) :: X1s(nVVs,nVVs)
double precision,intent(in) :: Y1s(nOOs,nVVs)
double precision,intent(in) :: Omega2s(nOOs)
double precision,intent(in) :: X2s(nVVs,nOOs)
double precision,intent(in) :: Y2s(nOOs,nOOs)
double precision,intent(in) :: rho1s(nBas,nBas,nVVs)
double precision,intent(in) :: rho2s(nBas,nBas,nOOs)
double precision,intent(in) :: Omega1t(nVVt)
double precision,intent(in) :: X1t(nVVt,nVVt)
double precision,intent(in) :: Y1t(nOOt,nVVt)
double precision,intent(in) :: Omega2t(nOOt)
double precision,intent(in) :: X2t(nVVt,nOOt)
double precision,intent(in) :: Y2t(nOOt,nOOt)
double precision,intent(in) :: rho1t(nBas,nBas,nVVt)
double precision,intent(in) :: rho2t(nBas,nBas,nOOt)
! Local variables
integer :: ispin
@ -62,6 +39,9 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
double precision :: lambda
double precision,allocatable :: Ec(:,:)
integer :: nOOs,nOOt
integer :: nVVs,nVVt
double precision :: EcRPA(nspin)
double precision,allocatable :: TA(:,:)
double precision,allocatable :: TB(:,:)
@ -69,14 +49,40 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
double precision,allocatable :: XpY(:,:,:)
double precision,allocatable :: XmY(:,:,:)
double precision,allocatable :: Omega1s(:),Omega1t(:)
double precision,allocatable :: X1s(:,:),X1t(:,:)
double precision,allocatable :: Y1s(:,:),Y1t(:,:)
double precision,allocatable :: rho1s(:,:,:),rho1t(:,:,:)
double precision,allocatable :: Omega2s(:),Omega2t(:)
double precision,allocatable :: X2s(:,:),X2t(:,:)
double precision,allocatable :: Y2s(:,:),Y2t(:,:)
double precision,allocatable :: rho2s(:,:,:),rho2t(:,:,:)
! Output variables
double precision,intent(out) :: EcAC(nspin)
! Useful quantities
! nOOs = nO*nO
! nVVs = nV*nV
nOOs = nO*(nO+1)/2
nVVs = nV*(nV+1)/2
nOOt = nO*(nO-1)/2
nVVt = nV*(nV-1)/2
! Memory allocation
allocate(Ec(nAC,nspin))
allocate(Omega1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), &
Omega2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), &
rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs), &
Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), &
Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), &
rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt))
allocate(TA(nS,nS),TB(nS,nS),Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin))
allocate(Ec(nAC,nspin))
! Antisymmetrized kernel version
@ -115,31 +121,31 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
TA(:,:) = 0d0
TB(:,:) = 0d0
if(doXBS) then
! if(doXBS) then
isp_T = 1
iblock = 3
! isp_T = 1
! iblock = 3
call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,lambda,eT,ERI, &
Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(isp_T))
! call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,lambda,eT,ERI, &
! Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(isp_T))
call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TA)
if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TB)
! call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TA)
! if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TB)
isp_T = 2
iblock = 4
! isp_T = 2
! iblock = 4
call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,lambda,eT,ERI, &
Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(isp_T))
! call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,lambda,eT,ERI, &
! Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(isp_T))
call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TA)
if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TB)
! call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TA)
! if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TB)
end if
! end if
call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TA,TB, &
EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
@ -180,31 +186,36 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
lambda = rAC(iAC)
if(doXBS) then
! Initialize T matrix
isp_T = 1
iblock = 3
TA(:,:) = 0d0
TB(:,:) = 0d0
! if(doXBS) then
! isp_T = 1
! iblock = 3
! call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,lambda,eT,ERI, &
! Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(isp_T))
call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TA)
if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TB)
! call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TA)
! if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TB)
isp_T = 2
iblock = 4
! isp_T = 2
! iblock = 4
! call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,lambda,eT,ERI, &
! Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(isp_T))
call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TA)
if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TB)
! call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TA)
! if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TB)
end if
! end if
call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TA,TB, &
EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))

173
src/RPA/ACFDT_cr.f90 Normal file
View File

@ -0,0 +1,173 @@
subroutine ACFDT_cr(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eW,e,EcAC)
! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem
! for the crossed-ring contribution
implicit none
include 'parameters.h'
include 'quadrature.h'
! Input variables
logical,intent(in) :: doXBS
logical,intent(in) :: exchange_kernel
logical,intent(in) :: dRPA
logical,intent(in) :: TDA_W
logical,intent(in) :: TDA
logical,intent(in) :: BSE
logical,intent(in) :: singlet
logical,intent(in) :: triplet
double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: eW(nBas)
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: ispin
integer :: isp_W
integer :: iAC
double precision :: lambda
double precision,allocatable :: Ec(:,:)
double precision :: EcRPA
double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rho_RPA(:,:,:)
double precision,allocatable :: Omega(:,:)
double precision,allocatable :: XpY(:,:,:)
double precision,allocatable :: XmY(:,:,:)
! Output variables
double precision,intent(out) :: EcAC(nspin)
! Memory allocation
allocate(Ec(nAC,nspin))
allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS))
allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin))
! Antisymmetrized kernel version
if(exchange_kernel) then
write(*,*)
write(*,*) '*** Exchange kernel version ***'
write(*,*)
end if
EcAC(:) = 0d0
Ec(:,:) = 0d0
! Compute (singlet) RPA screening
isp_W = 1
EcRPA = 0d0
call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,OmRPA, &
rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
! Singlet manifold
if(singlet) then
ispin = 1
write(*,*) '--------------'
write(*,*) 'Singlet states'
write(*,*) '--------------'
write(*,*)
write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','Ec(lambda)','Tr(K x P_lambda)'
write(*,*) '-----------------------------------------------------------------------------------'
do iAC=1,nAC
lambda = -rAC(iAC)
if(doXBS) then
call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, &
rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
! call print_excitation('W^lambda: ',isp_W,nS,OmRPA)
end if
call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,OmRPA, &
rho_RPA,EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS, &
ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin)
end do
EcAC(ispin) = -0.5d0*dot_product(wAC,Ec(:,ispin))
if(exchange_kernel) EcAC(ispin) = 0.5d0*EcAC(ispin)
write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin)
write(*,*) '-----------------------------------------------------------------------------------'
write(*,*)
end if
! Triplet manifold
if(triplet) then
ispin = 2
write(*,*) '--------------'
write(*,*) 'Triplet states'
write(*,*) '--------------'
write(*,*)
write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','Ec(lambda)','Tr(K x P_lambda)'
write(*,*) '-----------------------------------------------------------------------------------'
do iAC=1,nAC
lambda = -rAC(iAC)
if(doXBS) then
call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, &
rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
end if
call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,OmRPA, &
rho_RPA,EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin)
end do
EcAC(ispin) = -0.5d0*dot_product(wAC,Ec(:,ispin))
if(exchange_kernel) EcAC(ispin) = 1.5d0*EcAC(ispin)
write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin)
write(*,*) '-----------------------------------------------------------------------------------'
write(*,*)
end if
end subroutine ACFDT_cr

145
src/RPA/ACFDT_pp.f90 Normal file
View File

@ -0,0 +1,145 @@
subroutine ACFDT_pp(TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC)
! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem for pp sector
implicit none
include 'parameters.h'
include 'quadrature.h'
! Input variables
logical,intent(in) :: TDA
logical,intent(in) :: singlet
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) :: e(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: ispin
integer :: iAC
double precision :: lambda
double precision,allocatable :: Ec(:,:)
integer :: nOOs,nOOt
integer :: nVVs,nVVt
double precision,allocatable :: Omega1s(:),Omega1t(:)
double precision,allocatable :: X1s(:,:),X1t(:,:)
double precision,allocatable :: Y1s(:,:),Y1t(:,:)
double precision,allocatable :: rho1s(:,:,:),rho1t(:,:,:)
double precision,allocatable :: Omega2s(:),Omega2t(:)
double precision,allocatable :: X2s(:,:),X2t(:,:)
double precision,allocatable :: Y2s(:,:),Y2t(:,:)
double precision,allocatable :: rho2s(:,:,:),rho2t(:,:,:)
! Output variables
double precision,intent(out) :: EcAC(nspin)
! Useful quantities
nOOs = nO*(nO+1)/2
nVVs = nV*(nV+1)/2
nOOt = nO*(nO-1)/2
nVVt = nV*(nV-1)/2
! Memory allocation
allocate(Omega1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), &
Omega2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), &
rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs), &
Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), &
Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), &
rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt))
allocate(Ec(nAC,nspin))
! Antisymmetrized kernel version
EcAC(:) = 0d0
Ec(:,:) = 0d0
! Singlet manifold
if(singlet) then
ispin = 1
write(*,*) '--------------'
write(*,*) 'Singlet states'
write(*,*) '--------------'
write(*,*)
write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','Ec(lambda)','Tr(K x P_lambda)'
write(*,*) '-----------------------------------------------------------------------------------'
do iAC=1,nAC
lambda = rAC(iAC)
call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOOs,nVVs,lambda,e,ERI,Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcAC(ispin))
call ACFDT_pp_correlation_energy(ispin,nBas,nC,nO,nV,nR,nS,ERI,nOOs,nVVs,X1s,Y1s,X2s,Y2s,Ec(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin)
end do
EcAC(ispin) = 0.5d0*dot_product(wAC,Ec(:,ispin))
write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin)
write(*,*) '-----------------------------------------------------------------------------------'
write(*,*)
end if
! Triplet manifold
if(triplet) then
ispin = 2
write(*,*) '--------------'
write(*,*) 'Triplet states'
write(*,*) '--------------'
write(*,*)
write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','Ec(lambda)','Tr(K x P_lambda)'
write(*,*) '-----------------------------------------------------------------------------------'
do iAC=1,nAC
lambda = rAC(iAC)
! Initialize T matrix
call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOOt,nVVt,lambda,e,ERI,Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcAC(ispin))
call ACFDT_pp_correlation_energy(ispin,nBas,nC,nO,nV,nR,nS,ERI,nOOt,nVVt,X1t,Y1t,X2t,Y2t,Ec(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin)
end do
EcAC(ispin) = 1.5d0*dot_product(wAC,Ec(:,ispin))
write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin)
write(*,*) '-----------------------------------------------------------------------------------'
write(*,*)
end if
end subroutine ACFDT_pp

51
src/RPA/ACFDT_pp_correlation_energy.f90 Normal file
View File

@ -0,0 +1,51 @@
subroutine ACFDT_pp_correlation_energy(ispin,nBas,nC,nO,nV,nR,nS,ERI,nOO,nVV,X1,Y1,X2,Y2,EcAC)
! Compute the correlation energy via the adiabatic connection formula for the pp sector
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: ispin
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
integer,intent(in) :: nOO
integer,intent(in) :: nVV
double precision,intent(in) :: X1(nVV,nVV)
double precision,intent(in) :: Y1(nOO,nVV)
double precision,intent(in) :: X2(nVV,nOO)
double precision,intent(in) :: Y2(nOO,nOO)
! Local variables
double precision,allocatable :: B(:,:)
double precision,allocatable :: C(:,:)
double precision,allocatable :: D(:,:)
double precision,external :: trace_matrix
! Output variables
double precision,intent(out) :: EcAC
! Memory allocation
allocate(B(nVV,nOO),C(nVV,nVV),D(nOO,nOO))
! Build pp matrices
call linear_response_B_pp (ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,B)
call linear_response_C_pp_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,C)
call linear_response_D_pp_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,D)
! Compute Tr(K x P_lambda)
EcAC = trace_matrix(nVV,matmul(transpose(X1),matmul(B,Y1)) + matmul(transpose(Y1),matmul(transpose(B),X1))) &
+ trace_matrix(nVV,matmul(transpose(X1),matmul(C,X1)) + matmul(transpose(Y1),matmul(D,Y1))) &
+ trace_matrix(nOO,matmul(transpose(X2),matmul(B,Y2)) + matmul(transpose(Y2),matmul(transpose(B),X2))) &
+ trace_matrix(nOO,matmul(transpose(X2),matmul(C,X2)) + matmul(transpose(Y2),matmul(D,Y2))) &
- trace_matrix(nVV,C) - trace_matrix(nOO,D)
end subroutine ACFDT_pp_correlation_energy

8
src/RPA/RPA.f90
View File

@ -87,12 +87,12 @@ subroutine RPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,
endif
if(exchange_kernel) then
! if(exchange_kernel) then
EcRPA(1) = 0.5d0*EcRPA(1)
EcRPA(2) = 1.5d0*EcRPA(2)
! EcRPA(1) = 0.5d0*EcRPA(1)
! EcRPA(2) = 1.5d0*EcRPA(2)
end if
! end if
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'

6
src/RPA/RPAx.f90
View File

@ -13,8 +13,8 @@ subroutine RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR
logical,intent(in) :: doACFDT
logical,intent(in) :: exchange_kernel
logical,intent(in) :: singlet
double precision,intent(in) :: eta
logical,intent(in) :: triplet
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
@ -89,12 +89,12 @@ subroutine RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR
endif
if(exchange_kernel) then
! if(exchange_kernel) then
EcRPAx(1) = 0.5d0*EcRPAx(1)
EcRPAx(2) = 1.5d0*EcRPAx(2)
end if
! end if
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'

56
src/RPA/ehRPA.f90 → src/RPA/crRPA.f90
View File

@ -1,7 +1,7 @@
subroutine ehRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, &
subroutine crRPA(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)
! Crossed-ring channel of the random phase approximation
implicit none
include 'parameters.h'
@ -42,7 +42,7 @@ subroutine ehRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n
write(*,*)
write(*,*)'***********************************************************'
write(*,*)'| Random phase approximation calculation: eh channel |'
write(*,*)'| Random phase approximation calculation: cr channel |'
write(*,*)'***********************************************************'
write(*,*)
@ -70,7 +70,7 @@ subroutine ehRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n
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_excitation('crRPA@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
@ -83,48 +83,48 @@ subroutine ehRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n
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_excitation('crRPA@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
! if(exchange_kernel) then
EcRPAx(1) = 0.5d0*EcRPAx(1)
EcRPAx(2) = 1.5d0*EcRPAx(2)
end if
! 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(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy (singlet) =',EcRPAx(1)
write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy (triplet) =',EcRPAx(2)
write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy =',EcRPAx(1) + EcRPAx(2)
write(*,'(2X,A50,F20.10)') 'Tr@crRPA total energy =',ENuc + ERHF + EcRPAx(1) + EcRPAx(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
! Compute the correlation energy via the adiabatic connection
! if(doACFDT) then
if(doACFDT) then
! write(*,*) '-------------------------------------------------------'
! write(*,*) 'Adiabatic connection version of ehRPA correlation energy'
! write(*,*) '-------------------------------------------------------'
! write(*,*)
write(*,*) '-------------------------------------------------------'
write(*,*) 'Adiabatic connection version of crRPA 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)
call ACFDT_cr(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(*,*)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'AC@crRPA correlation energy (singlet) =',EcAC(1)
write(*,'(2X,A50,F20.10)') 'AC@crRPA correlation energy (triplet) =',EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@crRPA correlation energy =',EcAC(1) + EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@crRPA total energy =',ENuc + ERHF + EcAC(1) + EcAC(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
! end if
end if
end subroutine ehRPA
end subroutine crRPA

108
src/RPA/ppRPA.f90
View File

@ -1,4 +1,4 @@
subroutine ppRPA(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e)
subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e)
! Perform pp-RPA calculation
@ -8,6 +8,7 @@ subroutine ppRPA(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e)
! Input variables
logical,intent(in) :: TDA
logical,intent(in) :: doACFDT
logical,intent(in) :: singlet
logical,intent(in) :: triplet
integer,intent(in) :: nBas
@ -23,16 +24,18 @@ subroutine ppRPA(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e)
! Local variables
integer :: ispin
integer :: nOO
integer :: nVV
double precision,allocatable :: Omega1(:,:)
double precision,allocatable :: X1(:,:,:)
double precision,allocatable :: Y1(:,:,:)
double precision,allocatable :: Omega2(:,:)
double precision,allocatable :: X2(:,:,:)
double precision,allocatable :: Y2(:,:,:)
integer :: nS
integer :: nOOs,nOOt
integer :: nVVs,nVVt
double precision,allocatable :: Omega1s(:),Omega1t(:)
double precision,allocatable :: X1s(:,:),X1t(:,:)
double precision,allocatable :: Y1s(:,:),Y1t(:,:)
double precision,allocatable :: Omega2s(:),Omega2t(:)
double precision,allocatable :: X2s(:,:),X2t(:,:)
double precision,allocatable :: Y2s(:,:),Y2t(:,:)
double precision :: Ec_ppRPA(nspin)
double precision :: EcAC(nspin)
! Hello world
@ -45,6 +48,25 @@ subroutine ppRPA(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e)
! Initialization
Ec_ppRPA(:) = 0d0
EcAC(:) = 0d0
! Useful quantities
nS = nO*nV
nOOs = nO*(nO+1)/2
nVVs = nV*(nV+1)/2
nOOt = nO*(nO-1)/2
nVVt = nV*(nV-1)/2
! Memory allocation
allocate(Omega1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), &
Omega2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs))
allocate(Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), &
Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt))
! Singlet manifold
@ -52,25 +74,11 @@ subroutine ppRPA(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e)
ispin = 1
! Useful quantities
call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,e,ERI, &
Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,Ec_ppRPA(ispin))
nOO = nO*(nO+1)/2
nVV = nV*(nV+1)/2
! Memory allocation
allocate(Omega1(nVV,nspin),X1(nVV,nVV,nspin),Y1(nOO,nVV,nspin), &
Omega2(nOO,nspin),X2(nVV,nOO,nspin),Y2(nOO,nOO,nspin))
call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,1d0,e,ERI, &
Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), &
Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), &
Ec_ppRPA(ispin))
call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:,ispin))
call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2(:,ispin))
deallocate(Omega1,X1,Y1,Omega2,X2,Y2)
call print_excitation('pp-RPA (N+2)',ispin,nVVs,Omega1s)
call print_excitation('pp-RPA (N-2)',ispin,nOOs,Omega2s)
endif
@ -80,26 +88,11 @@ subroutine ppRPA(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e)
ispin = 2
! Useful quantities
call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,e,ERI, &
Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,Ec_ppRPA(ispin))
nOO = nO*(nO-1)/2
nVV = nV*(nV-1)/2
! Memory allocation
allocate(Omega1(nVV,nspin),X1(nVV,nVV,nspin),Y1(nOO,nVV,nspin), &
Omega2(nOO,nspin),X2(nVV,nOO,nspin),Y2(nOO,nOO,nspin))
call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,1d0,e,ERI, &
Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), &
Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), &
Ec_ppRPA(ispin))
call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:,ispin))
call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2(:,ispin))
deallocate(Omega1,X1,Y1,Omega2,X2,Y2)
call print_excitation('pp-RPA (N+2)',ispin,nVVt,Omega1t)
call print_excitation('pp-RPA (N-2)',ispin,nOOt,Omega2t)
endif
@ -112,4 +105,27 @@ subroutine ppRPA(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
! Compute the correlation energy via the adiabatic connection
if(doACFDT) then
write(*,*) '---------------------------------------------------------'
write(*,*) 'Adiabatic connection version of pp-RPA correlation energy'
write(*,*) '---------------------------------------------------------'
write(*,*)
call ACFDT_pp(TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy (singlet) =',EcAC(1),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy (triplet) =',EcAC(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy =',EcAC(1) + EcAC(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA total energy =',ENuc + ERHF + EcAC(1) + EcAC(2),' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
end if
end subroutine ppRPA