PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-22 20:34:46 +01:00
Code Releases Activity

ACFDT for ppRPA

Browse Source
This commit is contained in:
Pierre-Francois Loos 2021-11-18 13:34:08 +01:00
parent 6f886a61a6
commit c379b80eac
18 changed files with 446 additions and 61 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
input/methods
View File

@ -9,7 +9,7 @@
# CIS* CIS(D) CID CISD FCI
F F F F F
# RPA* RPAx* crRPA ppRPA
F T T F
F F F T
# G0F2* evGF2* qsGF2* G0F3 evGF3
F F F F F
# G0W0* evGW* qsGW* ufG0W0 ufGW

2
input/options
View File

@ -11,7 +11,7 @@
# 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
T T T
T T F
# BSE: BSE dBSE dTDA evDyn
T F T F
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift

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

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

2
src/MBPT/G0T0.f90
View File

@ -245,7 +245,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
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

2
src/MBPT/evGT.f90
View File

@ -287,7 +287,7 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
if(exchange_kernel) then
EcAC(1) = 0.5d0*EcAC(1)
EcAC(2) = 1.5d0*EcAC(1)
EcAC(2) = 1.5d0*EcAC(2)
end if

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

4
src/QuAcK/QuAcK.f90
View File

@ -442,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
@ -824,7 +824,7 @@ program QuAcK
if(doppRPA) then
call cpu_time(start_RPA)
call ppRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF)
call ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF)
call cpu_time(end_RPA)
t_RPA = end_RPA - start_RPA

76
src/RPA/ACFDT_Tmatrix.f90
View File

@ -64,8 +64,11 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
! Useful quantities
nOOs = nO*nO
nVVs = nV*nV
! 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
@ -118,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))
@ -183,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
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))
! if(doXBS) then
call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
! isp_T = 1
! iblock = 3
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 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))
isp_T = 2
iblock = 4
! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
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 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 excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
! isp_T = 2
! iblock = 4
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 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))
end if
! 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)
! 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))

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

14
src/RPA/ppRPA.f90
View File

@ -1,4 +1,4 @@
subroutine ppRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,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
@ -9,10 +9,8 @@ subroutine ppRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n
logical,intent(in) :: TDA
logical,intent(in) :: doACFDT
logical,intent(in) :: exchange_kernel
logical,intent(in) :: singlet
logical,intent(in) :: triplet
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
@ -116,15 +114,7 @@ subroutine ppRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n
write(*,*) '---------------------------------------------------------'
write(*,*)
call ACFDT_Tmatrix(exchange_kernel,.false.,.false.,.false.,TDA,.false.,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, &
ERI,e,e,EcAC)
if(exchange_kernel) then
EcAC(1) = 0.5d0*EcAC(1)
EcAC(2) = 1.5d0*EcAC(1)
end if
call ACFDT_pp(TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'