From 31e99117979ba912d09f1d2c3dfef2ea76456b05 Mon Sep 17 00:00:00 2001 From: EnzoMonino Date: Wed, 16 Feb 2022 22:24:43 +0100 Subject: [PATCH] evUGT not working yet --- input/options | 2 +- src/GT/evUGT.f90 | 54 +++++++++++++++++------------------------- src/GT/print_evUGT.f90 | 34 ++++++++++++++++---------- 3 files changed, 45 insertions(+), 45 deletions(-) diff --git a/input/options b/input/options index e9b8028..0cd766f 100644 --- a/input/options +++ b/input/options @@ -1,5 +1,5 @@ # 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: # CC: maxSCF thresh DIIS n_diis diff --git a/src/GT/evUGT.f90 b/src/GT/evUGT.f90 index 5e3866b..33c1702 100644 --- a/src/GT/evUGT.f90 +++ b/src/GT/evUGT.f90 @@ -81,8 +81,7 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & write(*,*) write(*,*)'************************************************' - write(*,*)'| One-shot G0T0 calculation |' - write(*,*)'| *** Unrestricted version *** |' + write(*,*)'| Self-consistent evUGT calculation |' write(*,*)'************************************************' write(*,*) @@ -161,10 +160,7 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & ERI_aabb,ERI_bbbb,Omega1ab,X1ab,Y1ab, & Omega2ab,X2ab,Y2ab,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(:)) +! EcRPA(ispin) = 1d0*EcRPA(ispin) !---------------------------------------------- ! alpha-alpha block @@ -181,10 +177,7 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & Omega2aa,X2aa,Y2aa,EcRPA(ispin)) ! EcRPA(ispin) = 2d0*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(:)) +! EcRPA(ispin) = 3d0*EcRPA(ispin) !---------------------------------------------- ! beta-beta block @@ -201,10 +194,7 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & Omega2bb,X2bb,Y2bb,EcRPA(ispin)) ! EcRPA(ispin) = 2d0*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(:)) +! EcRPA(ispin) = 3d0*EcRPA(ispin) !---------------------------------------------- ! 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 - ispin = 1 - iblock = 3 +! ispin = 1 +! iblock = 3 - call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nPab,nHaa,nHab,nHbb,nHab,1d0,eG0T0,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Omega1ab,X1ab,Y1ab, & - Omega2ab,X2ab,Y2ab,EcRPA(ispin)) +! call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & +! nPab,nHaa,nHab,nHbb,nHab,1d0,eG0T0,ERI_aaaa, & +! ERI_aabb,ERI_bbbb,Omega1ab,X1ab,Y1ab, & +! Omega2ab,X2ab,Y2ab,EcRPA(ispin)) !alpha-alpha block - ispin = 2 - iblock = 4 +! ispin = 2 +! iblock = 4 - call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nPaa,nHaa,nHab,nHbb,nHaa,1d0,eG0T0,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Omega1aa,X1aa,Y1aa, & - Omega2aa,X2aa,Y2aa,EcRPA(ispin)) +! call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & +! nPaa,nHaa,nHab,nHbb,nHaa,1d0,eG0T0,ERI_aaaa, & +! ERI_aabb,ERI_bbbb,Omega1aa,X1aa,Y1aa, & +! Omega2aa,X2aa,Y2aa,EcRPA(ispin)) Ecaa = EcRPA(2) !beta-beta block - iblock = 7 +! iblock = 7 - call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & - nPbb,nHaa,nHab,nHbb,nHbb,1d0,eG0T0,ERI_aaaa, & - ERI_aabb,ERI_bbbb,Omega1bb,X1bb,Y1bb, & - Omega2bb,X2bb,Y2bb,EcRPA(ispin)) +! call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & +! nPbb,nHaa,nHab,nHbb,nHbb,1d0,eG0T0,ERI_aaaa, & +! ERI_aabb,ERI_bbbb,Omega1bb,X1bb,Y1bb, & +! Omega2bb,X2bb,Y2bb,EcRPA(ispin)) Ecbb = EcRPA(2) EcRPA(2) = Ecaa + Ecbb EcRPA(1) = EcRPA(1) - 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 diff --git a/src/GT/print_evUGT.f90 b/src/GT/print_evUGT.f90 index 25f87f6..135566f 100644 --- a/src/GT/print_evUGT.f90 +++ b/src/GT/print_evUGT.f90 @@ -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 @@ -7,6 +7,8 @@ subroutine print_evUGT(nBas,nO,eHF,ENuc,EUHF,SigT,Z,eGT,EcGM,EcRPA) integer,intent(in) :: nBas integer,intent(in) :: nO(nspin) + integer,intent(in) :: nSCF + double precision,intent(in) :: Conv double precision,intent(in) :: ENuc double precision,intent(in) :: EUHF double precision,intent(in) :: EcGM @@ -36,9 +38,13 @@ subroutine print_evUGT(nBas,nO,eHF,ENuc,EUHF,SigT,Z,eGT,EcGM,EcRPA) end do ! Dump results - + 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(*,'(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)','|' @@ -51,16 +57,20 @@ subroutine print_evUGT(nBas,nO,eHF,ENuc,EUHF,SigT,Z,eGT,EcGM,EcRPA) enddo write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F15.6,A3)') 'UG0T0 HOMO energy (eV) =',maxval(HOMO(:))*HaToeV,' eV' - write(*,'(2X,A50,F15.6,A3)') 'UG0T0 LUMO energy (eV) =',minval(LUMO(:))*HaToeV,' eV' - write(*,'(2X,A50,F15.6,A3)') 'UG0T0 HOMO-LUMO gap (eV) =',(minval(LUMO(:))-maxval(HOMO(:)))*HaToeV,' eV' + write(*,'(2X,A10,I3)') 'Iteration ',nSCF + write(*,'(2X,A14,F15.5)')'Convergence = ',Conv write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@UG0T0 correlation energy (singlet) =',EcRPA(1),' au' - write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@UG0T0 correlation energy (triplet) =',EcRPA(2),' au' - write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@UG0T0 correlation energy =',EcRPA(1) + EcRPA(2),' au' - write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@UG0T0 total energy =',ENuc + EUHF + EcRPA(1) + EcRPA(2),' au' - write(*,'(2X,A50,F20.10,A3)') ' GM@UG0T0 correlation energy =',EcGM,' au' - write(*,'(2X,A50,F20.10,A3)') ' GM@UG0T0 total energy =',ENuc + EUHF + EcGM,' au' + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F15.6,A3)') 'evUGT HOMO energy (eV) =',maxval(HOMO(:))*HaToeV,' eV' + write(*,'(2X,A50,F15.6,A3)') 'evUGT LUMO energy (eV) =',minval(LUMO(:))*HaToeV,' eV' + write(*,'(2X,A50,F15.6,A3)') 'evUGT HOMO-LUMO gap (eV) =',(minval(LUMO(:))-maxval(HOMO(:)))*HaToeV,' eV' + write(*,*)'-------------------------------------------------------------------------------' + 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(*,*)