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

clean up GM again

Browse Source
This commit is contained in:
Pierre-Francois Loos 2021-03-06 23:08:43 +01:00
parent 96fa82931d
commit 7fd5abe76f
14 changed files with 148 additions and 103 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
input/methods
View File

@ -1,7 +1,7 @@
# RHF UHF KS MOM
T F F F
# MP2* MP3 MP2-F12
F F F
T F F
# CCD DCD CCSD CCSD(T)
F F F F
# drCCD rCCD lCCD pCCD
@ -11,9 +11,9 @@
# RPA* RPAx* ppRPA
F F F
# G0F2 evGF2 qsGF2 G0F3 evGF3
F F T F F
F F F F F
# G0W0* evGW* qsGW*
F F F
F F T
# G0T0 evGT qsGT
F F F
# MCMP2

4
input/options
View File

@ -7,12 +7,12 @@
# spin: TDA singlet triplet spin_conserved spin_flip
F T T T T
# GF: maxSCF thresh DIIS n_diis lin eta renorm
256 0.00001 T 5 T 0.0 3
256 0.00001 T 5 T 0.001 3
# GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0
256 0.00001 T 5 T 0.0 F F F F F
# ACFDT: AC Kx XBS
F F T
# BSE: BSE dBSE dTDA evDyn
F T T F
F F T F
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift
1000000 100000 10 0.3 10000 1234 T

5
src/MBPT/BSE2.f90
View File

@ -56,6 +56,9 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,
call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGF,ERI, &
OmBSE(:,ispin),rho,EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('BSE2 ',ispin,nS,OmBSE(:,ispin))
call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int, &
OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
! Compute dynamic correction for BSE via perturbation theory
@ -90,6 +93,8 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,
call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGF(:),ERI(:,:,:,:), &
OmBSE(:,ispin),rho,EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('BSE2 ',ispin,nS,OmBSE(:,ispin))
call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int, &
OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
! Compute dynamic correction for BSE via perturbation theory

5
src/MBPT/G0F2.f90
View File

@ -30,6 +30,7 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,nBas,nC,nO
! Local variables
double precision :: Ec
double precision :: EcBSE(nspin)
double precision,allocatable :: eGF2(:)
double precision,allocatable :: SigC(:)
@ -56,7 +57,7 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,nBas,nC,nO
! Frequency-dependent second-order contribution
call self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eHF,ERI,SigC,Z)
call self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eHF,ERI,SigC,Z,Ec)
if(linearize) then
@ -70,7 +71,7 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,nBas,nC,nO
! Print results
call print_G0F2(nBas,nO,eHF,SigC,eGF2,Z)
call print_G0F2(nBas,nO,eHF,SigC,eGF2,Z,ENuc,ERHF,Ec)
! Perform BSE2 calculation

5
src/MBPT/evGF2.f90
View File

@ -36,6 +36,7 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet,triplet,
integer :: nSCF
integer :: n_diis
double precision :: Ec
double precision :: EcBSE(nspin)
double precision :: Conv
double precision :: rcond
@ -76,7 +77,7 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet,triplet,
! Frequency-dependent second-order contribution
call self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z)
call self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z,Ec)
if(linearize) then
@ -92,7 +93,7 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet,triplet,
! Print results
call print_evGF2(nBas,nO,nSCF,Conv,eHF,SigC,Z,eGF2)
call print_evGF2(nBas,nO,nSCF,Conv,eHF,SigC,Z,eGF2,ENuc,ERHF,Ec)
! DIIS extrapolation

14
src/MBPT/print_G0F2.f90
View File

