PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2025-01-03 01:55:57 +01:00
Code Releases Activity

GM correlaction energy added and correction to COHSEX self-energy

Browse Source
This commit is contained in:
arjanberger10 2019-02-14 22:42:07 +01:00
parent 53cfa3a04d
commit 28ca539933
6 changed files with 52 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
src/MCQC/print_G0W0.f90
View File

@ -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 ! 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' include 'parameters.h'
integer,intent(in) :: nBas,nO 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) double precision,intent(in) :: e(nBas),SigmaC(nBas),Z(nBas),eGW(nBas)
integer :: x,HOMO,LUMO 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 LUMO energy (eV):',eGW(LUMO)*HaToeV
write(*,'(2X,A27,F15.6)') 'G0W0 HOMO-LUMO gap (eV):',Gap*HaToeV write(*,'(2X,A27,F15.6)') 'G0W0 HOMO-LUMO gap (eV):',Gap*HaToeV
write(*,*)'-------------------------------------------------------------------------------' 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)') '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(*,*)'-------------------------------------------------------------------------------'
write(*,*) write(*,*)

8
src/MCQC/print_evGW.f90
View File

@ -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 ! 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' include 'parameters.h'
integer,intent(in) :: nBas,nO,nSCF 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) double precision,intent(in) :: Conv,e(nBas),SigmaC(nBas),Z(nBas),eGW(nBas)
integer :: x,HOMO,LUMO 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 LUMO energy (eV):',eGW(LUMO)*HaToeV
write(*,'(2X,A27,F15.6)') 'evGW HOMO-LUMO gap (eV):',Gap*HaToeV write(*,'(2X,A27,F15.6)') 'evGW HOMO-LUMO gap (eV):',Gap*HaToeV
write(*,*)'-------------------------------------------------------------------------------' 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)') '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(*,*)'-------------------------------------------------------------------------------'
write(*,*) write(*,*)

6
src/MCQC/print_qsGW.f90
View File

@ -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 ! 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 ! Input variables
integer,intent(in) :: nBas,nO,nSCF 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) :: 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) :: T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas)
double precision,intent(in) :: J(nBas,nBas),K(nBas,nBas),F(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(*,'(2X,A27,F15.6)') 'qsGW HOMO-LUMO gap (eV):',Gap*HaToeV
write(*,*)'-------------------------------------------' write(*,*)'-------------------------------------------'
write(*,'(2X,A27,F15.6)') 'qsGW total energy =',EqsGW + ENuc 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 exchange energy =',Ex
write(*,'(2X,A27,F15.6)') 'qsGW correlation energy =',Ec write(*,'(2X,A27,F15.6)') 'qsGW correlation energy =',Ec
write(*,'(2X,A27,F15.6)') 'RPA correlation energy =',EcRPA write(*,'(2X,A27,F15.6)') 'RPA correlation energy =',EcRPA
write(*,'(2X,A27,F15.6)') 'GM correlation energy =',EcGM
write(*,*)'-------------------------------------------' write(*,*)'-------------------------------------------'
write(*,*) write(*,*)

8
src/MCQC/qsGW.f90
View File

@ -21,7 +21,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
logical :: dRPA logical :: dRPA
integer :: nSCF,nBasSq,ispin,n_diis integer :: nSCF,nBasSq,ispin,n_diis
double precision :: EcRPA,Conv double precision :: EcRPA,EcGM,Conv
double precision,external :: trace_matrix double precision,external :: trace_matrix
double precision,allocatable :: error_diis(:,:),F_diis(:,:) double precision,allocatable :: error_diis(:,:),F_diis(:,:)
double precision,allocatable :: Omega(:,:),XpY(:,:,:),rho(:,:,:,:),rhox(:,:,:,:) 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 if(G0W) then
call self_energy_correlation(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,eHF, & 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) call renormalization_factor(SOSEX,nBas,nC,nO,nV,nR,nS,eHF,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z)
else else
call self_energy_correlation(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e, & 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) call renormalization_factor(SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z)
endif endif
@ -152,7 +152,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
! Print results ! Print results
call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) 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 ! Increment

11
src/MCQC/self_energy_correlation.f90
View File

@ -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 ! 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 ! Output variables
double precision,intent(out) :: SigC(nBas,nBas) double precision,intent(out) :: SigC(nBas,nBas)
double precision,intent(out) :: EcGM
! Initialize ! 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 j=nC+1,nO
do b=nO+1,nBas-nR do b=nO+1,nBas-nR
jb = jb + 1 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 enddo
enddo enddo
@ -64,6 +66,11 @@ subroutine self_energy_correlation(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega,rho,
enddo enddo
enddo enddo
EcGM=0d0
do i=nC+1,nO
EcGM = EcGM + SigC(i,i)
enddo
else else
! Occupied part of the correlation self-energy ! Occupied part of the correlation self-energy

27
src/MCQC/self_energy_correlation_diag.f90
View File

@ -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 ! 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 ! Output variables
double precision,intent(out) :: SigC(nBas) double precision,intent(out) :: SigC(nBas)
double precision,intent(out) :: EcGM
! Initialize ! 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 j=nC+1,nO
do b=nO+1,nBas-nR do b=nO+1,nBas-nR
jb = jb + 1 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 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 enddo
enddo enddo
EcGM=0d0
do i=nC+1,nO
EcGM = EcGM + SigC(i)
enddo
else else
@ -98,6 +105,22 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega
enddo enddo
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 if(SOSEX) then
! SOSEX: occupied part of the correlation self-energy ! SOSEX: occupied part of the correlation self-energy