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

GW AC

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-01-14 16:36:11 +01:00
parent 51016fd5c1
commit 3ebcd29b19
9 changed files with 176 additions and 350 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
input/basis
View File

@ -1,9 +1,16 @@
1 3
S 3 1.00
38.3600000 0.0238090
5.7700000 0.1548910
1.2400000 0.4699870
1 6
S 4 1.00
234.0000000 0.0025870
35.1600000 0.0195330
7.9890000 0.0909980
2.2120000 0.2720500
S 1 1.00
0.2976000 1.0000000
0.6669000 1.0000000
S 1 1.00
0.2089000 1.0000000
P 1 1.00
1.2750000 1.0000000
3.0440000 1.0000000
P 1 1.00
0.7580000 1.0000000
D 1 1.00
1.9650000 1.0000000

4
input/methods
View File

@ -3,9 +3,9 @@
# MP2 MP3 MP2-F12
F F F
# CCD CCSD CCSD(T) ringCCD ladderCCD
F F F T T
F F F F F
# CIS RPA TDHF ppRPA ADC
F T T T F
F T T F F
# GF2 GF3
F F
# G0W0 evGW qsGW

21
input/weight
View File

@ -1,9 +1,16 @@
1 3
S 3 1.00
38.3600000 0.0238090
5.7700000 0.1548910
1.2400000 0.4699870
1 6
S 4 1.00
234.0000000 0.0025870
35.1600000 0.0195330
7.9890000 0.0909980
2.2120000 0.2720500
S 1 1.00
0.2976000 1.0000000
0.6669000 1.0000000
S 1 1.00
0.2089000 1.0000000
P 1 1.00
1.2750000 1.0000000
3.0440000 1.0000000
P 1 1.00
0.7580000 1.0000000
D 1 1.00
1.9650000 1.0000000

114
src/QuAcK/G0W0.f90
View File