@ -1,4 +1,4 @@
subroutine print_G0F2(nBas,nO,eHF,Sig,eGF2,Z)
subroutine print_G0F2(nBas,nO,eHF,Sig,eGF2,Z,ENuc,ERHF,Ec)
! Print one-electron energies and other stuff for G0F2
@ -11,6 +11,9 @@ subroutine print_G0F2(nBas,nO,eHF,Sig,eGF2,Z)
double precision,intent(in) :: Sig(nBas)
double precision,intent(in) :: eGF2(nBas)
double precision,intent(in) :: Z(nBas)
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: Ec
integer :: p
integer :: HOMO
@ -38,10 +41,13 @@ subroutine print_G0F2(nBas,nO,eHF,Sig,eGF2,Z)
enddo
write(*,*)'--------------------------------------------------------------------------'
write(*,'(2X,A27,F15.6)') 'G0F2 HOMO energy (eV):',eGF2(HOMO)*HaToeV
write(*,'(2X,A27,F15.6)') 'G0F2 LUMO energy (eV):',eGF2(LUMO)*HaToeV
write(*,'(2X,A27,F15.6)') 'G0F2 HOMO-LUMO gap (eV):',Gap*HaToeV
write(*,'(2X,A30,F15.6)') 'G0F2 HOMO energy (eV):',eGF2(HOMO)*HaToeV
write(*,'(2X,A30,F15.6)') 'G0F2 LUMO energy (eV):',eGF2(LUMO)*HaToeV
write(*,'(2X,A30,F15.6)') 'G0F2 HOMO-LUMO gap (eV):',Gap*HaToeV
write(*,*)'--------------------------------------------------------------------------'
write(*,'(2X,A30,F15.6,A3)') 'G0F2 total energy :',ENuc + ERHF + Ec,' au'
write(*,'(2X,A30,F15.6,A3)') 'G0F2 correlation energy:',Ec,' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
end subroutine print_G0F2

14
src/MBPT/print_evGF2.f90
View File

@ -1,4 +1,4 @@
subroutine print_evGF2(nBas,nO,nSCF,Conv,eHF,Sig,Z,eGF2)
subroutine print_evGF2(nBas,nO,nSCF,Conv,eHF,Sig,Z,eGF2,ENuc,ERHF,Ec)
! Print one-electron energies and other stuff for G0F2
@ -13,6 +13,9 @@ subroutine print_evGF2(nBas,nO,nSCF,Conv,eHF,Sig,Z,eGF2)
double precision,intent(in) :: Sig(nBas)
double precision,intent(in) :: eGF2(nBas)
double precision,intent(in) :: Z(nBas)
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: Ec
integer :: p
integer :: HOMO
@ -44,10 +47,13 @@ subroutine print_evGF2(nBas,nO,nSCF,Conv,eHF,Sig,Z,eGF2)
write(*,'(2X,A14,F15.5)')'Convergence = ',Conv
write(*,*)'--------------------------------------------------------------------------'
write(*,'(2X,A27,F15.6)') 'evGF2 HOMO energy (eV):',eGF2(HOMO)*HaToeV
write(*,'(2X,A27,F15.6)') 'evGF2 LUMO energy (eV):',eGF2(LUMO)*HaToeV
write(*,'(2X,A27,F15.6)') 'evGF2 HOMO-LUMO gap (eV):',Gap*HaToeV
write(*,'(2X,A30,F15.6)') 'evGF2 HOMO energy (eV):',eGF2(HOMO)*HaToeV
write(*,'(2X,A30,F15.6)') 'evGF2 LUMO energy (eV):',eGF2(LUMO)*HaToeV
write(*,'(2X,A30,F15.6)') 'evGF2 HOMO-LUMO gap (eV):',Gap*HaToeV
write(*,*)'--------------------------------------------------------------------------'
write(*,'(2X,A30,F15.6,A3)') 'evGF2 total energy :',ENuc + ERHF + Ec,' au'
write(*,'(2X,A30,F15.6,A3)') 'evGF2 correlation energy:',Ec,' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
end subroutine print_evGF2

12
src/MBPT/print_qsGF2.f90
View File

