LCPQ
/
quack
mirror of https://github.com/pfloos/quack
4
1
Fork 0
Code Releases Activity

working on ACFDT for T-matrix

This commit is contained in:
Pierre-Francois Loos 2021-11-10 22:41:40 +01:00
parent 3dfc2c3526
commit 23d7d83091
11 changed files with 151 additions and 101 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
input/methods
View File

@ -5,11 +5,11 @@
# CCD pCCD DCD CCSD CCSD(T)
F F F F F
# drCCD rCCD crCCD lCCD
T T T T
F F F F
# CIS* CIS(D) CID CISD FCI
F F F F F
# RPA* RPAx* crRPA ppRPA
T T T T
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
F F T
T T T
# BSE: BSE dBSE dTDA evDyn
T T T F
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift

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

9
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,8 +239,8 @@ 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

4
src/MBPT/evGT.f90
View File

@ -281,8 +281,8 @@ 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

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/QuAcK/QuAcK.f90
View File

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

61
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,37 @@ 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
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
@ -185,8 +188,8 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
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)
@ -196,8 +199,8 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
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)

2
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

32
src/RPA/crRPA.f90
View File

@ -106,25 +106,25 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n
! Compute the correlation energy via the adiabatic connection
! if(doACFDT) then
if(doACFDT) then
! write(*,*) '-------------------------------------------------------'
! write(*,*) 'Adiabatic connection version of crRPA 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(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 crRPA

118
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,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e)
! Perform pp-RPA calculation
@ -8,8 +8,11 @@ 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) :: 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
@ -23,16 +26,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 +50,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 +76,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 +90,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 +107,35 @@ 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_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
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