@ -29,8 +29,7 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
! Local variables
logical :: dRPA
integer :: ispin,jspin
integer :: ispin
double precision :: EcRPA(nspin)
double precision :: EcBSE(nspin)
double precision :: EcGM
@ -42,11 +41,8 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
double precision,allocatable :: rho(:,:,:,:)
double precision,allocatable :: rhox(:,:,:,:)
logical :: AC
integer :: iAC
double precision :: lambda
double precision,allocatable :: EcACBSE(:,:)
double precision,allocatable :: EcAC(:,:)
logical :: adiabatic_connection
logical :: scaled_screening
! Output variables
@ -70,10 +66,6 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
if(COHSEX) write(*,*) 'COHSEX approximation activated!'
write(*,*)
! Switch off exchange for G0W0
dRPA = .true.
! Spin manifold
ispin = 1
@ -83,12 +75,9 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
allocate(SigC(nBas),Z(nBas),Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin), &
rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin))
AC = .true.
allocate(EcACBSE(nAC,nspin),EcAC(nAC,nspin))
! Compute linear response
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, &
call linear_response(ispin,.true.,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
! Compute correlation part of the self-energy
@ -129,11 +118,11 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
ispin = 1
EcBSE(ispin) = 0d0
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, &
call linear_response(ispin,.true.,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eG0W0,ERI, &
call linear_response(ispin,.true.,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eG0W0,ERI, &
rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('BSE ',ispin,nS,Omega(:,ispin))
@ -146,11 +135,11 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
ispin = 2
EcBSE(ispin) = 0d0
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, &
call linear_response(ispin,.true.,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eG0W0,ERI, &
call linear_response(ispin,.true.,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eG0W0,ERI, &
rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('BSE ',ispin,nS,Omega(:,ispin))
@ -158,98 +147,35 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A40,F15.6)') 'BSE@G0W0 correlation energy (singlet) =',EcBSE(1)
write(*,'(2X,A40,F15.6)') 'BSE@G0W0 correlation energy (triplet) =',EcBSE(2)
write(*,'(2X,A40,F15.6)') 'BSE@G0W0 correlation energy =',EcBSE(1) + EcBSE(2)
write(*,'(2X,A40,F15.6)') 'BSE@G0W0 total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2)
write(*,'(2X,A40,F15.6)') 'Tr@BSE@G0W0 correlation energy (singlet) =',EcBSE(1)
write(*,'(2X,A40,F15.6)') 'Tr@BSE@G0W0 correlation energy (triplet) =',EcBSE(2)
write(*,'(2X,A40,F15.6)') 'Tr@BSE@G0W0 correlation energy =',EcBSE(1) + EcBSE(2)
write(*,'(2X,A40,F15.6)') 'Tr@BSE@G0W0 total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
! Compute the BSE correlation energy via the adiabatic connection
if(AC) then
adiabatic_connection = .true.
scaled_screening = .true.
if(adiabatic_connection) then
write(*,*) '------------------------------------------------------'
write(*,*) 'Adiabatic connection version of BSE correlation energy'
write(*,*) '------------------------------------------------------'
write(*,*)
if(singlet_manifold) then
ispin = 1
EcACBSE(:,ispin) = 0d0
if(scaled_screening) then
write(*,*) '--------------'
write(*,*) 'Singlet states'
write(*,*) '--------------'
write(*,*)
write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','EcBSE(lambda)','Tr(V x P_lambda)'
write(*,*) '-----------------------------------------------------------------------------------'
do iAC=1,nAC
lambda = rAC(iAC)
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,eG0W0,ERI, &
rho(:,:,:,ispin),EcACBSE(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),EcAC(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACBSE(iAC,ispin),EcAC(iAC,ispin)
end do
write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A50,1X,F15.6)') ' Ec(BSE) via Gauss-Legendre quadrature:',0.5d0*dot_product(wAC,EcAC(:,ispin))
write(*,*) '-----------------------------------------------------------------------------------'
write(*,*) '*** scaled screening version (extended BSE) ***'
write(*,*)
end if
if(triplet_manifold) then
ispin = 2
EcACBSE(:,ispin) = 0d0
call ACDFT(scaled_screening,.true.,TDA,BSE,singlet_manifold,triplet_manifold, &
nBas,nC,nO,nV,nR,nS,ERI,eG0W0,Omega,XpY,XmY,rho)
write(*,*) '--------------'
write(*,*) 'Triplet states'
write(*,*) '--------------'
write(*,*)
write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','EcBSE(lambda)','Tr(V x P_lambda)'
write(*,*) '-----------------------------------------------------------------------------------'
do iAC=1,nAC
lambda = rAC(iAC)
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,eG0W0,ERI, &
rho(:,:,:,ispin),EcACBSE(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),EcAC(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACBSE(iAC,ispin),EcAC(iAC,ispin)
end do
write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A50,1X,F15.6)') ' Ec(BSE) via Gauss-Legendre quadrature:',0.5d0*dot_product(wAC,EcAC(:,ispin))
write(*,*) '-----------------------------------------------------------------------------------'
write(*,*)
end if
end if
end if

121
src/QuAcK/RPA.f90
View File

@ -23,9 +23,6 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E
! Local variables
logical :: dRPA
logical :: TDA
logical :: BSE
integer :: ispin
double precision,allocatable :: Omega(:,:)
double precision,allocatable :: XpY(:,:,:)
@ -34,11 +31,7 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E
double precision :: rho
double precision :: EcRPA(nspin)
logical :: AC
integer :: iAC
double precision :: lambda
double precision,allocatable :: EcACRPA(:,:)
double precision,allocatable :: EcAC(:,:)
logical :: adiabatic_connection
! Hello world
@ -52,32 +45,17 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E
EcRPA(:) = 0d0
! Switch off exchange for RPA
dRPA = .true.
! Switch off Tamm-Dancoff approximation for RPA
TDA = .false.
! Switch off Bethe-Salpeter equation for RPA
BSE = .false.
! Memory allocation
allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin))
AC = .true.
allocate(EcACRPA(nAC,nspin),EcAC(nAC,nspin))
! Singlet manifold
if(singlet_manifold) then
ispin = 1
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
call linear_response(ispin,.true.,.false.,.false.,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
@ -89,7 +67,7 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E
ispin = 2
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
call linear_response(ispin,.true.,.false.,.false.,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
@ -97,90 +75,27 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A40,F15.6)') 'RPA@RPA correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A40,F15.6)') 'RPA@RPA correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A40,F15.6)') 'RPA@RPA correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A40,F15.6)') 'RPA@RPA total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
write(*,'(2X,A40,F15.6)') 'Tr@RPA correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A40,F15.6)') 'Tr@RPA correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A40,F15.6)') 'Tr@RPA correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A40,F15.6)') 'Tr@RPA total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
! Compute the correlation energy via the adiabatic connection
! Compute the correlation energy via the adiabatic connection
if(AC) then
adiabatic_connection = .true.
write(*,*) '------------------------------------------------------'
write(*,*) 'Adiabatic connection version of RPA correlation energy'
write(*,*) '------------------------------------------------------'
write(*,*)
if(adiabatic_connection) then
write(*,*) '------------------------------------------------------'
write(*,*) 'Adiabatic connection version of RPA correlation energy'
write(*,*) '------------------------------------------------------'
write(*,*)
if(singlet_manifold) then
call ACDFT(.false.,.true.,.false.,.false.,singlet_manifold,triplet_manifold, &
nBas,nC,nO,nV,nR,nS,ERI,e,Omega,XpY,XmY,rho)
ispin = 1
EcACRPA(:,ispin) = 0d0
write(*,*) '--------------'
write(*,*) 'Singlet states'
write(*,*) '--------------'
write(*,*)
write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','EcRPA(lambda)','Tr(V x P_lambda)'
write(*,*) '-----------------------------------------------------------------------------------'
do iAC=1,nAC
lambda = rAC(iAC)
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,rho, &
EcACRPA(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),EcAC(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACRPA(iAC,ispin),EcAC(iAC,ispin)
end do
write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A50,1X,F15.6)') ' Ec(RPA) via Gauss-Legendre quadrature:',0.5d0*dot_product(wAC,EcAC(:,ispin))
write(*,*) '-----------------------------------------------------------------------------------'
write(*,*)
end if
if(triplet_manifold) then
ispin = 2
EcACRPA(:,ispin) = 0d0
write(*,*) '--------------'
write(*,*) 'Triplet states'
write(*,*) '--------------'
write(*,*)
write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','EcRPA(lambda)','Tr(V x P_lambda)'
write(*,*) '-----------------------------------------------------------------------------------'
do iAC=1,nAC
lambda = rAC(iAC)
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,rho, &
EcACRPA(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),EcAC(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACRPA(iAC,ispin),EcAC(iAC,ispin)
end do
write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A50,1X,F15.6)') ' Ec(RPA) via Gauss-Legendre quadrature:',0.5d0*dot_product(wAC,EcAC(:,ispin))
write(*,*) '-----------------------------------------------------------------------------------'
write(*,*)
end if
end if
end if
end subroutine RPA

