4
1
mirror of https://github.com/pfloos/quack synced 2025-05-06 23:34:42 +02:00

fix cRHF print

This commit is contained in:
Loris Burth 2025-03-21 11:15:59 +01:00
parent ef50c07685
commit 8c7b3ad994
3 changed files with 15 additions and 8 deletions

View File

@ -27,7 +27,7 @@ subroutine print_cRG0W0(nBas,nO,eHF,ENuc,ERHF,Re_SigC,Im_SigC,Re_Z,Im_Z,Re_eGW,I
index_homo = maxloc(Re_eGW(1:nO),1)
Re_eHOMO = Re_eGW(index_homo)
Im_eHOMO = Im_eGW(index_homo)
index_lumo = minloc(Re_eGW(1:nO),1)
index_lumo = minloc(Re_eGW(nO+1:nBas),1)
Re_eLUMO = Re_eGW(index_lumo)
Im_eLUMO = Im_eGW(index_lumo)
Re_Gap = Re_eLUMO-Re_eHOMO
@ -56,8 +56,8 @@ subroutine print_cRG0W0(nBas,nO,eHF,ENuc,ERHF,Re_SigC,Im_SigC,Re_Z,Im_Z,Re_eGW,I
end do
write(*,*)'---------------------------------------------------------------------------------------------------'
write(*,'(2X,A60,F15.6,A5,F15.6,A3)') 'cG0W0@RHF HOMO energy = ',Re_eHOMO*HaToeV,' + i*',Im_eHOMO*HaToeV,' eV'
write(*,'(2X,A60,F15.6,A5,F15.6,A3)') 'cG0W0@RHF LUMO energy = ',Re_eLUMO*HaToeV,' + i*',Im_eHOMO*HaToeV,' eV'
write(*,'(2X,A60,F15.6,A5,F15.6,A3)') 'cG0W0@RHF HOMO-LUMO gap = ',Re_Gap*HaToeV,' + i*',Im_eHOMO*HaToeV,' eV'
write(*,'(2X,A60,F15.6,A5,F15.6,A3)') 'cG0W0@RHF LUMO energy = ',Re_eLUMO*HaToeV,' + i*',Im_eLUMO*HaToeV,' eV'
write(*,'(2X,A60,F15.6,A5,F15.6,A3)') 'cG0W0@RHF HOMO-LUMO gap = ',Re_Gap*HaToeV,' + i*',Im_Gap*HaToeV,' eV'
write(*,*)'---------------------------------------------------------------------------------------------------'
write(*,'(2X,A60,F15.6,A3)') 'phRPA@cG0W0@RHF total energy = ',ENuc + ERHF + EcRPA,' au'
write(*,'(2X,A60,F15.6,A3)') 'phRPA@cG0W0@RHF correlation energy = ',EcRPA,' au'

View File

@ -34,6 +34,7 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,ENuc, &
complex*16 :: EV
complex*16 :: EJ
complex*16 :: EK
complex*16 :: EW
double precision :: Conv
double precision :: rcond
@ -101,7 +102,7 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,ENuc, &
write(*,*)
write(*,*)'-------------------------------------------------------------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A36,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X)') &
'|','#','|','E(RHF)','|','EJ(RHF)','|','EK(RHF)','|','Conv','|'
'|','#','|','E(RHF)','|','RE(EJ(RHF))','|','Re(EK(RHF))','|','Conv','|'
write(*,*)'-------------------------------------------------------------------------------------------------'
do while(Conv > thresh .and. nSCF < maxSCF)
@ -124,10 +125,14 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,ENuc, &
! Kinetic energy
ET = cmplx(trace_matrix(nBas,real(matmul(P,T))),trace_matrix(nBas,aimag(matmul(P,T))),kind=8)
! Potential energy
EV = cmplx(trace_matrix(nBas,real(matmul(P,V))),trace_matrix(nBas,aimag(matmul(P,V))),kind=8)
! CAP energy
EW = cmplx(trace_matrix(nBas,real(matmul(P,CAP))),trace_matrix(nBas,aimag(matmul(P,CAP))),kind=8)
! Hartree energy
@ -141,7 +146,7 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,ENuc, &
! Total energy
ERHF = ET + EV + EJ + EK
ERHF = ET + EV+ EW + EJ + EK
! DIIS extrapolation !
@ -188,7 +193,7 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,ENuc, &
end if
call print_cRHF(nBas,nBas,nO,eHF,C,ENuc,ET,EV,EJ,EK,ERHF)
call print_cRHF(nBas,nBas,nO,eHF,C,ENuc,ET,EV,EW,EJ,EK,ERHF)
! Testing zone

View File

@ -1,7 +1,7 @@
! ---
subroutine print_cRHF(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, ERHF)
subroutine print_cRHF(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EW,EJ, EK, ERHF)
! Print one-electron energies and other stuff for G0W0
@ -19,6 +19,7 @@ subroutine print_cRHF(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, ERHF)
complex*16,intent(in) :: EJ
complex*16,intent(in) :: EK
complex*16,intent(in) :: EV
complex*16,intent(in) :: EW
complex*16,intent(in) :: ERHF
! Local variables
@ -47,7 +48,8 @@ subroutine print_cRHF(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, ERHF)
write(*,'(A69)') '---------------------------------------------------------'
write(*,'(A33,1X,F16.10,1X,A1,F16.10,A1,A3)') ' One-electron energy = ',real(ET + EV),'+',aimag(ET+EV),'i',' au'
write(*,'(A33,1X,F16.10,1X,A1,F16.10,A1,A3)') ' Kinetic energy = ',real(ET),'+',aimag(ET),'i',' au'
write(*,'(A33,1X,F16.10,1X,A1,F16.10,A1,A3)') ' Potential energy = ',real(EV),'+',aimag(Ev),'i',' au'
write(*,'(A33,1X,F16.10,1X,A1,F16.10,A1,A3)') ' Potential energy = ',real(EV),'+',aimag(EV),'i',' au'
write(*,'(A33,1X,F16.10,1X,A1,F16.10,A1,A3)') ' CAP energy = ',real(EW),'+',aimag(EW),'i',' au'
write(*,'(A69)') '---------------------------------------------------'
write(*,'(A33,1X,F16.10,1X,A1,F16.10,A1,A3)') ' Two-electron energy = ',real(EJ + EK),'+',aimag(EJ+EK),'i',' au'
write(*,'(A33,1X,F16.10,1X,A1,F16.10,A1,A3)') ' Hartree energy = ',real(EJ),'+',aimag(EJ),'i',' au'