@ -1,4 +1,4 @@
subroutine print_qsGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,c,ENuc,P,T,V,J,K,F,SigC,Z,EqsGF2,dipole)
subroutine print_qsGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,c,ENuc,P,T,V,J,K,F,SigC,Z,EqsGF2,Ec,dipole)
! Print one-electron energies and other stuff for qsGF2
@ -44,8 +44,7 @@ subroutine print_qsGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,c,ENuc,P,T,V,J,K,F,SigC
EV = trace_matrix(nBas,matmul(P,V))
EJ = 0.5d0*trace_matrix(nBas,matmul(P,J))
Ex = 0.25d0*trace_matrix(nBas,matmul(P,K))
Ec = 0.50d0*trace_matrix(nBas,matmul(P,SigC))
EqsGF2 = ET + EV + EJ + Ex + Ec
EqsGF2 = ET + EV + EJ + Ex
! Dump results
@ -73,8 +72,9 @@ subroutine print_qsGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,c,ENuc,P,T,V,J,K,F,SigC
write(*,'(2X,A30,F15.6,A3)') 'qsGF2 LUMO energy:',eGF2(LUMO)*HaToeV,' eV'
write(*,'(2X,A30,F15.6,A3)') 'qsGF2 HOMO-LUMO gap :',Gap*HaToeV,' eV'
write(*,*)'-------------------------------------------'
write(*,'(2X,A30,F15.6,A3)') ' qsGF2 total energy:',EqsGF2 + ENuc,' au'
write(*,'(2X,A30,F15.6,A3)') ' qsGF2 total energy:',ENuc + EqsGF2 + Ec,' au'
write(*,'(2X,A30,F15.6,A3)') ' qsGF2 exchange energy:',Ex,' au'
write(*,'(2X,A30,F15.6,A3)') ' qsGF2 correlation energy:',Ec,' au'
write(*,*)'-------------------------------------------'
write(*,*)
@ -95,9 +95,9 @@ subroutine print_qsGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,c,ENuc,P,T,V,J,K,F,SigC
write(*,'(A32,1X,F16.10,A3)') ' Exchange energy: ',Ex,' au'
write(*,'(A32,1X,F16.10,A3)') ' Correlation energy: ',Ec,' au'
write(*,'(A50)') '---------------------------------------'
write(*,'(A32,1X,F16.10,A3)') ' Electronic energy: ',EqsGF2,' au'
write(*,'(A32,1X,F16.10,A3)') ' Electronic energy: ',EqsGF2 + Ec,' au'
write(*,'(A32,1X,F16.10,A3)') ' Nuclear repulsion: ',ENuc,' au'
write(*,'(A32,1X,F16.10,A3)') ' qsGF2 energy: ',ENuc + EqsGF2,' au'
write(*,'(A32,1X,F16.10,A3)') ' qsGF2 energy: ',ENuc + EqsGF2 + Ec,' au'
write(*,'(A50)') '---------------------------------------'
write(*,'(A35)') ' Dipole moment (Debye) '
write(*,'(10X,4A10)') 'X','Y','Z','Tot.'

11
src/MBPT/print_qsGW.f90
View File

@ -27,7 +27,7 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z
! Local variables
integer :: x,ixyz,HOMO,LUMO
double precision :: Gap,ET,EV,EJ,Ex,Ec
double precision :: Gap,ET,EV,EJ,Ex
double precision,external :: trace_matrix
! Output variables
@ -46,8 +46,7 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z
EV = trace_matrix(nBas,matmul(P,V))
EJ = 0.5d0*trace_matrix(nBas,matmul(P,J))
Ex = 0.25d0*trace_matrix(nBas,matmul(P,K))
! Ec = -0.50d0*trace_matrix(nBas,matmul(P,SigC))
EqsGW = ET + EV + EJ + Ex + EcGM
EqsGW = ET + EV + EJ + Ex
! Dump results
@ -75,9 +74,9 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z
write(*,'(2X,A30,F15.6,A3)') 'qsGW LUMO energy:',eGW(LUMO)*HaToeV,' eV'
write(*,'(2X,A30,F15.6,A3)') 'qsGW HOMO-LUMO gap :',Gap*HaToeV,' eV'
write(*,*)'-------------------------------------------'
write(*,'(2X,A30,F15.6,A3)') ' qsGW total energy:',EqsGW + ENuc,' au'
write(*,'(2X,A30,F15.6,A3)') ' qsGW total energy:',ENuc + EqsGW + EcGM,' au'
write(*,'(2X,A30,F15.6,A3)') ' qsGW exchange energy:',Ex,' au'
write(*,'(2X,A30,F15.6,A3)') ' qsGW correlation energy:',EcGM,' au'
write(*,'(2X,A30,F15.6,A3)') ' GM@qsGW correlation energy:',EcGM,' au'
write(*,'(2X,A30,F15.6,A3)') 'RPA@qsGW correlation energy:',EcRPA,' au'
write(*,*)'-------------------------------------------'
write(*,*)
@ -101,7 +100,7 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z
write(*,'(A50)') '---------------------------------------'
write(*,'(A32,1X,F16.10,A3)') ' Electronic energy: ',EqsGW,' au'
write(*,'(A32,1X,F16.10,A3)') ' Nuclear repulsion: ',ENuc,' au'
write(*,'(A32,1X,F16.10,A3)') ' qsGW energy: ',ENuc + EqsGW,' au'
write(*,'(A32,1X,F16.10,A3)') ' qsGW energy: ',ENuc + EqsGW + EcGM,' au'
write(*,'(A50)') '---------------------------------------'
write(*,'(A35)') ' Dipole moment (Debye) '
write(*,'(10X,4A10)') 'X','Y','Z','Tot.'