119
src/QuAcK/TDHF.f90
View File

@ -23,9 +23,6 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,
! Local variables
logical :: dRPA
logical :: TDA
logical :: BSE
integer :: ispin
double precision,allocatable :: Omega(:,:)
double precision,allocatable :: XpY(:,:,:)
@ -34,11 +31,7 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,
double precision :: rho
double precision :: EcRPA(nspin)
logical :: AC
integer :: iAC
double precision :: lambda
double precision,allocatable :: EcACRPA(:,:)
double precision,allocatable :: EcAC(:,:)
logical :: adiabatic_connection
! Hello world
@ -52,32 +45,17 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,
EcRPA(:) = 0d0
! Switch on exchange for TDHF
dRPA = .false.
! Switch off Tamm-Dancoff approximation for TDHF
TDA = .false.
! Switch off Bethe-Salpeter equation for TDHF
BSE = .false.
! Memory allocation
allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin))
AC = .true.
allocate(EcACRPA(nAC,nspin),EcAC(nAC,nspin))
! Singlet manifold
if(singlet_manifold) then
ispin = 1
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
call linear_response(ispin,.false.,.false.,.false.,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('TDHF ',ispin,nS,Omega(:,ispin))
@ -89,7 +67,7 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,
ispin = 2
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
call linear_response(ispin,.false.,.false.,.false.,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('TDHF ',ispin,nS,Omega(:,ispin))
@ -97,90 +75,27 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A40,F15.6)') 'RPA@TDHF correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A40,F15.6)') 'RPA@TDHF correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A40,F15.6)') 'RPA@TDHF correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A40,F15.6)') 'RPA@TDHF total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
write(*,'(2X,A40,F15.6)') 'Tr@TDHF correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A40,F15.6)') 'Tr@TDHF correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A40,F15.6)') 'Tr@TDHF correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A40,F15.6)') 'Tr@TDHF total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
! Compute the correlation energy via the adiabatic connection
! Compute the correlation energy via the adiabatic connection
if(AC) then
adiabatic_connection = .true.
write(*,*) '------------------------------------------------------'
write(*,*) 'Adiabatic connection version of RPA correlation energy'
write(*,*) '------------------------------------------------------'
write(*,*)
if(singlet_manifold) then
if(adiabatic_connection) then
ispin = 1
EcACRPA(:,ispin) = 0d0
write(*,*) '-------------------------------------------------------'
write(*,*) 'Adiabatic connection version of TDHF correlation energy'
write(*,*) '-------------------------------------------------------'
write(*,*)
write(*,*) '--------------'
write(*,*) 'Singlet states'
write(*,*) '--------------'
write(*,*)
call ACDFT(.false.,.false.,.false.,.false.,singlet_manifold,triplet_manifold, &
nBas,nC,nO,nV,nR,nS,ERI,e,Omega,XpY,XmY,rho)
write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','EcRPA(lambda)','Tr(V x P_lambda)'
write(*,*) '-----------------------------------------------------------------------------------'
do iAC=1,nAC
lambda = rAC(iAC)
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
rho,EcACRPA(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),EcAC(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACRPA(iAC,ispin),EcAC(iAC,ispin)
end do
write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A50,1X,F15.6)') ' Ec(RPA) via Gauss-Legendre quadrature:',0.5d0*dot_product(wAC,EcAC(:,ispin))
write(*,*) '-----------------------------------------------------------------------------------'
write(*,*)
end if
if(triplet_manifold) then
ispin = 2
EcACRPA(:,ispin) = 0d0
write(*,*) '--------------'
write(*,*) 'Triplet states'
write(*,*) '--------------'
write(*,*)
write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','EcRPA(lambda)','Tr(V x P_lambda)'
write(*,*) '-----------------------------------------------------------------------------------'
do iAC=1,nAC
lambda = rAC(iAC)
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
rho,EcACRPA(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),EcAC(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACRPA(iAC,ispin),EcAC(iAC,ispin)
end do
write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A50,1X,F15.6)') ' Ec(RPA) via Gauss-Legendre quadrature:',0.5d0*dot_product(wAC,EcAC(:,ispin))
write(*,*) '-----------------------------------------------------------------------------------'
write(*,*)
end if
end if
end if
end subroutine TDHF

