From e1f67c94d42094ca5df3066346e339c5722495bc Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Thu, 6 Jan 2022 13:48:15 +0100 Subject: [PATCH] GM in GT --- input/dft | 4 ++-- input/methods | 4 ++-- input/options | 4 ++-- mol/h2.xyz | 2 +- src/GT/G0T0.f90 | 8 ++++--- src/GT/evGT.f90 | 34 +++++++---------------------- src/GT/print_G0T0.f90 | 16 +++++++++----- src/GT/print_evGT.f90 | 13 ++++++++++- src/GT/print_qsGT.f90 | 12 +++++----- src/GT/qsGT.f90 | 30 +++++-------------------- src/GT/self_energy_Tmatrix.f90 | 31 ++++++++++++++++++++++---- src/GT/self_energy_Tmatrix_diag.f90 | 27 +++++++++++++++++++++-- src/GW/self_energy_correlation.f90 | 4 ++-- 13 files changed, 108 insertions(+), 81 deletions(-) diff --git a/input/dft b/input/dft index bf9a080..18dbd67 100644 --- a/input/dft +++ b/input/dft @@ -13,7 +13,7 @@ # GGA = 2: LYP,PBE # MGGA = 3: # Hybrid = 4: HF,B3LYP,PBE -0 H +1 VWN5 # quadrature grid SG-n 1 # Number of states in ensemble (nEns) @@ -31,7 +31,7 @@ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 # Ensemble weights: wEns(1),...,wEns(nEns-1) - 1.0 0.0 0.0 + 0.6 0.0 0.0 # Ncentered ? F # Parameters for CC weight-dependent exchange functional diff --git a/input/methods b/input/methods index 8124aba..e19f51c 100644 --- a/input/methods +++ b/input/methods @@ -1,5 +1,5 @@ # RHF UHF KS MOM - F F T F + T F F F # MP2* MP3 MP2-F12 F F F # CCD pCCD DCD CCSD CCSD(T) @@ -15,7 +15,7 @@ # G0W0* evGW* qsGW* ufG0W0 ufGW F F F F F # G0T0 evGT qsGT - F F F + F F T # MCMP2 F # * unrestricted version available diff --git a/input/options b/input/options index 2d0e99d..a3a8628 100644 --- a/input/options +++ b/input/options @@ -1,5 +1,5 @@ # HF: maxSCF thresh DIIS n_diis guess_type ortho_type mix_guess stability - 1024 0.00001 T 5 1 1 T F + 1024 0.00001 T 5 1 1 F F # MP: # CC: maxSCF thresh DIIS n_diis @@ -15,6 +15,6 @@ # ACFDT: AC Kx XBS F F F # BSE: BSE dBSE dTDA evDyn - T T T F + F T T F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/mol/h2.xyz b/mol/h2.xyz index a4e936a..d955cc4 100644 --- a/mol/h2.xyz +++ b/mol/h2.xyz @@ -1,4 +1,4 @@ 2 H 0. 0. 0. -H 0. 0. 2.000000 +H 0. 0. 0.7 diff --git a/src/GT/G0T0.f90 b/src/GT/G0T0.f90 index 94f886c..16225b3 100644 --- a/src/GT/G0T0.f90 +++ b/src/GT/G0T0.f90 @@ -48,6 +48,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing double precision :: EcRPA(nspin) double precision :: EcBSE(nspin) double precision :: EcAC(nspin) + double precision :: EcGM double precision,allocatable :: Omega1s(:),Omega1t(:) double precision,allocatable :: X1s(:,:),X1t(:,:) double precision,allocatable :: Y1s(:,:),Y1t(:,:) @@ -134,6 +135,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing ! Compute T-matrix version of the self-energy !---------------------------------------------- + EcGM = 0d0 SigT(:) = 0d0 Z(:) = 0d0 @@ -141,7 +143,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI_MO,X1s,Y1s,rho1s,X2s,Y2s,rho2s) - call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF,Omega1s,rho1s,Omega2s,rho2s,SigT) + call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF,Omega1s,rho1s,Omega2s,rho2s,EcGM,SigT) call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF,Omega1s,rho1s,Omega2s,rho2s,Z) @@ -149,7 +151,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI_MO,X1t,Y1t,rho1t,X2t,Y2t,rho2t) - call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF,Omega1t,rho1t,Omega2t,rho2t,SigT) + call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF,Omega1t,rho1t,Omega2t,rho2t,EcGM,SigT) call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF,Omega1t,rho1t,Omega2t,rho2t,Z) @@ -194,7 +196,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing EcRPA(1) = EcRPA(1) - EcRPA(2) EcRPA(2) = 3d0*EcRPA(2) - call print_G0T0(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eG0T0,EcRPA) + call print_G0T0(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eG0T0,EcGM,EcRPA) ! Perform BSE calculation diff --git a/src/GT/evGT.f90 b/src/GT/evGT.f90 index 65dcfa4..016cf32 100644 --- a/src/GT/evGT.f90 +++ b/src/GT/evGT.f90 @@ -55,6 +55,7 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & integer :: iblock integer :: nOOs,nOOt integer :: nVVs,nVVt + double precision :: EcGM double precision :: EcRPA(nspin) double precision :: EcBSE(nspin) double precision :: EcAC(nspin) @@ -148,10 +149,14 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT,ERI_MO, & Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin)) + EcRPA(1) = EcRPA(1) - EcRPA(2) + EcRPA(2) = 3d0*EcRPA(2) + !---------------------------------------------- ! Compute T-matrix version of the self-energy !---------------------------------------------- + EcGM = 0d0 SigT(:) = 0d0 Z(:) = 0d0 @@ -161,7 +166,7 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & X1s,Y1s,rho1s,X2s,Y2s,rho2s) call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, & - Omega1s,rho1s,Omega2s,rho2s,SigT) + Omega1s,rho1s,Omega2s,rho2s,EcGM,SigT) call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, & Omega1s,rho1s,Omega2s,rho2s,Z) @@ -172,7 +177,7 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & X1t,Y1t,rho1t,X2t,Y2t,rho2t) call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, & - Omega1t,rho1t,Omega2t,rho2t,SigT) + Omega1t,rho1t,Omega2t,rho2t,EcGM,SigT) call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, & Omega1t,rho1t,Omega2t,rho2t,Z) @@ -195,7 +200,7 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & ! Dump results !---------------------------------------------- - call print_evGT(nBas,nO,nSCF,Conv,eHF,SigT,Z,eGT) + call print_evGT(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigT,Z,eGT,EcGM,EcRPA) ! DIIS extrapolation @@ -219,29 +224,6 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & ! End main loop !------------------------------------------------------------------------ -! Compute the ppRPA correlation energy - - ispin = 1 - iblock = 3 - call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT,ERI_MO, & - Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin)) - ispin = 2 - iblock = 4 - call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT,ERI_MO, & - Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin)) - EcRPA(1) = EcRPA(1) - EcRPA(2) - EcRPA(2) = 3d0*EcRPA(2) - - write(*,*) - write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@evGT correlation energy (singlet) =',EcRPA(1) - write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@evGT correlation energy (triplet) =',EcRPA(2) - write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@evGT correlation energy =',EcRPA(1) + EcRPA(2) - write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@evGT total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2) - write(*,*)'-------------------------------------------------------------------------------' - write(*,*) - - ! Perform BSE calculation if(BSE) then diff --git a/src/GT/print_G0T0.f90 b/src/GT/print_G0T0.f90 index 2150cd4..95b4949 100644 --- a/src/GT/print_G0T0.f90 +++ b/src/GT/print_G0T0.f90 @@ -1,13 +1,15 @@ -subroutine print_G0T0(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eGT,EcRPA) +subroutine print_G0T0(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eGT,EcGM,EcRPA) ! Print one-electron energies and other stuff for G0T0 implicit none include 'parameters.h' - integer,intent(in) :: nBas,nO + integer,intent(in) :: nBas + integer,intent(in) :: nO double precision,intent(in) :: ENuc double precision,intent(in) :: ERHF + double precision,intent(in) :: EcGM double precision,intent(in) :: EcRPA(nspin) double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: SigT(nBas) @@ -50,10 +52,12 @@ subroutine print_G0T0(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eGT,EcRPA) write(*,'(2X,A50,F15.6,A3)') 'G0T0 LUMO energy (eV) =',eGT(LUMO)*HaToeV,' eV' write(*,'(2X,A50,F15.6,A3)') 'G0T0 HOMO-LUMO gap (eV) =',Gap*HaToeV,' eV' write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10,A3)') 'Tr@RPA@G0T0 correlation energy (singlet) =',EcRPA(1),' au' - write(*,'(2X,A50,F20.10,A3)') 'Tr@RPA@G0T0 correlation energy (triplet) =',EcRPA(2),' au' - write(*,'(2X,A50,F20.10,A3)') 'Tr@RPA@G0T0 correlation energy =',EcRPA(1) + EcRPA(2),' au' - write(*,'(2X,A50,F20.10,A3)') 'Tr@RPA@G0T0 total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2),' au' + write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@G0T0 correlation energy (singlet) =',EcRPA(1),' au' + write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@G0T0 correlation energy (triplet) =',EcRPA(2),' au' + write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@G0T0 correlation energy =',EcRPA(1) + EcRPA(2),' au' + write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@G0T0 total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2),' au' + write(*,'(2X,A50,F20.10,A3)') ' GM@G0T0 correlation energy =',EcGM,' au' + write(*,'(2X,A50,F20.10,A3)') ' GM@G0T0 total energy =',ENuc + ERHF + EcGM,' au' write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/GT/print_evGT.f90 b/src/GT/print_evGT.f90 index 9b15356..7a11544 100644 --- a/src/GT/print_evGT.f90 +++ b/src/GT/print_evGT.f90 @@ -1,4 +1,4 @@ -subroutine print_evGT(nBas,nO,nSCF,Conv,eHF,SigT,Z,eGT) +subroutine print_evGT(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigT,Z,eGT,EcGM,EcRPA) ! Print one-electron energies and other stuff for evGT @@ -10,6 +10,10 @@ subroutine print_evGT(nBas,nO,nSCF,Conv,eHF,SigT,Z,eGT) integer,intent(in) :: nSCF double precision,intent(in) :: Conv double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: ENuc + double precision,intent(in) :: ERHF + double precision,intent(in) :: EcGM + double precision,intent(in) :: EcRPA(nspin) double precision,intent(in) :: SigT(nBas) double precision,intent(in) :: Z(nBas) double precision,intent(in) :: eGT(nBas) @@ -49,6 +53,13 @@ subroutine print_evGT(nBas,nO,nSCF,Conv,eHF,SigT,Z,eGT) write(*,'(2X,A30,F15.6)') 'evGT LUMO energy (eV):',eGT(LUMO)*HaToeV write(*,'(2X,A30,F15.6)') 'evGT HOMO-LUMO gap (eV):',Gap*HaToeV write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@G0T0 correlation energy (singlet) =',EcRPA(1),' au' + write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@G0T0 correlation energy (triplet) =',EcRPA(2),' au' + write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@G0T0 correlation energy =',EcRPA(1) + EcRPA(2),' au' + write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@G0T0 total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2),' au' + write(*,'(2X,A50,F20.10,A3)') ' GM@G0T0 correlation energy =',EcGM,' au' + write(*,'(2X,A50,F20.10,A3)') ' GM@G0T0 total energy =',ENuc + ERHF + EcGM,' au' + write(*,*)'-------------------------------------------------------------------------------' write(*,*) end subroutine print_evGT diff --git a/src/GT/print_qsGT.f90 b/src/GT/print_qsGT.f90 index 5883b62..e383c6e 100644 --- a/src/GT/print_qsGT.f90 +++ b/src/GT/print_qsGT.f90 @@ -48,9 +48,9 @@ subroutine print_qsGT(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ,Ex write(*,*)'-------------------------------------------------------------------------------' if(nSCF < 10) then - write(*,'(1X,A21,I1,A1,I1,A12)')' Self-consistent qsG',nSCF,'W',nSCF,' calculation' + write(*,'(1X,A21,I1,A1,I1,A12)')' Self-consistent qsG',nSCF,'T',nSCF,' calculation' else - write(*,'(1X,A21,I2,A1,I2,A12)')' Self-consistent qsG',nSCF,'W',nSCF,' calculation' + write(*,'(1X,A21,I2,A1,I2,A12)')' Self-consistent qsG',nSCF,'T',nSCF,' calculation' endif 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)') & @@ -70,10 +70,10 @@ subroutine print_qsGT(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ,Ex write(*,'(2X,A30,F15.6,A3)') 'qsGT LUMO energy:',eGT(LUMO)*HaToeV,' eV' write(*,'(2X,A30,F15.6,A3)') 'qsGT HOMO-LUMO gap :',Gap*HaToeV,' eV' write(*,*)'-------------------------------------------' - write(*,'(2X,A30,F15.6,A3)') ' qsGT total energy:',ENuc + EqsGT,' au' - write(*,'(2X,A30,F15.6,A3)') ' qsGT exchange energy:',Ex,' au' - write(*,'(2X,A30,F15.6,A3)') ' GM@qsGT correlation energy:',EcGM,' au' - write(*,'(2X,A30,F15.6,A3)') 'RPA@qsGT correlation energy:',sum(EcRPA(:)),' au' + write(*,'(2X,A30,F15.6,A3)') ' qsGT total energy:',ENuc + EqsGT,' au' + write(*,'(2X,A30,F15.6,A3)') ' qsGT exchange energy:',Ex,' au' + write(*,'(2X,A30,F15.6,A3)') ' GM@qsGT correlation energy:',EcGM,' au' + write(*,'(2X,A30,F15.6,A3)') 'ppRPA@qsGT correlation energy:',sum(EcRPA(:)),' au' write(*,*)'-------------------------------------------' write(*,*) diff --git a/src/GT/qsGT.f90 b/src/GT/qsGT.f90 index 6a5aa5a..7b7e474 100644 --- a/src/GT/qsGT.f90 +++ b/src/GT/qsGT.f90 @@ -198,8 +198,12 @@ subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,T call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT,ERI_MO, & Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin)) + EcRPA(1) = EcRPA(1) - EcRPA(2) + EcRPA(2) = 3d0*EcRPA(2) + ! Compute correlation part of the self-energy + EcGM = 0d0 SigT(:,:) = 0d0 Z(:) = 0d0 @@ -209,7 +213,7 @@ subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,T X1s,Y1s,rho1s,X2s,Y2s,rho2s) call self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, & - Omega1s,rho1s,Omega2s,rho2s,SigT) + Omega1s,rho1s,Omega2s,rho2s,EcGM,SigT) call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, & Omega1s,rho1s,Omega2s,rho2s,Z) @@ -220,7 +224,7 @@ subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,T X1t,Y1t,rho1t,X2t,Y2t,rho2t) call self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, & - Omega1t,rho1t,Omega2t,rho2t,SigT) + Omega1t,rho1t,Omega2t,rho2t,EcGM,SigT) call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, & Omega1t,rho1t,Omega2t,rho2t,Z) @@ -302,28 +306,6 @@ subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,T ! End main loop !------------------------------------------------------------------------ -! Compute the ppRPA correlation energy - - ispin = 1 - iblock = 3 - call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT,ERI_MO, & - Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin)) - ispin = 2 - iblock = 4 - call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT,ERI_MO, & - Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin)) - EcRPA(1) = EcRPA(1) - EcRPA(2) - EcRPA(2) = 3d0*EcRPA(2) - - write(*,*) - write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@qsGT correlation energy (singlet) =',EcRPA(1) - write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@qsGT correlation energy (triplet) =',EcRPA(2) - write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@qsGT correlation energy =',EcRPA(1) + EcRPA(2) - write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@qsGT total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2) - write(*,*)'-------------------------------------------------------------------------------' - write(*,*) - ! Did it actually converge? if(nSCF == maxSCF+1) then diff --git a/src/GT/self_energy_Tmatrix.f90 b/src/GT/self_energy_Tmatrix.f90 index 3a6c515..0840419 100644 --- a/src/GT/self_energy_Tmatrix.f90 +++ b/src/GT/self_energy_Tmatrix.f90 @@ -1,4 +1,4 @@ -subroutine self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,SigT) +subroutine self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,EcGM,SigT) ! Compute the correlation part of the T-matrix self-energy @@ -23,12 +23,13 @@ subroutine self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2 ! Local variables - integer :: i,a,p,q,cd,kl + integer :: i,j,a,b,p,q,cd,kl double precision :: eps ! Output variables - double precision,intent(inout) :: SigT(nBas,nBas) + double precision,intent(inout):: EcGM + double precision,intent(inout):: SigT(nBas,nBas) !---------------------------------------------- ! Occupied part of the T-matrix self-energy @@ -46,7 +47,7 @@ subroutine self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2 enddo !---------------------------------------------- - ! Virtual part of the T-matrix self-energy +! Virtual part of the T-matrix self-energy !---------------------------------------------- do p=nC+1,nBas-nR @@ -60,4 +61,26 @@ subroutine self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2 enddo enddo +!---------------------------------------------- +! Galitskii-Migdal correlation energy +!---------------------------------------------- + + do i=nC+1,nO + do j=nC+1,nO + do cd=1,nVV + eps = e(i) + e(j) - Omega2(cd) + EcGM = EcGM - rho1(i,j,cd)*rho1(i,j,cd)*eps/(eps**2 + eta**2) + enddo + enddo + enddo + + do a=nO+1,nBas-nR + do b=nO+1,nBas-nR + do kl=1,nOO + eps = e(a) + e(b) - Omega2(kl) + EcGM = EcGM + rho2(a,b,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2) + enddo + enddo + enddo + end subroutine self_energy_Tmatrix diff --git a/src/GT/self_energy_Tmatrix_diag.f90 b/src/GT/self_energy_Tmatrix_diag.f90 index 76884d0..d6d1b48 100644 --- a/src/GT/self_energy_Tmatrix_diag.f90 +++ b/src/GT/self_energy_Tmatrix_diag.f90 @@ -1,4 +1,4 @@ -subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,SigT) +subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,EcGM,SigT) ! Compute diagonal of the correlation part of the T-matrix self-energy @@ -23,11 +23,12 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,O ! Local variables - integer :: i,a,p,cd,kl + integer :: i,j,a,b,p,cd,kl double precision :: eps ! Output variables + double precision,intent(inout) :: EcGM double precision,intent(inout) :: SigT(nBas) !---------------------------------------------- @@ -56,4 +57,26 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,O enddo enddo +!---------------------------------------------- +! Galitskii-Migdal correlation energy +!---------------------------------------------- + + do i=nC+1,nO + do j=nC+1,nO + do cd=1,nVV + eps = e(i) + e(j) - Omega2(cd) + EcGM = EcGM - rho1(i,j,cd)*rho1(i,j,cd)*eps/(eps**2 + eta**2) + enddo + enddo + enddo + + do a=nO+1,nBas-nR + do b=nO+1,nBas-nR + do kl=1,nOO + eps = e(a) + e(b) - Omega2(kl) + EcGM = EcGM + rho2(a,b,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2) + enddo + enddo + enddo + end subroutine self_energy_Tmatrix_diag diff --git a/src/GW/self_energy_correlation.f90 b/src/GW/self_energy_correlation.f90 index 4d62b43..ebe0b67 100644 --- a/src/GW/self_energy_correlation.f90 +++ b/src/GW/self_energy_correlation.f90 @@ -28,8 +28,8 @@ subroutine self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Ec ! Output variables - double precision,intent(out) :: SigC(nBas,nBas) double precision,intent(out) :: EcGM + double precision,intent(out) :: SigC(nBas,nBas) ! Initialize @@ -102,7 +102,7 @@ subroutine self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Ec end do end do - ! GM correlation energy + ! Galitskii-Migdal correlation energy EcGM = 0d0 do i=nC+1,nO