18
src/MBPT/print_qsUGW.f90
View File

@ -1,5 +1,5 @@
subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,PGW,Ov,T,V,J,K, &
ENuc,ET,EV,EJ,Ex,Ec,EcGM,EcRPA,EqsGW,SigC,Z,dipole)
ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGW,SigC,Z,dipole)
! Print one-electron energies and other stuff for qsUGW
@ -16,7 +16,6 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,PGW,Ov,T,V,J,K, &
double precision,intent(in) :: EV(nspin)
double precision,intent(in) :: EJ(nsp)
double precision,intent(in) :: Ex(nspin)
double precision,intent(in) :: Ec(nsp)
double precision,intent(in) :: EcGM(nspin)
double precision,intent(in) :: EcRPA
double precision,intent(in) :: EqsGW
@ -101,9 +100,8 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,PGW,Ov,T,V,J,K, &
write(*,'(2X,A30,F15.6,A3)') 'qsUGW HOMO-LUMO gap :',(minval(LUMO(:))-maxval(HOMO(:)))*HaToeV,' eV'
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
write(*,'(2X,A30,F15.6,A3)') ' qsUGW total energy:',EqsGW + ENuc,' au'
write(*,'(2X,A30,F15.6,A3)') ' qsUGW total energy:',ENuc + EqsGW + sum(EcGM(:)),' au'
write(*,'(2X,A30,F15.6,A3)') ' qsUGW exchange energy:',sum(Ex(:)),' au'
! write(*,'(2X,A30,F15.6,A3)') ' qsUGW correlation energy:',sum(Ec(:)),' au'
write(*,'(2X,A30,F15.6,A3)') ' GM@qsUGW correlation energy:',sum(EcGM(:)),' au'
write(*,'(2X,A30,F15.6,A3)') 'RPA@qsUGW correlation energy:',EcRPA,' au'
write(*,*)'-------------------------------------------------------------------------------&
@ -113,7 +111,6 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,PGW,Ov,T,V,J,K, &
! Dump results for final iteration
if(Conv < thresh) then
! if(.true.) then
write(*,*)
write(*,'(A60)') '-------------------------------------------------'
@ -145,14 +142,13 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,PGW,Ov,T,V,J,K, &
write(*,'(A40,1X,F16.10,A3)') ' Exchange a energy: ',Ex(1),' au'
write(*,'(A40,1X,F16.10,A3)') ' Exchange b energy: ',Ex(2),' au'
write(*,*)
! write(*,'(A40,1X,F16.10,A3)') ' Correlation energy: ',sum(Ec(:)),' au'
! write(*,'(A40,1X,F16.10,A3)') ' Correlation aa energy: ',Ec(1),' au'
! write(*,'(A40,1X,F16.10,A3)') ' Correlation ab energy: ',Ec(2),' au'
! write(*,'(A40,1X,F16.10,A3)') ' Correlation bb energy: ',Ec(3),' au'
write(*,'(A40,1X,F16.10,A3)') ' Correlation energy: ',sum(EcGM(:)),' au'
write(*,'(A40,1X,F16.10,A3)') ' Correlation aa energy: ',EcGM(1),' au'
write(*,'(A40,1X,F16.10,A3)') ' Correlation bb energy: ',EcGM(2),' au'
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A40,1X,F16.10,A3)') ' Electronic energy: ',EqsGW,' au'
write(*,'(A40,1X,F16.10,A3)') ' Electronic energy: ',EqsGW + sum(EcGM(:)),' au'
write(*,'(A40,1X,F16.10,A3)') ' Nuclear repulsion: ',ENuc,' au'
write(*,'(A40,1X,F16.10,A3)') ' qsUGW energy: ',EqsGW + ENuc,' au'
write(*,'(A40,1X,F16.10,A3)') ' qsUGW energy: ',ENuc + EqsGW + sum(EcGM(:)),' au'
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A40,F13.6)') ' S (exact) :',2d0*S_exact + 1d0
write(*,'(A40,F13.6)') ' S :',2d0*S + 1d0

