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

evUGT not working yet

Browse Source
This commit is contained in:
EnzoMonino 2022-02-16 22:24:43 +01:00
parent a878fc6e3b
commit 31e9911797
3 changed files with 45 additions and 45 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
input/options
View File

@ -1,5 +1,5 @@
# HF: maxSCF thresh DIIS n_diis guess_type ortho_type mix_guess level_shift stability # HF: maxSCF thresh DIIS n_diis guess_type ortho_type mix_guess level_shift stability
256 0.0000001 T 5 2 1 F 0.0 F 256 0.0000001 T 5 2 1 T 0.0 F
# MP: # MP:
# CC: maxSCF thresh DIIS n_diis # CC: maxSCF thresh DIIS n_diis

54
src/GT/evUGT.f90
View File

@ -81,8 +81,7 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
write(*,*) write(*,*)
write(*,*)'************************************************' write(*,*)'************************************************'
write(*,*)'| One-shot G0T0 calculation |' write(*,*)'| Self-consistent evUGT calculation |'
write(*,*)'| *** Unrestricted version *** |'
write(*,*)'************************************************' write(*,*)'************************************************'
write(*,*) write(*,*)
@ -161,10 +160,7 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
ERI_aabb,ERI_bbbb,Omega1ab,X1ab,Y1ab, & ERI_aabb,ERI_bbbb,Omega1ab,X1ab,Y1ab, &
Omega2ab,X2ab,Y2ab,EcRPA(ispin)) Omega2ab,X2ab,Y2ab,EcRPA(ispin))
! EcRPA(ispin) = 1d0*EcRPA(ispin) ! EcRPA(ispin) = 1d0*EcRPA(ispin)
call print_excitation('pp-RPA (N+2)',iblock,nPab,Omega1ab(:))
call print_excitation('pp-RPA (N-2)',iblock,nHab,Omega2ab(:))
!---------------------------------------------- !----------------------------------------------
! alpha-alpha block ! alpha-alpha block
@ -181,10 +177,7 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
Omega2aa,X2aa,Y2aa,EcRPA(ispin)) Omega2aa,X2aa,Y2aa,EcRPA(ispin))
! EcRPA(ispin) = 2d0*EcRPA(ispin) ! EcRPA(ispin) = 2d0*EcRPA(ispin)
! EcRPA(ispin) = 3d0*EcRPA(ispin) ! EcRPA(ispin) = 3d0*EcRPA(ispin)
call print_excitation('pp-RPA (N+2)',iblock,nPaa,Omega1aa(:))
call print_excitation('pp-RPA (N-2)',iblock,nHaa,Omega2aa(:))
!---------------------------------------------- !----------------------------------------------
! beta-beta block ! beta-beta block
@ -201,10 +194,7 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
Omega2bb,X2bb,Y2bb,EcRPA(ispin)) Omega2bb,X2bb,Y2bb,EcRPA(ispin))
! EcRPA(ispin) = 2d0*EcRPA(ispin) ! EcRPA(ispin) = 2d0*EcRPA(ispin)
! EcRPA(ispin) = 3d0*EcRPA(ispin) ! EcRPA(ispin) = 3d0*EcRPA(ispin)
call print_excitation('pp-RPA (N+2)',iblock,nPbb,Omega1bb(:))
call print_excitation('pp-RPA (N-2)',iblock,nHbb,Omega2bb(:))
!---------------------------------------------- !----------------------------------------------
! Compute T-matrix version of the self-energy ! Compute T-matrix version of the self-energy
@ -269,41 +259,41 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
!alpha-beta block !alpha-beta block
ispin = 1 ! ispin = 1
iblock = 3 ! iblock = 3
call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & ! call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, &
nPab,nHaa,nHab,nHbb,nHab,1d0,eG0T0,ERI_aaaa, & ! nPab,nHaa,nHab,nHbb,nHab,1d0,eG0T0,ERI_aaaa, &
ERI_aabb,ERI_bbbb,Omega1ab,X1ab,Y1ab, & ! ERI_aabb,ERI_bbbb,Omega1ab,X1ab,Y1ab, &
Omega2ab,X2ab,Y2ab,EcRPA(ispin)) ! Omega2ab,X2ab,Y2ab,EcRPA(ispin))
!alpha-alpha block !alpha-alpha block
ispin = 2 ! ispin = 2
iblock = 4 ! iblock = 4
call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & ! call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, &
nPaa,nHaa,nHab,nHbb,nHaa,1d0,eG0T0,ERI_aaaa, & ! nPaa,nHaa,nHab,nHbb,nHaa,1d0,eG0T0,ERI_aaaa, &
ERI_aabb,ERI_bbbb,Omega1aa,X1aa,Y1aa, & ! ERI_aabb,ERI_bbbb,Omega1aa,X1aa,Y1aa, &
Omega2aa,X2aa,Y2aa,EcRPA(ispin)) ! Omega2aa,X2aa,Y2aa,EcRPA(ispin))
Ecaa = EcRPA(2) Ecaa = EcRPA(2)
!beta-beta block !beta-beta block
iblock = 7 ! iblock = 7
call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & ! call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, &
nPbb,nHaa,nHab,nHbb,nHbb,1d0,eG0T0,ERI_aaaa, & ! nPbb,nHaa,nHab,nHbb,nHbb,1d0,eG0T0,ERI_aaaa, &
ERI_aabb,ERI_bbbb,Omega1bb,X1bb,Y1bb, & ! ERI_aabb,ERI_bbbb,Omega1bb,X1bb,Y1bb, &
Omega2bb,X2bb,Y2bb,EcRPA(ispin)) ! Omega2bb,X2bb,Y2bb,EcRPA(ispin))
Ecbb = EcRPA(2) Ecbb = EcRPA(2)
EcRPA(2) = Ecaa + Ecbb EcRPA(2) = Ecaa + Ecbb
EcRPA(1) = EcRPA(1) - EcRPA(2) EcRPA(1) = EcRPA(1) - EcRPA(2)
EcRPA(2) = 3d0*EcRPA(2) EcRPA(2) = 3d0*EcRPA(2)
call print_evUGT(nBas,nO,eHF,ENuc,EUHF,SigT,Z,eGT,EcGM,EcRPA) call print_evUGT(nBas,nO,nSCF,Conv,eHF,ENuc,EUHF,SigT,Z,eGT,EcGM,EcRPA)
! DIIS extrapolation ! DIIS extrapolation