48
src/QuAcK/evGW.f90
View File

@ -53,9 +53,13 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
double precision,allocatable :: SigC(:)
double precision,allocatable :: Omega(:,:)
double precision,allocatable :: XpY(:,:,:)
double precision,allocatable :: XmY(:,:,:)
double precision,allocatable :: rho(:,:,:,:)
double precision,allocatable :: rhox(:,:,:,:)
logical :: adiabatic_connection
logical :: scaled_screening
! Hello world
write(*,*)
@ -86,7 +90,7 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
! Memory allocation
allocate(eGW(nBas),eOld(nBas),Z(nBas),SigC(nBas),Omega(nS,nspin), &
XpY(nS,nS,nspin),rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin), &
XpY(nS,nS,nspin),XmY(nS,nS,nspin),rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin), &
error_diis(nBas,max_diis),e_diis(nBas,max_diis))
! Initialization
@ -112,7 +116,7 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
if(.not. GW0 .or. nSCF == 0) then
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eGW,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
endif
@ -221,11 +225,11 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
EcBSE(ispin) = 0d0
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin))
rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('BSE ',ispin,nS,Omega(:,ispin))
endif
@ -238,24 +242,48 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
EcBSE(ispin) = 0d0
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin))
rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('BSE ',ispin,nS,Omega(:,ispin))
endif
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A40,F15.6)') 'BSE@evGW correlation energy (singlet) =',EcBSE(1)
write(*,'(2X,A40,F15.6)') 'BSE@evGW correlation energy (triplet) =',EcBSE(2)
write(*,'(2X,A40,F15.6)') 'BSE@evGW correlation energy =',EcBSE(1) + EcBSE(2)
write(*,'(2X,A40,F15.6)') 'BSE@evGW total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2)
write(*,'(2X,A40,F15.6)') 'Tr@BSE@evGW correlation energy (singlet) =',EcBSE(1)
write(*,'(2X,A40,F15.6)') 'Tr@BSE@evGW correlation energy (triplet) =',EcBSE(2)
write(*,'(2X,A40,F15.6)') 'Tr@BSE@evGW correlation energy =',EcBSE(1) + EcBSE(2)
write(*,'(2X,A40,F15.6)') 'Tr@BSE@evGW total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
! Compute the BSE correlation energy via the adiabatic connection
adiabatic_connection = .true.
scaled_screening = .true.
if(adiabatic_connection) then
write(*,*) '------------------------------------------------------'
write(*,*) 'Adiabatic connection version of BSE correlation energy'
write(*,*) '------------------------------------------------------'
write(*,*)
if(scaled_screening) then
write(*,*) '*** scaled screening version (extended BSE) ***'
write(*,*)
end if
call ACDFT(scaled_screening,.true.,TDA,BSE,singlet_manifold,triplet_manifold, &
nBas,nC,nO,nV,nR,nS,ERI,eGW,Omega,XpY,XmY,rho)
end if
endif
end subroutine evGW

