From 7fd5abe76f4e1d2f549c74d6e04ce5f74f6399f7 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Sat, 6 Mar 2021 23:08:43 +0100 Subject: [PATCH] clean up GM again --- input/methods | 6 +-- input/options | 4 +- src/MBPT/BSE2.f90 | 5 +++ src/MBPT/G0F2.f90 | 5 ++- src/MBPT/evGF2.f90 | 5 ++- src/MBPT/print_G0F2.f90 | 14 +++++-- src/MBPT/print_evGF2.f90 | 14 +++++-- src/MBPT/print_qsGF2.f90 | 12 +++--- src/MBPT/print_qsGW.f90 | 11 +++-- src/MBPT/print_qsUGW.f90 | 18 ++++---- src/MBPT/qsGF2.f90 | 5 ++- src/MBPT/qsUGW.f90 | 18 ++------ src/MBPT/self_energy_GF2.f90 | 68 +++++++++++++++++++++---------- src/MBPT/self_energy_GF2_diag.f90 | 66 ++++++++++++++++++------------ 14 files changed, 148 insertions(+), 103 deletions(-) diff --git a/input/methods b/input/methods index 89fb8f0..d790e69 100644 --- a/input/methods +++ b/input/methods @@ -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 diff --git a/input/options b/input/options index 347efea..f64a274 100644 --- a/input/options +++ b/input/options @@ -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 diff --git a/src/MBPT/BSE2.f90 b/src/MBPT/BSE2.f90 index 123cfea..25bc56c 100644 --- a/src/MBPT/BSE2.f90 +++ b/src/MBPT/BSE2.f90 @@ -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 diff --git a/src/MBPT/G0F2.f90 b/src/MBPT/G0F2.f90 index 303e5ea..bc1c7ab 100644 --- a/src/MBPT/G0F2.f90 +++ b/src/MBPT/G0F2.f90 @@ -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 diff --git a/src/MBPT/evGF2.f90 b/src/MBPT/evGF2.f90 index ff39619..eb5d775 100644 --- a/src/MBPT/evGF2.f90 +++ b/src/MBPT/evGF2.f90 @@ -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 diff --git a/src/MBPT/print_G0F2.f90 b/src/MBPT/print_G0F2.f90 index e560db8..87b545b 100644 --- a/src/MBPT/print_G0F2.f90 +++ b/src/MBPT/print_G0F2.f90 @@ -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 diff --git a/src/MBPT/print_evGF2.f90 b/src/MBPT/print_evGF2.f90 index fdab254..07c8b42 100644 --- a/src/MBPT/print_evGF2.f90 +++ b/src/MBPT/print_evGF2.f90 @@ -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 diff --git a/src/MBPT/print_qsGF2.f90 b/src/MBPT/print_qsGF2.f90 index c78100b..558e8b8 100644 --- a/src/MBPT/print_qsGF2.f90 +++ b/src/MBPT/print_qsGF2.f90 @@ -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.' diff --git a/src/MBPT/print_qsGW.f90 b/src/MBPT/print_qsGW.f90 index d6be17b..97b65be 100644 --- a/src/MBPT/print_qsGW.f90 +++ b/src/MBPT/print_qsGW.f90 @@ -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.' diff --git a/src/MBPT/print_qsUGW.f90 b/src/MBPT/print_qsUGW.f90 index 32dc633..d776764 100644 --- a/src/MBPT/print_qsUGW.f90 +++ b/src/MBPT/print_qsUGW.f90 @@ -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 diff --git a/src/MBPT/qsGF2.f90 b/src/MBPT/qsGF2.f90 index 0248685..7ae192b 100644 --- a/src/MBPT/qsGF2.f90 +++ b/src/MBPT/qsGF2.f90 @@ -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 !------------------------------------------------------------------------ diff --git a/src/MBPT/qsUGW.f90 b/src/MBPT/qsUGW.f90 index d5b84d2..de27201 100644 --- a/src/MBPT/qsUGW.f90 +++ b/src/MBPT/qsUGW.f90 @@ -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(*,*) diff --git a/src/MBPT/self_energy_GF2.f90 b/src/MBPT/self_energy_GF2.f90 index c64d507..21c1520 100644 --- a/src/MBPT/self_energy_GF2.f90 +++ b/src/MBPT/self_energy_GF2.f90 @@ -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 diff --git a/src/MBPT/self_energy_GF2_diag.f90 b/src/MBPT/self_energy_GF2_diag.f90 index e4dcdd1..e0611a0 100644 --- a/src/MBPT/self_energy_GF2_diag.f90 +++ b/src/MBPT/self_energy_GF2_diag.f90 @@ -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