5
src/MBPT/qsGF2.f90
View File

@ -52,6 +52,7 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,
double precision :: rcond
double precision,external :: trace_matrix
double precision :: dipole(ncart)
double precision :: Ec
double precision :: EcBSE(nspin)
double precision,allocatable :: error_diis(:,:)
@ -136,7 +137,7 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,
! Compute self-energy and renormalization factor
call self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI_MO,SigC,Z)
call self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI_MO,SigC,Z,Ec)
! Make correlation self-energy Hermitian and transform it back to AO basis
@ -177,7 +178,7 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,
! Print results
call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
call print_qsGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,c,ENuc,P,T,V,J,K,F,SigCp,Z,EqsGF2,dipole)
call print_qsGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,c,ENuc,P,T,V,J,K,F,SigCp,Z,EqsGF2,Ec,dipole)
enddo
!------------------------------------------------------------------------

18
src/MBPT/qsUGW.f90
View File

@ -75,7 +75,6 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS
double precision :: EV(nspin)
double precision :: EJ(nsp)
double precision :: Ex(nspin)
double precision :: Ec(nsp)
double precision :: EcRPA
double precision :: EcGM(nspin)
double precision :: EqsGW
@ -332,25 +331,16 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS
Ex(is) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,is),K(:,:,is)))
end do
! Correlation energy
Ec(:) = 0d0
! Ec(1) = - 0.25d0*trace_matrix(nBas,matmul(P(:,:,1),SigCp(:,:,1)))
! Ec(2) = - 0.25d0*trace_matrix(nBas,matmul(P(:,:,1),SigCp(:,:,2))) &
! - 0.25d0*trace_matrix(nBas,matmul(P(:,:,2),SigCp(:,:,1)))
! Ec(3) = - 0.25d0*trace_matrix(nBas,matmul(P(:,:,2),SigCp(:,:,2)))
! Total energy
EqsGW = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(Ex(:)) + sum(Ec(:))
EqsGW = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(Ex(:))
!------------------------------------------------------------------------
! Print results
!------------------------------------------------------------------------
call dipole_moment(nBas,P(:,:,1)+P(:,:,2),nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
call print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,P,S,T,V,J,K,ENuc,ET,EV,EJ,Ex,Ec,EcGM,EcRPA,EqsGW,SigCp,Z,dipole)
call print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,P,S,T,V,J,K,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGW,SigCp,Z,dipole)
enddo
!------------------------------------------------------------------------
@ -398,7 +388,7 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS
write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsUGW correlation energy (spin-conserved) =',EcBSE(1)
write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsUGW correlation energy (spin-flip) =',EcBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsUGW correlation energy =',EcBSE(1) + EcBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsUGW total energy =',ENuc + EUHF + EcBSE(1) + EcBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsUGW total energy =',ENuc + EqsGW + EcBSE(1) + EcBSE(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
@ -426,7 +416,7 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS
write(*,'(2X,A50,F20.10)') 'AC@BSE@qsUGW correlation energy (spin-conserved) =',EcAC(1)
write(*,'(2X,A50,F20.10)') 'AC@BSE@qsUGW correlation energy (spin-flip) =',EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@BSE@qsUGW correlation energy =',EcAC(1) + EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@BSE@qsUGW total energy =',ENuc + EUHF + EcAC(1) + EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@BSE@qsUGW total energy =',ENuc + EqsGW + EcAC(1) + EcAC(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)

68
src/MBPT/self_energy_GF2.f90
View File

@ -1,6 +1,6 @@
subroutine self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z)
subroutine self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z,Ec)
! Compute diagonal part of the GF2 self-energy and its renormalization factor
! Compute GF2 self-energy and its renormalization factor
implicit none
include 'parameters.h'
@ -16,49 +16,54 @@ subroutine self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z)
! Local variables
integer :: i,j,a,b
integer :: p
integer :: p,q
double precision :: eps
double precision :: num
! Output variables
double precision,intent(out) :: SigC(nBas)
double precision,intent(out) :: SigC(nBas,nBas)
double precision,intent(out) :: Z(nBas)
double precision,intent(out) :: Ec
! Initialize
SigC(:) = 0d0
Z(:) = 0d0
SigC(:,:) = 0d0
Z(:) = 0d0
! Compute GF2 self-energy
! Compute GF2 self-energy and renormalization factor
do p=nC+1,nBas-nR
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do q=nC+1,nBas-nR
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
eps = eGF2(p) + eHF(a) - eHF(i) - eHF(j)
num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)
eps = eGF2(p) + eHF(a) - eHF(i) - eHF(j)
num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(q,a,i,j)
SigC(p) = SigC(p) + num*eps/(eps**2 + eta**2)
Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
SigC(p,q) = SigC(p,q) + num*eps/(eps**2 + eta**2)
if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
end do
end do
end do
end do
do p=nC+1,nBas-nR
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do q=nC+1,nBas-nR
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
eps = eGF2(p) + eHF(i) - eHF(a) - eHF(b)
num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)
eps = eGF2(p) + eHF(i) - eHF(a) - eHF(b)
num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(q,i,a,b)
SigC(p) = SigC(p) + num*eps/(eps**2 + eta**2)
Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
SigC(p,q) = SigC(p,q) + num*eps/(eps**2 + eta**2)
if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
end do
end do
end do
@ -66,4 +71,23 @@ subroutine self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z)
Z(:) = 1d0/(1d0 - Z(:))
end subroutine self_energy_GF2_diag
! Compute correlaiton energy
Ec = 0d0
do j=nC+1,nO
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
eps = eGF2(j) + eHF(i) - eHF(a) - eHF(b)
num = (2d0*ERI(j,i,a,b) - ERI(j,i,b,a))*ERI(j,i,a,b)
Ec = Ec + num*eps/(eps**2 + eta**2)
end do
end do
end do
end do
end subroutine self_energy_GF2