2
src/QuAcK/ladder_CCD.f90
View File

@ -187,7 +187,7 @@ subroutine ladder_CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
write(*,*)' ladder-CCD energy '
write(*,*)'----------------------------------------------------'
write(*,'(1X,A30,1X,F15.10)')' E(ladder-CCD) = ',ECCD
write(*,'(1X,A30,1X,F10.6)')' Ec(ladder-CCSD) = ',EcCCD
write(*,'(1X,A30,1X,F15.10)')' Ec(ladder-CCSD) = ',EcCCD
write(*,*)'----------------------------------------------------'
write(*,*)

76
src/QuAcK/qsGW.f90
View File

@ -52,11 +52,12 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
double precision,allocatable :: F_diis(:,:)
double precision,allocatable :: Omega(:,:)
double precision,allocatable :: XpY(:,:,:)
double precision,allocatable :: XmY(:,:,:)
double precision,allocatable :: rho(:,:,:,:)
double precision,allocatable :: rhox(:,:,:,:)
double precision,allocatable :: c(:,:)
double precision,allocatable :: cp(:,:)
double precision,allocatable :: e(:)
double precision,allocatable :: eGW(:)
double precision,allocatable :: P(:,:)
double precision,allocatable :: F(:,:)
double precision,allocatable :: Fp(:,:)
@ -68,6 +69,9 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
double precision,allocatable :: Z(:)
double precision,allocatable :: error(:,:)
logical :: adiabatic_connection
logical :: scaled_screening
! Hello world
write(*,*)
@ -101,9 +105,9 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
! Memory allocation
allocate(e(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), &
allocate(eGW(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), &
J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),SigCm(nBas,nBas),Z(nBas), &
Omega(nS,nspin),XpY(nS,nS,nspin),rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin), &
Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin),rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin), &
error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
! Initialization
@ -113,7 +117,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
ispin = 1
Conv = 1d0
P(:,:) = PHF(:,:)
e(:) = eHF(:)
eGW(:) = eHF(:)
c(:,:) = cHF(:,:)
F_diis(:,:) = 0d0
error_diis(:,:) = 0d0
@ -140,8 +144,8 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
if(.not. GW0 .or. nSCF == 0) then
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,e,ERI_MO_basis, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO_basis, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
endif
@ -160,9 +164,9 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
else
call self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e, &
call self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eGW, &
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC)
call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e, &
call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eGW, &
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z)
endif
@ -196,7 +200,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
Fp = matmul(transpose(X),matmul(F,X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,e)
call diagonalize_matrix(nBas,cp,eGW)
c = matmul(X,cp)
! Compute new density matrix in the AO basis
@ -206,7 +210,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
! Print results
call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,e,c,ENuc,P,T,V,Hc,J,K,F,SigCp,Z,EcRPA(ispin),EcGM,EqsGW)
call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,SigCp,Z,EcRPA(ispin),EcGM,EqsGW)
! Increment
@ -219,7 +223,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
! Compute second-order correction of the Hermitization error
call qsGW_PT(nBas,nC,nO,nV,nR,nS,e,SigCm)
call qsGW_PT(nBas,nC,nO,nV,nR,nS,eGW,SigCm)
! Compute the overlap between HF and GW orbitals
@ -253,12 +257,12 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
ispin = 1
EcBSE(ispin) = 0d0
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,e,ERI_MO_basis, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO_basis, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,e,ERI_MO_basis, &
rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin))
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO_basis, &
rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('BSE ',ispin,nS,Omega(:,ispin))
endif
@ -269,25 +273,49 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
ispin = 2
EcBSE(ispin) = 0d0
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,e,ERI_MO_basis, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO_basis, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,e,ERI_MO_basis, &
rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin))
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO_basis, &
rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('BSE ',ispin,nS,Omega(:,ispin))
endif
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A40,F15.6)') 'BSE@qsGW correlation energy (singlet) =',EcBSE(1)
write(*,'(2X,A40,F15.6)') 'BSE@qsGW correlation energy (triplet) =',EcBSE(2)
write(*,'(2X,A40,F15.6)') 'BSE@qsGW correlation energy =',EcBSE(1) + EcBSE(2)
write(*,'(2X,A40,F15.6)') 'BSE@qsGW total energy =',ENuc + EqsGW + EcBSE(1) + EcBSE(2)
write(*,'(2X,A40,F15.6)') 'Tr@BSE@qsGW correlation energy (singlet) =',EcBSE(1)
write(*,'(2X,A40,F15.6)') 'Tr@BSE@qsGW correlation energy (triplet) =',EcBSE(2)
write(*,'(2X,A40,F15.6)') 'Tr@BSE@qsGW correlation energy =',EcBSE(1) + EcBSE(2)
write(*,'(2X,A40,F15.6)') 'Tr@BSE@qsGW total energy =',ENuc + EqsGW + EcBSE(1) + EcBSE(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
endif
! Compute the BSE correlation energy via the adiabatic connection
adiabatic_connection = .true.
scaled_screening = .true.
if(adiabatic_connection) then
write(*,*) '------------------------------------------------------'
write(*,*) 'Adiabatic connection version of BSE correlation energy'
write(*,*) '------------------------------------------------------'
write(*,*)
if(scaled_screening) then
write(*,*) '*** scaled screening version (extended BSE) ***'
write(*,*)
end if
call ACDFT(scaled_screening,.true.,TDA,BSE,singlet_manifold,triplet_manifold, &
nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,Omega,XpY,XmY,rho)
end if
end if
end subroutine qsGW