From 907bb2aa782e51a908902f796b5f57b03e593761 Mon Sep 17 00:00:00 2001 From: EnzoMonino Date: Wed, 16 Feb 2022 13:29:22 +0100 Subject: [PATCH] evUGT --- src/GT/evUGT.f90 | 342 +++++++++++++++++++++++++++++++++++++++++ src/GT/print_evUGT.f90 | 69 +++++++++ src/QuAcK/QuAcK.f90 | 7 +- 3 files changed, 417 insertions(+), 1 deletion(-) create mode 100644 src/GT/evUGT.f90 create mode 100644 src/GT/print_evUGT.f90 diff --git a/src/GT/evUGT.f90 b/src/GT/evUGT.f90 new file mode 100644 index 0000000..5e3866b --- /dev/null +++ b/src/GT/evUGT.f90 @@ -0,0 +1,342 @@ +subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & + TDA_T,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,& + eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,ERI_AO,ERI_aaaa, & + ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF, & + Vxc,eG0T0) + +! Perform one-shot calculation with a T-matrix self-energy (G0T0) + + implicit none + include 'parameters.h' + +! Input variables + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + logical,intent(in) :: doACFDT + logical,intent(in) :: exchange_kernel + logical,intent(in) :: doXBS + logical,intent(in) :: BSE + logical,intent(in) :: TDA_T + logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA + logical,intent(in) :: evDyn + logical,intent(in) :: spin_conserved + logical,intent(in) :: spin_flip + double precision,intent(in) :: eta + logical,intent(in) :: regularize + + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nV(nspin) + integer,intent(in) :: nR(nspin) + integer,intent(in) :: nS(nspin) + double precision,intent(in) :: ENuc + double precision,intent(in) :: EUHF + double precision,intent(in) :: Vxc(nBas,nspin) + double precision,intent(in) :: eHF(nBas,nspin) + double precision,intent(in) :: cHF(nBas,nBas,nspin) + double precision,intent(in) :: PHF(nBas,nBas,nspin) + double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart) + double precision,intent(in) :: eG0T0(nBas,nspin) + +! Local variables + integer :: nSCF + integer :: n_diis + double precision :: rcond(nspin) + double precision :: Conv + integer :: ispin,is + integer :: iblock + integer :: nH_sc,nH_sf,nHaa,nHab,nHbb + integer :: nP_sc,nP_sf,nPaa,nPab,nPbb + double precision :: EcRPA(nspin),Ecaa,Ecbb + double precision :: EcBSE(nspin) + double precision :: EcAC(nspin) + double precision :: EcGM + double precision,allocatable :: Omega1ab(:),Omega1aa(:),Omega1bb(:) + double precision,allocatable :: X1ab(:,:),X1aa(:,:),X1bb(:,:) + double precision,allocatable :: Y1ab(:,:),Y1aa(:,:),Y1bb(:,:) + double precision,allocatable :: rho1ab(:,:,:),rho1aa(:,:,:),rho1bb(:,:,:) + double precision,allocatable :: Omega2ab(:),Omega2aa(:),Omega2bb(:) + double precision,allocatable :: X2ab(:,:),X2aa(:,:),X2bb(:,:) + double precision,allocatable :: Y2ab(:,:),Y2aa(:,:),Y2bb(:,:) + double precision,allocatable :: rho2ab(:,:,:),rho2aa(:,:,:),rho2bb(:,:,:) + double precision,allocatable :: SigX(:,:) + double precision,allocatable :: SigT(:,:) + double precision,allocatable :: Z(:,:) + double precision,allocatable :: eGT(:,:) + double precision,allocatable :: eOld(:,:) + double precision,allocatable :: error_diis(:,:,:) + double precision,allocatable :: e_diis(:,:,:) + + +! Hello world + + write(*,*) + write(*,*)'************************************************' + write(*,*)'| One-shot G0T0 calculation |' + write(*,*)'| *** Unrestricted version *** |' + write(*,*)'************************************************' + write(*,*) + +! Dimensions of the pp-URPA linear reponse matrices + + nPaa = nV(1)*(nV(1)-1)/2 + nPbb = nV(2)*(nV(2)-1)/2 + + nHaa = nO(1)*(nO(1)-1)/2; + nHbb = nO(2)*(nO(2)-1)/2; + + nPab = nV(1)*nV(2) + nHab = nO(1)*nO(2) + + nP_sc = nPab + nH_sc = nHab + + nP_sf = nPaa + nPbb + nH_sf = nHaa + nHbb + + +! Memory allocation + + allocate(Omega1ab(nPab),X1ab(nPab,nPab),Y1ab(nHab,nPab), & + Omega2ab(nHab),X2ab(nPab,nHab),Y2ab(nHab,nHab), & + rho1ab(nBas,nBas,nPab),rho2ab(nBas,nBas,nHab), & + Omega1aa(nPaa),X1aa(nPaa,nPaa),Y1aa(nHaa,nPaa), & + Omega2aa(nHaa),X2aa(nPaa,nHaa),Y2aa(nHaa,nHaa), & + rho1aa(nBas,nBas,nPaa),rho2aa(nBas,nBas,nHaa), & + Omega1bb(nPbb),X1bb(nPbb,nPbb),Y1bb(nHbb,nPbb), & + Omega2bb(nPbb),X2bb(nPbb,nPbb),Y2bb(nHbb,nPbb), & + rho1bb(nBas,nBas,nPbb),rho2bb(nBas,nBas,nHbb), & + SigX(nBas,nspin),SigT(nBas,nspin),Z(nBas,nspin), & + eGT(nBas,nspin),eOld(nBas,nspin),error_diis(nBas,max_diis,nspin), & + e_diis(nBas,max_diis,nspin)) + +!---------------------------------------------- +! Compute the exchange part of the self-energy +!---------------------------------------------- + + do is=1,nspin + call self_energy_exchange_diag(nBas,cHF(:,:,is),PHF(:,:,is),ERI_AO, & + SigX(:,is)) + end do + +!Initialization + + nSCF = 0 + n_diis = 0 + Conv = 1d0 + e_diis(:,:,:) = 0d0 + error_diis(:,:,:) = 0d0 + eGT(:,:) = eG0T0(:,:) + eOld(:,:) = eGT(:,:) + Z(:,:) = 1d0 + rcond(:) = 0d0 + +!------------------------------------------------------------------------ +! Main loop +!------------------------------------------------------------------------ + + do while(Conv > thresh .and. nSCF <= maxSCF) + +!---------------------------------------------- +! alpha-beta block +!---------------------------------------------- + + ispin = 1 + iblock = 3 +! iblock = 1 + +! Compute linear response + + call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & + nPab,nHaa,nHab,nHbb,nHab,1d0,eHF,ERI_aaaa, & + 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(:)) + +!---------------------------------------------- +! alpha-alpha block +!---------------------------------------------- + + ispin = 2 + iblock = 4 + +! Compute linear response + + call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & + nPaa,nHaa,nHab,nHbb,nHaa,1d0,eHF,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Omega1aa,X1aa,Y1aa, & + 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(:)) + +!---------------------------------------------- +! beta-beta block +!---------------------------------------------- + + ispin = 2 + iblock = 7 + +! Compute linear response + + call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & + nPbb,nHaa,nHab,nHbb,nHbb,1d0,eHF,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Omega1bb,X1bb,Y1bb, & + 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(:)) + +!---------------------------------------------- +! Compute T-matrix version of the self-energy +!---------------------------------------------- + + EcGM = 0d0 + SigT(:,:) = 0d0 + Z(:,:) = 0d0 + +!alpha-beta block + + iblock = 3 + + call unrestricted_excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nHab,nPab, & + ERI_aaaa,ERI_aabb,ERI_bbbb,X1ab,Y1ab, & + rho1ab,X2ab,Y2ab,rho2ab) +!alpha-alpha block + + iblock = 4 + + call unrestricted_excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nHaa,nPaa, & + ERI_aaaa,ERI_aabb,ERI_bbbb,X1aa,Y1aa, & + rho1aa,X2aa,Y2aa,rho2aa) + +!beta-beta block + + iblock = 7 + + call unrestricted_excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nHbb,nPbb, & + ERI_aaaa,ERI_aabb,ERI_bbbb,X1bb,Y1bb, & + rho1bb,X2bb,Y2bb,rho2bb) + + call unrestricted_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nPaa,& + nPab,nPbb,eHF,Omega1aa,Omega1ab,Omega1bb,& + rho1aa,rho1ab,rho1bb,Omega2aa,Omega2ab,& + Omega2bb,rho2aa,rho2ab,rho2bb,EcGM,SigT) + + call unrestricted_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,& + nPaa,nPab,nPbb,eHF,Omega1aa,Omega1ab,& + Omega1bb,rho1aa,rho1ab,rho1bb, & + Omega2aa,Omega2ab,Omega2bb,rho2aa, & + rho2ab,rho2bb,Z) + + + Z(:,:) = 1d0/(1d0 - Z(:,:)) + +!---------------------------------------------- +! Solve the quasi-particle equation +!---------------------------------------------- + + eGT(:,:) = eHF(:,:) + SigX(:,:) + SigT(:,:) - Vxc(:,:) + +! Convergence criteria + + Conv = maxval(abs(eGT(:,:) - eOld(:,:))) + +!---------------------------------------------- +! Dump results +!---------------------------------------------- + +! Compute the ppRPA correlation energy + +!alpha-beta block + + 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)) + +!alpha-alpha block + + 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)) + + Ecaa = EcRPA(2) + +!beta-beta block + + 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)) + + 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) + + ! DIIS extrapolation + + n_diis = min(n_diis+1,max_diis) + do is=1,nspin + call DIIS_extrapolation(rcond(ispin),nBas,nBas,n_diis, & + error_diis(:,1:n_diis,is), & + e_diis(:,1:n_diis,is),eGT(:,is)-eOld(:,is), & + eGT(:,is)) + end do + + ! Reset DIIS if required + + if(minval(rcond(:)) < 1d-15) n_diis = 0 + +! Save quasiparticles energy for next cycle + + eOld(:,:) = eGT(:,:) + + ! Increment + + nSCF = nSCF + 1 + + enddo + +!------------------------------------------------------------------------ +! End main loop +!------------------------------------------------------------------------ + +! Free memory + + deallocate(Omega1ab,X1ab,Y1ab,Omega2ab,X2ab,Y2ab,rho1ab,rho2ab, & + Omega1aa,X1aa,Y1aa,Omega2aa,X2aa,Y2aa,rho1aa,rho2aa, & + Omega1bb,X1bb,Y1bb,Omega2bb,X2bb,Y2bb,rho1bb,rho2bb) + +end subroutine evUGT diff --git a/src/GT/print_evUGT.f90 b/src/GT/print_evUGT.f90 new file mode 100644 index 0000000..25f87f6 --- /dev/null +++ b/src/GT/print_evUGT.f90 @@ -0,0 +1,69 @@ +subroutine print_evUGT(nBas,nO,eHF,ENuc,EUHF,SigT,Z,eGT,EcGM,EcRPA) + +! Print one-electron energies and other stuff for UG0T0 + + implicit none + include 'parameters.h' + + integer,intent(in) :: nBas + integer,intent(in) :: nO(nspin) + double precision,intent(in) :: ENuc + double precision,intent(in) :: EUHF + double precision,intent(in) :: EcGM + double precision,intent(in) :: EcRPA(nspin) + double precision,intent(in) :: eHF(nBas,nspin) + double precision,intent(in) :: SigT(nBas,nspin) + double precision,intent(in) :: Z(nBas,nspin) + double precision,intent(in) :: eGT(nBas,nspin) + + integer :: p + integer :: ispin + double precision :: HOMO(nspin) + double precision :: LUMO(nspin) + double precision :: Gap(nspin) + +! HOMO and LUMO + do ispin=1,nspin + if(nO(ispin) > 0) then + HOMO(ispin) = eGT(nO(ispin),ispin) + LUMO(ispin) = eGT(nO(ispin)+1,ispin) + Gap(ispin) = LUMO(ispin) - HOMO(ispin) + else + HOMO(ispin) = 0d0 + LUMO(ispin) = eGT(1,ispin) + Gap(ispin) = 0d0 + end if + end do + +! Dump results + + write(*,*)'-------------------------------------------------------------------------------' + write(*,*)' Unrestricted one-shot G0T0 calculation (T-matrix self-energy) ' + 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)','|' + write(*,*)'-------------------------------------------------------------------------------' + + do p=1,nBas + write(*,'(A1,I3,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1)') & + '|',p,'|',eHF(p,1)*HaToeV,eHF(p,2)*HaToeV,'|',SigT(p,1)*HaToeV,SigT(p,2)*HaToeV,'|', & + Z(p,1),Z(p,2),'|',eGT(p,1)*HaToeV,eGT(p,2)*HaToeV,'|' + 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(*,*)'-------------------------------------------------------------------------------' + 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(*,*) + +end subroutine print_evUGT + + diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 903fcbf..ff07b53 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -1188,7 +1188,12 @@ program QuAcK if(unrestricted) then - print*,'!!! evGT NYI at the unrestricted level !!!' + !print*,'!!! evGT NYI at the unrestricted level !!!' + call evUGT(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS, & + BSE,TDA_T,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,& + eta_GT,regGT,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,ERI_AO, & + ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_aa, & + dipole_int_bb,PHF,cHF,eHF,Vxc,eG0T0) else