34
src/GT/print_evUGT.f90
View File

@ -1,4 +1,4 @@
subroutine print_evUGT(nBas,nO,eHF,ENuc,EUHF,SigT,Z,eGT,EcGM,EcRPA) subroutine print_evUGT(nBas,nO,nSCF,Conv,eHF,ENuc,EUHF,SigT,Z,eGT,EcGM,EcRPA)
! Print one-electron energies and other stuff for UG0T0 ! Print one-electron energies and other stuff for UG0T0
@ -7,6 +7,8 @@ subroutine print_evUGT(nBas,nO,eHF,ENuc,EUHF,SigT,Z,eGT,EcGM,EcRPA)
integer,intent(in) :: nBas integer,intent(in) :: nBas
integer,intent(in) :: nO(nspin) integer,intent(in) :: nO(nspin)
integer,intent(in) :: nSCF
double precision,intent(in) :: Conv
double precision,intent(in) :: ENuc double precision,intent(in) :: ENuc
double precision,intent(in) :: EUHF double precision,intent(in) :: EUHF
double precision,intent(in) :: EcGM double precision,intent(in) :: EcGM
@ -36,9 +38,13 @@ subroutine print_evUGT(nBas,nO,eHF,ENuc,EUHF,SigT,Z,eGT,EcGM,EcRPA)
end do end do
! Dump results ! Dump results
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,*)' Unrestricted one-shot G0T0 calculation (T-matrix self-energy) ' if(nSCF < 10) then
write(*,'(1X,A21,I1,A1,I1,A12)')' Self-consistent evG',nSCF,'T',nSCF,' calculation'
else
write(*,'(1X,A21,I2,A1,I2,A12)')' Self-consistent evG',nSCF,'T',nSCF,' calculation'
endif
write(*,*)'-------------------------------------------------------------------------------' 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)') & write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
'|','#','|','e_HF (eV)','|','Sigma_T (eV)','|','Z','|','e_QP (eV)','|' '|','#','|','e_HF (eV)','|','Sigma_T (eV)','|','Z','|','e_QP (eV)','|'
@ -51,16 +57,20 @@ subroutine print_evUGT(nBas,nO,eHF,ENuc,EUHF,SigT,Z,eGT,EcGM,EcRPA)
enddo enddo
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F15.6,A3)') 'UG0T0 HOMO energy (eV) =',maxval(HOMO(:))*HaToeV,' eV' write(*,'(2X,A10,I3)') 'Iteration ',nSCF
write(*,'(2X,A50,F15.6,A3)') 'UG0T0 LUMO energy (eV) =',minval(LUMO(:))*HaToeV,' eV' write(*,'(2X,A14,F15.5)')'Convergence = ',Conv
write(*,'(2X,A50,F15.6,A3)') 'UG0T0 HOMO-LUMO gap (eV) =',(minval(LUMO(:))-maxval(HOMO(:)))*HaToeV,' eV'
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@UG0T0 correlation energy (singlet) =',EcRPA(1),' au' write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@UG0T0 correlation energy (triplet) =',EcRPA(2),' au' write(*,'(2X,A50,F15.6,A3)') 'evUGT HOMO energy (eV) =',maxval(HOMO(:))*HaToeV,' eV'
write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@UG0T0 correlation energy =',EcRPA(1) + EcRPA(2),' au' write(*,'(2X,A50,F15.6,A3)') 'evUGT LUMO energy (eV) =',minval(LUMO(:))*HaToeV,' eV'
write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@UG0T0 total energy =',ENuc + EUHF + EcRPA(1) + EcRPA(2),' au' write(*,'(2X,A50,F15.6,A3)') 'evUGT HOMO-LUMO gap (eV) =',(minval(LUMO(:))-maxval(HOMO(:)))*HaToeV,' eV'
write(*,'(2X,A50,F20.10,A3)') ' GM@UG0T0 correlation energy =',EcGM,' au' write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10,A3)') ' GM@UG0T0 total energy =',ENuc + EUHF + EcGM,' au' write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@evUGT correlation energy (singlet) =',EcRPA(1),' au'
write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@evUGT correlation energy (triplet) =',EcRPA(2),' au'
write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@evUGT correlation energy =',EcRPA(1) + EcRPA(2),' au'
write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@evUGT total energy =',ENuc + EUHF + EcRPA(1) + EcRPA(2),' au'
write(*,'(2X,A50,F20.10,A3)') ' GM@evUGT correlation energy =',EcGM,' au'
write(*,'(2X,A50,F20.10,A3)') ' GM@evUGT total energy =',ENuc + EUHF + EcGM,' au'
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,*) write(*,*)