From 28ca539933969a0761f897c54eb3aa41d4975898 Mon Sep 17 00:00:00 2001 From: arjanberger10 Date: Thu, 14 Feb 2019 22:42:07 +0100 Subject: [PATCH] GM correlaction energy added and correction to COHSEX self-energy --- src/MCQC/print_G0W0.f90 | 8 ++++--- src/MCQC/print_evGW.f90 | 8 ++++--- src/MCQC/print_qsGW.f90 | 6 +++-- src/MCQC/qsGW.f90 | 8 +++---- src/MCQC/self_energy_correlation.f90 | 11 +++++++-- src/MCQC/self_energy_correlation_diag.f90 | 27 +++++++++++++++++++++-- 6 files changed, 52 insertions(+), 16 deletions(-) diff --git a/src/MCQC/print_G0W0.f90 b/src/MCQC/print_G0W0.f90 index e3f96a5..e069f4c 100644 --- a/src/MCQC/print_G0W0.f90 +++ b/src/MCQC/print_G0W0.f90 @@ -1,4 +1,4 @@ -subroutine print_G0W0(nBas,nO,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA) +subroutine print_G0W0(nBas,nO,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA,EcGM) ! Print one-electron energies and other stuff for G0W0 @@ -6,7 +6,7 @@ subroutine print_G0W0(nBas,nO,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA) include 'parameters.h' integer,intent(in) :: nBas,nO - double precision,intent(in) :: ENuc,EHF,EcRPA + double precision,intent(in) :: ENuc,EHF,EcRPA,EcGM double precision,intent(in) :: e(nBas),SigmaC(nBas),Z(nBas),eGW(nBas) integer :: x,HOMO,LUMO @@ -37,8 +37,10 @@ subroutine print_G0W0(nBas,nO,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA) write(*,'(2X,A27,F15.6)') 'G0W0 LUMO energy (eV):',eGW(LUMO)*HaToeV write(*,'(2X,A27,F15.6)') 'G0W0 HOMO-LUMO gap (eV):',Gap*HaToeV write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A27,F15.6)') 'G0W0 total energy =',ENuc + EHF + EcRPA + write(*,'(2X,A27,F15.6)') 'G0W0 RPA total energy =',ENuc + EHF + EcRPA write(*,'(2X,A27,F15.6)') 'RPA correlation energy =',EcRPA + write(*,'(2X,A27,F15.6)') 'G0W0 GM total energy =',ENuc + EHF + EcGM + write(*,'(2X,A27,F15.6)') 'GM correlation energy =',EcGM write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/MCQC/print_evGW.f90 b/src/MCQC/print_evGW.f90 index a3b143c..54acbcb 100644 --- a/src/MCQC/print_evGW.f90 +++ b/src/MCQC/print_evGW.f90 @@ -1,4 +1,4 @@ -subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA) +subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA,EcGM) ! Print one-electron energies and other stuff for evGW @@ -6,7 +6,7 @@ subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA) include 'parameters.h' integer,intent(in) :: nBas,nO,nSCF - double precision,intent(in) :: ENuc,EHF,EcRPA + double precision,intent(in) :: ENuc,EHF,EcRPA,EcGM double precision,intent(in) :: Conv,e(nBas),SigmaC(nBas),Z(nBas),eGW(nBas) integer :: x,HOMO,LUMO @@ -44,8 +44,10 @@ subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA) write(*,'(2X,A27,F15.6)') 'evGW LUMO energy (eV):',eGW(LUMO)*HaToeV write(*,'(2X,A27,F15.6)') 'evGW HOMO-LUMO gap (eV):',Gap*HaToeV write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A27,F15.6)') 'evGW total energy =',ENuc + EHF + EcRPA + write(*,'(2X,A27,F15.6)') 'evGW RPA total energy =',ENuc + EHF + EcRPA write(*,'(2X,A27,F15.6)') 'RPA correlation energy =',EcRPA + write(*,'(2X,A27,F15.6)') 'evGW GM total energy =',ENuc + EHF + EcGM + write(*,'(2X,A27,F15.6)') 'GM correlation energy =',EcGM write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/MCQC/print_qsGW.f90 b/src/MCQC/print_qsGW.f90 index bf7eb87..f92be3d 100644 --- a/src/MCQC/print_qsGW.f90 +++ b/src/MCQC/print_qsGW.f90 @@ -1,4 +1,4 @@ -subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,SigmaC,Z,EcRPA) +subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,SigmaC,Z,EcRPA,EcGM) ! Print one-electron energies and other stuff for qsGW @@ -9,7 +9,7 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,Sig ! Input variables integer,intent(in) :: nBas,nO,nSCF - double precision,intent(in) :: ENuc,EcRPA,Conv,thresh + double precision,intent(in) :: ENuc,EcRPA,EcGM,Conv,thresh double precision,intent(in) :: eHF(nBas),eGW(nBas),c(nBas),P(nBas,nBas) double precision,intent(in) :: T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas) double precision,intent(in) :: J(nBas,nBas),K(nBas,nBas),F(nBas,nBas) @@ -63,9 +63,11 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,Sig write(*,'(2X,A27,F15.6)') 'qsGW HOMO-LUMO gap (eV):',Gap*HaToeV write(*,*)'-------------------------------------------' write(*,'(2X,A27,F15.6)') 'qsGW total energy =',EqsGW + ENuc + write(*,'(2X,A27,F15.6)') 'qsGW GM total energy =',EqsGW + ENuc + EcGM write(*,'(2X,A27,F15.6)') 'qsGW exchange energy =',Ex write(*,'(2X,A27,F15.6)') 'qsGW correlation energy =',Ec write(*,'(2X,A27,F15.6)') 'RPA correlation energy =',EcRPA + write(*,'(2X,A27,F15.6)') 'GM correlation energy =',EcGM write(*,*)'-------------------------------------------' write(*,*) diff --git a/src/MCQC/qsGW.f90 b/src/MCQC/qsGW.f90 index 364ad66..1b34714 100644 --- a/src/MCQC/qsGW.f90 +++ b/src/MCQC/qsGW.f90 @@ -21,7 +21,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani logical :: dRPA integer :: nSCF,nBasSq,ispin,n_diis - double precision :: EcRPA,Conv + double precision :: EcRPA,EcGM,Conv double precision,external :: trace_matrix double precision,allocatable :: error_diis(:,:),F_diis(:,:) double precision,allocatable :: Omega(:,:),XpY(:,:,:),rho(:,:,:,:),rhox(:,:,:,:) @@ -106,13 +106,13 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani if(G0W) then call self_energy_correlation(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,eHF, & - Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),SigC) + Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC) call renormalization_factor(SOSEX,nBas,nC,nO,nV,nR,nS,eHF,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z) else call self_energy_correlation(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e, & - Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),SigC) + Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC) call renormalization_factor(SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z) endif @@ -152,7 +152,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) + call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,e,c,ENuc,P,T,V,Hc,J,K,F,SigCp,Z,EcRPA,EcGM) ! Increment diff --git a/src/MCQC/self_energy_correlation.f90 b/src/MCQC/self_energy_correlation.f90 index 07f6436..d1d1ad2 100644 --- a/src/MCQC/self_energy_correlation.f90 +++ b/src/MCQC/self_energy_correlation.f90 @@ -1,4 +1,4 @@ -subroutine self_energy_correlation(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega,rho,rhox,SigC) +subroutine self_energy_correlation(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega,rho,rhox,EcGM,SigC) ! Compute correlation part of the self-energy @@ -19,6 +19,7 @@ subroutine self_energy_correlation(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega,rho, ! Output variables double precision,intent(out) :: SigC(nBas,nBas) + double precision,intent(out) :: EcGM ! Initialize @@ -41,7 +42,8 @@ subroutine self_energy_correlation(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega,rho, do j=nC+1,nO do b=nO+1,nBas-nR jb = jb + 1 - SigC(x,y) = SigC(x,y) + 4d0*rho(x,i,jb)*rho(y,i,jb)/Omega(jb) +! SigC(x,y) = SigC(x,y) + 4d0*rho(x,i,jb)*rho(y,i,jb)/Omega(jb) + SigC(x,y) = SigC(x,y) + 2d0*rho(x,i,jb)*rho(y,i,jb)/Omega(jb) enddo enddo enddo @@ -64,6 +66,11 @@ subroutine self_energy_correlation(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega,rho, enddo enddo + EcGM=0d0 + do i=nC+1,nO + EcGM = EcGM + SigC(i,i) + enddo + else ! Occupied part of the correlation self-energy diff --git a/src/MCQC/self_energy_correlation_diag.f90 b/src/MCQC/self_energy_correlation_diag.f90 index 9116561..f9b426e 100644 --- a/src/MCQC/self_energy_correlation_diag.f90 +++ b/src/MCQC/self_energy_correlation_diag.f90 @@ -1,4 +1,4 @@ -subroutine self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega,rho,rhox,SigC) +subroutine self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega,rho,rhox,EcGM,SigC) ! Compute diagonal of the correlation part of the self-energy @@ -20,6 +20,7 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega ! Output variables double precision,intent(out) :: SigC(nBas) + double precision,intent(out) :: EcGM ! Initialize @@ -42,7 +43,8 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega do j=nC+1,nO do b=nO+1,nBas-nR jb = jb + 1 - SigC(x) = SigC(x) + 4d0*rho(x,i,jb)**2/Omega(jb) +! SigC(x) = SigC(x) + 4d0*rho(x,i,jb)**2/Omega(jb) + SigC(x) = SigC(x) + 2d0*rho(x,i,jb)**2/Omega(jb) enddo enddo enddo @@ -61,6 +63,11 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega enddo enddo enddo + + EcGM=0d0 + do i=nC+1,nO + EcGM = EcGM + SigC(i) + enddo else @@ -98,6 +105,22 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega enddo enddo +! GM correlation energy + + EcGM=0d0 + do i=nC+1,nO + do a=nO+1,nBas-nR + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + eps = e(a) - e(i) + Omega(jb) + EcGM = EcGM - 4d0*rho(a,i,jb)*rho(a,i,jb)/eps + enddo + enddo + enddo + enddo + if(SOSEX) then ! SOSEX: occupied part of the correlation self-energy