diff --git a/input/methods b/input/methods index bdc5b27..2daf6f2 100644 --- a/input/methods +++ b/input/methods @@ -14,6 +14,6 @@ T F F F F # G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW T F F F F F -# G0T0 evGT qsGT ehG0T0 - T F F T +# G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh + T F F T T F # * unrestricted version available diff --git a/src/GT/ehG0T0.f90 b/src/GT/G0T0eh.f90 similarity index 96% rename from src/GT/ehG0T0.f90 rename to src/GT/G0T0eh.f90 index d1eddf3..de6e17f 100644 --- a/src/GT/ehG0T0.f90 +++ b/src/GT/G0T0eh.f90 @@ -1,4 +1,4 @@ -subroutine ehG0T0(doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_T,TDA,dBSE,dTDA,evDyn,ppBSE, & +subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_T,TDA,dBSE,dTDA,evDyn,ppBSE, & singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, & ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eGT) @@ -115,7 +115,7 @@ subroutine ehG0T0(doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_T,TDA,dBSE,dTDA,evD ! Compute spectral weights ! !--------------------------! - call ehGT_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,XmY_RPA,rhoL_RPA,rhoR_RPA) + call GTeh_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,XmY_RPA,rhoL_RPA,rhoR_RPA) !------------------------! ! Compute GW self-energy ! @@ -130,8 +130,8 @@ subroutine ehG0T0(doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_T,TDA,dBSE,dTDA,evD else - call ehGT_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rhoL_RPA,rhoR_RPA,EcGM,SigC) - call ehGT_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rhoL_RPA,rhoR_RPA,Z) + call GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rhoL_RPA,rhoR_RPA,EcGM,SigC) + call GTeh_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rhoL_RPA,rhoR_RPA,Z) end if @@ -172,7 +172,7 @@ subroutine ehG0T0(doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_T,TDA,dBSE,dTDA,evD ! Dump results ! !--------------! - call print_ehG0T0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGT,EcRPA,EcGM) + call print_G0T0eh(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGT,EcRPA,EcGM) ! Deallocate memory diff --git a/src/GT/ehGT_excitation_density.f90 b/src/GT/GTeh_excitation_density.f90 similarity index 95% rename from src/GT/ehGT_excitation_density.f90 rename to src/GT/GTeh_excitation_density.f90 index cc26650..22072bf 100644 --- a/src/GT/ehGT_excitation_density.f90 +++ b/src/GT/GTeh_excitation_density.f90 @@ -1,4 +1,4 @@ -subroutine ehGT_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,XmY,rhoL,rhoR) +subroutine GTeh_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,XmY,rhoL,rhoR) ! Compute excitation densities diff --git a/src/GT/ehGT_renormalization_factor.f90 b/src/GT/GTeh_renormalization_factor.f90 similarity index 95% rename from src/GT/ehGT_renormalization_factor.f90 rename to src/GT/GTeh_renormalization_factor.f90 index c9ccd2b..75ad573 100644 --- a/src/GT/ehGT_renormalization_factor.f90 +++ b/src/GT/GTeh_renormalization_factor.f90 @@ -1,4 +1,4 @@ -subroutine ehGT_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,Z) +subroutine GTeh_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,Z) ! Compute renormalization factor for GW diff --git a/src/GT/ehGT_self_energy_diag.f90 b/src/GT/GTeh_self_energy_diag.f90 similarity index 96% rename from src/GT/ehGT_self_energy_diag.f90 rename to src/GT/GTeh_self_energy_diag.f90 index 842fa4d..4dd666f 100644 --- a/src/GT/ehGT_self_energy_diag.f90 +++ b/src/GT/GTeh_self_energy_diag.f90 @@ -1,4 +1,4 @@ -subroutine ehGT_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,SigC) +subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,SigC) ! Compute diagonal of the correlation part of the self-energy diff --git a/src/GT/evGTeh.f90 b/src/GT/evGTeh.f90 new file mode 100644 index 0000000..7d69c8c --- /dev/null +++ b/src/GT/evGTeh.f90 @@ -0,0 +1,321 @@ +subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_T,TDA,dBSE,dTDA,evDyn,ppBSE, & + singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF, & + cHF,eHF,Vxc,eG0T0) + +! Perform self-consistent eigenvalue-only ehGT calculation + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + double precision,intent(in) :: ENuc + double precision,intent(in) :: ERHF + logical,intent(in) :: doACFDT + logical,intent(in) :: exchange_kernel + logical,intent(in) :: doXBS + logical,intent(in) :: BSE + logical,intent(in) :: BSE2 + logical,intent(in) :: TDA_T + logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA + logical,intent(in) :: evDyn + logical,intent(in) :: ppBSE + 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) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + double precision,intent(in) :: PHF(nBas,nBas) + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: cHF(nBas,nBas) + double precision,intent(in) :: Vxc(nBas) + double precision,intent(in) :: eG0T0(nBas) + double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + +! Local variables + + logical :: linear_mixing + integer :: ispin + integer :: nSCF + integer :: n_diis + integer :: i,a,jb,p + double precision :: rcond + double precision :: Conv + double precision :: EcRPA + double precision :: EcBSE(nspin) + double precision :: EcAC(nspin) + double precision :: EcppBSE(nspin) + double precision :: EcGM + double precision :: alpha + double precision :: Dpijb,Dpajb + double precision,allocatable :: error_diis(:,:) + double precision,allocatable :: e_diis(:,:) + double precision,allocatable :: eGT(:) + double precision,allocatable :: eOld(:) + double precision,allocatable :: Z(:) + double precision,allocatable :: SigX(:) + double precision,allocatable :: SigC(:) + double precision,allocatable :: OmRPA(:) + double precision,allocatable :: XpY_RPA(:,:) + double precision,allocatable :: XmY_RPA(:,:) + double precision,allocatable :: rhoL_RPA(:,:,:) + double precision,allocatable :: rhoR_RPA(:,:,:) + + double precision,allocatable :: eGTlin(:) + +! Hello world + + write(*,*) + write(*,*)'************************************************' + write(*,*)'| Self-consistent evGTeh calculation |' + write(*,*)'************************************************' + write(*,*) + +! TDA for T + + if(TDA_T) then + write(*,*) 'Tamm-Dancoff approximation for eh T-matrix!' + write(*,*) + end if + +! TDA + + if(TDA) then + write(*,*) 'Tamm-Dancoff approximation activated!' + write(*,*) + end if + +! Linear mixing + + linear_mixing = .false. + alpha = 0.2d0 + +! Memory allocation + + allocate(eGT(nBas),eOld(nBas),Z(nBas),SigX(nBas),SigC(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS), & + rhoL_RPA(nBas,nBas,nS),rhoR_RPA(nBas,nBas,nS),error_diis(nBas,max_diis),e_diis(nBas,max_diis),eGTlin(nBas)) + +! Compute the exchange part of the self-energy + + call self_energy_exchange_diag(nBas,cHF,PHF,ERI_AO,SigX) + +! Initialization + + nSCF = 0 + ispin = 3 + 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) + + ! Compute screening + + call linear_response(ispin,.false.,TDA_T,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI_MO, & + EcRPA,OmRPA,XpY_RPA,XmY_RPA) + + ! Compute spectral weights + + call GTeh_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,XmY_RPA,rhoL_RPA,rhoR_RPA) + + ! Compute correlation part of the self-energy + + if(regularize) then + +! call regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC) +! call renormalization_factor_SRG(eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z) + + else + + call GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eGT,OmRPA,rhoL_RPA,rhoR_RPA,EcGM,SigC) + call GTeh_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,eGT,OmRPA,rhoL_RPA,rhoR_RPA,Z) + + end if + + ! Solve the quasi-particle equation + + eGTlin(:) = eHF(:) + SigX(:) + SigC(:) - Vxc(:) + + ! Linearized or graphical solution? + + if(linearize) then + + write(*,*) ' *** Quasiparticle energies obtained by linearization *** ' + write(*,*) + + eGT(:) = eGTlin(:) + + else + +! write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** ' +! write(*,*) +! +! call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,OmRPA,rho_RPA,eGWlin,eGW,regularize) + + end if + + ! Convergence criteria + + Conv = maxval(abs(eGT - eOld)) + + ! Print results + + call print_evGTeh(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGT,EcRPA,EcGM) + + ! Linear mixing or DIIS extrapolation + + if(linear_mixing) then + + eGT(:) = alpha*eGT(:) + (1d0 - alpha)*eOld(:) + + else + + n_diis = min(n_diis+1,max_diis) + if(abs(rcond) > 1d-7) then + call DIIS_extrapolation(rcond,nBas,nBas,n_diis,error_diis,e_diis,eGT-eOld,eGT) + else + n_diis = 0 + end if + + end if + + ! Save quasiparticles energy for next cycle + + eOld(:) = eGT(:) + + ! Increment + + nSCF = nSCF + 1 + + end do +!------------------------------------------------------------------------ +! End main loop +!------------------------------------------------------------------------ + +! Did it actually converge? + + if(nSCF == maxSCF+1) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)' Convergence failed ' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*) + + stop + + end if + +! Deallocate memory + + deallocate(eOld,Z,SigC,OmRPA,XpY_RPA,XmY_RPA,rhoL_RPA,rhoR_RPA,error_diis,e_diis) + +! Perform BSE calculation + +! if(BSE) then + +! call Bethe_Salpeter(BSE2,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eGW,eGW,EcBSE) + +! if(exchange_kernel) then + +! EcBSE(1) = 0.5d0*EcBSE(1) +! EcBSE(2) = 1.5d0*EcBSE(2) + +! end if + +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW correlation energy (singlet) =',EcBSE(1) +! write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW correlation energy (triplet) =',EcBSE(2) +! write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW correlation energy =',EcBSE(1) + EcBSE(2) +! write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) + +! Compute the BSE correlation energy via the adiabatic connection + +! if(doACFDT) then + +! write(*,*) '------------------------------------------------------' +! write(*,*) 'Adiabatic connection version of BSE correlation energy' +! write(*,*) '------------------------------------------------------' +! write(*,*) + +! if(doXBS) then + +! write(*,*) '*** scaled screening version (XBS) ***' +! write(*,*) + +! end if + +! call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC) + +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy (singlet) =',EcAC(1) +! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy (triplet) =',EcAC(2) +! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy =',EcAC(1) + EcAC(2) +! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW total energy =',ENuc + ERHF + EcAC(1) + EcAC(2) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) + +! end if + +! end if + +! if(ppBSE) then + +! call Bethe_Salpeter_pp(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eGW,EcppBSE) + +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@G0W0 correlation energy (singlet) =',EcppBSE(1) +! write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@G0W0 correlation energy (triplet) =',3d0*EcppBSE(2) +! write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@G0W0 correlation energy =',EcppBSE(1) + 3d0*EcppBSE(2) +! write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@G0W0 total energy =',ENuc + ERHF + EcppBSE(1) + 3d0*EcppBSE(2) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) + +! nBas2 = 2*nBas +! nO2 = 2*nO +! nV2 = 2*nV +! nC2 = 2*nC +! nR2 = 2*nR +! nS2 = nO2*nV2 +! +! allocate(seHF(nBas2),seGW(nBas2),sERI(nBas2,nBas2,nBas2,nBas2)) +! +! call spatial_to_spin_MO_energy(nBas,eHF,nBas2,seHF) +! call spatial_to_spin_MO_energy(nBas,eGW,nBas2,seGW) +! call spatial_to_spin_ERI(nBas,ERI_MO,nBas2,sERI) +! +! call Bethe_Salpeter_pp_so(TDA_W,TDA,singlet,triplet,eta,nBas2,nC2,nO2,nV2,nR2,nS2,sERI,dipole_int,seHF,seGW,EcppBSE) + +! end if + +end subroutine diff --git a/src/GT/print_ehG0T0.f90 b/src/GT/print_G0T0eh.f90 similarity index 76% rename from src/GT/print_ehG0T0.f90 rename to src/GT/print_G0T0eh.f90 index 82f4db5..238cd41 100644 --- a/src/GT/print_ehG0T0.f90 +++ b/src/GT/print_G0T0eh.f90 @@ -1,4 +1,4 @@ -subroutine print_ehG0T0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGT,EcRPA,EcGM) +subroutine print_G0T0eh(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGT,EcRPA,EcGM) ! Print one-electron energies and other stuff for G0W0 @@ -27,7 +27,7 @@ subroutine print_ehG0T0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGT,EcRPA,EcGM) ! Dump results write(*,*)'-------------------------------------------------------------------------------' - write(*,*)' One-shot eh G0T0 calculation' + write(*,*)' One-shot G0T0eh 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)','|','Sig_c (eV)','|','Z','|','e_QP (eV)','|' @@ -39,14 +39,14 @@ subroutine print_ehG0T0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGT,EcRPA,EcGM) enddo write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A30,F15.6,A3)') 'ehG0T0 HOMO energy:',eGT(HOMO)*HaToeV,' eV' - write(*,'(2X,A30,F15.6,A3)') 'ehG0T0 LUMO energy:',eGT(LUMO)*HaToeV,' eV' - write(*,'(2X,A30,F15.6,A3)') 'ehG0T0 HOMO-LUMO gap :',Gap*HaToeV,' eV' + write(*,'(2X,A30,F15.6,A3)') 'G0T0eh HOMO energy:',eGT(HOMO)*HaToeV,' eV' + write(*,'(2X,A30,F15.6,A3)') 'G0T0eh LUMO energy:',eGT(LUMO)*HaToeV,' eV' + write(*,'(2X,A30,F15.6,A3)') 'G0T0eh HOMO-LUMO gap :',Gap*HaToeV,' eV' write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A30,F15.6,A3)') 'RPA@ehG0T0 total energy :',ENuc + ERHF + EcRPA,' au' - write(*,'(2X,A30,F15.6,A3)') 'RPA@ehG0T0 correlation energy:',EcRPA,' au' - write(*,'(2X,A30,F15.6,A3)') 'GM@ehG0T0 total energy :',ENuc + ERHF + EcGM,' au' - write(*,'(2X,A30,F15.6,A3)') 'GM@ehG0T0 correlation energy:',EcGM,' au' + write(*,'(2X,A30,F15.6,A3)') 'RPA@G0T0eh total energy :',ENuc + ERHF + EcRPA,' au' + write(*,'(2X,A30,F15.6,A3)') 'RPA@G0T0eh correlation energy:',EcRPA,' au' + write(*,'(2X,A30,F15.6,A3)') 'GM@G0T0eh total energy :',ENuc + ERHF + EcGM,' au' + write(*,'(2X,A30,F15.6,A3)') 'GM@G0T0eh correlation energy:',EcGM,' au' write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/GT/print_evGTeh.f90 b/src/GT/print_evGTeh.f90 new file mode 100644 index 0000000..0107962 --- /dev/null +++ b/src/GT/print_evGTeh.f90 @@ -0,0 +1,61 @@ +subroutine print_evGTeh(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGT,EcRPA,EcGM) + +! Print one-electron energies and other stuff for evGTeh + + implicit none + include 'parameters.h' + + integer,intent(in) :: nBas,nO,nSCF + double precision,intent(in) :: ENuc + double precision,intent(in) :: ERHF + double precision,intent(in) :: Conv + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: SigC(nBas) + double precision,intent(in) :: Z(nBas) + double precision,intent(in) :: eGT(nBas) + double precision,intent(in) :: EcRPA + double precision,intent(in) :: EcGM + + integer :: p,HOMO,LUMO + double precision :: Gap + +! HOMO and LUMO + + HOMO = nO + LUMO = HOMO + 1 + Gap = eGT(LUMO)-eGT(HOMO) + +! Dump results + + write(*,*)'-------------------------------------------------------------------------------' + if(nSCF < 10) then + write(*,'(1X,A21,I1,A3,I1,A12)')' Self-consistent evG',nSCF,'Teh',nSCF,' calculation' + else + write(*,'(1X,A21,I2,A3,I2,A12)')' Self-consistent evG',nSCF,'Teh',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_c (eV)','|','Z','|','e_QP (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,'|',SigC(p)*HaToeV,'|',Z(p),'|',eGT(p)*HaToeV,'|' + enddo + + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A10,I3)') 'Iteration ',nSCF + write(*,'(2X,A14,F15.5)')'Convergence = ',Conv + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A30,F15.6,A3)') 'evGTeh HOMO energy:',eGT(HOMO)*HaToeV,' eV' + write(*,'(2X,A30,F15.6,A3)') 'evGTeh LUMO energy:',eGT(LUMO)*HaToeV,' eV' + write(*,'(2X,A30,F15.6,A3)') 'evGTeh HOMO-LUMO gap :',Gap*HaToeV,' eV' + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A30,F15.6,A3)') 'RPA@evGTeh total energy :',ENuc + ERHF + EcRPA,' au' + write(*,'(2X,A30,F15.6,A3)') 'RPA@evGTeh correlation energy:',EcRPA,' au' + write(*,'(2X,A30,F15.6,A3)') 'GM@evGTeh total energy :',ENuc + ERHF + EcGM,' au' + write(*,'(2X,A30,F15.6,A3)') 'GM@evGTeh correlation energy:',EcGM,' au' + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + +end subroutine diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 5f0a6ed..9043513 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -16,8 +16,8 @@ program QuAcK logical :: doADC logical :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3 logical :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW - logical :: doG0T0,doevGT,doqsGT - logical :: doehG0T0 + logical :: doG0T0pp,doevGTpp,doqsGTpp + logical :: doG0T0eh,doevGTeh,doqsGTeh integer :: nNuc,nBas,nBasCABS integer :: nEl(nspin) @@ -170,8 +170,8 @@ program QuAcK doG0F3,doevGF3, & doG0W0,doevGW,doqsGW,doSRGqsGW, & doufG0W0,doufGW, & - doG0T0,doevGT,doqsGT, & - doehG0T0) + doG0T0pp,doevGTpp,doqsGTpp, & + doG0T0eh,doevGTeh,doqsGTeh) ! Read options for methods @@ -1131,7 +1131,7 @@ program QuAcK eG0T0(:,:) = eHF(:,:) - if(doG0T0) then + if(doG0T0pp) then call cpu_time(start_G0T0) @@ -1164,7 +1164,7 @@ program QuAcK ! Perform evGT calculatiom !------------------------------------------------------------------------ - if(doevGT) then + if(doevGTpp) then call cpu_time(start_evGT) @@ -1197,7 +1197,7 @@ program QuAcK ! Perform qsGT calculation !------------------------------------------------------------------------ - if(doqsGT) then + if(doqsGTpp) then call cpu_time(start_qsGT) @@ -1225,22 +1225,22 @@ program QuAcK end if !------------------------------------------------------------------------ -! Perform ehG0T0 calculatiom +! Perform G0T0eh calculatiom !------------------------------------------------------------------------ eG0T0(:,:) = eHF(:,:) - if(doehG0T0) then + if(doG0T0eh) then call cpu_time(start_G0T0) if(unrestricted) then - print*,'!!! ehG0T0 NYI at the unrestricted level !!!' + print*,'!!! eh G0T0 NYI at the unrestricted level !!!' else - call ehG0T0(doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE,singlet,triplet, & + call G0T0eh(doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE,singlet,triplet, & linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0T0) end if @@ -1253,6 +1253,29 @@ program QuAcK end if +!------------------------------------------------------------------------ +! Perform evGTeh calculation +!------------------------------------------------------------------------ + + if(doevGTeh) then + + call cpu_time(start_evGT) + if(unrestricted) then + + else + + call evGTeh(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS, & + BSE,BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE,singlet,triplet,linGT,eta_GT,regGT, & + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0T0) + end if + call cpu_time(end_evGT) + + t_evGT = end_evGT - start_evGT + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGT = ',t_evGT,' seconds' + write(*,*) + + end if + !------------------------------------------------------------------------ ! Compute FCI !------------------------------------------------------------------------ diff --git a/src/QuAcK/read_methods.f90 b/src/QuAcK/read_methods.f90 index 459516a..451d349 100644 --- a/src/QuAcK/read_methods.f90 +++ b/src/QuAcK/read_methods.f90 @@ -8,8 +8,8 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & doG0F3,doevGF3, & doG0W0,doevGW,doqsGW,doSRGqsGW, & doufG0W0,doufGW, & - doG0T0,doevGT,doqsGT, & - doehG0T0) + doG0T0pp,doevGTpp,doqsGTpp, & + doG0T0eh,doevGTeh,doqsGTeh) ! Read desired methods @@ -25,8 +25,8 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & logical,intent(out) :: doRPA,doRPAx,docrRPA,doppRPA logical,intent(out) :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3 logical,intent(out) :: doG0W0,doevGW,doqsGW,doSRGqsGW,doufG0W0,doufGW - logical,intent(out) :: doG0T0,doevGT,doqsGT - logical,intent(out) :: doehG0T0 + logical,intent(out) :: doG0T0pp,doevGTpp,doqsGTpp + logical,intent(out) :: doG0T0eh,doevGTeh,doqsGTeh ! Local variables @@ -81,9 +81,12 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & doufG0W0 = .false. doufGW = .false. - doG0T0 = .false. - doevGT = .false. - doqsGT = .false. + doG0T0pp = .false. + doevGTpp = .false. + doqsGTpp = .false. + doG0T0eh = .false. + doevGTeh = .false. + doqsGTeh = .false. ! Read mean-field methods @@ -161,11 +164,13 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & ! Read GT methods read(1,*) - read(1,*) answer1,answer2,answer3,answer4 - if(answer1 == 'T') doG0T0 = .true. - if(answer2 == 'T') doevGT = .true. - if(answer3 == 'T') doqsGT = .true. - if(answer4 == 'T') doehG0T0 = .true. + read(1,*) answer1,answer2,answer3,answer4,answer5,answer6 + if(answer1 == 'T') doG0T0pp = .true. + if(answer2 == 'T') doevGTpp = .true. + if(answer3 == 'T') doqsGTpp = .true. + if(answer4 == 'T') doG0T0eh = .true. + if(answer5 == 'T') doevGTeh = .true. + if(answer6 == 'T') doqsGTeh = .true. ! Close file with geometry specification