mirror of
https://github.com/pfloos/quack
synced 2024-06-18 19:25:31 +02:00
starting G0T0
This commit is contained in:
parent
3c848e9244
commit
6e3abd8927
|
|
@ -2,4 +2,4 @@
|
|||
2 2 2 0 0
|
||||
# Znuc x y z
|
||||
Li 0. 0. 0.
|
||||
H 0. 0. 3.015
|
||||
H 0. 0. 3.061356
|
||||
|
|
|
|||
|
|
@ -5,10 +5,12 @@
|
|||
# CCD CCSD CCSD(T)
|
||||
F F F
|
||||
# CIS TDHF ppRPA ADC
|
||||
F F T F
|
||||
F F F F
|
||||
# GF2 GF3
|
||||
F F
|
||||
# G0W0 evGW qsGW
|
||||
# G0W0 evGW qsGW
|
||||
F F F
|
||||
# G0T0 evGT qsGT
|
||||
T F F
|
||||
# MCMP2
|
||||
F
|
||||
|
|
|
|||
153
src/QuAcK/G0T0.f90
Normal file
153
src/QuAcK/G0T0.f90
Normal file
|
|
@ -0,0 +1,153 @@
|
|||
subroutine G0T0(BSE,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nOO,nVV,ENuc,ERHF,Hc,H,ERI,PHF,cHF,eHF,eG0T0)
|
||||
|
||||
! Perform G0W0 calculation with a T-matrix self-energy (G0T0)
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
logical,intent(in) :: BSE
|
||||
logical,intent(in) :: singlet_manifold
|
||||
logical,intent(in) :: triplet_manifold
|
||||
double precision,intent(in) :: eta
|
||||
|
||||
integer,intent(in) :: nBas,nC,nO,nV,nR,nOO,nVV
|
||||
double precision,intent(in) :: ENuc
|
||||
double precision,intent(in) :: ERHF
|
||||
double precision,intent(in) :: eHF(nBas)
|
||||
double precision,intent(in) :: cHF(nBas,nBas)
|
||||
double precision,intent(in) :: PHF(nBas,nBas)
|
||||
double precision,intent(in) :: Hc(nBas,nBas)
|
||||
double precision,intent(in) :: H(nBas,nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
logical :: dRPA
|
||||
integer :: ispin,jspin
|
||||
double precision :: EcRPA(nspin)
|
||||
double precision :: EcBSE(nspin)
|
||||
double precision :: EcGM
|
||||
double precision,allocatable :: SigT(:)
|
||||
double precision,allocatable :: Z(:)
|
||||
double precision,allocatable :: Omega1(:,:)
|
||||
double precision,allocatable :: X1(:,:,:)
|
||||
double precision,allocatable :: Y1(:,:,:)
|
||||
double precision,allocatable :: rho1(:,:,:,:)
|
||||
double precision,allocatable :: Omega2(:,:)
|
||||
double precision,allocatable :: X2(:,:,:)
|
||||
double precision,allocatable :: Y2(:,:,:)
|
||||
double precision,allocatable :: rho2(:,:,:,:)
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision :: eG0T0(nBas)
|
||||
|
||||
! Hello world
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'************************************************'
|
||||
write(*,*)'| One-shot G0T0 calculation |'
|
||||
write(*,*)'************************************************'
|
||||
write(*,*)
|
||||
|
||||
! Spin manifold
|
||||
|
||||
ispin = 1
|
||||
|
||||
! 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), &
|
||||
rho1(nBas,nBas,nVV,nspin),rho2(nBas,nBas,nOO,nspin), &
|
||||
SigT(nBas),Z(nBas))
|
||||
|
||||
! Compute linear response
|
||||
|
||||
call linear_response_pp(ispin,.false.,nBas,nC,nO,nV,nR,nOO,nVV,eHF,ERI, &
|
||||
Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), &
|
||||
Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), &
|
||||
EcRPA(ispin))
|
||||
|
||||
! Compute excitation densities for the T-matrix
|
||||
|
||||
call excitation_density_Tmatrix(nBas,nC,nO,nR,nOO,nVV,ERI, &
|
||||
X1(:,:,ispin),Y1(:,:,ispin),rho1(:,:,:,ispin), &
|
||||
X2(:,:,ispin),Y2(:,:,ispin),rho2(:,:,:,ispin))
|
||||
|
||||
! Compute T-matrix version of the self-energy
|
||||
|
||||
call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF, &
|
||||
Omega1(:,ispin),rho1(:,:,:,ispin), &
|
||||
Omega2(:,ispin),rho2(:,:,:,ispin), &
|
||||
SigT)
|
||||
|
||||
! Compute renormalization factor for T-matrix self-energy
|
||||
|
||||
call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF, &
|
||||
Omega1(:,ispin),rho1(:,:,:,ispin), &
|
||||
Omega2(:,ispin),rho2(:,:,:,ispin), &
|
||||
Z(:))
|
||||
|
||||
! Solve the quasi-particle equation
|
||||
|
||||
eG0T0(:) = eHF(:) + Z(:)*SigT(:)
|
||||
|
||||
! Dump results
|
||||
|
||||
call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:,ispin))
|
||||
call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2(:,ispin))
|
||||
|
||||
call print_G0T0(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eG0T0,EcRPA(ispin))
|
||||
|
||||
! Perform BSE calculation
|
||||
|
||||
! if(BSE) then
|
||||
|
||||
! ! Singlet manifold
|
||||
|
||||
! if(singlet_manifold) then
|
||||
|
||||
! ispin = 1
|
||||
! EcBSE(ispin) = 0d0
|
||||
|
||||
! call linear_response(ispin,.false.,.false.,.false.,nBas,nC,nO,nV,nR,nS,eHF,ERI, &
|
||||
! rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
|
||||
! call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
|
||||
|
||||
! call linear_response(ispin,.false.,.false.,BSE,nBas,nC,nO,nV,nR,nS,eG0T0,ERI, &
|
||||
! rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin))
|
||||
! call print_excitation('BSE ',ispin,nS,Omega(:,ispin))
|
||||
|
||||
! endif
|
||||
|
||||
! ! Triplet manifold
|
||||
|
||||
! if(triplet_manifold) then
|
||||
|
||||
! ispin = 2
|
||||
! EcBSE(ispin) = 0d0
|
||||
|
||||
! call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eHF,ERI, &
|
||||
! rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,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,eG0T0,ERI, &
|
||||
! rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin))
|
||||
! call print_excitation('BSE ',ispin,nS,Omega(:,ispin))
|
||||
|
||||
! endif
|
||||
|
||||
! write(*,*)
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,'(2X,A40,F15.6)') 'BSE@G0T0 correlation energy (singlet) =',EcBSE(1)
|
||||
! write(*,'(2X,A40,F15.6)') 'BSE@G0T0 correlation energy (triplet) =',EcBSE(2)
|
||||
! write(*,'(2X,A40,F15.6)') 'BSE@G0T0 correlation energy =',EcBSE(1) + EcBSE(2)
|
||||
! write(*,'(2X,A40,F15.6)') 'BSE@G0T0 total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2)
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,*)
|
||||
|
||||
! endif
|
||||
|
||||
end subroutine G0T0
|
||||
|
|
@ -10,18 +10,20 @@ program QuAcK
|
|||
logical :: doCIS,doTDHF,doppRPA,doADC
|
||||
logical :: doGF2,doGF3
|
||||
logical :: doG0W0,doevGW,doqsGW
|
||||
logical :: doG0T0,doevGT,doqsGT
|
||||
logical :: doMCMP2,doMinMCMP2
|
||||
logical :: doeNcusp
|
||||
logical :: doBas
|
||||
|
||||
integer :: nNuc,nBas,nBasCABS
|
||||
integer :: nEl(nspin),nC(nspin),nO(nspin),nV(nspin),nR(nspin),nS(nspin)
|
||||
integer :: nEl(nspin),nC(nspin),nO(nspin),nV(nspin),nR(nspin)
|
||||
integer :: nS(nspin),nOO(nspin),nVV(nspin)
|
||||
double precision :: ENuc,ERHF,EUHF,Norm
|
||||
double precision :: EcMP2(3),EcMP3,EcMP2F12(3),EcMCMP2(3),Err_EcMCMP2(3),Var_EcMCMP2(3)
|
||||
|
||||
double precision,allocatable :: ZNuc(:),rNuc(:,:)
|
||||
double precision,allocatable :: cHF(:,:,:),eHF(:,:),PHF(:,:,:)
|
||||
double precision,allocatable :: eG0W0(:)
|
||||
double precision,allocatable :: eG0T0(:)
|
||||
|
||||
integer :: nShell
|
||||
integer,allocatable :: TotAngMomShell(:),KShell(:)
|
||||
|
|
@ -50,7 +52,9 @@ program QuAcK
|
|||
double precision :: start_G0W0 ,end_G0W0 ,t_G0W0
|
||||
double precision :: start_evGW ,end_evGW ,t_evGW
|
||||
double precision :: start_qsGW ,end_qsGW ,t_qsGW
|
||||
double precision :: start_eNcusp ,end_eNcusp ,t_eNcusp
|
||||
double precision :: start_G0T0 ,end_G0T0 ,t_G0T0
|
||||
double precision :: start_evGT ,end_evGT ,t_evGT
|
||||
double precision :: start_qsGT ,end_qsGT ,t_qsGT
|
||||
double precision :: start_MP2 ,end_MP2 ,t_MP2
|
||||
double precision :: start_MP3 ,end_MP3 ,t_MP3
|
||||
double precision :: start_MP2F12 ,end_MP2F12 ,t_MP2F12
|
||||
|
|
@ -108,6 +112,7 @@ program QuAcK
|
|||
doCIS,doTDHF,doppRPA,doADC, &
|
||||
doGF2,doGF3, &
|
||||
doG0W0,doevGW,doqsGW, &
|
||||
doG0T0,doevGT,doqsGT, &
|
||||
doMCMP2)
|
||||
|
||||
! Read options for methods
|
||||
|
|
@ -121,7 +126,6 @@ program QuAcK
|
|||
|
||||
! Weird stuff
|
||||
|
||||
doeNCusp = .false.
|
||||
doMinMCMP2 = .false.
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
|
|
@ -175,7 +179,7 @@ program QuAcK
|
|||
|
||||
! Memory allocation for one- and two-electron integrals
|
||||
|
||||
allocate(cHF(nBas,nBas,nspin),eHF(nBas,nspin),eG0W0(nBas),PHF(nBas,nBas,nspin), &
|
||||
allocate(cHF(nBas,nBas,nspin),eHF(nBas,nspin),eG0W0(nBas),eG0T0(nBas),PHF(nBas,nBas,nspin), &
|
||||
S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),H(nBas,nBas),X(nBas,nBas), &
|
||||
ERI_AO_basis(nBas,nBas,nBas,nBas),ERI_MO_basis(nBas,nBas,nBas,nBas))
|
||||
|
||||
|
|
@ -522,16 +526,23 @@ program QuAcK
|
|||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Compute e-N cusp dressing
|
||||
! Perform G0T0 calculatiom
|
||||
!------------------------------------------------------------------------
|
||||
if(doeNcusp) then
|
||||
|
||||
call cpu_time(start_eNcusp)
|
||||
! call eNcusp()
|
||||
call cpu_time(end_eNcusp)
|
||||
eG0T0(:) = eHF(:,1)
|
||||
|
||||
t_eNcusp = end_eNcusp - start_eNcusp
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for e-N cusp dressing = ',t_eNcusp,' seconds'
|
||||
if(doG0T0) then
|
||||
|
||||
nOO(:) = nO(:)*(nO(:)+1)/2
|
||||
nVV(:) = nV(:)*(nV(:)+1)/2
|
||||
|
||||
call cpu_time(start_G0T0)
|
||||
call G0T0(BSE,singlet_manifold,triplet_manifold,eta, &
|
||||
nBas,nC(1),nO(1),nV(1),nR(1),nOO(1),nVV(1),ENuc,ERHF,Hc,H,ERI_MO_basis,PHF,cHF,eHF,eG0T0)
|
||||
call cpu_time(end_G0T0)
|
||||
|
||||
t_G0T0 = end_G0T0 - start_G0T0
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0 = ',t_G0T0,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
|
|
|||
|
|
@ -1,4 +1,4 @@
|
|||
subroutine excitation_density_Tmatrix(nBas,nC,nO,nR,nOO,nVV,ERI,X1,Y1,X2,Y2,rho1,rho2)
|
||||
subroutine excitation_density_Tmatrix(nBas,nC,nO,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2)
|
||||
|
||||
! Compute excitation densities for T-matrix self-energy
|
||||
|
||||
|
|
|
|||
49
src/QuAcK/print_G0T0.f90
Normal file
49
src/QuAcK/print_G0T0.f90
Normal file
|
|
@ -0,0 +1,49 @@
|
|||
subroutine print_G0T0(nBas,nO,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA)
|
||||
|
||||
! Print one-electron energies and other stuff for G0T0
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
integer,intent(in) :: nBas,nO
|
||||
double precision,intent(in) :: ENuc
|
||||
double precision,intent(in) :: EHF
|
||||
double precision,intent(in) :: EcRPA
|
||||
double precision,intent(in) :: e(nBas),SigmaC(nBas),Z(nBas),eGW(nBas)
|
||||
|
||||
integer :: x,HOMO,LUMO
|
||||
double precision :: Gap
|
||||
|
||||
! HOMO and LUMO
|
||||
|
||||
HOMO = nO
|
||||
LUMO = HOMO + 1
|
||||
Gap = eGW(LUMO)-eGW(HOMO)
|
||||
|
||||
! Dump results
|
||||
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,*)' One-shot G0T0 calculation (T-matrix self-energy) '
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
|
||||
'|','#','|','e_HF (eV)','|','Sigma_c (eV)','|','Z','|','e_QP (eV)','|'
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
|
||||
do x=1,nBas
|
||||
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
|
||||
'|',x,'|',e(x)*HaToeV,'|',SigmaC(x)*HaToeV,'|',Z(x),'|',eGW(x)*HaToeV,'|'
|
||||
enddo
|
||||
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,'(2X,A30,F15.6)') 'G0T0 HOMO energy (eV):',eGW(HOMO)*HaToeV
|
||||
write(*,'(2X,A30,F15.6)') 'G0T0 LUMO energy (eV):',eGW(LUMO)*HaToeV
|
||||
write(*,'(2X,A30,F15.6)') 'G0T0 HOMO-LUMO gap (eV):',Gap*HaToeV
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,'(2X,A30,F15.6)') 'RPA@G0T0 total energy =',ENuc + EHF + EcRPA
|
||||
write(*,'(2X,A30,F15.6)') 'RPA@G0T0 correlation energy =',EcRPA
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
end subroutine print_G0T0
|
||||
|
||||
|
||||
|
|
@ -4,6 +4,7 @@ subroutine read_methods(doRHF,doUHF,doMOM, &
|
|||
doCIS,doTDHF,doppRPA,doADC, &
|
||||
doGF2,doGF3, &
|
||||
doG0W0,doevGW,doqsGW, &
|
||||
doG0T0,doevGT,doqsGT, &
|
||||
doMCMP2)
|
||||
|
||||
! Read desired methods
|
||||
|
|
@ -18,6 +19,7 @@ subroutine read_methods(doRHF,doUHF,doMOM, &
|
|||
logical,intent(out) :: doCIS,doTDHF,doppRPA,doADC
|
||||
logical,intent(out) :: doGF2,doGF3
|
||||
logical,intent(out) :: doG0W0,doevGW,doqsGW
|
||||
logical,intent(out) :: doG0T0,doevGT,doqsGT
|
||||
logical,intent(out) :: doMCMP2
|
||||
|
||||
! Local variables
|
||||
|
|
@ -51,6 +53,10 @@ subroutine read_methods(doRHF,doUHF,doMOM, &
|
|||
doGF3 = .false.
|
||||
|
||||
doG0W0 = .false.
|
||||
doevGT = .false.
|
||||
doqsGT = .false.
|
||||
|
||||
doG0T0 = .false.
|
||||
doevGW = .false.
|
||||
doqsGW = .false.
|
||||
|
||||
|
|
@ -104,6 +110,14 @@ subroutine read_methods(doRHF,doUHF,doMOM, &
|
|||
if(answer2 == 'T') doevGW = .true.
|
||||
if(answer3 == 'T') doqsGW = .true.
|
||||
|
||||
! Read GT methods
|
||||
|
||||
read(1,*)
|
||||
read(1,*) answer1,answer2,answer3
|
||||
if(answer1 == 'T') doG0T0 = .true.
|
||||
if(answer2 == 'T') doevGT = .true.
|
||||
if(answer3 == 'T') doqsGT = .true.
|
||||
|
||||
! Read stochastic methods
|
||||
|
||||
read(1,*)
|
||||
|
|
|
|||
|
|
@ -1,4 +1,4 @@
|
|||
subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,Omega2,rho1,rho2,Z)
|
||||
subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,Z)
|
||||
|
||||
! Compute renormalization factor of the T-matrix self-energy
|
||||
|
||||
|
|
|
|||
|
|
@ -1,4 +1,4 @@
|
|||
subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,Omega2,rho1,rho2,SigT)
|
||||
subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,SigT)
|
||||
|
||||
! Compute diagonal of the correlation part of the T-matrix self-energy
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue
Block a user