From 0ed01dc905aa60697527eff3d3630ff692446907 Mon Sep 17 00:00:00 2001 From: Loris Burth Date: Tue, 25 Mar 2025 15:36:04 +0100 Subject: [PATCH] added complex G0F2 --- src/GF/cRG0F2.f90 | 153 ++++++++++++++++++++++++++++++ src/GF/cRGF2_self_energy_diag.f90 | 88 +++++++++++++++++ src/GF/print_cRG0F2.f90 | 58 +++++++++++ 3 files changed, 299 insertions(+) create mode 100644 src/GF/cRG0F2.f90 create mode 100644 src/GF/cRGF2_self_energy_diag.f90 create mode 100644 src/GF/print_cRG0F2.f90 diff --git a/src/GF/cRG0F2.f90 b/src/GF/cRG0F2.f90 new file mode 100644 index 0000000..2cd2c29 --- /dev/null +++ b/src/GF/cRG0F2.f90 @@ -0,0 +1,153 @@ +subroutine cRG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, & + nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,CAP,dipole_int,eHF) + +! Perform a one-shot second-order Green function calculation + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: dotest + + logical,intent(in) :: dophBSE + logical,intent(in) :: doppBSE + logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA + logical,intent(in) :: singlet + logical,intent(in) :: triplet + logical,intent(in) :: linearize + double precision,intent(in) :: eta + logical,intent(in) :: regularize + + integer,intent(in) :: nBas + integer,intent(in) :: nOrb + integer,intent(in) :: nO + integer,intent(in) :: nC + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + double precision,intent(in) :: ENuc + double precision,intent(in) :: ERHF + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart) + double precision,intent(in) :: CAP(nOrb,nOrb) + +! Local variables + + integer :: p + double precision :: Ec + double precision :: EcBSE(nspin) + double precision,allocatable :: Re_SigC(:) + double precision,allocatable :: Im_SigC(:) + double precision,allocatable :: Re_Z(:) + double precision,allocatable :: Im_Z(:) + double precision,allocatable :: Re_eGFlin(:) + double precision, allocatable :: Im_eGFlin(:) + double precision,allocatable :: Re_eGF(:) + double precision,allocatable :: Im_eGF(:) + double precision, allocatable :: e_CAP(:) + + ! Hello world + + write(*,*) + write(*,*)'*******************************' + write(*,*)'* Restricted G0F2 Calculation *' + write(*,*)'*******************************' + write(*,*) + +! Memory allocation + + allocate(Re_SigC(nOrb),Im_SigC(nOrb), Re_Z(nOrb),Im_Z(nOrb),& + Re_eGFlin(nOrb),Im_eGFlin(nOrb), Re_eGF(nOrb),Im_eGF(nOrb),e_CAP(nOrb)) + do p = 1, nOrb + e_CAP(p) = CAP(p,p) + end do + +! Frequency-dependent second-order contribution + + if(regularize) then + write(*,*) "Regularisation not implemented (yet)" + ! call RGF2_reg_self_energy_diag(eta,nOrb,nC,nO,nV,nR,eHF,ERI,SigC,Z) + + else + + call RGF2_self_energy_diag(eta,nOrb,nC,nO,nV,nR,eHF,ERI,Re_SigC,Im_SigC,Re_Z,Im_Z,e_CAP) + + end if + + Re_eGFlin(:) = eHF(:) + Re_Z(:)*Re_SigC(:) - Im_Z(:)*Im_SigC(:) + Im_eGFlin(:) = e_CAP(:) + Re_Z(:)*Im_SigC(:) + Im_Z(:)*Re_SigC(:) + + if(linearize) then + + write(*,*) '*** Quasiparticle energies obtained by linearization ***' + + Re_eGF(:) = Re_eGFlin(:) + Im_eGF(:) = Im_eGFlin(:) + + else + + write(*,*) ' *** Quasiparticle energies obtained by root search *** NOT IMPLEMEMTED YET' + write(*,*) + + !call RGF2_QP_graph(eta,nOrb,nC,nO,nV,nR,eHF,ERI,eGFlin,eHF,eGF,Z) + + end if + + ! Print results + +! call RMP2(.false.,regularize,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eGF,Ec) + call print_cRG0F2(nOrb,nO,eHF,e_CAP,Re_SigC,Im_SigC,Re_eGF,Im_eGF,Re_Z,Im_Z,ENuc,ERHF,Ec) + +! Perform BSE@GF2 calculation +! +! if(dophBSE) then +! +! call RGF2_phBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE) +! +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy (singlet) =',EcBSE(1) +! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy (triplet) =',EcBSE(2) +! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy =',sum(EcBSE) +! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 total energy =',ENuc + ERHF + sum(EcBSE) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) +! +! end if +! +!! Perform ppBSE@GF2 calculation +! +! if(doppBSE) then +! +! call RGF2_ppBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE) +! +! EcBSE(2) = 3d0*EcBSE(2) +! +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 correlation energy (singlet) =',EcBSE(1),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 correlation energy (triplet) =',EcBSE(2),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 correlation energy =',sum(EcBSE),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 total energy =',ENuc + ERHF + sum(EcBSE),' au' +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) +! +! end if +! +!! Testing zone +! +! if(dotest) then +! +! call dump_test_value('R','G0F2 correlation energy',Ec) +! call dump_test_value('R','G0F2 HOMO energy',eGF(nO)) +! call dump_test_value('R','G0F2 LUMO energy',eGF(nO+1)) +! +! end if +! +! deallocate(SigC, Z, eGFlin, eGF) +! +end subroutine diff --git a/src/GF/cRGF2_self_energy_diag.f90 b/src/GF/cRGF2_self_energy_diag.f90 new file mode 100644 index 0000000..1d8b677 --- /dev/null +++ b/src/GF/cRGF2_self_energy_diag.f90 @@ -0,0 +1,88 @@ +subroutine cRGF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,e,ERI,Re_SigC,Im_SigC,Re_Z,Im_Z,e_cap) + +! Compute diagonal part of the GF2 self-energy and its renormalization factor + + implicit none + include 'parameters.h' + +! Input variables + + double precision,intent(in) :: eta + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: e_cap(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: i,j,a,b + integer :: p + double precision :: eps + double precision :: eta_tilde + double precision :: num + double precision,allocatable :: Re_DS(:) + double precision,allocatable :: Im_DS(:) + +! Output variables + + double precision,intent(out) :: Re_SigC(nBas) + double precision,intent(out) :: Im_SigC(nBas) + double precision,intent(out) :: Re_Z(nBas) + double precision,intent(out) :: Im_Z(nBas) + +! Initialize + allocate(Re_DS(nBas),Im_DS(nBas)) + Re_SigC(:) = 0d0 + Im_SigC(:) = 0d0 + Re_DS(:) = 0d0 + Im_DS(:) = 0d0 + + +! Compute GF2 self-energy + + do p=nC+1,nBas-nR + do i=nC+1,nO + do j=nC+1,nO + do a=nO+1,nBas-nR + + eps = e(p) + e(a) - e(i) - e(j) + eta_tilde = eta - e_cap(p) + e_cap(i) + e_cap(a) - e_cap(j) + num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j) + + Re_SigC(p) = Re_SigC(p) + num*eps/(eps**2 + eta_tilde**2) + Im_SigC(p) = Im_SigC(p) - num*eta_tilde/(eps**2 + eta_tilde**2) + Re_DS(p) = Re_DS(p) - num*(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2 + Im_DS(p) = Im_DS(p) + 2*num*eta_tilde*eps/(eps**2 + eta_tilde**2)**2 + + end do + end do + end do + end do + + do p=nC+1,nBas-nR + do i=nC+1,nO + do a=nO+1,nBas-nR + do b=nO+1,nBas-nR + + eps = e(p) + e(i) - e(a) - e(b) + eta_tilde = eta + e_cap(p) - e_cap(a) - e_cap(b) + e_cap(i) + num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b) + + Re_SigC(p) = Re_SigC(p) + num*eps/(eps**2 + eta_tilde**2) + Im_SigC(p) = Im_SigC(p) - num*eta_tilde/(eps**2 + eta_tilde**2) + Re_DS(p) = Re_DS(p) - num*(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2 + Im_DS(p) = Im_DS(p) + 2*num*eta_tilde*eps/(eps**2 + eta_tilde**2)**2 + + end do + end do + end do + end do + + Re_Z(:) = (1d0-Re_DS(:))/((1d0 - Re_DS(:))**2 + Im_DS(:)**2) + Im_Z(:) = Im_DS(:)/((1d0 - Re_DS(:))**2 + Im_DS(:)**2) + deallocate(Re_DS,Im_DS) +end subroutine diff --git a/src/GF/print_cRG0F2.f90 b/src/GF/print_cRG0F2.f90 new file mode 100644 index 0000000..278d400 --- /dev/null +++ b/src/GF/print_cRG0F2.f90 @@ -0,0 +1,58 @@ +subroutine print_cRG0F2(nBas,nO,eHF,e_cap,Re_Sig,Im_Sig,Re_eGF,Im_eGF,Re_Z,Im_Z,ENuc,ERHF,Ec) + +! Print one-electron energies and other stuff for G0F2 + + implicit none + include 'parameters.h' + + integer,intent(in) :: nBas + integer,intent(in) :: nO + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: e_cap(nBas) + double precision,intent(in) :: Re_Sig(nBas) + double precision,intent(in) :: Im_Sig(nBas) + double precision,intent(in) :: Re_eGF(nBas) + double precision,intent(in) :: Im_eGF(nBas) + double precision,intent(in) :: Re_Z(nBas) + double precision,intent(in) :: Im_Z(nBas) + double precision,intent(in) :: ENuc + double precision,intent(in) :: ERHF + double precision,intent(in) :: Ec + + integer :: p + integer :: HOMO + integer :: LUMO + double precision :: Gap + + + +! Dump results + + write(*,*)'-------------------------------------------------------------------------------' + write(*,*)' One-shot G0F2 calculation' + 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)','|','Re(Sig_GF2) (eV)','|','Re(Z)','|','Re(e_GF2) (eV)','|' + write(*,*)'-------------------------------------------------------------------------------' + + do p=1,nBas + write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & + '|',p,'|',eHF(p)*HaToeV,'|',Re_Sig(p)*HaToeV,'|',Re_Z(p),'|',Re_eGF(p)*HaToeV,'|' + end do + + 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)') & + '|','#','|','Im(e_HF) (eV)','|','Im(Sig_GF2) (eV)','|','Im(Z)','|','Im(e_GF2) (eV)','|' + write(*,*)'-------------------------------------------------------------------------------' + + do p=1,nBas + write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & + '|',p,'|',e_cap(p)*HaToeV,'|',Im_Sig(p)*HaToeV,'|',Im_Z(p),'|',Im_eGF(p)*HaToeV,'|' + end do + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A60,F15.6,A3)') 'G0F2 total energy =',ENuc + ERHF + Ec,' au' + write(*,'(2X,A60,F15.6,A3)') 'G0F2 correlation energy =',Ec,' au' + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + +end subroutine