66
src/MBPT/self_energy_GF2_diag.f90
View File

@ -1,6 +1,6 @@
subroutine self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z)
subroutine self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z,Ec)
! Compute GF2 self-energy and its renormalization factor
! Compute diagonal part of the GF2 self-energy and its renormalization factor
implicit none
include 'parameters.h'
@ -16,53 +16,50 @@ subroutine self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z)
! Local variables
integer :: i,j,a,b
integer :: p,q
integer :: p
double precision :: eps
double precision :: num
! Output variables
double precision,intent(out) :: SigC(nBas,nBas)
double precision,intent(out) :: SigC(nBas)
double precision,intent(out) :: Z(nBas)
double precision,intent(out) :: Ec
! Initialize
SigC(:,:) = 0d0
Z(:) = 0d0
SigC(:) = 0d0
Z(:) = 0d0
! Compute GF2 self-energy
do p=nC+1,nBas-nR
do q=nC+1,nBas-nR
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
eps = eGF2(p) + eHF(a) - eHF(i) - eHF(j)
num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(q,a,i,j)
eps = eGF2(p) + eHF(a) - eHF(i) - eHF(j)
num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)
SigC(p,q) = SigC(p,q) + num*eps/(eps**2 + eta**2)
if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
SigC(p) = SigC(p) + num*eps/(eps**2 + eta**2)
Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
end do
end do
end do
end do
do p=nC+1,nBas-nR
do q=nC+1,nBas-nR
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
eps = eGF2(p) + eHF(i) - eHF(a) - eHF(b)
num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(q,i,a,b)
eps = eGF2(p) + eHF(i) - eHF(a) - eHF(b)
num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)
SigC(p,q) = SigC(p,q) + num*eps/(eps**2 + eta**2)
if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
SigC(p) = SigC(p) + num*eps/(eps**2 + eta**2)
Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
end do
end do
end do
@ -70,4 +67,23 @@ subroutine self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z)
Z(:) = 1d0/(1d0 - Z(:))
end subroutine self_energy_GF2
! Compute correlaiton energy
Ec = 0d0
do j=nC+1,nO
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
eps = eGF2(j) + eHF(i) - eHF(a) - eHF(b)
num = (2d0*ERI(j,i,a,b) - ERI(j,i,b,a))*ERI(j,i,a,b)
Ec = Ec + num*eps/(eps**2 + eta**2)
end do
end do
end do
end do
end subroutine self_energy_GF2_diag